← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

Viewing changes to src/bin/cints/DFT/calc_den.c

  • 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/bin/cints/DFT/calc_den.c
1
 
/*----------------------------------------------------
2
 
  
3
 
  get_density.c
4
 
  
5
 
  By Shawn Brown
6
 
  
7
 
  The code contains functions that will retrieve 
8
 
  and process the density at a given x,y,z coord
9
 
  in space
10
 
  
11
 
  --------------------------------------------------*/
12
 
#include <math.h>
13
 
#include <stdio.h>
14
 
#include <stdlib.h>
15
 
#include <libipv1/ip_lib.h>
16
 
#include <libciomr/libciomr.h>
17
 
#include <libqt/qt.h>
18
 
#include <libint/libint.h>
19
 
 
20
 
#include"defines.h"
21
 
#define EXTERN
22
 
#include"global.h"
23
 
#include"bas_comp_functions.h"
24
 
 
25
 
struct den_info_s calc_density(struct coordinates geom){
26
 
    
27
 
    int i,j,k;
28
 
    int shell_type;
29
 
    int shell_start;
30
 
    int shell_end;
31
 
    int n_shells;
32
 
    int num_ao,ndocc;
33
 
    int shell_center;
34
 
    double x,y,z;
35
 
    double xa,ya,za;
36
 
    double rr;
37
 
    double rrtmp;
38
 
    double bastmp;
39
 
    double bastmp1;
40
 
    double den_sum;
41
 
    double coeff;
42
 
    double expon;
43
 
    double *norm_ptr;
44
 
    double *dist_atom;
45
 
    double *temp_arr;
46
 
    
47
 
    struct coordinates *dist_coord;
48
 
    struct den_info_s den_info;
49
 
    struct shell_pair *sp;
50
 
    
51
 
    x = geom.x;
52
 
    y = geom.y;
53
 
    z = geom.z;
54
 
    
55
 
    num_ao = BasisSet.num_ao;
56
 
    ndocc = MOInfo.ndocc;
57
 
    temp_arr = init_array(num_ao);
58
 
    dist_atom = init_array(Molecule.num_atoms);
59
 
    dist_coord = (struct coordinates *)malloc(sizeof(struct coordinates)*Molecule.num_atoms);
60
 
    timer_on("distance");
61
 
    for(i=0;i<Molecule.num_atoms;i++){
62
 
        dist_coord[i].x = x-Molecule.centers[i].x;
63
 
        dist_coord[i].y = y-Molecule.centers[i].y;
64
 
        dist_coord[i].z = z-Molecule.centers[i].z;
65
 
        dist_atom[i] = dist_coord[i].x*dist_coord[i].x
66
 
            +dist_coord[i].y*dist_coord[i].y
67
 
            +dist_coord[i].z*dist_coord[i].z;
68
 
    }
69
 
    n_shells = BasisSet.num_shells;
70
 
timer_off("distance");
71
 
timer_on("basis");
72
 
    for(i=k=0;i<n_shells;i++){
73
 
        
74
 
        shell_type = BasisSet.shells[i].am;
75
 
        shell_center = BasisSet.shells[i].center-1;
76
 
        xa = dist_coord[shell_center].x;
77
 
        ya = dist_coord[shell_center].y;
78
 
        za = dist_coord[shell_center].z;
79
 
        rr = dist_atom[shell_center];
80
 
 
81
 
        
82
 
        shell_start = BasisSet.shells[i].fprim-1;
83
 
        shell_end = shell_start+BasisSet.shells[i].n_prims;
84
 
        
85
 
        norm_ptr = GTOs.bf_norm[shell_type-1];
86
 
        timer_on("exponent");
87
 
        /*bastmp1 = calc_exp_basis(i,rr);*/
88
 
        
89
 
        bastmp = 0;
90
 
        for(j=shell_start;j<shell_end;j++){
91
 
            expon = -BasisSet.cgtos[j].exp;
92
 
            coeff = BasisSet.cgtos[j].ccoeff[shell_type-1];
93
 
            bastmp += coeff*exp(expon*rr);
94
 
        }
95
 
        /*if(bastmp != bastmp1)
96
 
            fprintf(outfile,"\nbastmp1 = %10.15lf bastmp2 = %10.15lf for shell %d at rr = %e",bastmp1,bastmp,i,rr);*/
97
 
timer_off("exponent");  
98
 
        /*----------------------------------
99
 
          Compute values of basis functions
100
 
 
101
 
          NOTE: using Psi 3 ordering of
102
 
                functions within shells
103
 
         ----------------------------------*/
104
 
        switch (shell_type) {
105
 
            
106
 
        case 1:
107
 
            DFT_options.basis[k] = norm_ptr[0]*bastmp;
108
 
            k++;
109
 
            break;
110
 
            
111
 
        case 2:
112
 
            DFT_options.basis[k] = norm_ptr[0]*bastmp*xa;
113
 
            k++;
114
 
            
115
 
            DFT_options.basis[k] = norm_ptr[1]*bastmp*ya;
116
 
            k++;
117
 
            
118
 
            DFT_options.basis[k] = norm_ptr[2]*bastmp*za;
119
 
            k++;
120
 
            break;
121
 
        case 3:
122
 
            DFT_options.basis[k] = norm_ptr[0]*bastmp*xa*xa;
123
 
            k++;
124
 
            
125
 
            DFT_options.basis[k] = norm_ptr[1]*bastmp*xa*ya;
126
 
            k++;
127
 
            
128
 
            DFT_options.basis[k] = norm_ptr[2]*bastmp*xa*za;
129
 
            k++;
130
 
            
131
 
            DFT_options.basis[k] = norm_ptr[3]*bastmp*ya*ya;
132
 
            k++;
133
 
            
134
 
            DFT_options.basis[k] = norm_ptr[4]*bastmp*ya*za;
135
 
            k++;
136
 
            
137
 
            DFT_options.basis[k] = norm_ptr[5]*bastmp*za*za;
138
 
            k++;
139
 
            break;
140
 
        case 4:
141
 
            DFT_options.basis[k] = norm_ptr[0]*bastmp*xa*xa*xa;
142
 
            k++;
143
 
            
144
 
            DFT_options.basis[k] = norm_ptr[1]*bastmp*xa*xa*ya;
145
 
            k++;
146
 
            
147
 
            DFT_options.basis[k] = norm_ptr[2]*bastmp*xa*xa*za;
148
 
            k++;
149
 
            
150
 
            DFT_options.basis[k] = norm_ptr[3]*bastmp*xa*ya*ya;
151
 
            k++;
152
 
            
153
 
            DFT_options.basis[k] = norm_ptr[4]*bastmp*xa*ya*za;
154
 
            k++;
155
 
 
156
 
            DFT_options.basis[k] = norm_ptr[5]*bastmp*xa*za*za;
157
 
            k++;
158
 
            
159
 
            DFT_options.basis[k] = norm_ptr[6]*bastmp*ya*ya*ya;
160
 
            k++;
161
 
            
162
 
            DFT_options.basis[k] = norm_ptr[7]*bastmp*ya*ya*za;
163
 
            k++;
164
 
            
165
 
            DFT_options.basis[k] = norm_ptr[8]*bastmp*ya*za*za;
166
 
            k++;
167
 
            
168
 
            DFT_options.basis[k] = norm_ptr[9]*bastmp*za*za*za;
169
 
            k++;
170
 
            break;
171
 
        case 5:
172
 
            DFT_options.basis[k] = norm_ptr[0]*bastmp*xa*xa*xa*xa;
173
 
            k++;
174
 
            
175
 
            DFT_options.basis[k] = norm_ptr[1]*bastmp*xa*xa*xa*ya;
176
 
            k++;
177
 
            
178
 
            DFT_options.basis[k] = norm_ptr[2]*bastmp*xa*xa*xa*za;
179
 
            k++;
180
 
            
181
 
            DFT_options.basis[k] = norm_ptr[3]*bastmp*xa*xa*ya*ya;
182
 
            k++;
183
 
            
184
 
            DFT_options.basis[k] = norm_ptr[4]*bastmp*xa*xa*ya*za;
185
 
            k++;
186
 
            
187
 
            DFT_options.basis[k] = norm_ptr[5]*bastmp*xa*xa*za*za;
188
 
            k++;
189
 
            
190
 
            DFT_options.basis[k] = norm_ptr[6]*bastmp*xa*ya*ya*ya;
191
 
            k++;
192
 
            
193
 
            DFT_options.basis[k] = norm_ptr[7]*bastmp*xa*ya*ya*za;
194
 
            k++;
195
 
            
196
 
            DFT_options.basis[k] = norm_ptr[8]*bastmp*xa*ya*za*za;
197
 
            k++;
198
 
 
199
 
            DFT_options.basis[k] = norm_ptr[9]*bastmp*xa*za*za*za;
200
 
            k++;
201
 
            
202
 
            DFT_options.basis[k] = norm_ptr[10]*bastmp*ya*ya*ya*ya;
203
 
            k++;
204
 
            
205
 
            DFT_options.basis[k] = norm_ptr[11]*bastmp*ya*ya*ya*za;
206
 
            k++;
207
 
            
208
 
            DFT_options.basis[k] = norm_ptr[12]*bastmp*ya*ya*za*za;
209
 
            k++;
210
 
            
211
 
            DFT_options.basis[k] = norm_ptr[13]*bastmp*ya*za*za*za;
212
 
            k++;
213
 
            
214
 
            DFT_options.basis[k] = norm_ptr[14]*bastmp*za*za*za*za;
215
 
            k++;
216
 
            break;
217
 
        default:
218
 
            fprintf(stderr,"Basis Functions of Angular Momentum %d not implemented in get_density function",shell_type);
219
 
            fprintf(outfile,"Basis Functions of Angular Momentum %d not implemented in get_density function",shell_type);
220
 
            punt("");
221
 
        }
222
 
    }
223
 
 
224
 
    timer_off("basis");
225
 
    /* Now contract the basis functions with the AO density matrix elements */
226
 
    timer_on("density"); 
227
 
    
228
 
   if(UserOptions.reftype == rhf){
229
 
       den_sum = 0.0;
230
 
#if USE_BLAS
231
 
       C_DGEMV('t',num_ao,ndocc,1.0,Cocc[0],ndocc,
232
 
               DFT_options.basis,1,0.0,temp_arr,1);
233
 
       for(i=0;i<ndocc;i++)
234
 
       den_sum = C_DDOT(ndocc,temp_arr,1,temp_arr,1);
235
 
#else
236
 
       for(i=0;i<ndocc;i++){
237
 
           for(j=0;j<num_ao;j++){
238
 
               temp_arr[i] += Cocc[j][i]*DFT_options.basis[j];
239
 
           }
240
 
           
241
 
       }
242
 
       dot_arr(temp_arr,temp_arr,MOInfo.ndocc,&den_sum);
243
 
#endif
244
 
       den_info.den = den_sum;
245
 
        
246
 
    }
247
 
   free(temp_arr);
248
 
   timer_off("density");
249
 
   free(dist_coord);
250
 
   free(dist_atom);
251
 
   return den_info;
252
 
}
253
 
 
254
 
        
255