← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/lib/libmints/gshell.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/gshell.cc
 
1
#include <cstdlib>
 
2
#include <cmath>
 
3
#include <libmints/gshell.h>
 
4
 
 
5
#include <libmints/wavefunction.h>
 
6
 
 
7
#include <psiconfig.h>
 
8
 
 
9
using namespace psi;
 
10
 
 
11
extern FILE *outfile;
 
12
 
 
13
double norm(int x1, int x2,double c,double ss)
 
14
{
 
15
    if (x1 < x2) 
 
16
        return norm(x2,x1,c,ss);
 
17
    if (x1 == 1) {
 
18
        if (x2 == 1)
 
19
            return c * ss;
 
20
        else
 
21
            return 0.0;
 
22
    }
 
23
    if (x1 == 0)
 
24
        return ss;
 
25
    return c * ( (x1-1) * norm(x1-2,x2,c,ss) + (x2 * norm(x1-1,x2-1,c,ss)));
 
26
}
 
27
 
 
28
 
 
29
GaussianShell::GaussianShell(int ncn, int nprm, double* e, int* am, GaussianType pure,
 
30
    double** c, int nc, Vector3& center, int start, PrimitiveType pt):
 
31
    nprimitives_(nprm), ncontractions_(ncn), nc_(nc), center_(center), start_(start), sym_transfrom_(0)
 
32
{
 
33
    puream_ = new int[ncontraction()];
 
34
    for (int i=0; i<ncontraction(); ++i) {
 
35
        puream_[i] = (pure == Pure);
 
36
    }
 
37
 
 
38
    copy_data(am, e, c);
 
39
    // Compute the number of basis functions in this shell
 
40
    init_data();
 
41
    
 
42
    // Convert the coefficients to coefficients for unnormalized primitives, if needed
 
43
    // if (pt == Normalized)
 
44
    //     convert_coefficients();
 
45
        
 
46
    // Compute the normalization constants
 
47
    if (pt == Unnormalized)
 
48
        normalize_shell();
 
49
}
 
50
 
 
51
GaussianShell::~GaussianShell()
 
52
{
 
53
    delete[] l_;
 
54
    delete[] puream_;
 
55
    delete[] exp_;
 
56
    
 
57
    for (int i=0; i<ncontractions_; ++i)
 
58
        delete[] coef_[i];
 
59
    
 
60
    delete[] coef_;
 
61
    
 
62
    if (sym_transfrom_)
 
63
        delete[] sym_transfrom_;
 
64
}
 
65
 
 
66
// expects coef to be in primitive x contraction format. the data is transposed here
 
67
void GaussianShell::copy_data(int *l, double *exp, double **coef)
 
68
{
 
69
    l_ = new int[ncontraction()];
 
70
    for (int c=0; c<ncontraction(); ++c)
 
71
        l_[c] = l[c];
 
72
    
 
73
    exp_ = new double[nprimitive()];
 
74
    coef_ = new double*[ncontraction()];
 
75
    for (int p=0; p<nprimitive(); ++p) {
 
76
        exp_[p] = exp[p];
 
77
    }
 
78
    for (int c=0; c<ncontraction(); ++c) {
 
79
        coef_[c] = new double[nprimitive()];
 
80
        for (int p=0; p<nprimitive(); ++p) {
 
81
            // Yes, I want it stored c x p, but it came in from chkpt as p x c
 
82
            coef_[c][p] = coef[p][c];
 
83
        }
 
84
    }
 
85
}
 
86
 
 
87
void GaussianShell::convert_coefficients()
 
88
{
 
89
    int i,gc;
 
90
    double c,ss;
 
91
 
 
92
    // Convert the contraction coefficients from coefficients over
 
93
    // normalized primitives to coefficients over unnormalized primitives
 
94
    for (gc=0; gc<ncontractions_; gc++) {
 
95
        for (i=0; i<nprimitives_; i++) {
 
96
            c = 0.25/exp_[i];
 
97
            ss = pow(M_PI/(exp_[i]+exp_[i]),1.5);
 
98
            coef_[gc][i] *= 1.0/sqrt(::norm(l_[gc],l_[gc],c,ss));
 
99
        }
 
100
    }
 
101
}
 
102
 
 
103
double GaussianShell::shell_normalization(int gs)
 
104
{
 
105
    int i,j;
 
106
    double result,c,ss;
 
107
 
 
108
    result = 0.0;
 
109
    for (i=0; i<nprimitives_; i++) {
 
110
        for (j=0; j<nprimitives_; j++) {
 
111
            c = 0.50/(exp_[i] + exp_[j]);
 
112
            ss = pow(M_PI/(exp_[i]+exp_[j]),1.5);
 
113
            result += coef_[gs][i] * coef_[gs][j] *
 
114
                ::norm(l_[gs],l_[gs],c,ss);
 
115
        }
 
116
    }
 
117
 
 
118
    return 1.0/sqrt(result);
 
119
}
 
120
 
 
121
void GaussianShell::normalize_shell()
 
122
{
 
123
    int i, gs;
 
124
    
 
125
    for (gs = 0; gs < ncontractions_; ++gs) {
 
126
        double normalization = shell_normalization(gs);
 
127
        for (i = 0; i < nprimitives_; ++i) {
 
128
            coef_[gs][i] *= normalization;
 
129
            #ifdef DEBUG
 
130
            fprintf(outfile, "New c = %20.16f for %20.16f\n", coef_[gs][i], exp_[i]);
 
131
            #endif
 
132
        }
 
133
    }
 
134
}
 
135
 
 
136
int GaussianShell::nfunction(int c) const
 
137
{
 
138
    return INT_NFUNC(puream_[c], l_[c]);
 
139
}
 
140
 
 
141
void GaussianShell::init_data()
 
142
{
 
143
    int max = 0;
 
144
    int min = 0;
 
145
    int nc = 0;
 
146
    int nf = 0;
 
147
    has_pure_ = false;
 
148
    
 
149
    for (int i=0; i<ncontraction(); ++i) {
 
150
        int maxi = l_[i];
 
151
        if (max < maxi)
 
152
            max = maxi;
 
153
            
 
154
        int mini = l_[i];
 
155
        if (min > mini || i == 0)
 
156
            min = mini;
 
157
            
 
158
        nc += ncartesian(i);
 
159
        nf += nfunction(i);
 
160
        
 
161
        if (is_pure(i))
 
162
            has_pure_ = true;
 
163
    }
 
164
    
 
165
    max_am_ = max;
 
166
    min_am_ = min;
 
167
    ncartesians_ = nc;
 
168
    nfunctions_ = nf;
 
169
}
 
170
 
 
171
void GaussianShell::print(FILE *out) const
 
172
{
 
173
    fprintf(out, "      Number of contractions: %d\n", ncontraction());
 
174
    fprintf(out, "      Number of primitives: %d\n", nprimitive());
 
175
    fprintf(out, "      Number of Cartesian Gaussians: %d\n", ncartesian());
 
176
    fprintf(out, "      Spherical Harmonics?: %s\n", has_pure() ? "true" : "false");
 
177
    if (max_am() == min_am())
 
178
        fprintf(out, "      Angular momentum: %d\n", max_am());
 
179
    else {
 
180
        fprintf(out, "      Max angular momentum: %d\n", max_am());
 
181
        fprintf(out, "      Min angular momentum: %d\n", min_am());
 
182
    }
 
183
    fprintf(out, "      Exponent       ");
 
184
    for(int c=0; c<ncontraction(); c++)
 
185
        fprintf(out, " Contr. %3d  ",c);
 
186
    fprintf(out, "\n");
 
187
    for(int p=0; p<nprimitive(); p++) {
 
188
        fprintf(out, "      %15.10f",exp_[p]);
 
189
        for(int c=0; c<ncontraction(); c++)
 
190
            fprintf(out, " %12.9f",coef_[c][p]);
 
191
        fprintf(out, "\n");
 
192
    }
 
193
    fprintf(out, "\n");
 
194
}
 
195
 
 
196
double GaussianShell::normalize(int l, int m, int n)
 
197
{
 
198
    static int use_cca_integrals_standard = (PSI_INTEGRALS_STANDARD == 1);
 
199
    if (use_cca_integrals_standard) {
 
200
        return 1.0;
 
201
    } else {
 
202
        double numer = df[2*l_[0]];
 
203
        double denom = df[2*l] * df[2*m] * df[2*n];
 
204
        // printf("l_=%d, l=%d, m=%d, n=%d, norm=%15.10f\n", l_[0], l, m, n, sqrt(numer/denom));
 
205
        return sqrt(numer/denom);        
 
206
    }
 
207
}
 
208
 
 
209
void GaussianShell::set_sym_transform(int nirreps, int *vec)
 
210
{
 
211
    if (sym_transfrom_)
 
212
        delete[] sym_transfrom_;
 
213
    
 
214
    sym_transfrom_ = new int[nirreps];
 
215
    for (int i=0; i<nirreps; ++i)
 
216
        sym_transfrom_[i] = vec[0];
 
217
}