← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/cints/Fock/hf_fock_thread.cc

  • 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/Fock/hf_fock_thread.cc
 
1
/*! \file
 
2
    \ingroup CINTS
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include<cmath>
 
6
#include<cstdio>
 
7
#include<cstring>
 
8
#include<cstdlib>
 
9
#include<memory.h>
 
10
#include<pthread.h>
 
11
#include<psitypes.h>
 
12
#include<libipv1/ip_lib.h>
 
13
#include<libciomr/libciomr.h>
 
14
#include<libpsio/psio.h>
 
15
#include<libint/libint.h>
 
16
 
 
17
#include"defines.h"
 
18
#define EXTERN
 
19
#include"global.h"
 
20
#include <stdexcept>
 
21
 
 
22
#include"quartet_data.h"      /* From Default_Ints */
 
23
#include"norm_quartet.h"
 
24
#include"hash.h"
 
25
#include"transmat.h"
 
26
#include"read_scf_opdm.h"
 
27
#ifdef USE_TAYLOR_FM
 
28
  #include"taylor_fm_eval.h"
 
29
#else
 
30
  #include"int_fjt.h"
 
31
  #include"fjt.h"
 
32
#endif
 
33
#include"schwartz.h"
 
34
#include"shell_block_matrix.h"
 
35
 
 
36
using namespace std;
 
37
      
 
38
extern pthread_mutex_t fock_mutex;    /* Used to lock Fock matrices during update */
 
39
 
 
40
 
 
41
namespace psi {
 
42
  namespace CINTS {
 
43
 
 
44
    extern double ****Gskel, ****Gskel_o;   /* Global skeleton Fock matrices, updated in critical sections */
 
45
void *hf_fock_thread(void *tnum_ptr)
 
46
{
 
47
  const long int thread_num = (long int) tnum_ptr;
 
48
  /*--- Various data structures ---*/
 
49
  struct tebuf *tot_data;             /* buffer for non-zero integrals */
 
50
  struct shell_pair *sp_ij, *sp_kl, *sp_ik, *sp_il, *sp_jk, *sp_jl;
 
51
  Libint_t Libint;                    /* Integrals library object */
 
52
  htable_t htable;                    /* hashing table */
 
53
#ifndef USE_TAYLOR_FM
 
54
  double_array_t fjt_table;           /* table of auxiliary function F_m(u) for each primitive combination */
 
55
#endif
 
56
  
 
57
  int total_te_count = 0;
 
58
  int ij, kl, ik, jl, ijkl;
 
59
  int ioffset, joffset, koffset, loffset;
 
60
  int count ;
 
61
  int dum;
 
62
  int n, num;
 
63
  int total_am, am;
 
64
  int orig_am[4];
 
65
  int pkblock_end_index = -1;
 
66
  int g, i, j, k, l, m, ii, jj, kk, ll;
 
67
  int a, b, c, d;
 
68
  int si;                                  /* GCC compiler screwes up if static is left out */
 
69
  int sj, sk, sl, si_g, sj_g;
 
70
  int sii, sjj, skk, sll , slll;
 
71
  int sij, skl, sijkl;
 
72
  int pi, pj, pk, pl ;
 
73
  int max_pj, max_pl;
 
74
  int pii, pjj, pkk, pll;
 
75
  int upk, num_unique_pk;
 
76
  int usi_arr[3], usj_arr[3], usk_arr[3], usl_arr[3];
 
77
  int *si_arr, *sj_arr, *sk_arr, *sl_arr, *key_arr;
 
78
  int usii,usjj,uskk,usll,usi,usj,usk,usl;
 
79
  int usi_eq_usj, usi_eq_usk, usi_eq_usl, usj_eq_usl, usk_eq_usj, usk_eq_usl;
 
80
  int usij_eq_uskl, usik_eq_usjl, usil_eq_uskj;
 
81
  int stab_i,stab_j,stab_k,stab_l,stab_ij,stab_kl;
 
82
  int *R_list, *S_list, *T_list;
 
83
  int R,S,T;
 
84
  int dcr_ij, dcr_kl, dcr_ijkl;
 
85
  int num_unique_quartets;
 
86
  int plquartet;
 
87
  int max_num_unique_quartets;
 
88
  int max_num_prim_comb;
 
89
 
 
90
  int class_size;
 
91
  int max_cart_class_size;
 
92
 
 
93
  int bf_i, bf_j, bf_k, bf_l, so_i, so_j, so_k, so_l, s;
 
94
  int np_i, np_j, np_k, np_l;
 
95
  int ni, nj, nk, nl, li, lj;
 
96
 
 
97
  int index;
 
98
  int iimax, jjmax, kkmax, llmax;
 
99
  int irrep, npi_ij, npi_kl, npi_ik, npi_jl, ind_offset;
 
100
 
 
101
  int num_prim_comb, p;
 
102
  PSI_INT_LEAST64 key, key1, key2, key3;
 
103
  int new_quartet, htable_ptr, nstri;
 
104
  PSI_INT_LEAST64 quartet_index;
 
105
 
 
106
  double so_int;
 
107
  double lambda_T = 0.5/Symmetry.nirreps;
 
108
  double AB2, CD2;
 
109
  double *data;
 
110
  double pkblock_end_value = 0.0;
 
111
  double temp;
 
112
  double **tmpmat1;
 
113
  double *qijkl_arr, *qikjl_arr, *qiljk_arr;
 
114
  double q4ijkl;
 
115
  double fac1, fac2, fac3;
 
116
  double ffac1, ffac2, ffac3;
 
117
  double c1, c2, c3, c4;
 
118
  double dmax;
 
119
  double ****G, ****G_o;       /* Shell-blocked skeleton G matrices from this thread */
 
120
 
 
121
  /*---
 
122
    init hashing table to store and retrieve quartet data
 
123
    init table for Fj(T)
 
124
    ---*/
 
125
  if (Symmetry.nirreps > 1)
 
126
    init_htable( &htable, Symmetry.max_stab_index );
 
127
 
 
128
#ifndef USE_TAYLOR_FM
 
129
  init_fjt_table(&fjt_table);
 
130
#endif
 
131
 
 
132
  max_cart_class_size = (ioff[BasisSet.max_am])*
 
133
                        (ioff[BasisSet.max_am])*
 
134
                        (ioff[BasisSet.max_am])*
 
135
                        (ioff[BasisSet.max_am]);
 
136
  max_num_unique_quartets = Symmetry.max_stab_index*
 
137
                            Symmetry.max_stab_index*
 
138
                            Symmetry.max_stab_index;
 
139
  max_num_prim_comb = (BasisSet.max_num_prims*BasisSet.max_num_prims)*
 
140
                      (BasisSet.max_num_prims*BasisSet.max_num_prims);
 
141
  /*--- init a LIBINT object ---*/
 
142
  pthread_mutex_lock(&fock_mutex);
 
143
  UserOptions.memory -= init_libint(&Libint,BasisSet.max_am-1,max_num_prim_comb);
 
144
  pthread_mutex_unlock(&fock_mutex);
 
145
 
 
146
  /*--- Allocate this thread's shell-blocked skeleton G's ---*/
 
147
  G = init_shell_block_matrix();
 
148
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf)
 
149
    G_o = init_shell_block_matrix();
 
150
  
 
151
  tot_data = (struct tebuf*) malloc(max_cart_class_size*sizeof(struct tebuf));
 
152
  if (Symmetry.nirreps == 1) {
 
153
    si_arr = (int *)malloc(sizeof(int)*3*max_num_unique_quartets);
 
154
    sj_arr = (int *)malloc(sizeof(int)*3*max_num_unique_quartets);
 
155
    sk_arr = (int *)malloc(sizeof(int)*3*max_num_unique_quartets);
 
156
    sl_arr = (int *)malloc(sizeof(int)*3*max_num_unique_quartets);
 
157
  }
 
158
  key_arr = (int *)malloc(sizeof(int)*3*max_num_unique_quartets);
 
159
 
 
160
  c1 = 1.0 - 0.5*UserOptions.hf_exch;
 
161
  c2 = 1.0 - 0.25*UserOptions.hf_exch;
 
162
  c3 = -0.5*UserOptions.hf_exch;
 
163
  c4 = -0.25*UserOptions.hf_exch;
 
164
  
 
165
/*-------------------------------------------------
 
166
  generate all unique shell quartets with ordering
 
167
  suitable for building the PK-matrix
 
168
 -------------------------------------------------*/
 
169
  quartet_index = 0;
 
170
  for (usii=0; usii<Symmetry.num_unique_shells; usii++)
 
171
    for (usjj=0; usjj<=usii; usjj++)
 
172
      for (uskk=0; uskk<=usjj; uskk++)
 
173
        for (usll=0; usll<=uskk; usll++, quartet_index++){
 
174
 
 
175
          /*--- Decide if this thread will do this ---*/
 
176
          if ( quartet_index%UserOptions.num_threads != thread_num )
 
177
            continue;
 
178
            
 
179
  /*--- Decide what unique shell quartets out of (ij|kl), (ik|jl), and (il|jk) are unique
 
180
        with respect to permutations ---*/
 
181
          usi_arr[0] = usii; usj_arr[0] = usjj; usk_arr[0] = uskk; usl_arr[0] = usll;
 
182
          if (usii == usjj && usii == uskk || usjj == uskk && usjj == usll)
 
183
            num_unique_pk = 1;
 
184
          else if (usii == uskk || usjj == usll) {
 
185
            num_unique_pk = 2;
 
186
            usi_arr[1] = usii; usj_arr[1] = uskk; usk_arr[1] = usjj; usl_arr[1] = usll;
 
187
          }
 
188
          else if (usjj == uskk) {
 
189
            num_unique_pk = 2;
 
190
            usi_arr[1] = usii; usj_arr[1] = usll; usk_arr[1] = usjj; usl_arr[1] = uskk;
 
191
          }
 
192
          else if (usii == usjj || uskk == usll) {
 
193
            num_unique_pk = 2;
 
194
            usi_arr[1] = usii; usj_arr[1] = uskk; usk_arr[1] = usjj; usl_arr[1] = usll;
 
195
          }
 
196
          else {
 
197
            num_unique_pk = 3;
 
198
            usi_arr[1] = usii; usj_arr[1] = uskk; usk_arr[1] = usjj; usl_arr[1] = usll;
 
199
            usi_arr[2] = usii; usj_arr[2] = usll; usk_arr[2] = usjj; usl_arr[2] = uskk;
 
200
          }
 
201
 
 
202
          count = 0;
 
203
          for(upk=0;upk<num_unique_pk;upk++) {
 
204
            /*--- For each combination of unique shells generate "petit list" of shells ---*/
 
205
            usi = usi_arr[upk]; usj = usj_arr[upk]; usk = usk_arr[upk]; usl = usl_arr[upk];
 
206
 
 
207
            si = Symmetry.us2s[usi];
 
208
            sjj = Symmetry.us2s[usj];
 
209
            skk = Symmetry.us2s[usk];
 
210
            sll = Symmetry.us2s[usl];
 
211
            total_am = BasisSet.shells[si].am+BasisSet.shells[sjj].am+BasisSet.shells[skk].am+BasisSet.shells[sll].am;
 
212
 
 
213
            if (Symmetry.nirreps > 1) { /*--- Non-C1 symmetry case ---*/
 
214
              /*--- Generate the petite list of shell quadruplets using DCD approach of Davidson ---*/
 
215
              stab_i = Symmetry.atom_positions[BasisSet.shells[si].center-1];
 
216
              stab_j = Symmetry.atom_positions[BasisSet.shells[sjj].center-1];
 
217
              stab_k = Symmetry.atom_positions[BasisSet.shells[skk].center-1];
 
218
              stab_l = Symmetry.atom_positions[BasisSet.shells[sll].center-1];
 
219
              stab_ij = Symmetry.GnG[stab_i][stab_j];
 
220
              stab_kl = Symmetry.GnG[stab_k][stab_l];
 
221
              R_list = Symmetry.dcr[stab_i][stab_j];
 
222
              S_list = Symmetry.dcr[stab_k][stab_l];
 
223
              T_list = Symmetry.dcr[stab_ij][stab_kl];
 
224
              fac1 = Symmetry.nirreps/
 
225
                     Symmetry.dcr_deg[Symmetry.GnG[stab_i][stab_j]][Symmetry.GnG[stab_k][stab_l]];
 
226
              usi_eq_usj = (usi == usj);
 
227
              usi_eq_usk = (usi == usk);
 
228
              usi_eq_usl = (usi == usl);
 
229
              usj_eq_usl = (usj == usl);
 
230
              usk_eq_usj = (usk == usj);
 
231
              usk_eq_usl = (usk == usl);
 
232
              usij_eq_uskl = (INDEX(usi,usj) == INDEX(usk,usl));
 
233
              usik_eq_usjl = (INDEX(usi,usk) == INDEX(usj,usl));
 
234
              usil_eq_uskj = (INDEX(usi,usl) == INDEX(usk,usj));
 
235
              if (!usi_eq_usj)
 
236
                fac1 *= 2.0;
 
237
              if (!usk_eq_usl)
 
238
                fac1 *= 2.0;
 
239
              if (usij_eq_uskl)
 
240
                fac1 *= 0.5;
 
241
                
 
242
 
 
243
              for(dcr_ij=0;dcr_ij<Symmetry.dcr_dim[stab_i][stab_j];dcr_ij++){
 
244
                R = R_list[dcr_ij];
 
245
                sj = BasisSet.shells[sjj].trans_vec[R]-1;
 
246
                for(dcr_ijkl=0;dcr_ijkl<Symmetry.dcr_dim[stab_ij][stab_kl];dcr_ijkl++){
 
247
                  T = T_list[dcr_ijkl];
 
248
                  sk = BasisSet.shells[skk].trans_vec[T]-1;
 
249
                  slll = BasisSet.shells[sll].trans_vec[T]-1;
 
250
                  for(dcr_kl=0;dcr_kl<Symmetry.dcr_dim[stab_k][stab_l];dcr_kl++) {
 
251
                    S = S_list[dcr_kl];
 
252
                    sl = BasisSet.shells[slll].trans_vec[S]-1;
 
253
 
 
254
                    /*----------------------------------------------------------
 
255
                      Obviously redundant or zero cases are be eliminated here!
 
256
                      Redundancies arise in DCD approach when usk == usl etc.
 
257
                     -------------------------------------------------------------*/
 
258
                      
 
259
                    if(!(total_am%2)||
 
260
                       (BasisSet.shells[si].center!=BasisSet.shells[sj].center)||
 
261
                       (BasisSet.shells[sj].center!=BasisSet.shells[sk].center)||
 
262
                       (BasisSet.shells[sk].center!=BasisSet.shells[sl].center)) {
 
263
 
 
264
                      q4ijkl = fac1;
 
265
 
 
266
                      key1 = compute_key(si,sj,sk,sl);
 
267
                      key2 = compute_key(si,sk,sj,sl);
 
268
                      key3 = compute_key(si,sl,sk,sj);
 
269
                      new_quartet = put_entry(&htable,key1,si,sj,sk,sl,q4ijkl,0,0);
 
270
                      if (new_quartet > -1) {
 
271
                        key_arr[count] = new_quartet;
 
272
                        count++;
 
273
                      }
 
274
 
 
275
                      if ( (key1 == key3 && key3 != key2) ||
 
276
                           (key2 == key3 && key3 != key1) ||
 
277
                           (key2 != key3 && key1 != key3 && key2 != key1)) {
 
278
                        new_quartet = put_entry(&htable,key2,si,sk,sj,sl,0,q4ijkl,0);
 
279
                        if (new_quartet > -1) {
 
280
                          key_arr[count] = new_quartet;
 
281
                          count++;
 
282
                        }
 
283
                      }
 
284
 
 
285
                      if ( (key1 == key2 && key3 != key1) ||
 
286
                           (key2 != key3 && key1 != key3 && key1 != key2)) {
 
287
                        new_quartet = put_entry(&htable,key3,si,sl,sk,sj,0,0,q4ijkl);
 
288
                        if (new_quartet > -1) {
 
289
                          key_arr[count] = new_quartet;
 
290
                          count++;
 
291
                        }
 
292
                      }
 
293
                    }
 
294
                  }
 
295
                }
 
296
              } /* petite list is ready to be used */
 
297
            }
 
298
            else {
 
299
              if(!(total_am%2)||
 
300
                 (BasisSet.shells[si].center!=BasisSet.shells[sjj].center)||
 
301
                 (BasisSet.shells[sjj].center!=BasisSet.shells[skk].center)||
 
302
                 (BasisSet.shells[skk].center!=BasisSet.shells[sll].center)) {
 
303
                si_arr[count] = si;
 
304
                sj_arr[count] = sjj;
 
305
                sk_arr[count] = skk;
 
306
                sl_arr[count] = sll;
 
307
                count++;
 
308
              }
 
309
            }
 
310
          }
 
311
          num_unique_quartets = count;
 
312
          if (count > 3*max_num_unique_quartets)
 
313
            throw std::domain_error("Problem with hashing?");
 
314
 
 
315
            /*----------------------------------
 
316
              Compute the nonredundant quartets
 
317
             ----------------------------------*/
 
318
            for(plquartet=0;plquartet<num_unique_quartets;plquartet++) {
 
319
              if (Symmetry.nirreps == 1) {
 
320
                si = si_arr[plquartet];
 
321
                sj = sj_arr[plquartet];
 
322
                sk = sk_arr[plquartet];
 
323
                sl = sl_arr[plquartet];
 
324
                dmax = MAX(BasisSet.shell_pairs[si][sj].Dmax, BasisSet.shell_pairs[sk][sl].Dmax);
 
325
                dmax = MAX(dmax, 0.25*BasisSet.shell_pairs[si][sk].Dmax);
 
326
                dmax = MAX(dmax, 0.25*BasisSet.shell_pairs[si][sl].Dmax);
 
327
                dmax = MAX(dmax, 0.25*BasisSet.shell_pairs[sj][sk].Dmax);
 
328
                dmax = MAX(dmax, 0.25*BasisSet.shell_pairs[sj][sl].Dmax);
 
329
                if (BasisSet.schwartz_eri[si][sj]*BasisSet.schwartz_eri[sk][sl]*dmax < UserOptions.cutoff)
 
330
                  continue;
 
331
                
 
332
                fac1 = fac2 = fac3 = 1.0;
 
333
                if (si != sj)
 
334
                  fac1 *= 2.0;
 
335
                if (sk != sl)
 
336
                  fac1 *= 2.0;
 
337
                if (si != sk)
 
338
                  fac2 *= 2.0;
 
339
                if (sj != sl)
 
340
                  fac2 *= 2.0;
 
341
                if (si != sl)
 
342
                  fac3 *= 2.0;
 
343
                if (sj != sk)
 
344
                  fac3 *= 2.0;
 
345
                if (INDEX(si,sj) == INDEX(sk,sl))
 
346
                  fac1 *= 0.5;
 
347
                if (INDEX(si,sk) == INDEX(sj,sl))
 
348
                  fac2 *= 0.5;
 
349
                if (INDEX(si,sl) == INDEX(sj,sk))
 
350
                  fac3 *= 0.5;
 
351
              }
 
352
              else {
 
353
                htable_ptr = key_arr[plquartet];
 
354
                if (htable_ptr >= htable.size || htable_ptr < 0)
 
355
                    throw std::domain_error("Problem with hashing?");
 
356
                htable.table[htable_ptr].key = EMPTY_KEY;
 
357
                si = htable.table[htable_ptr].si;
 
358
                sj = htable.table[htable_ptr].sj;
 
359
                sk = htable.table[htable_ptr].sk;
 
360
                sl = htable.table[htable_ptr].sl;
 
361
                dmax = MAX(BasisSet.shell_pairs[si][sj].Dmax, BasisSet.shell_pairs[sk][sl].Dmax);
 
362
                dmax = MAX(dmax, 0.25*BasisSet.shell_pairs[si][sk].Dmax);
 
363
                dmax = MAX(dmax, 0.25*BasisSet.shell_pairs[si][sl].Dmax);
 
364
                dmax = MAX(dmax, 0.25*BasisSet.shell_pairs[sj][sk].Dmax);
 
365
                dmax = MAX(dmax, 0.25*BasisSet.shell_pairs[sj][sl].Dmax);
 
366
                if (BasisSet.schwartz_eri[si][sj]*BasisSet.schwartz_eri[sk][sl]*dmax < UserOptions.cutoff)
 
367
                  continue;
 
368
 
 
369
                fac1 = htable.table[htable_ptr].q4ijkl;
 
370
                fac2 = htable.table[htable_ptr].q4ikjl;
 
371
                fac3 = htable.table[htable_ptr].q4ilkj;
 
372
                if (si == sj && si == sk ||
 
373
                    sj == sk && sj == sl ||
 
374
                    si == sj && si == sl ||
 
375
                    si == sk && si == sl) {
 
376
                  fac2 = fac3 = fac1;
 
377
                }
 
378
                else if (si == sk || sj == sl) {
 
379
                  fac1 = fac3 = (fac1 + fac3);
 
380
                }
 
381
                else if (sj == sk || si == sl) {
 
382
                  fac1 = fac2 = (fac1 + fac2);
 
383
                }
 
384
                else if (si == sj || sk == sl) {
 
385
                  fac3 = fac2 = (fac3 + fac2);
 
386
                }
 
387
              }
 
388
 
 
389
#if DEBUG
 
390
              if (si < 0 || si >= BasisSet.num_shells)
 
391
                  throw std::domain_error("Problem with shell indices");
 
392
              if (sj < 0 || sj >= BasisSet.num_shells)
 
393
                  throw std::domain_error("Problem with shell indices");
 
394
              if (sk < 0 || sk >= BasisSet.num_shells)
 
395
                  throw std::domain_error("Problem with shell indices");
 
396
              if (sl < 0 || sl >= BasisSet.num_shells)
 
397
                  throw std::domain_error("Problem with shell indices");
 
398
#endif
 
399
              
 
400
              /* place in "ascending" angular mom-
 
401
               my simple way of optimizing PHG recursion (VRR) */
 
402
              /* these first two are good for the HRR */
 
403
              if(BasisSet.shells[si].am < BasisSet.shells[sj].am){
 
404
                dum = si;
 
405
                si = sj;
 
406
                sj = dum;
 
407
                temp = fac2;
 
408
                fac2 = fac3;
 
409
                fac3 = temp;
 
410
              }
 
411
              if(BasisSet.shells[sk].am < BasisSet.shells[sl].am){
 
412
                dum = sk;
 
413
                sk = sl;
 
414
                sl = dum;
 
415
                temp = fac2;
 
416
                fac2 = fac3;
 
417
                fac3 = temp;
 
418
              }
 
419
              /* this should be /good/ for the VRR */
 
420
              if(BasisSet.shells[si].am + BasisSet.shells[sj].am >
 
421
                 BasisSet.shells[sk].am + BasisSet.shells[sl].am){
 
422
                dum = si;
 
423
                si = sk;
 
424
                sk = dum;
 
425
                dum = sj;
 
426
                sj = sl;
 
427
                sl = dum;
 
428
              }
 
429
              ioffset = BasisSet.shells[si].fao - 1;
 
430
              joffset = BasisSet.shells[sj].fao - 1;
 
431
              koffset = BasisSet.shells[sk].fao - 1;
 
432
              loffset = BasisSet.shells[sl].fao - 1;
 
433
              ni = ioff[BasisSet.shells[si].am];
 
434
              nj = ioff[BasisSet.shells[sj].am];
 
435
              nk = ioff[BasisSet.shells[sk].am];
 
436
              nl = ioff[BasisSet.shells[sl].am];
 
437
              np_i = BasisSet.shells[si].n_prims;
 
438
              np_j = BasisSet.shells[sj].n_prims;
 
439
              np_k = BasisSet.shells[sk].n_prims;
 
440
              np_l = BasisSet.shells[sl].n_prims;
 
441
              orig_am[0] = BasisSet.shells[si].am-1;
 
442
              orig_am[1] = BasisSet.shells[sj].am-1;
 
443
              orig_am[2] = BasisSet.shells[sk].am-1;
 
444
              orig_am[3] = BasisSet.shells[sl].am-1;
 
445
              am = orig_am[0] + orig_am[1] + orig_am[2] + orig_am[3];
 
446
            
 
447
              sp_ij = &(BasisSet.shell_pairs[si][sj]);
 
448
              sp_kl = &(BasisSet.shell_pairs[sk][sl]);
 
449
              sp_ik = &(BasisSet.shell_pairs[si][sk]);
 
450
              sp_il = &(BasisSet.shell_pairs[si][sl]);
 
451
              sp_jk = &(BasisSet.shell_pairs[sj][sk]);
 
452
              sp_jl = &(BasisSet.shell_pairs[sj][sl]);
 
453
              Libint.AB[0] = sp_ij->AB[0];
 
454
              Libint.AB[1] = sp_ij->AB[1];
 
455
              Libint.AB[2] = sp_ij->AB[2];
 
456
              Libint.CD[0] = sp_kl->AB[0];
 
457
              Libint.CD[1] = sp_kl->AB[1];
 
458
              Libint.CD[2] = sp_kl->AB[2];
 
459
                  
 
460
              AB2 = Libint.AB[0]*Libint.AB[0]+
 
461
                    Libint.AB[1]*Libint.AB[1]+
 
462
                    Libint.AB[2]*Libint.AB[2];
 
463
              CD2 = Libint.CD[0]*Libint.CD[0]+
 
464
                    Libint.CD[1]*Libint.CD[1]+
 
465
                    Libint.CD[2]*Libint.CD[2];
 
466
 
 
467
              /*--- Compute data for primitive quartets here ---*/
 
468
              num_prim_comb = 0;
 
469
              for (pi = 0; pi < np_i; pi++) {
 
470
                max_pj = (si == sj) ? pi+1 : np_j;
 
471
                for (pj = 0; pj < max_pj; pj++) {
 
472
                  m = (1 + (si == sj && pi != pj));
 
473
                  for (pk = 0; pk < np_k; pk++) {
 
474
                    max_pl = (sk == sl) ? pk+1 : np_l;
 
475
                    for (pl = 0; pl < max_pl; pl++){
 
476
                      n = m * (1 + (sk == sl && pk != pl));
 
477
#ifdef USE_TAYLOR_FM
 
478
                      quartet_data(&(Libint.PrimQuartet[num_prim_comb++]), NULL, AB2, CD2,
 
479
                                   sp_ij, sp_kl, am, pi, pj, pk, pl, n*lambda_T);
 
480
#else
 
481
                      quartet_data(&(Libint.PrimQuartet[num_prim_comb++]), &fjt_table, AB2, CD2,
 
482
                                   sp_ij, sp_kl, am, pi, pj, pk, pl, n*lambda_T);
 
483
#endif
 
484
                    }
 
485
                  }
 
486
                }
 
487
              }
 
488
 
 
489
              /*-----------------------------------------------------
 
490
                Compute the quartet and compute its contributions to
 
491
                appropriate blocks of Gsh and Gsh_o
 
492
               -----------------------------------------------------*/
 
493
              if (am) {
 
494
                data = build_eri[orig_am[0]][orig_am[1]][orig_am[2]][orig_am[3]](&Libint,num_prim_comb);
 
495
                /*                    zero here means no transformation to puream basis */
 
496
                /*                                    |                                 */
 
497
                data = norm_quartet(data, NULL, orig_am, 0);
 
498
 
 
499
                /*--- Here just put non-redundant integrals to tot_data ---*/
 
500
                num = 0;
 
501
                if(si==sj && sk==sl && si==sk) {
 
502
                  iimax = ni - 1;
 
503
                  for(ii=0; ii <= iimax; ii++){
 
504
                    jjmax = ii;
 
505
                    for(jj=0; jj <= jjmax; jj++){
 
506
                      kkmax = ii;
 
507
                      for(kk=0; kk <= kkmax; kk++){
 
508
                        llmax = (kk==ii)? jj : kk ;
 
509
                        for(ll=0; ll <= llmax; ll++){
 
510
                          ijkl = ll+nl*(kk+nk*(jj+nj*ii));
 
511
/*                        if (fabs(data[ijkl])>UserOptions.cutoff) {*/
 
512
                            tot_data[num].i = (short int) ii;
 
513
                            tot_data[num].j = (short int) jj;
 
514
                            tot_data[num].k = (short int) kk;
 
515
                            tot_data[num].l = (short int) ll;
 
516
                            tot_data[num].val = data[ijkl];
 
517
                            num++;
 
518
/*                        }*/
 
519
                        }
 
520
                      }
 
521
                    }
 
522
                  }
 
523
                }
 
524
                else if(si==sk && sj==sl){
 
525
                  iimax = ni - 1;
 
526
                  for(ii=0; ii <= iimax; ii++){
 
527
                    jjmax = nj - 1;
 
528
                    for(jj=0; jj <= jjmax; jj++){
 
529
                      kkmax = ii;
 
530
                      for(kk=0; kk <= kkmax; kk++){
 
531
                        llmax = (kk==ii)? jj : nl - 1;
 
532
                        for(ll=0; ll <= llmax; ll++){
 
533
                          ijkl = ll+nl*(kk+nk*(jj+nj*ii));
 
534
/*                        if(fabs(data[ijkl])>UserOptions.cutoff){*/
 
535
                            tot_data[num].i = (short int) ii;
 
536
                            tot_data[num].j = (short int) jj;
 
537
                            tot_data[num].k = (short int) kk;
 
538
                            tot_data[num].l = (short int) ll;
 
539
                            tot_data[num].val = data[ijkl];
 
540
                            num++;
 
541
/*                        }*/
 
542
                        }
 
543
                      }
 
544
                    }
 
545
                  }
 
546
                }
 
547
                else {
 
548
                  iimax = ni - 1;
 
549
                  kkmax = nk - 1;
 
550
                  for(ii=0; ii <= iimax; ii++){
 
551
                    jjmax = (si == sj) ? ii : nj - 1;
 
552
                    for(jj=0; jj <= jjmax; jj++){
 
553
                      for(kk=0; kk <= kkmax; kk++){
 
554
                        llmax = (sk == sl) ? kk : nl - 1;
 
555
                        for(ll=0; ll <= llmax; ll++){
 
556
                          ijkl = ll+nl*(kk+nk*(jj+nj*ii));
 
557
/*                        if(fabs(data[ijkl])>UserOptions.cutoff){*/
 
558
                            tot_data[num].i = (short int) ii;
 
559
                            tot_data[num].j = (short int) jj;
 
560
                            tot_data[num].k = (short int) kk;
 
561
                            tot_data[num].l = (short int) ll;
 
562
                            tot_data[num].val = data[ijkl];
 
563
                            num++;
 
564
/*                        }*/
 
565
                        }
 
566
                      }
 
567
                    }
 
568
                  }
 
569
                }
 
570
 
 
571
              }
 
572
              else {
 
573
                temp = 0.0;
 
574
                for(p=0;p<num_prim_comb;p++)
 
575
                  temp += Libint.PrimQuartet[p].F[0];
 
576
 
 
577
                tot_data[0].i = tot_data[0].j = tot_data[0].k = tot_data[0].l = 0;
 
578
                tot_data[0].val = temp;
 
579
                num = 0;
 
580
/*              if (fabs(temp) > UserOptions.cutoff)*/
 
581
                  num = 1;
 
582
 
 
583
              }
 
584
 
 
585
              total_te_count += num;
 
586
 
 
587
              if (UserOptions.reftype == rhf) {
 
588
                for(n=0;n<num;n++) {
 
589
                  i = tot_data[n].i;
 
590
                  j = tot_data[n].j;
 
591
                  k = tot_data[n].k;
 
592
                  l = tot_data[n].l;
 
593
                  ii = i + ioffset;
 
594
                  jj = j + joffset;
 
595
                  kk = k + koffset;
 
596
                  ll = l + loffset;
 
597
                  temp = tot_data[n].val;
 
598
                  
 
599
                  ffac1 = fac1*temp;
 
600
                  ffac2 = fac2*temp;
 
601
                  ffac3 = fac3*temp;
 
602
                  if (ii != jj && si == sj)
 
603
                      ffac1 *= 2.0;
 
604
                  if (kk != ll && sk == sl)
 
605
                      ffac1 *= 2.0;
 
606
                  if (ii != kk && si == sk)
 
607
                      ffac2 *= 2.0;
 
608
                  if (jj != ll && sj == sl)
 
609
                      ffac2 *= 2.0;
 
610
                  if (ii != ll && si == sl)
 
611
                      ffac3 *= 2.0;
 
612
                  if (jj != kk && sj == sk)
 
613
                      ffac3 *= 2.0;
 
614
                  if (INDEX(ii,jj) != INDEX(kk,ll) &&
 
615
                      INDEX(si,sj) == INDEX(sk,sl))
 
616
                      ffac1 *= 2.0;
 
617
                  if (INDEX(ii,kk) != INDEX(jj,ll) &&
 
618
                      INDEX(si,sk) == INDEX(sj,sl))
 
619
                      ffac2 *= 2.0;
 
620
                  if (INDEX(ii,ll) != INDEX(jj,kk) &&
 
621
                      INDEX(si,sl) == INDEX(sj,sk))
 
622
                      ffac3 *= 2.0;
 
623
                  
 
624
                  if (ii == jj && ii == kk ||
 
625
                      jj == kk && jj == ll ||
 
626
                      ii == jj && ii == ll ||
 
627
                      ii == kk && ii == ll) {
 
628
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*(c1*ffac1);
 
629
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*(c1*ffac1);
 
630
                  }
 
631
                  else if (ii == kk || jj == ll) {
 
632
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*(c2*ffac1);
 
633
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*(c2*ffac1);
 
634
                      G[si][sk][i][k] += sp_jl->dmat[j][l]*(c3*ffac2);
 
635
                      G[sj][sl][j][l] += sp_ik->dmat[i][k]*(c3*ffac2);
 
636
                  }
 
637
                  else if (jj == kk || ii == ll) {
 
638
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*(c2*ffac1);
 
639
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*(c2*ffac1);
 
640
                      G[si][sl][i][l] += sp_jk->dmat[j][k]*(c3*ffac3);
 
641
                      G[sj][sk][j][k] += sp_il->dmat[i][l]*(c3*ffac3);
 
642
                  }
 
643
                  else if (ii == jj || kk == ll) {
 
644
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*(ffac1);
 
645
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*(ffac1);
 
646
                      G[si][sk][i][k] += sp_jl->dmat[j][l]*(c4*ffac2);
 
647
                      G[sj][sl][j][l] += sp_ik->dmat[i][k]*(c4*ffac2);
 
648
                  }
 
649
                  else {
 
650
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*(ffac1);
 
651
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*(ffac1);
 
652
                      G[si][sk][i][k] += sp_jl->dmat[j][l]*(c4*ffac2);
 
653
                      G[sj][sl][j][l] += sp_ik->dmat[i][k]*(c4*ffac2);
 
654
                      G[si][sl][i][l] += sp_jk->dmat[j][k]*(c4*ffac3);
 
655
                      G[sj][sk][j][k] += sp_il->dmat[i][l]*(c4*ffac3);
 
656
                  }
 
657
                }
 
658
              }
 
659
              else if (UserOptions.reftype == uhf || UserOptions.reftype == rohf) {
 
660
                for(n=0;n<num;n++) {
 
661
                  i = tot_data[n].i;
 
662
                  j = tot_data[n].j;
 
663
                  k = tot_data[n].k;
 
664
                  l = tot_data[n].l;
 
665
                  ii = i + ioffset;
 
666
                  jj = j + joffset;
 
667
                  kk = k + koffset;
 
668
                  ll = l + loffset;
 
669
                  temp = tot_data[n].val;
 
670
                  
 
671
                  ffac1 = fac1*temp;
 
672
                  ffac2 = fac2*temp;
 
673
                  ffac3 = fac3*temp;
 
674
                  if (ii != jj && si == sj)
 
675
                      ffac1 *= 2.0;
 
676
                  if (kk != ll && sk == sl)
 
677
                      ffac1 *= 2.0;
 
678
                  if (ii != kk && si == sk)
 
679
                      ffac2 *= 2.0;
 
680
                  if (jj != ll && sj == sl)
 
681
                      ffac2 *= 2.0;
 
682
                  if (ii != ll && si == sl)
 
683
                      ffac3 *= 2.0;
 
684
                  if (jj != kk && sj == sk)
 
685
                      ffac3 *= 2.0;
 
686
                  if (INDEX(ii,jj) != INDEX(kk,ll) &&
 
687
                      INDEX(si,sj) == INDEX(sk,sl))
 
688
                      ffac1 *= 2.0;
 
689
                  if (INDEX(ii,kk) != INDEX(jj,ll) &&
 
690
                      INDEX(si,sk) == INDEX(sj,sl))
 
691
                      ffac2 *= 2.0;
 
692
                  if (INDEX(ii,ll) != INDEX(jj,kk) &&
 
693
                      INDEX(si,sl) == INDEX(sj,sk))
 
694
                      ffac3 *= 2.0;
 
695
                  
 
696
                  if (ii == jj && ii == kk ||
 
697
                      jj == kk && jj == ll ||
 
698
                      ii == jj && ii == ll ||
 
699
                      ii == kk && ii == ll) {
 
700
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*(c1*ffac1);
 
701
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*(c1*ffac1);
 
702
                      G_o[si][sj][i][j] -= sp_kl->dmato[k][l]*(c3*ffac1);
 
703
                      G_o[sk][sl][k][l] -= sp_ij->dmato[i][j]*(c3*ffac1);
 
704
                  }
 
705
                  else if (ii == kk || jj == ll) {
 
706
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*(c2*ffac1);
 
707
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*(c2*ffac1);
 
708
                      G[si][sk][i][k] += sp_jl->dmat[j][l]*(c3*ffac2);
 
709
                      G[sj][sl][j][l] += sp_ik->dmat[i][k]*(c3*ffac2);
 
710
                      G_o[si][sj][i][j] -= sp_kl->dmato[k][l]*(c4*ffac1);
 
711
                      G_o[sk][sl][k][l] -= sp_ij->dmato[i][j]*(c4*ffac1);
 
712
                      G_o[si][sk][i][k] -= sp_jl->dmato[j][l]*(c3*ffac2);
 
713
                      G_o[sj][sl][j][l] -= sp_ik->dmato[i][k]*(c3*ffac2);
 
714
                  }
 
715
                  else if (jj == kk || ii == ll) {
 
716
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*(c2*ffac1);
 
717
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*(c2*ffac1);
 
718
                      G[si][sl][i][l] += sp_jk->dmat[j][k]*(c3*ffac3);
 
719
                      G[sj][sk][j][k] += sp_il->dmat[i][l]*(c3*ffac3);
 
720
                      G_o[si][sj][i][j] -= sp_kl->dmato[k][l]*(c4*ffac1);
 
721
                      G_o[sk][sl][k][l] -= sp_ij->dmato[i][j]*(c4*ffac1);
 
722
                      G_o[si][sl][i][l] -= sp_jk->dmato[j][k]*(c3*ffac3);
 
723
                      G_o[sj][sk][j][k] -= sp_il->dmato[i][l]*(c3*ffac3);
 
724
                  }
 
725
                  else if (ii == jj || kk == ll) {
 
726
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*ffac1;
 
727
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*ffac1;
 
728
                      G[si][sk][i][k] += sp_jl->dmat[j][l]*(c4*ffac2);
 
729
                      G[sj][sl][j][l] += sp_ik->dmat[i][k]*(c4*ffac2);
 
730
                      G_o[si][sk][i][k] -= sp_jl->dmato[j][l]*(c4*ffac2);
 
731
                      G_o[sj][sl][j][l] -= sp_ik->dmato[i][k]*(c4*ffac2);
 
732
                  }
 
733
                  else {
 
734
                      G[si][sj][i][j] += sp_kl->dmat[k][l]*ffac1;
 
735
                      G[sk][sl][k][l] += sp_ij->dmat[i][j]*ffac1;
 
736
                      G[si][sk][i][k] += sp_jl->dmat[j][l]*(c4*ffac2);
 
737
                      G[sj][sl][j][l] += sp_ik->dmat[i][k]*(c4*ffac2);
 
738
                      G[si][sl][i][l] += sp_jk->dmat[j][k]*(c4*ffac3);
 
739
                      G[sj][sk][j][k] += sp_il->dmat[i][l]*(c4*ffac3);
 
740
                      G_o[si][sk][i][k] -= sp_jl->dmato[j][l]*(c4*ffac2);
 
741
                      G_o[sj][sl][j][l] -= sp_ik->dmato[i][k]*(c4*ffac2);
 
742
                      G_o[si][sl][i][l] -= sp_jk->dmato[j][k]*(c4*ffac3);
 
743
                      G_o[sj][sk][j][k] -= sp_il->dmato[i][l]*(c4*ffac3);
 
744
                  }
 
745
                }
 
746
              }
 
747
              
 
748
            } /* end of computing "petit" list */
 
749
        } /* end getting unique shell combination */
 
750
 
 
751
  pthread_mutex_lock(&fock_mutex);
 
752
  {
 
753
      for(si=0;si<BasisSet.num_shells;si++) {
 
754
        ni = ioff[BasisSet.shells[si].am];
 
755
        for(sj=0;sj<BasisSet.num_shells;sj++) {
 
756
          nj = ioff[BasisSet.shells[sj].am];
 
757
          for(i=0;i<ni;i++)
 
758
            for(j=0;j<nj;j++)
 
759
              Gskel[si][sj][i][j] += G[si][sj][i][j];
 
760
        }
 
761
      }
 
762
      if (UserOptions.reftype == rohf || UserOptions.reftype == uhf)
 
763
        for(si=0;si<BasisSet.num_shells;si++) {
 
764
          ni = ioff[BasisSet.shells[si].am];
 
765
          for(sj=0;sj<BasisSet.num_shells;sj++) {
 
766
            nj = ioff[BasisSet.shells[sj].am];
 
767
            for(i=0;i<ni;i++)
 
768
              for(j=0;j<nj;j++)
 
769
                Gskel_o[si][sj][i][j] += G_o[si][sj][i][j];
 
770
          }
 
771
        }
 
772
  }
 
773
   pthread_mutex_unlock(&fock_mutex);
 
774
  
 
775
 
 
776
  /*---------
 
777
    Clean-up
 
778
   ---------*/
 
779
  free_shell_block_matrix(G);
 
780
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
781
    free_shell_block_matrix(G_o);
 
782
  }
 
783
  if (Symmetry.nirreps > 1) {
 
784
    free_htable(&htable);
 
785
  }
 
786
  free_libint(&Libint);
 
787
  free(tot_data);
 
788
  if (Symmetry.nirreps == 1) {
 
789
    free(si_arr);
 
790
    free(sj_arr);
 
791
    free(sk_arr);
 
792
    free(sl_arr);
 
793
  }
 
794
  free(key_arr);
 
795
 
 
796
#ifndef USE_TAYLOR_FM
 
797
  free_fjt_table(&fjt_table);
 
798
#endif
 
799
 
 
800
  return NULL;
 
801
}
 
802
 
 
803
};};