← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/Default_Deriv1/te_deriv1_corr.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/Default_Deriv1/te_deriv1_corr.c
 
1
#include <math.h>
 
2
#include <stdio.h>
 
3
#include <string.h>
 
4
#include <memory.h>
 
5
#include <stdlib.h>
 
6
#include <libipv1/ip_lib.h>
 
7
#include <libiwl/iwl.h>
 
8
#include <libciomr/libciomr.h>
 
9
#include <libint/libint.h>
 
10
#include <libderiv/libderiv.h>
 
11
#include "defines.h"
 
12
#define EXTERN
 
13
#include "global.h"
 
14
#ifdef USE_TAYLOR_FM
 
15
  #include"taylor_fm_eval.h"
 
16
#else
 
17
  #include"int_fjt.h"
 
18
#endif
 
19
#include "deriv1_quartet_data.h"
 
20
#include "small_fns.h"
 
21
 
 
22
 
 
23
void te_deriv1_corr()
 
24
{
 
25
  /*--- Various data structures ---*/
 
26
  struct iwlbuf TPDM;                 /* IWL buffer for two-pdm matrix elements */
 
27
  struct shell_pair *sp_ij, *sp_kl;
 
28
  Libderiv_t Libderiv;                    /* Integrals library object */
 
29
#ifndef USE_TAYLOR_FM
 
30
  double_array_t fjt_table;               /* table of auxiliary function F_m(u) for each primitive combination */
 
31
#endif
 
32
 
 
33
  int ij, kl, ik, jl, ijkl;
 
34
  int ioffset, joffset, koffset, loffset;
 
35
  int count ;
 
36
  int dum;
 
37
  int n, num;
 
38
  int total_am, am;
 
39
  int orig_am[4];
 
40
  register int i, j, k, l, m, ii, jj, kk, ll;
 
41
  register int si, sj, sk, sl ;
 
42
  register int sii, sjj, skk, sll, slll;
 
43
  register int pi, pj, pk, pl ;
 
44
  int max_pj, max_pl;
 
45
  register int pii, pjj, pkk, pll ;
 
46
  int switch_ij, switch_kl, switch_ijkl;
 
47
  int center_i, center_j, center_k, center_l;
 
48
 
 
49
  int class_size;
 
50
  int max_class_size;
 
51
  int max_cart_class_size;
 
52
 
 
53
  int bf_i, bf_j, bf_k, bf_l, so_i, so_j, so_k, so_l, s;
 
54
  int np_i, np_j, np_k, np_l;
 
55
  int ni, nj, nk, nl, quartet_size;
 
56
 
 
57
  int si_fao, sj_fao, sk_fao, sl_fao;
 
58
  int sii_fao, sjj_fao, skk_fao, sll_fao;
 
59
  int ao_i, imax, ao_j, jmax, ao_k, kmax, ao_l, lmax;
 
60
 
 
61
  int index;
 
62
  int iimax, jjmax, kkmax, llmax;
 
63
  int irrep, npi_ij, npi_kl, npi_ik, npi_jl, ind_offset;
 
64
 
 
65
  int num_prim_comb, p, max_num_prim_comb;
 
66
 
 
67
  int buf_offset, buf_4offset, buf_size, last_buf;
 
68
  int quartet_done, offset;
 
69
 
 
70
  int mosh_i, mosh_j;
 
71
 
 
72
  double AB2, CD2;
 
73
  double *FourInd;
 
74
  double **grad_te_local;
 
75
  double pfac;
 
76
  double temp;
 
77
  double alpha, beta;
 
78
  double **dens_i, **dens_j;
 
79
  double ddax, dday, ddaz, ddbx, ddby, ddbz,
 
80
         ddcx, ddcy, ddcz, dddx, dddy, dddz;
 
81
 
 
82
  /*---------------
 
83
    Initialization
 
84
   ---------------*/
 
85
  iwl_buf_init(&TPDM, IOUnits.itapG, 0.0, 1, 1);
 
86
  buf_offset = 0;
 
87
  buf_4offset = 0;
 
88
  buf_size = TPDM.inbuf;
 
89
#ifdef USE_TAYLOR_FM
 
90
  init_Taylor_Fm_Eval(BasisSet.max_am*4-4+DERIV_LVL,UserOptions.cutoff);
 
91
#else
 
92
  init_fjt(BasisSet.max_am*4+DERIV_LVL);
 
93
  init_fjt_table(&fjt_table);
 
94
#endif
 
95
  init_libderiv_base();
 
96
 
 
97
  max_cart_class_size = ioff[BasisSet.max_am]*ioff[BasisSet.max_am]*ioff[BasisSet.max_am]*ioff[BasisSet.max_am];
 
98
  max_class_size = max_cart_class_size;
 
99
  max_num_prim_comb = (BasisSet.max_num_prims*BasisSet.max_num_prims)*
 
100
                      (BasisSet.max_num_prims*BasisSet.max_num_prims);
 
101
  init_libderiv1(&Libderiv,BasisSet.max_am-1,max_num_prim_comb,max_cart_class_size);
 
102
  FourInd = init_array(max_cart_class_size);
 
103
 
 
104
  grad_te_local = block_matrix(Molecule.num_atoms,3);
 
105
 
 
106
/*-------------------------------------------------
 
107
  generate all unique shell quartets with ordering
 
108
  suitable for building the PK-matrix
 
109
 -------------------------------------------------*/
 
110
  for (sii=0; sii<BasisSet.num_shells; sii++)
 
111
    for (sjj=0; sjj<=sii; sjj++)
 
112
      for (skk=0; skk<=sii; skk++)
 
113
        for (sll=0; sll<= ((sii == skk) ? sjj : skk); sll++){
 
114
 
 
115
            si = sii; sj = sjj; sk = skk; sl = sll;
 
116
 
 
117
            /*--- Skip this quartet if all four centers are the same ---*/
 
118
            if (BasisSet.shells[si].center == BasisSet.shells[sj].center &&
 
119
                BasisSet.shells[si].center == BasisSet.shells[sk].center &&
 
120
                BasisSet.shells[si].center == BasisSet.shells[sl].center) {
 
121
              /*--- If reading in density - need to skip the appropriate block of that too ---*/
 
122
/*            if (read_dens) {
 
123
                last_buf = TPDM.lastbuf;
 
124
                quartet_done = 0;
 
125
                do {
 
126
                  if (buf_offset < buf_size) {
 
127
                    i = TPDM.labels[buf_4offset]   - si_fao;
 
128
                    if (i < 0)
 
129
                      quartet_done = 1;
 
130
                    buf_offset++;
 
131
                    buf_4offset += 4;
 
132
                  }
 
133
                  else if (!last_buf) {
 
134
                    iwl_buf_fetch(&TPDM);
 
135
                    buf_offset = 0;
 
136
                    buf_4offset = 0;
 
137
                    last_buf = TPDM.lastbuf;
 
138
                  }
 
139
                  else {
 
140
                    punt(fpo,"  The last TPDM quartet not marked\n");
 
141
                  }
 
142
                } while (!quartet_done);
 
143
              }*/
 
144
 
 
145
              continue;
 
146
            }
 
147
 
 
148
            switch_ij = 0;
 
149
            switch_kl = 0;
 
150
            switch_ijkl = 0;
 
151
            /* place in "ascending" angular mom-
 
152
               my simple way of optimizing PHG recursion (VRR) */
 
153
            /* these first two are good for the HRR */
 
154
            if(BasisSet.shells[si].am < BasisSet.shells[sj].am){
 
155
              dum = si;
 
156
              si = sj;
 
157
              sj = dum;
 
158
              switch_ij = 1;
 
159
            }
 
160
            if(BasisSet.shells[sk].am < BasisSet.shells[sl].am){
 
161
              dum = sk;
 
162
              sk = sl;
 
163
              sl = dum;
 
164
              switch_kl = 1;
 
165
            }
 
166
            /* this should be /good/ for the VRR */
 
167
            if(BasisSet.shells[si].am + BasisSet.shells[sj].am > BasisSet.shells[sk].am + BasisSet.shells[sl].am){
 
168
              dum = si;
 
169
              si = sk;
 
170
              sk = dum;
 
171
              dum = sj;
 
172
              sj = sl;
 
173
              sl = dum;
 
174
              switch_ijkl = 1;
 
175
            }
 
176
 
 
177
            ni = ioff[BasisSet.shells[si].am];
 
178
            nj = ioff[BasisSet.shells[sj].am];
 
179
            nk = ioff[BasisSet.shells[sk].am];
 
180
            nl = ioff[BasisSet.shells[sl].am];
 
181
            quartet_size = ni*nj*nk*nl;
 
182
 
 
183
            np_i = BasisSet.shells[si].n_prims;
 
184
            np_j = BasisSet.shells[sj].n_prims;
 
185
            np_k = BasisSet.shells[sk].n_prims;
 
186
            np_l = BasisSet.shells[sl].n_prims;
 
187
            
 
188
            orig_am[0] = BasisSet.shells[si].am-1;
 
189
            orig_am[1] = BasisSet.shells[sj].am-1;
 
190
            orig_am[2] = BasisSet.shells[sk].am-1;
 
191
            orig_am[3] = BasisSet.shells[sl].am-1;
 
192
            am = orig_am[0] + orig_am[1] + orig_am[2] + orig_am[3];
 
193
 
 
194
            sp_ij = &(BasisSet.shell_pairs[si][sj]);
 
195
            sp_kl = &(BasisSet.shell_pairs[sk][sl]);
 
196
 
 
197
            Libderiv.AB[0] = sp_ij->AB[0];
 
198
            Libderiv.AB[1] = sp_ij->AB[1];
 
199
            Libderiv.AB[2] = sp_ij->AB[2];
 
200
            Libderiv.CD[0] = sp_kl->AB[0];
 
201
            Libderiv.CD[1] = sp_kl->AB[1];
 
202
            Libderiv.CD[2] = sp_kl->AB[2];
 
203
            
 
204
            AB2 = Libderiv.AB[0]*Libderiv.AB[0]+
 
205
                  Libderiv.AB[1]*Libderiv.AB[1]+
 
206
                  Libderiv.AB[2]*Libderiv.AB[2];
 
207
            CD2 = Libderiv.CD[0]*Libderiv.CD[0]+
 
208
                  Libderiv.CD[1]*Libderiv.CD[1]+
 
209
                  Libderiv.CD[2]*Libderiv.CD[2];
 
210
 
 
211
            /*-------------------------
 
212
              Figure out the prefactor
 
213
             -------------------------*/
 
214
            pfac = 1.0;
 
215
            if (si == sj)
 
216
              pfac *= 0.5;
 
217
            if (sk == sl)
 
218
              pfac *= 0.5;
 
219
            if (si == sk && sj == sl || si == sl && sj == sk)
 
220
              pfac *= 0.5;
 
221
            pfac *= 8.0;   /*--- The factor of 8 needed for correlated densities ---*/
 
222
 
 
223
            /*--- Compute data for primitive quartets here ---*/
 
224
            num_prim_comb = 0;
 
225
            for (pi = 0; pi < np_i; pi++) {
 
226
              max_pj = (si == sj) ? pi+1 : np_j;
 
227
              for (pj = 0; pj < max_pj; pj++) {
 
228
                m = (1 + (si == sj && pi != pj));
 
229
                for (pk = 0; pk < np_k; pk++) {
 
230
                  max_pl = (sk == sl) ? pk+1 : np_l;
 
231
                  for (pl = 0; pl < max_pl; pl++){
 
232
                    n = m * (1 + (sk == sl && pk != pl));
 
233
#ifdef USE_TAYLOR_FM
 
234
                    deriv1_quartet_data(&(Libderiv.PrimQuartet[num_prim_comb++]),
 
235
                                        NULL, AB2, CD2,
 
236
                                        sp_ij, sp_kl, am, pi, pj, pk, pl, n*pfac);
 
237
#else
 
238
                    deriv1_quartet_data(&(Libderiv.PrimQuartet[num_prim_comb++]),
 
239
                                        &fjt_table, AB2, CD2,
 
240
                                        sp_ij, sp_kl, am, pi, pj, pk, pl, n*pfac);
 
241
#endif    
 
242
                    
 
243
                  }
 
244
                }
 
245
              }
 
246
            }
 
247
 
 
248
            /*--- Read in a shell quartet from disk ---*/
 
249
              memset(FourInd,0,sizeof(double)*quartet_size);
 
250
              last_buf = TPDM.lastbuf;
 
251
              quartet_done = 0;
 
252
              sii_fao = BasisSet.shells[sii].fao-1;
 
253
              sjj_fao = BasisSet.shells[sjj].fao-1;
 
254
              skk_fao = BasisSet.shells[skk].fao-1;
 
255
              sll_fao = BasisSet.shells[sll].fao-1;
 
256
              do {
 
257
                if (buf_offset < buf_size) {
 
258
                  i = TPDM.labels[buf_4offset]   - sii_fao;
 
259
                  if (i >= 0) {
 
260
                    j = TPDM.labels[buf_4offset+1] - sjj_fao;
 
261
                    k = TPDM.labels[buf_4offset+2] - skk_fao;
 
262
                    l = TPDM.labels[buf_4offset+3] - sll_fao;
 
263
                    if (switch_ij) {
 
264
                      dum = i;
 
265
                      i = j;
 
266
                      j = dum;
 
267
                    }
 
268
                    if (switch_kl) {
 
269
                      dum = k;
 
270
                      k = l;
 
271
                      l = dum;
 
272
                    }
 
273
                    if (switch_ijkl) {
 
274
                      dum = i;
 
275
                      i = k;
 
276
                      k = dum;
 
277
                      dum = j;
 
278
                      j = l;
 
279
                      l = dum;
 
280
                    }
 
281
                    offset = ((i*nj+j)*nk+k)*nl+l;
 
282
                    FourInd[offset] += TPDM.values[buf_offset]*
 
283
                                      GTOs.bf_norm[orig_am[0]][i]*
 
284
                                      GTOs.bf_norm[orig_am[1]][j]*
 
285
                                      GTOs.bf_norm[orig_am[2]][k]*
 
286
                                      GTOs.bf_norm[orig_am[3]][l];  
 
287
                  }
 
288
                  else
 
289
                    quartet_done = 1;
 
290
                  buf_offset++;
 
291
                  buf_4offset += 4;
 
292
                }
 
293
                else if (!last_buf) {
 
294
                  iwl_buf_fetch(&TPDM);
 
295
                  buf_offset = 0;
 
296
                  buf_4offset = 0;
 
297
                  last_buf = TPDM.lastbuf;
 
298
                }
 
299
                else {
 
300
                  punt("The last TPDM quartet not marked");
 
301
                }
 
302
              } while (!quartet_done);
 
303
              
 
304
            build_deriv1_eri[orig_am[0]][orig_am[1]][orig_am[2]][orig_am[3]](&Libderiv,num_prim_comb);
 
305
 
 
306
            center_i = BasisSet.shells[si].center-1;
 
307
            center_j = BasisSet.shells[sj].center-1;
 
308
            center_k = BasisSet.shells[sk].center-1;
 
309
            center_l = BasisSet.shells[sl].center-1;
 
310
            
 
311
            ddax = 0.0;
 
312
            for(k=0;k<quartet_size;k++)
 
313
              ddax += Libderiv.ABCD[0][k]*FourInd[k];
 
314
            grad_te_local[center_i][0] += ddax;
 
315
 
 
316
            dday = 0.0;
 
317
            for(k=0;k<quartet_size;k++)
 
318
              dday += Libderiv.ABCD[1][k]*FourInd[k];
 
319
            grad_te_local[center_i][1] += dday;
 
320
 
 
321
            ddaz = 0.0;
 
322
            for(k=0;k<quartet_size;k++)
 
323
              ddaz += Libderiv.ABCD[2][k]*FourInd[k];
 
324
            grad_te_local[center_i][2] += ddaz;
 
325
 
 
326
            /*ddbx = 0.0;
 
327
            for(k=0;k<quartet_size;k++)
 
328
              ddbx += Libderiv.ABCD[3][k]*FourInd[k];
 
329
            grad_te_local[center_j][0] += ddbx;
 
330
 
 
331
            ddby = 0.0;
 
332
            for(k=0;k<quartet_size;k++)
 
333
              ddby += Libderiv.ABCD[4][k]*FourInd[k];
 
334
            grad_te_local[center_j][1] += ddby;
 
335
 
 
336
            ddbz = 0.0;
 
337
            for(k=0;k<quartet_size;k++)
 
338
              ddbz += Libderiv.ABCD[5][k]*FourInd[k];
 
339
              grad_te_local[center_j][2] += ddbz;*/
 
340
 
 
341
            ddcx = 0.0;
 
342
            for(k=0;k<quartet_size;k++)
 
343
              ddcx += Libderiv.ABCD[6][k]*FourInd[k];
 
344
            grad_te_local[center_k][0] += ddcx;
 
345
 
 
346
            ddcy = 0.0;
 
347
            for(k=0;k<quartet_size;k++)
 
348
              ddcy += Libderiv.ABCD[7][k]*FourInd[k];
 
349
            grad_te_local[center_k][1] += ddcy;
 
350
 
 
351
            ddcz = 0.0;
 
352
            for(k=0;k<quartet_size;k++)
 
353
              ddcz += Libderiv.ABCD[8][k]*FourInd[k];
 
354
            grad_te_local[center_k][2] += ddcz;
 
355
 
 
356
            dddx = 0.0;
 
357
            for(k=0;k<quartet_size;k++)
 
358
              dddx += Libderiv.ABCD[9][k]*FourInd[k];
 
359
            grad_te_local[center_l][0] += dddx;
 
360
 
 
361
            dddy = 0.0;
 
362
            for(k=0;k<quartet_size;k++)
 
363
              dddy += Libderiv.ABCD[10][k]*FourInd[k];
 
364
            grad_te_local[center_l][1] += dddy;
 
365
 
 
366
            dddz = 0.0;
 
367
            for(k=0;k<quartet_size;k++)
 
368
              dddz += Libderiv.ABCD[11][k]*FourInd[k];
 
369
            grad_te_local[center_l][2] += dddz;
 
370
 
 
371
            grad_te_local[center_j][0] -= ddax + ddcx + dddx;
 
372
            grad_te_local[center_j][1] -= dday + ddcy + dddy;
 
373
            grad_te_local[center_j][2] -= ddaz + ddcz + dddz;
 
374
        }
 
375
 
 
376
  if (UserOptions.print_lvl >= PRINT_TEDERIV)
 
377
    print_atomvec("Two-electron contribution to the forces (a.u.)",grad_te_local);
 
378
 
 
379
  for(i=0;i<Molecule.num_atoms;i++) {
 
380
    Grad[i][0] += grad_te_local[i][0];
 
381
    Grad[i][1] += grad_te_local[i][1];
 
382
    Grad[i][2] += grad_te_local[i][2];
 
383
  }
 
384
  
 
385
  /*---------
 
386
    Clean-up
 
387
   ---------*/
 
388
  free(FourInd);
 
389
  free_block(grad_te_local);
 
390
  free_libderiv(&Libderiv);
 
391
#ifdef USE_TAYLOR_FM
 
392
  free_Taylor_Fm_Eval();
 
393
#else
 
394
  free_fjt_table(&fjt_table);
 
395
  free_fjt();
 
396
#endif
 
397
  iwl_buf_close(&TPDM,1);
 
398
 
 
399
  return;
 
400
}
 
401
 
 
402
 
 
403