← Back to branch summary

~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/lib/libmints/transform.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/libmints/transform.cc
 
1
#include <cmath>
 
2
 
 
3
#include <libmints/wavefunction.h>
 
4
 
 
5
#include <libmints/basisset.h>
 
6
#include <libmints/onebody.h>
 
7
#include <libmints/integral.h>
 
8
 
 
9
#include <psiconfig.h>
 
10
 
 
11
using namespace psi;
 
12
 
 
13
extern FILE *outfile;
 
14
 
 
15
#define parity(m) ((m)%2 ? -1 : 1) // returns (-1)^m
 
16
#define MAX(a, b) (((a) > (b)) ? (a) : (b))
 
17
#define MIN(a, b) (((a) < (b)) ? (a) : (b))
 
18
 
 
19
double xyz2lm_coeff(int l, int m, int lx, int ly, int lz)
 
20
{
 
21
    static int use_cca_integrals_standard = (PSI_INTEGRALS_STANDARD == 1);
 
22
    int i, j, k, i_max;
 
23
    int k_min, k_max;
 
24
    int abs_m;
 
25
    int comp;
 
26
    double pfac, pfac1, sum, sum1;
 
27
 
 
28
    abs_m = abs(m);
 
29
    if ((lx + ly - abs(m))%2)
 
30
        return 0.0;
 
31
    else
 
32
        j = (lx + ly - abs(m))/2;
 
33
 
 
34
    if (j < 0)
 
35
        return 0.0;
 
36
 
 
37
    comp = (m >= 0) ? 1 : -1;
 
38
    i = abs_m-lx;
 
39
    if (comp != parity(i))
 
40
        return 0.0;
 
41
 
 
42
    pfac = sqrt(fac[2*lx]*fac[2*ly]*fac[2*lz]*fac[l-abs_m]/(fac[2*l]*fac[l]
 
43
        *fac[lx]*fac[ly]*fac[lz]*fac[l+abs_m]));
 
44
    pfac /= (1 << l);
 
45
 
 
46
    if (m < 0)
 
47
        pfac *= parity((i-1)/2);
 
48
    else
 
49
        pfac *= parity(i/2);
 
50
 
 
51
    i_max = (l-abs_m)/2;
 
52
    sum = 0.0;
 
53
    for (i=0; i<=i_max; i++) {
 
54
        pfac1 = bc[l][i]*bc[i][j];
 
55
        if (pfac1 == 0.0)
 
56
            continue;
 
57
        else
 
58
            pfac1 *= (parity(i)*fac[2*(l-i)]/fac[l-abs_m-2*i]);
 
59
        sum1 = 0.0;
 
60
        k_min = MAX((lx-abs_m)/2,0);
 
61
        k_max = MIN(j,lx/2);
 
62
        for (k=k_min; k<=k_max; k++)
 
63
            sum1 += bc[j][k]*bc[abs_m][lx-2*k]*parity(k);
 
64
        sum += pfac1*sum1;
 
65
    }
 
66
 
 
67
    if (use_cca_integrals_standard)
 
68
        sum *= sqrt(df[2*l]/(df[2*lx]*df[2*ly]*df[2*lz]));
 
69
 
 
70
    if (m == 0)
 
71
        return pfac*sum;
 
72
    else
 
73
        return M_SQRT2*pfac*sum;
 
74
}
 
75
 
 
76
void SphericalTransformComponent::init(int a, int b, int c, double coef, int cartindex, int pureindex)
 
77
{
 
78
    // printf("a = %d, b = %d, c = %d, coef = %f, cartindex = %d, pureindex = %d\n", a, b, c, coef, cartindex, pureindex);
 
79
    a_ = a;
 
80
    b_ = b;
 
81
    c_ = c;
 
82
    coef_ = coef;
 
83
    cartindex_ = cartindex;
 
84
    pureindex_ = pureindex;
 
85
}
 
86
 
 
87
SphericalTransform::SphericalTransform(int l)
 
88
{    
 
89
    int pureindex = 0;
 
90
    int cartindex = 0;
 
91
 
 
92
    l_ = l;
 
93
    
 
94
    // Go through all the spherical harmonic possibilities
 
95
    for (int m = -l_; m <= l_; ++m) {
 
96
        // Generate all the possible cartesian functions
 
97
        cartindex = 0;
 
98
        for(int ii = 0; ii <= l_; ii++) {
 
99
            int l1 = l_ - ii;
 
100
            for(int jj = 0; jj <= ii; jj++) {
 
101
                int m1 = ii - jj;
 
102
                int n1 = jj;
 
103
                
 
104
                // Compute the coefficient needed to convert from cartesian to spherical
 
105
                double coef = xyz2lm_coeff(l, m, l1, m1, n1);
 
106
                
 
107
                // Save the information if it is nonzero
 
108
                if (fabs(coef) > 1.0e-16) {
 
109
                    SphericalTransformComponent component;
 
110
                    component.init(l1, m1, n1, coef, cartindex, pureindex);
 
111
                    components_.push_back(component);
 
112
                }
 
113
                
 
114
                cartindex++;
 
115
            }
 
116
        }
 
117
        pureindex++;
 
118
    }
 
119
}