← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/DFT/calc_den_fast.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_fast.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_fast(struct coordinates geom, int atom_num){
26
 
    
27
 
    int i,j,k,l,m;
28
 
    int am2shell;
29
 
    int shell_type;
30
 
    int shell_start;
31
 
    int shell_end;
32
 
    int n_shells;
33
 
    int num_ao,ndocc;
34
 
    int shell_center;
35
 
    double x,y,z;
36
 
    double xa,ya,za;
37
 
    double rr;
38
 
    double rrtmp;
39
 
    double bastmp;
40
 
    double den_sum;
41
 
    double coeff;
42
 
    double expon;
43
 
    double *norm_ptr;
44
 
    double *dist_atom;
45
 
    double *temp_arr;
46
 
    double *coord_array[14];
47
 
    double *expon_arr;
48
 
    double norm_ptr0,norm_ptr1,norm_ptr2;
49
 
    double norm_ptr3,norm_ptr4,norm_ptr5;
50
 
    double norm_ptr6,norm_ptr7,norm_ptr8;
51
 
    double norm_ptr9,norm_ptr10,norm_ptr11;
52
 
    double norm_ptr12,norm_ptr13,norm_ptr14;
53
 
    double xx,yy,zz,xy,xz,yz;
54
 
    double bxx,byy,bzz,bxy,bxz,byz;
55
 
 
56
 
    
57
 
    struct coordinates *dist_coord;
58
 
    struct leb_chunk_s *chunk;
59
 
    struct den_info_s den_info;
60
 
    struct close_shell_info_s *close;
61
 
    
62
 
    x = geom.x;
63
 
    y = geom.y;
64
 
    z = geom.z;
65
 
    
66
 
    num_ao = BasisSet.num_ao;
67
 
    ndocc = MOInfo.ndocc;
68
 
    temp_arr = init_array(ndocc);
69
 
    close = &(DFT_options.close_shell_info);
70
 
    /* ---------------------------------
71
 
       Compute distances from atom that 
72
 
       the basis function is centered on
73
 
       to the grid point
74
 
       --------------------------------*/
75
 
    zero_arr(DFT_options.basis,num_ao);
76
 
    dist_atom = init_array(Molecule.num_atoms);
77
 
    dist_coord = (struct coordinates *)
78
 
        malloc(sizeof(struct coordinates)*Molecule.num_atoms);
79
 
    timer_on("distance");
80
 
    for(i=0;i<Molecule.num_atoms;i++){
81
 
        dist_coord[i].x = x-Molecule.centers[i].x;
82
 
        dist_coord[i].y = y-Molecule.centers[i].y;
83
 
        dist_coord[i].z = z-Molecule.centers[i].z;
84
 
        dist_atom[i] = dist_coord[i].x*dist_coord[i].x
85
 
            +dist_coord[i].y*dist_coord[i].y
86
 
            +dist_coord[i].z*dist_coord[i].z;
87
 
    }
88
 
    
89
 
    n_shells = BasisSet.num_shells;
90
 
    timer_off("distance");
91
 
    timer_on("basis");
92
 
    
93
 
    l=0;
94
 
    m=0;
95
 
    for(i=0;i<BasisSet.max_am;i++){
96
 
        
97
 
        norm_ptr = GTOs.bf_norm[i];
98
 
        shell_type = i;
99
 
        switch(i){
100
 
            
101
 
            /* ------------------------------------------------------*/
102
 
            /* S-type functions */
103
 
            
104
 
        case 0:
105
 
            norm_ptr0 = norm_ptr[0];
106
 
            for(j=0;j<close->close_shells_per_am[i];j++){
107
 
                
108
 
                am2shell = close->shells_close_to_chunk[m];
109
 
                shell_center = BasisSet.shells[am2shell].center-1;
110
 
                xa = dist_coord[shell_center].x;
111
 
                ya = dist_coord[shell_center].y;
112
 
                za = dist_coord[shell_center].z;
113
 
                rr = dist_atom[shell_center];
114
 
                
115
 
                shell_start = BasisSet.shells[am2shell].fprim-1;
116
 
                shell_end = shell_start
117
 
                    +BasisSet.shells[am2shell].n_prims;
118
 
                
119
 
                
120
 
                timer_on("exponent");
121
 
                bastmp = 0.0;
122
 
                
123
 
                /*bastmp = calc_exp_basis(am2shell,rr);*/
124
 
                for(k=shell_start;k<shell_end;k++){
125
 
                    expon = -BasisSet.cgtos[k].exp;
126
 
                    coeff = BasisSet.cgtos[k].ccoeff[shell_type];
127
 
                    bastmp += coeff*exp(expon*rr);
128
 
      
129
 
                }
130
 
                
131
 
                timer_off("exponent");
132
 
                
133
 
                DFT_options.basis[l] = norm_ptr0*bastmp;
134
 
                l++;
135
 
                
136
 
                m++;
137
 
            }
138
 
 
139
 
            break;
140
 
            /* ------------------------------------------------------*/
141
 
            
142
 
            
143
 
            /* ------------------------------------------------------*/
144
 
            /* P-type functions */
145
 
            
146
 
        case 1:
147
 
            norm_ptr0 = norm_ptr[0];
148
 
            norm_ptr1 = norm_ptr[1];
149
 
            norm_ptr2 = norm_ptr[2];
150
 
            
151
 
            for(j=0;j<close->close_shells_per_am[i];j++){
152
 
                
153
 
                am2shell = close->shells_close_to_chunk[m];
154
 
                shell_center = BasisSet.shells[am2shell].center-1;
155
 
 
156
 
                xa = dist_coord[shell_center].x;
157
 
                ya = dist_coord[shell_center].y;
158
 
                za = dist_coord[shell_center].z;
159
 
                rr = dist_atom[shell_center];
160
 
                
161
 
                shell_start = BasisSet.shells[am2shell].fprim-1;
162
 
                shell_end = shell_start
163
 
                    +BasisSet.shells[am2shell].n_prims;
164
 
                
165
 
                
166
 
                timer_on("exponent");
167
 
                bastmp = 0.0;
168
 
                for(k=shell_start;k<shell_end;k++){
169
 
                    expon = -BasisSet.cgtos[k].exp*rr;
170
 
                    coeff = BasisSet.cgtos[k].ccoeff[shell_type];
171
 
                    bastmp += coeff*exp(expon);
172
 
                    
173
 
                }
174
 
 
175
 
                timer_off("exponent");
176
 
                
177
 
                DFT_options.basis[l] = norm_ptr0*bastmp*xa;
178
 
                l++;
179
 
                
180
 
                DFT_options.basis[l] = norm_ptr1*bastmp*ya;
181
 
                l++;
182
 
                
183
 
                DFT_options.basis[l] = norm_ptr2*bastmp*za;
184
 
                l++;
185
 
                m++;
186
 
            }
187
 
            break;
188
 
            /* ------------------------------------------------------*/
189
 
            
190
 
            /* ------------------------------------------------------*/
191
 
            /* D-type functions */ 
192
 
        
193
 
        case 2:
194
 
            norm_ptr0 = norm_ptr[0];
195
 
            norm_ptr1 = norm_ptr[1];
196
 
            norm_ptr2 = norm_ptr[2];
197
 
            norm_ptr3 = norm_ptr[3];
198
 
            norm_ptr4 = norm_ptr[4];
199
 
            norm_ptr5 = norm_ptr[5];
200
 
            for(j=0;j<close->close_shells_per_am[i];j++){
201
 
                am2shell = close->shells_close_to_chunk[m];
202
 
                shell_center = BasisSet.shells[am2shell].center-1;
203
 
                xa = dist_coord[shell_center].x;
204
 
                ya = dist_coord[shell_center].y;
205
 
                za = dist_coord[shell_center].z;
206
 
                rr = dist_atom[shell_center];
207
 
                
208
 
                shell_start = BasisSet.shells[am2shell].fprim-1;
209
 
                shell_end = shell_start
210
 
                    +BasisSet.shells[am2shell].n_prims;
211
 
                
212
 
                
213
 
                timer_on("exponent");
214
 
                bastmp = 0.0;
215
 
                for(k=shell_start;k<shell_end;k++){
216
 
                    expon = -BasisSet.cgtos[k].exp*rr;
217
 
                    coeff = BasisSet.cgtos[k].ccoeff[shell_type];
218
 
                    bastmp += coeff*exp(expon);
219
 
                   
220
 
                }
221
 
                timer_off("exponent");
222
 
                
223
 
                DFT_options.basis[l] = norm_ptr0*bastmp*xa*xa;
224
 
                l++;
225
 
                
226
 
                DFT_options.basis[l] = norm_ptr1*bastmp*xa*ya;
227
 
                l++;
228
 
                
229
 
                DFT_options.basis[l] = norm_ptr2*bastmp*xa*za;
230
 
                l++;
231
 
                
232
 
                DFT_options.basis[l] = norm_ptr3*bastmp*ya*ya;
233
 
                l++;
234
 
                
235
 
                DFT_options.basis[l] = norm_ptr4*bastmp*ya*za;
236
 
                l++;
237
 
                
238
 
                DFT_options.basis[l] = norm_ptr5*bastmp*za*za;
239
 
                l++;
240
 
                m++;
241
 
            }
242
 
            break;
243
 
            /* ------------------------------------------------------*/
244
 
            
245
 
            /* ------------------------------------------------------*/
246
 
            /* F-type functions */
247
 
            
248
 
        case 3:
249
 
            norm_ptr0 = norm_ptr[0];
250
 
            norm_ptr1 = norm_ptr[1];
251
 
            norm_ptr2 = norm_ptr[2];
252
 
            norm_ptr3 = norm_ptr[3];
253
 
            norm_ptr4 = norm_ptr[4];
254
 
            norm_ptr5 = norm_ptr[5];
255
 
            norm_ptr6 = norm_ptr[6];
256
 
            norm_ptr7 = norm_ptr[7];
257
 
            norm_ptr8 = norm_ptr[8];
258
 
            norm_ptr9 = norm_ptr[9];
259
 
            
260
 
            for(j=0;j<close->close_shells_per_am[i];j++){
261
 
                am2shell = close->shells_close_to_chunk[m];
262
 
                shell_center = BasisSet.shells[am2shell].center-1;
263
 
                xa = dist_coord[shell_center].x;
264
 
                ya = dist_coord[shell_center].y;
265
 
                za = dist_coord[shell_center].z;
266
 
                rr = dist_atom[shell_center];
267
 
                
268
 
                xx = xa*xa;
269
 
                xy = xa*ya;
270
 
                yy = ya*ya;
271
 
                zz = za*za;
272
 
                
273
 
                bxx = bastmp*xx;
274
 
                byy = bastmp*yy;
275
 
                bxy = bastmp*xy;
276
 
                bzz = bastmp*zz;
277
 
                
278
 
                shell_start = BasisSet.shells[am2shell].fprim-1;
279
 
                shell_end = shell_start
280
 
                    +BasisSet.shells[am2shell].n_prims;
281
 
                
282
 
                
283
 
                timer_on("exponent");
284
 
                bastmp = 0.0;
285
 
                for(k=shell_start;k<shell_end;k++){
286
 
                    expon = -BasisSet.cgtos[k].exp*rr;
287
 
                    coeff = BasisSet.cgtos[k].ccoeff[shell_type];
288
 
                    bastmp += coeff*exp(expon);
289
 
                }
290
 
                timer_off("exponent");
291
 
                
292
 
                DFT_options.basis[l] = norm_ptr0*bxx*xa;
293
 
                l++;
294
 
                
295
 
                DFT_options.basis[l] = norm_ptr1*bxx*ya;
296
 
                l++;
297
 
                
298
 
                DFT_options.basis[l] = norm_ptr2*bxx*za;
299
 
                l++;
300
 
                
301
 
                DFT_options.basis[l] = norm_ptr3*byy*xa;
302
 
                l++;
303
 
                
304
 
                DFT_options.basis[l] = norm_ptr4*bxy*za;
305
 
                l++;
306
 
                
307
 
                DFT_options.basis[l] = norm_ptr5*bzz*xa;
308
 
                l++;
309
 
                
310
 
                DFT_options.basis[l] = norm_ptr6*byy*ya;
311
 
                l++;
312
 
                
313
 
                DFT_options.basis[l] = norm_ptr7*byy*za;
314
 
                l++;
315
 
                
316
 
                DFT_options.basis[l] = norm_ptr8*bzz*ya;
317
 
                l++;
318
 
                
319
 
                DFT_options.basis[l] = norm_ptr9*bzz*za;
320
 
                l++;
321
 
                m++;
322
 
            }
323
 
            
324
 
            break;
325
 
            
326
 
            /* ------------------------------------------------------*/
327
 
            
328
 
            /* ------------------------------------------------------*/
329
 
            /* G-type functions */
330
 
        case 4:
331
 
            norm_ptr0 = norm_ptr[0];
332
 
            norm_ptr1 = norm_ptr[1];
333
 
            norm_ptr2 = norm_ptr[2];
334
 
            norm_ptr3 = norm_ptr[3];
335
 
            norm_ptr4 = norm_ptr[4];
336
 
            norm_ptr5 = norm_ptr[5];
337
 
            norm_ptr6 = norm_ptr[6];
338
 
            norm_ptr7 = norm_ptr[7];
339
 
            norm_ptr8 = norm_ptr[8];
340
 
            norm_ptr9 = norm_ptr[9];
341
 
            norm_ptr10 = norm_ptr[10];
342
 
            norm_ptr11 = norm_ptr[11];
343
 
            norm_ptr12 = norm_ptr[12];
344
 
            norm_ptr13 = norm_ptr[13];
345
 
            norm_ptr14 = norm_ptr[14];
346
 
            for(j=0;j<close->close_shells_per_am[i];j++){
347
 
                am2shell = close->shells_close_to_chunk[m];
348
 
                shell_center = BasisSet.shells[am2shell].center-1;
349
 
                xa = dist_coord[shell_center].x;
350
 
                ya = dist_coord[shell_center].y;
351
 
                za = dist_coord[shell_center].z;
352
 
                rr = dist_atom[shell_center];
353
 
                
354
 
                shell_start = BasisSet.shells[am2shell].fprim-1;
355
 
                shell_end = shell_start
356
 
                    +BasisSet.shells[am2shell].n_prims;
357
 
                
358
 
                
359
 
                timer_on("exponent");
360
 
                bastmp = 0.0;
361
 
                for(k=shell_start;k<shell_end;k++){
362
 
                    expon = -BasisSet.cgtos[k].exp;
363
 
                    coeff = BasisSet.cgtos[k].ccoeff[shell_type];
364
 
                    bastmp += coeff*exp(expon*rr);
365
 
                }
366
 
                timer_off("exponent");
367
 
                DFT_options.basis[l] = norm_ptr0*bastmp*xa*xa*xa*xa;
368
 
                l++;
369
 
                
370
 
                DFT_options.basis[l] = norm_ptr1*bastmp*xa*xa*xa*ya;
371
 
                l++;
372
 
                
373
 
                DFT_options.basis[l] = norm_ptr2*bastmp*xa*xa*xa*za;
374
 
                l++;
375
 
                
376
 
                DFT_options.basis[l] = norm_ptr3*bastmp*xa*xa*ya*ya;
377
 
                l++;
378
 
                
379
 
                DFT_options.basis[l] = norm_ptr4*bastmp*xa*xa*ya*za;
380
 
                l++;
381
 
                
382
 
                DFT_options.basis[l] = norm_ptr5*bastmp*xa*xa*za*za;
383
 
                l++;
384
 
                
385
 
                DFT_options.basis[l] = norm_ptr6*bastmp*xa*ya*ya*ya;
386
 
                l++;
387
 
                
388
 
                DFT_options.basis[l] = norm_ptr7*bastmp*xa*ya*ya*za;
389
 
                l++;
390
 
                
391
 
                DFT_options.basis[l] = norm_ptr8*bastmp*xa*ya*za*za;
392
 
                l++;
393
 
                
394
 
                DFT_options.basis[l] = norm_ptr9*bastmp*xa*za*za*za;
395
 
                l++;
396
 
                
397
 
                DFT_options.basis[l] = norm_ptr10*bastmp*ya*ya*ya*ya;
398
 
                l++;
399
 
                
400
 
                DFT_options.basis[l] = norm_ptr11*bastmp*ya*ya*ya*za;
401
 
                l++;
402
 
                
403
 
                DFT_options.basis[l] = norm_ptr12*bastmp*ya*ya*za*za;
404
 
                l++;
405
 
                
406
 
                DFT_options.basis[l] = norm_ptr13*bastmp*ya*za*za*za;
407
 
                l++;
408
 
                
409
 
                DFT_options.basis[l] = norm_ptr14*bastmp*za*za*za*za;
410
 
                l++;
411
 
                m++;
412
 
            }
413
 
 
414
 
            break;
415
 
            
416
 
        default:
417
 
            fprintf(stderr,"Basis Functions of Angular Momentum \n \
418
 
%d not implemented in get_density function",shell_type);
419
 
            fprintf(outfile,"Basis Functions of Angular Momentum \n \
420
 
%d not implemented in get_density function",shell_type);
421
 
            punt("");
422
 
        }
423
 
        
424
 
    }
425
 
    timer_off("basis"); 
426
 
/* Now contract the basis functions with the AO density matrix elements */
427
 
    timer_on("density"); 
428
 
    /*for(i=0;i<close->num_close_aos;i++)
429
 
        fprintf(outfile,"\nBasis[%d] = %10.10lf",i,DFT_options.basis[i]);*/
430
 
   
431
 
    den_sum = 0.0; 
432
 
#if USE_BLAS
433
 
    C_DGEMV('t',close->num_close_aos,ndocc,1.0,close->close_COCC[0],ndocc,
434
 
            DFT_options.basis,1,0.0,temp_arr,1);
435
 
    den_sum = C_DDOT(ndocc,temp_arr,1,temp_arr,1);
436
 
#else
437
 
    for(i=0;i<ndocc;i++){
438
 
        for(j=0;j<close->num_close_aos;j++){
439
 
            temp_arr[i] += close->close_COCC[j][i]*DFT_options.basis[j];
440
 
        }
441
 
        /*fprintf(outfile,"\ntemp_arr[%d] = %10.10lf",i,temp_arr[i]);*/
442
 
    }
443
 
    dot_arr(temp_arr,temp_arr,MOInfo.ndocc,&den_sum);
444
 
#endif
445
 
    den_info.den = den_sum;
446
 
    /*fprintf(outfile,"\nden = %10.10lf",den_info.den);*/
447
 
    
448
 
    
449
 
    free(temp_arr);
450
 
    timer_off("density");
451
 
    free(dist_coord);
452
 
    free(dist_atom);
453
 
    return den_info;
454
 
}
455
 
 
456
 
        
457