← 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_grad_fast.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

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