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