← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/R12_Ints/r12_te_ints.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/R12_Ints/r12_te_ints.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<memory.h>
 
9
#include<cstdlib>
 
10
#include<libiwl/iwl.h>
 
11
#include<libciomr/libciomr.h>
 
12
#include<libint/libint.h>
 
13
#include<libr12/libr12.h>
 
14
 
 
15
#include"defines.h"
 
16
#define EXTERN
 
17
#include"global.h"
 
18
#include <stdexcept>
 
19
#include"r12_quartet_data.h"
 
20
#include"iwl_tebuf.h"
 
21
#include"norm_quartet.h"
 
22
#include"int_fjt.h"
 
23
 
 
24
 
 
25
namespace psi { namespace CINTS {
 
26
void r12_te_ints()
 
27
{
 
28
  const double toler = UserOptions.cutoff;
 
29
  
 
30
  /*--- ASCII file to print integrals ---*/
 
31
  FILE *teout[NUM_TE_TYPES];
 
32
  char teout_filename[] = "teout0.dat";
 
33
 
 
34
  /*--- Various data structures ---*/
 
35
  struct iwlbuf TEOUT[NUM_TE_TYPES-1];             /* IWL buffer for target integrals */
 
36
  int itapTE[NUM_TE_TYPES-1];
 
37
  struct coordinates ericent[4];       /* coordinates of centers A, B, C, and D */
 
38
  struct tebuf *tot_data[NUM_TE_TYPES];              /* accum. for contracted integrals */
 
39
  struct shell_pair *sp_ij, *sp_kl;
 
40
  struct unique_shell_pair *usp_ij,*usp_kl;
 
41
  Libr12_t Libr12;
 
42
  double_array_t fjt_table;
 
43
 
 
44
  static char *te_operator[] = { "1/r12", "r12", "[r12,T1]" };
 
45
  int total_te_count[NUM_TE_TYPES] = {0, 0, 0, 0};
 
46
  int ij, kl, ik, jl, ijkl;
 
47
  int ioffset, joffset, koffset, loffset;
 
48
  int count ;
 
49
  int dum;
 
50
  int te_type;
 
51
  int n, num[NUM_TE_TYPES];
 
52
  int total_am, am;
 
53
  int orig_am[4];
 
54
  int pkblock_end_index = -1;
 
55
  register int i, j, k, l, ii, jj, kk, ll;
 
56
  register int si, sj, sk, sl ;
 
57
  register int sii, sjj, skk, sll , slll;
 
58
  register int pi, pj, pk, pl ;
 
59
  register int pii, pjj, pkk, pll ;
 
60
  int upk, num_unique_pk;
 
61
  int usi_arr[3], usj_arr[3], usk_arr[3], usl_arr[3];
 
62
  int *sj_arr, *sk_arr, *sl_arr;
 
63
  int *sj_fbf_arr, *sk_fbf_arr, *sl_fbf_arr;
 
64
  int usii,usjj,uskk,usll,usi,usj,usk,usl;
 
65
  int stab_i,stab_j,stab_k,stab_l,stab_ij,stab_kl;
 
66
  int *R_list, *S_list, *T_list;
 
67
  int R,S,T;
 
68
  int dcr_ij, dcr_kl, dcr_ijkl;
 
69
  int lambda_T = 1;
 
70
  int num_unique_quartets;
 
71
  int plquartet;
 
72
  int max_num_unique_quartets;
 
73
  int max_num_prim_comb;
 
74
 
 
75
  int class_size;
 
76
  int max_class_size;
 
77
  int max_cart_class_size;
 
78
 
 
79
  int bf_i, bf_j, bf_k, bf_l, so_i, so_j, so_k, so_l, s;
 
80
  int so_ii, so_jj, so_kk, so_ll;
 
81
  int np_i, np_j, np_k, np_l;
 
82
  int ni, nj, nk, nl;
 
83
 
 
84
  int index;
 
85
  int iimax, jjmax, kkmax, llmax;
 
86
  int irrep, npi_ij, npi_kl, npi_ik, npi_jl, ind_offset;
 
87
 
 
88
  int num_prim_comb, p;
 
89
 
 
90
  double so_int[NUM_TE_TYPES];
 
91
  double AB2, CD2;
 
92
  double *data[NUM_TE_TYPES];
 
93
  double *puream_data[NUM_TE_TYPES];
 
94
  double **plist_data[NUM_TE_TYPES];
 
95
  double pkblock_end_value = 0.0;
 
96
  double value;
 
97
  double temp;
 
98
  double ssss, ss_r12_ss;
 
99
 
 
100
  /*---------------
 
101
    Initialization
 
102
   ---------------*/
 
103
#if PRINT
 
104
  for(te_type=0;te_type<NUM_TE_TYPES;te_type++) {
 
105
    teout_filename[5] = te_type + '0';
 
106
    teout[te_type] = fopen(teout_filename,"w");
 
107
  }
 
108
#endif
 
109
  itapTE[0] = IOUnits.itap33;
 
110
  itapTE[1] = IOUnits.itapR12;
 
111
  itapTE[2] = IOUnits.itapT1;
 
112
  for(te_type=0;te_type<NUM_TE_TYPES-1;te_type++) {
 
113
    iwl_buf_init(&TEOUT[te_type], itapTE[te_type], toler, 0, 0);
 
114
  }
 
115
  init_fjt(BasisSet.max_am*4);
 
116
  init_libr12_base();
 
117
 
 
118
  
 
119
  /*-------------------------
 
120
    Allocate data structures
 
121
   -------------------------*/
 
122
  max_cart_class_size = ioff[BasisSet.max_am]*
 
123
                        ioff[BasisSet.max_am]*
 
124
                        ioff[BasisSet.max_am]*
 
125
                        ioff[BasisSet.max_am];
 
126
  max_num_unique_quartets = Symmetry.max_stab_index*Symmetry.max_stab_index*Symmetry.max_stab_index;
 
127
  for(te_type=0;te_type<NUM_TE_TYPES;te_type++) {
 
128
    tot_data[te_type] = (struct tebuf*) malloc(max_num_unique_quartets*max_cart_class_size*sizeof(struct tebuf));
 
129
    memset(tot_data[te_type], 0, (max_num_unique_quartets*max_cart_class_size)*sizeof(struct tebuf));
 
130
  }
 
131
  if (BasisSet.puream) {
 
132
    if (Symmetry.nirreps == 1) { /*--- allocate NUM_TE_TYPES scratch arrays for cart->puream transformation ---*/
 
133
      for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
134
        puream_data[te_type] = (double *) malloc(sizeof(double)*
 
135
                                                 (BasisSet.max_am*2-1)*
 
136
                                                 ioff[BasisSet.max_am]*
 
137
                                                 ioff[BasisSet.max_am]*
 
138
                                                 ioff[BasisSet.max_am]);
 
139
    }
 
140
    else {
 
141
      puream_data[0] = (double *) malloc(sizeof(double)*
 
142
                                         (BasisSet.max_am*2-1)*
 
143
                                         ioff[BasisSet.max_am]*
 
144
                                         ioff[BasisSet.max_am]*
 
145
                                         ioff[BasisSet.max_am]);
 
146
      for(te_type=1;te_type<NUM_TE_TYPES;te_type++)
 
147
        puream_data[te_type] = puream_data[0];
 
148
    }
 
149
  }
 
150
  sj_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
151
  sk_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
152
  sl_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
153
  if (Symmetry.nirreps > 1) {
 
154
    if (BasisSet.puream)
 
155
      max_class_size = (2*BasisSet.max_am-1)*(2*BasisSet.max_am-1)*(2*BasisSet.max_am-1)*(2*BasisSet.max_am-1);
 
156
    else
 
157
      max_class_size = max_cart_class_size;
 
158
    for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
159
      plist_data[te_type] = block_matrix(max_num_unique_quartets,max_class_size);
 
160
    sj_fbf_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
161
    sk_fbf_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
162
    sl_fbf_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
163
  }
 
164
 
 
165
  max_num_prim_comb = (BasisSet.max_num_prims*
 
166
                       BasisSet.max_num_prims)*
 
167
                      (BasisSet.max_num_prims*
 
168
                       BasisSet.max_num_prims);
 
169
  UserOptions.memory -= init_libr12(&Libr12,BasisSet.max_am-1,max_num_prim_comb);
 
170
  init_fjt_table(&fjt_table);
 
171
 
 
172
/*-------------------------------------------------
 
173
  generate all unique shell quartets with ordering
 
174
  suitable for building the PK-matrix
 
175
 -------------------------------------------------*/
 
176
  for (usii=0; usii<Symmetry.num_unique_shells; usii++)
 
177
    for (usjj=0; usjj<=usii; usjj++)
 
178
      for (uskk=0; uskk<=usjj; uskk++)
 
179
        for (usll=0; usll<=uskk; usll++){
 
180
 
 
181
          /*--- Decide what shell quartets out of (ij|kl), (ik|jl), and (il|jk) are unique ---*/
 
182
          usi_arr[0] = usii; usj_arr[0] = usjj; usk_arr[0] = uskk; usl_arr[0] = usll;
 
183
          if (usii == usjj && usii == uskk || usjj == uskk && usjj == usll)
 
184
            num_unique_pk = 1;
 
185
          else if (usii == uskk || usjj == usll) {
 
186
            num_unique_pk = 2;
 
187
            usi_arr[1] = usii; usj_arr[1] = uskk; usk_arr[1] = usjj; usl_arr[1] = usll;
 
188
          }
 
189
          else if (usjj == uskk) {
 
190
            num_unique_pk = 2;
 
191
            usi_arr[1] = usii; usj_arr[1] = usll; usk_arr[1] = usjj; usl_arr[1] = uskk;
 
192
          }
 
193
          else if (usii == usjj || uskk == usll) {
 
194
            num_unique_pk = 2;
 
195
            usi_arr[1] = usii; usj_arr[1] = uskk; usk_arr[1] = usjj; usl_arr[1] = usll;
 
196
          }
 
197
          else {
 
198
            num_unique_pk = 3;
 
199
            usi_arr[1] = usii; usj_arr[1] = uskk; usk_arr[1] = usjj; usl_arr[1] = usll;
 
200
            usi_arr[2] = usii; usj_arr[2] = usll; usk_arr[2] = usjj; usl_arr[2] = uskk;
 
201
          }
 
202
 
 
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
            /* place in "ascending" angular momentum order-
 
208
               remember that [r12,T1] do not like it, I will have to worry about that later */
 
209
            /* these first two are good for the HRR */
 
210
            if(BasisSet.shells[Symmetry.us2s[usi]].am < BasisSet.shells[Symmetry.us2s[usj]].am){
 
211
              dum = usi;
 
212
              usi = usj;
 
213
              usj = dum;
 
214
            }
 
215
            if(BasisSet.shells[Symmetry.us2s[usk]].am < BasisSet.shells[Symmetry.us2s[usl]].am){
 
216
              dum = usk;
 
217
              usk = usl;
 
218
              usl = dum;
 
219
            }
 
220
            /* this should be /good/ for the VRR */
 
221
            if(BasisSet.shells[Symmetry.us2s[usi]].am + BasisSet.shells[Symmetry.us2s[usj]].am >
 
222
               BasisSet.shells[Symmetry.us2s[usk]].am + BasisSet.shells[Symmetry.us2s[usl]].am){
 
223
              dum = usi;
 
224
              usi = usk;
 
225
              usk = dum;
 
226
              dum = usj;
 
227
              usj = usl;
 
228
              usl = dum;
 
229
            }
 
230
 
 
231
            si = Symmetry.us2s[usi];
 
232
            sjj = Symmetry.us2s[usj];
 
233
            skk = Symmetry.us2s[usk];
 
234
            sll = Symmetry.us2s[usl];
 
235
            if (Symmetry.nirreps > 1) { /*--- Non-C1 symmetry case ---*/
 
236
              /*--- Generate the petite list of shell quadruplets using DCD approach of Davidson ---*/
 
237
              usp_ij = &(Symmetry.us_pairs[usi][usj]);
 
238
              usp_kl = &(Symmetry.us_pairs[usk][usl]);
 
239
              stab_i = Symmetry.atom_positions[BasisSet.shells[si].center-1];
 
240
              stab_j = Symmetry.atom_positions[BasisSet.shells[sjj].center-1];
 
241
              stab_k = Symmetry.atom_positions[BasisSet.shells[skk].center-1];
 
242
              stab_l = Symmetry.atom_positions[BasisSet.shells[sll].center-1];
 
243
              stab_ij = Symmetry.GnG[stab_i][stab_j];
 
244
              stab_kl = Symmetry.GnG[stab_k][stab_l];
 
245
              R_list = Symmetry.dcr[stab_i][stab_j];
 
246
              S_list = Symmetry.dcr[stab_k][stab_l];
 
247
              T_list = Symmetry.dcr[stab_ij][stab_kl];
 
248
              lambda_T = Symmetry.nirreps/Symmetry.dcr_deg[stab_ij][stab_kl];
 
249
              ni = (BasisSet.puream ? 2*BasisSet.shells[si].am - 1 : ioff[BasisSet.shells[si].am]);
 
250
              nj = (BasisSet.puream ? 2*BasisSet.shells[sjj].am - 1 : ioff[BasisSet.shells[sjj].am]);
 
251
              nk = (BasisSet.puream ? 2*BasisSet.shells[skk].am - 1 : ioff[BasisSet.shells[skk].am]);
 
252
              nl = (BasisSet.puream ? 2*BasisSet.shells[sll].am - 1 : ioff[BasisSet.shells[sll].am]);
 
253
              class_size = ni*nj*nk*nl;
 
254
 
 
255
              memset(sj_arr,0,sizeof(int)*max_num_unique_quartets);
 
256
              memset(sk_arr,0,sizeof(int)*max_num_unique_quartets);
 
257
              memset(sl_arr,0,sizeof(int)*max_num_unique_quartets);
 
258
              memset(sj_fbf_arr,0,sizeof(int)*max_num_unique_quartets);
 
259
              memset(sk_fbf_arr,0,sizeof(int)*max_num_unique_quartets);
 
260
              memset(sl_fbf_arr,0,sizeof(int)*max_num_unique_quartets);
 
261
              count = 0;
 
262
              for(dcr_ij=0;dcr_ij<Symmetry.dcr_dim[stab_i][stab_j];dcr_ij++){
 
263
                R = R_list[dcr_ij];
 
264
                sj = BasisSet.shells[sjj].trans_vec[R]-1;
 
265
                for(dcr_ijkl=0;dcr_ijkl<Symmetry.dcr_dim[stab_ij][stab_kl];dcr_ijkl++){
 
266
                  T = T_list[dcr_ijkl];
 
267
                  sk = BasisSet.shells[skk].trans_vec[T]-1;
 
268
                  slll = BasisSet.shells[sll].trans_vec[T]-1;
 
269
                  for(dcr_kl=0;dcr_kl<Symmetry.dcr_dim[stab_k][stab_l];dcr_kl++) {
 
270
                    S = S_list[dcr_kl];
 
271
                    sl = BasisSet.shells[slll].trans_vec[S]-1;
 
272
 
 
273
                    total_am = BasisSet.shells[si].am+BasisSet.shells[sj].am+BasisSet.shells[sk].am+BasisSet.shells[sl].am;
 
274
                    /*-------------------------------------------------------------
 
275
                      Obviously redundant or zero cases should be eliminated here!
 
276
                      Right now only zero case is eliminated. Redundancies arising
 
277
                      in DCD approach when usi == usj etc. may be eliminated too
 
278
                      but lambda_T will have to be replaced by an array (it won't
 
279
                      the same for every shell quartet in petite list anymore).
 
280
                     -------------------------------------------------------------*/
 
281
                    if(!(total_am%2)||
 
282
                       (BasisSet.shells[si].center!=BasisSet.shells[sj].center)||
 
283
                       (BasisSet.shells[sj].center!=BasisSet.shells[sk].center)||
 
284
                       (BasisSet.shells[sk].center!=BasisSet.shells[sl].center)) {
 
285
                      sj_arr[count] = sj;
 
286
                      sk_arr[count] = sk;
 
287
                      sl_arr[count] = sl;
 
288
                      sj_fbf_arr[count] = BasisSet.shells[sj].fbf-1;
 
289
                      sk_fbf_arr[count] = BasisSet.shells[sk].fbf-1;
 
290
                      sl_fbf_arr[count] = BasisSet.shells[sl].fbf-1;
 
291
                      count++;
 
292
                    }
 
293
                  }
 
294
                }
 
295
              } /* petite list is ready to be used */
 
296
              num_unique_quartets = count;
 
297
            }
 
298
            else { /*--- C1 symmetry case ---*/
 
299
              num_unique_quartets = 1;
 
300
              sj_arr[0] = usj;
 
301
              sk_arr[0] = usk;
 
302
              sl_arr[0] = usl;
 
303
              ni = (BasisSet.puream ? 2*BasisSet.shells[usi].am - 1 : ioff[BasisSet.shells[usi].am]);
 
304
              nj = (BasisSet.puream ? 2*BasisSet.shells[usj].am - 1 : ioff[BasisSet.shells[usj].am]);
 
305
              nk = (BasisSet.puream ? 2*BasisSet.shells[usk].am - 1 : ioff[BasisSet.shells[usk].am]);
 
306
              nl = (BasisSet.puream ? 2*BasisSet.shells[usl].am - 1 : ioff[BasisSet.shells[usl].am]);
 
307
              ioffset = BasisSet.shells[usi].fbf - 1;
 
308
              joffset = BasisSet.shells[usj].fbf - 1;
 
309
              koffset = BasisSet.shells[usk].fbf - 1;
 
310
              loffset = BasisSet.shells[usl].fbf - 1;
 
311
            }
 
312
 
 
313
            np_i = BasisSet.shells[si].n_prims;
 
314
            np_j = BasisSet.shells[sjj].n_prims;
 
315
            np_k = BasisSet.shells[skk].n_prims;
 
316
            np_l = BasisSet.shells[sll].n_prims;
 
317
            
 
318
            orig_am[0] = BasisSet.shells[si].am-1;
 
319
            orig_am[1] = BasisSet.shells[sjj].am-1;
 
320
            orig_am[2] = BasisSet.shells[skk].am-1;
 
321
            orig_am[3] = BasisSet.shells[sll].am-1;
 
322
            am = orig_am[0] + orig_am[1] + orig_am[2] + orig_am[3];
 
323
 
 
324
 
 
325
            /*----------------------------------
 
326
              Compute the nonredundant quartets
 
327
             ----------------------------------*/
 
328
            for(plquartet=0;plquartet<num_unique_quartets;plquartet++) {
 
329
              sj = sj_arr[plquartet];
 
330
              sk = sk_arr[plquartet];
 
331
              sl = sl_arr[plquartet];
 
332
 
 
333
              sp_ij = &(BasisSet.shell_pairs[si][sj]);
 
334
              sp_kl = &(BasisSet.shell_pairs[sk][sl]);
 
335
 
 
336
              Libr12.ShellQuartet.AB[0] = sp_ij->AB[0];
 
337
                Libr12.ShellQuartet.AB[1] = sp_ij->AB[1];
 
338
                Libr12.ShellQuartet.AB[2] = sp_ij->AB[2];
 
339
                Libr12.ShellQuartet.CD[0] = sp_kl->AB[0];
 
340
                Libr12.ShellQuartet.CD[1] = sp_kl->AB[1];
 
341
                Libr12.ShellQuartet.CD[2] = sp_kl->AB[2];
 
342
                Libr12.ShellQuartet.AC[0] = Molecule.centers[BasisSet.shells[si].center-1].x-
 
343
                        Molecule.centers[BasisSet.shells[sk].center-1].x;
 
344
                Libr12.ShellQuartet.AC[1] = Molecule.centers[BasisSet.shells[si].center-1].y-
 
345
                        Molecule.centers[BasisSet.shells[sk].center-1].y;
 
346
                Libr12.ShellQuartet.AC[2] = Molecule.centers[BasisSet.shells[si].center-1].z-
 
347
                        Molecule.centers[BasisSet.shells[sk].center-1].z;
 
348
                Libr12.ShellQuartet.ABdotAC = Libr12.ShellQuartet.AB[0]*Libr12.ShellQuartet.AC[0]+
 
349
                                              Libr12.ShellQuartet.AB[1]*Libr12.ShellQuartet.AC[1]+
 
350
                                              Libr12.ShellQuartet.AB[2]*Libr12.ShellQuartet.AC[2];
 
351
                Libr12.ShellQuartet.CDdotCA = -1.0*(Libr12.ShellQuartet.CD[0]*Libr12.ShellQuartet.AC[0]+
 
352
                                                    Libr12.ShellQuartet.CD[1]*Libr12.ShellQuartet.AC[1]+
 
353
                                                    Libr12.ShellQuartet.CD[2]*Libr12.ShellQuartet.AC[2]);
 
354
                AB2 = Libr12.ShellQuartet.AB[0]*Libr12.ShellQuartet.AB[0]+
 
355
                      Libr12.ShellQuartet.AB[1]*Libr12.ShellQuartet.AB[1]+
 
356
                      Libr12.ShellQuartet.AB[2]*Libr12.ShellQuartet.AB[2];
 
357
                CD2 = Libr12.ShellQuartet.CD[0]*Libr12.ShellQuartet.CD[0]+
 
358
                      Libr12.ShellQuartet.CD[1]*Libr12.ShellQuartet.CD[1]+
 
359
                      Libr12.ShellQuartet.CD[2]*Libr12.ShellQuartet.CD[2];
 
360
 
 
361
              /*--------------------------------
 
362
                contract by primitives out here
 
363
               --------------------------------*/
 
364
              num_prim_comb = 0;
 
365
              for (pi = 0; pi < np_i; pi++)
 
366
                for (pj = 0; pj < np_j; pj++)
 
367
                  for (pk = 0; pk < np_k; pk++)
 
368
                    for (pl = 0; pl < np_l; pl++){
 
369
                      r12_quartet_data(&(Libr12.PrimQuartet[num_prim_comb++]), &fjt_table, AB2, CD2,
 
370
                                       sp_ij, sp_kl, am, pi, pj, pk, pl, lambda_T);
 
371
                    }
 
372
 
 
373
              
 
374
              if (am) {
 
375
                build_r12_grt[orig_am[0]][orig_am[1]][orig_am[2]][orig_am[3]](&Libr12, num_prim_comb);
 
376
                /*--- copy transformed data to plist_data to be used in the symmetrization step ---*/
 
377
                if (Symmetry.nirreps > 1)
 
378
                  for(te_type=0;te_type<NUM_TE_TYPES;te_type++) {
 
379
                    data[te_type] = norm_quartet(Libr12.te_ptr[te_type], puream_data[te_type], orig_am, BasisSet.puream);
 
380
                    memcpy(plist_data[te_type][plquartet],data[te_type],sizeof(double)*class_size);
 
381
                  }
 
382
                /*--- or just transform data ---*/
 
383
                else
 
384
                  for(te_type=0;te_type<NUM_TE_TYPES;te_type++) {
 
385
                    data[te_type] = norm_quartet(Libr12.te_ptr[te_type], puream_data[te_type], orig_am, BasisSet.puream);
 
386
                  }
 
387
 
 
388
              }
 
389
              else {
 
390
                ssss = 0.0;
 
391
                ss_r12_ss = 0.0;
 
392
                for(p=0;p<num_prim_comb;p++) {
 
393
                  ssss += Libr12.PrimQuartet[p].F[0];
 
394
                  ss_r12_ss += Libr12.PrimQuartet[p].ss_r12_ss;
 
395
                }
 
396
                build_r12_grt[0][0][0][0](&Libr12,num_prim_comb);
 
397
                if (Symmetry.nirreps > 1) {
 
398
                  plist_data[0][plquartet][0] = ssss;
 
399
                  plist_data[1][plquartet][0] = ss_r12_ss;
 
400
                  plist_data[2][plquartet][0] = Libr12.te_ptr[2][0];
 
401
                  plist_data[3][plquartet][0] = Libr12.te_ptr[3][0];
 
402
                }
 
403
                else {
 
404
                    /*--------------------------------------------------------
 
405
                      This was so freaking careless on my part. I forgot that
 
406
                      in hrr_grt_order_0000 te_ptr[2] and te_ptr[3] point to
 
407
                      int_stack[1] and int_stack[0] respectively. Therefore I
 
408
                      have to copy these into int_stack[2] and int_stack[3] first
 
409
                      before copying ssss and ss_r12_ss into those locations.
 
410
                     --------------------------------------------------------*/
 
411
                  Libr12.int_stack[2] = Libr12.te_ptr[2][0];
 
412
                  Libr12.int_stack[3] = Libr12.te_ptr[3][0];
 
413
                  Libr12.int_stack[0] = ssss;
 
414
                  Libr12.int_stack[1] = ss_r12_ss;
 
415
                  data[0] = Libr12.int_stack;
 
416
                  data[1] = Libr12.int_stack+1;
 
417
                  data[2] = Libr12.int_stack+2;
 
418
                  data[3] = Libr12.int_stack+3;
 
419
                }
 
420
              }
 
421
 
 
422
            } /* end of computing "petit" list */
 
423
 
 
424
            memset(num,0,NUM_TE_TYPES*sizeof(int));
 
425
            if (Symmetry.nirreps > 1) { /*--- Non-C1 case ---*/
 
426
            /*------------------------------------------------------------------------
 
427
              Now we have everything to build SO's. Need to distinguish several cases
 
428
              that are slightly different from each other. To avoid extra if's inside
 
429
              the loops I separated them:
 
430
              1) usi == usj == usk == usl
 
431
              2) usi == usj != usk == usl
 
432
              3) usi == usk != usj == usl
 
433
              4) usi == usj
 
434
              5) usk == usl
 
435
              6) general case
 
436
              "Symmetrization" is based on Pitzer's equal contribution theorem.
 
437
 
 
438
              The general procedure is as follows:
 
439
              1) loop over irreps, if there are any SO products that arise from
 
440
                 this usp_ij pair and belong to this irrep, loop over them (ij)
 
441
              2) from ij figure out product of which SOs this is, and what AOs
 
442
                 (or should I say, basis functions) contribute to these SO.
 
443
                 so_i and so_j are absolute indices, but i and j are relative
 
444
                 indices (i and j are indices of cart/puream components in usi
 
445
                 and usj)
 
446
              3) if there are any products of SOs from usp_kl that belong to irrep,
 
447
                 loop over those (kl). Obviously, the whole point of introducing
 
448
                 restrictions on different cases was to improve efficiency, but now
 
449
                 it seems like a stupid idea. I should have left just one loop case.
 
450
 
 
451
                 For more comments see Default_Ints/te_ints.c
 
452
             ------------------------------------------------------------------------*/
 
453
            bf_i = BasisSet.shells[si].fbf-1;
 
454
            if (usi == usj && usi == usk && usi == usl || usi == usk && usj == usl)
 
455
              for(irrep=0;irrep<Symmetry.nirreps;irrep++) {
 
456
                if (npi_ij = usp_ij->SOpair_npi[irrep])
 
457
                  for(ij=0;ij<npi_ij;ij++) {
 
458
                    so_i = usp_ij->SOpair_so_i[irrep][ij];
 
459
                    so_j = usp_ij->SOpair_so_j[irrep][ij];
 
460
                    i = usp_ij->SOpair_bf_i[irrep][ij];
 
461
                    j = usp_ij->SOpair_bf_j[irrep][ij];
 
462
                    ind_offset = (i*nj + j)*nk*nl;
 
463
                    for(kl=0;kl<=ij;kl++) {
 
464
                      so_k = usp_kl->SOpair_so_i[irrep][kl];
 
465
                      so_l = usp_kl->SOpair_so_j[irrep][kl];
 
466
                      k = usp_kl->SOpair_bf_i[irrep][kl];
 
467
                      l = usp_kl->SOpair_bf_j[irrep][kl];
 
468
                      index = ind_offset + k*nl + l;
 
469
                      memset(so_int,0,NUM_TE_TYPES*sizeof(double));
 
470
                      for(s=0;s<num_unique_quartets;s++){
 
471
                        ioffset = bf_i + i;
 
472
                        bf_j = sj_fbf_arr[s]+j;
 
473
                        bf_k = sk_fbf_arr[s]+k;
 
474
                        bf_l = sl_fbf_arr[s]+l;
 
475
                        temp = Symmetry.usotao[so_i][ioffset]*
 
476
                               Symmetry.usotao[so_j][bf_j]*
 
477
                               Symmetry.usotao[so_k][bf_k]*
 
478
                               Symmetry.usotao[so_l][bf_l];
 
479
                        for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
480
                          so_int[te_type] += temp*plist_data[te_type][s][index];
 
481
                      }
 
482
                      for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
483
                        if (fabs(so_int[te_type])>toler) {
 
484
                          so_ii = so_i;
 
485
                          so_jj = so_j;
 
486
                          so_kk = so_k;
 
487
                          so_ll = so_l;
 
488
                          if (so_ii < so_jj) {
 
489
                            SWAP(so_ii,so_jj);
 
490
                            if (te_type == 2) /* [r12,T1] is antisymmetric WRT i and j */
 
491
                              so_int[te_type] *= -1.0;
 
492
                          }
 
493
                          if (so_kk < so_ll) {
 
494
                            SWAP(so_kk,so_ll);
 
495
                            if (te_type == 3) /* [r12,T2] is antisymmetric WRT to k and l */
 
496
                              so_int[te_type] *= -1.0;
 
497
                          }
 
498
                          if ((so_ii < so_kk) || (so_ii == so_kk && so_jj < so_ll)) {
 
499
                            if (te_type < 2) { /* only eri and r12 ints have bra-ket symmetry */
 
500
                              SWAP(so_ii,so_kk);
 
501
                              SWAP(so_jj,so_ll);
 
502
                            }
 
503
                          }
 
504
                          tot_data[te_type][num[te_type]].i = (short int) so_ii;
 
505
                          tot_data[te_type][num[te_type]].j = (short int) so_jj;
 
506
                          tot_data[te_type][num[te_type]].k = (short int) so_kk;
 
507
                          tot_data[te_type][num[te_type]].l = (short int) so_ll;
 
508
                          tot_data[te_type][num[te_type]].val = so_int[te_type];
 
509
                          num[te_type]++;
 
510
                        }
 
511
                    }
 
512
                  }
 
513
              }
 
514
            else
 
515
              for(irrep=0;irrep<Symmetry.nirreps;irrep++) {
 
516
                if ((npi_ij = usp_ij->SOpair_npi[irrep]) && (npi_kl = usp_kl->SOpair_npi[irrep]))
 
517
                  for(ij=0;ij<npi_ij;ij++) {
 
518
                    i = usp_ij->SOpair_bf_i[irrep][ij];
 
519
                    j = usp_ij->SOpair_bf_j[irrep][ij];
 
520
                    so_i = usp_ij->SOpair_so_i[irrep][ij];
 
521
                    so_j = usp_ij->SOpair_so_j[irrep][ij];
 
522
                    ind_offset = (i*nj + j)*nk*nl;
 
523
                    for(kl=0;kl<npi_kl;kl++) {
 
524
                      k = usp_kl->SOpair_bf_i[irrep][kl];
 
525
                      l = usp_kl->SOpair_bf_j[irrep][kl];
 
526
                      so_k = usp_kl->SOpair_so_i[irrep][kl];
 
527
                      so_l = usp_kl->SOpair_so_j[irrep][kl];
 
528
                      index = ind_offset + k*nl + l;
 
529
                      memset(so_int,0,NUM_TE_TYPES*sizeof(double));
 
530
                      for(s=0;s<num_unique_quartets;s++){
 
531
                        ioffset = bf_i + i;
 
532
                        bf_j = sj_fbf_arr[s]+j;
 
533
                        bf_k = sk_fbf_arr[s]+k;
 
534
                        bf_l = sl_fbf_arr[s]+l;
 
535
                        temp = Symmetry.usotao[so_i][ioffset]*
 
536
                               Symmetry.usotao[so_j][bf_j]*
 
537
                               Symmetry.usotao[so_k][bf_k]*
 
538
                               Symmetry.usotao[so_l][bf_l];
 
539
                        for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
540
                          so_int[te_type] += temp*plist_data[te_type][s][index];
 
541
                      }
 
542
                      for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
543
                        if (fabs(so_int[te_type])>toler) {
 
544
                          so_ii = so_i;
 
545
                          so_jj = so_j;
 
546
                          so_kk = so_k;
 
547
                          so_ll = so_l;
 
548
                          if (so_ii < so_jj) {
 
549
                            SWAP(so_ii,so_jj);
 
550
                            if (te_type == 2) /* [r12,T1] is antisymmetric WRT i and j */
 
551
                              so_int[te_type] *= -1.0;
 
552
                          }
 
553
                          if (so_kk < so_ll) {
 
554
                            SWAP(so_kk,so_ll);
 
555
                            if (te_type == 3) /* [r12,T2] is antisymmetric WRT to k and l */
 
556
                              so_int[te_type] *= -1.0;
 
557
                          }
 
558
                          if ((so_ii < so_kk) || (so_ii == so_kk && so_jj < so_ll)) {
 
559
                            if (te_type < 2) { /* only eri and r12 ints have bra-ket symmetry */
 
560
                              SWAP(so_ii,so_kk);
 
561
                              SWAP(so_jj,so_ll);
 
562
                            }
 
563
                          }
 
564
                          tot_data[te_type][num[te_type]].i = (short int) so_ii;
 
565
                          tot_data[te_type][num[te_type]].j = (short int) so_jj;
 
566
                          tot_data[te_type][num[te_type]].k = (short int) so_kk;
 
567
                          tot_data[te_type][num[te_type]].l = (short int) so_ll;
 
568
                          tot_data[te_type][num[te_type]].val = so_int[te_type];
 
569
                          num[te_type]++;
 
570
                        }
 
571
                    }
 
572
                  }
 
573
              }
 
574
            }
 
575
            else { /*--- C1 symmetry ---*/
 
576
              /*--- Here just put non-redundant integrals to tot_data ---*/
 
577
              if(usi==usj&&usk==usl&&usi==usk) { /*--- All shells are the same - the (aa|aa) case ---*/
 
578
                iimax = ni - 1;
 
579
                for(ii=0; ii <= iimax; ii++){
 
580
                  jjmax = ii;
 
581
                  for(jj=0; jj <= jjmax; jj++){
 
582
                    kkmax = ii;
 
583
                    for(kk=0; kk <= kkmax; kk++){
 
584
                      llmax = (kk==ii)? jj : kk ;
 
585
                      for(ll=0; ll <= llmax; ll++){
 
586
                        index = ll+nl*(kk+nk*(jj+nj*ii));
 
587
                        for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
588
                          if(fabs(data[te_type][index])>toler){
 
589
                            tot_data[te_type][num[te_type]].i = (short int) (ii+ioffset);
 
590
                            tot_data[te_type][num[te_type]].j = (short int) (jj+joffset);
 
591
                            tot_data[te_type][num[te_type]].k = (short int) (kk+koffset);
 
592
                            tot_data[te_type][num[te_type]].l = (short int) (ll+loffset);
 
593
                            tot_data[te_type][num[te_type]].val = data[te_type][index];
 
594
                            num[te_type]++;
 
595
                          }
 
596
                      }
 
597
                    }
 
598
                  }
 
599
                }
 
600
              }
 
601
              else if(usi==usk && usj==usl){    /*--- The (ab|ab) case ---*/
 
602
                iimax = ni - 1;
 
603
                for(ii=0; ii <= iimax; ii++){
 
604
                  jjmax = nj - 1;
 
605
                  for(jj=0; jj <= jjmax; jj++){
 
606
                    kkmax = ii;
 
607
                    for(kk=0; kk <= kkmax; kk++){
 
608
                      llmax = (kk==ii)? jj : nl - 1;
 
609
                      for(ll=0; ll <= llmax; ll++){
 
610
                        index = ll+nl*(kk+nk*(jj+nj*ii));
 
611
                        for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
612
                          if(fabs(data[te_type][index])>toler){
 
613
                            i = ii + ioffset;
 
614
                            j = jj + joffset;
 
615
                            k = kk + koffset;
 
616
                            l = ll + loffset;
 
617
                            value = data[te_type][index];
 
618
                            if (i < j) {
 
619
                              SWAP(i,j);
 
620
                              SWAP(k,l);
 
621
                              if (te_type > 1) /* [r12,Ti] are antisymmetric WRT to i -> j, k -> l */
 
622
                                value *= -1.0;
 
623
                            }
 
624
                            if (i < k) {
 
625
                              if (te_type < 2) { /* only eri and r12 ints have bra-ket symmetry */
 
626
                                SWAP(i,k);
 
627
                                SWAP(j,l);
 
628
                              }
 
629
                            }
 
630
                            tot_data[te_type][num[te_type]].i = (short int) i;
 
631
                            tot_data[te_type][num[te_type]].j = (short int) j;
 
632
                            tot_data[te_type][num[te_type]].k = (short int) k;
 
633
                            tot_data[te_type][num[te_type]].l = (short int) l;
 
634
                            tot_data[te_type][num[te_type]].val = value;
 
635
                            num[te_type]++;
 
636
                          }
 
637
                      }
 
638
                    }
 
639
                  }
 
640
                }
 
641
              }
 
642
              else {   /*--- The (ab|cd) case ---*/
 
643
                iimax = ni - 1;
 
644
                kkmax = nk - 1;
 
645
                for(ii=0; ii <= iimax; ii++){
 
646
                  jjmax = (usi==usj) ? ii : nj - 1;
 
647
                  for(jj=0; jj <= jjmax; jj++){
 
648
                    for(kk=0; kk <= kkmax; kk++){
 
649
                      llmax = (usk==usl) ? kk : nl - 1;
 
650
                      for(ll=0; ll <= llmax; ll++){
 
651
                        index = ll+nl*(kk+nk*(jj+nj*ii));
 
652
                        for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
653
                          if(fabs(data[te_type][index])>toler){
 
654
                            i = ii + ioffset;
 
655
                            j = jj + joffset;
 
656
                            k = kk + koffset;
 
657
                            l = ll + loffset;
 
658
                            value = data[te_type][index];
 
659
                            if (i < j) {
 
660
                              SWAP(i,j);
 
661
                              if (te_type == 2) /* [r12,T1] is antisymmetric WRT i and j */
 
662
                                value *= -1.0;
 
663
                            }
 
664
                            if (k < l) {
 
665
                              SWAP(k,l);
 
666
                              if (te_type == 3) /* [r12,T2] is antisymmetric WRT k and l */
 
667
                                value *= -1.0;
 
668
                            }
 
669
                            if ((i < k) || (i == k && j < l)) {
 
670
                              if (te_type < 2) { /* only eri and r12 ints have bra-ket symmetry */
 
671
                                SWAP(i,k);
 
672
                                SWAP(j,l);
 
673
                              }
 
674
                            }
 
675
                            tot_data[te_type][num[te_type]].i = (short int) i;
 
676
                            tot_data[te_type][num[te_type]].j = (short int) j;
 
677
                            tot_data[te_type][num[te_type]].k = (short int) k;
 
678
                            tot_data[te_type][num[te_type]].l = (short int) l;
 
679
                            tot_data[te_type][num[te_type]].val = value;
 
680
                            num[te_type]++;
 
681
                          }
 
682
                      }
 
683
                    }
 
684
                  }
 
685
                }
 
686
              }
 
687
            }
 
688
 
 
689
 
 
690
            /*---------------
 
691
              Write out ERIs
 
692
             ---------------*/
 
693
            if (num[0]) { /* Let's see if we need to write out something */
 
694
              total_te_count[0] += num[0];
 
695
              if (upk == num_unique_pk - 1) /* if this is the last quartet needed for a pk-block - let CSCF know
 
696
                                               by setting index i of the last integral to negative of itself.
 
697
                                               The only guy where this trick won't work will be (00|00).
 
698
                                               But normally (00|00) is of (ss|ss) type and is enough for computing one
 
699
                                               pk-matrix element. PK-buffer won't get overfull because of this one guy */
 
700
                tot_data[0][num[0]-1].i = -tot_data[0][num[0]-1].i;
 
701
              iwl_buf_wrt_struct_nocut(&TEOUT[0], tot_data[0], num[0]);
 
702
            }
 
703
 
 
704
            /*----------------
 
705
              Write out r12's
 
706
             ----------------*/
 
707
            if (num[1]) { /* Let's see if we need to write out something */
 
708
              total_te_count[1] += num[1];
 
709
              iwl_buf_wrt_struct_nocut(&TEOUT[1], tot_data[1], num[1]);
 
710
            }
 
711
 
 
712
            if (NUM_TE_TYPES > 2) {
 
713
              /*---------------------
 
714
                Write out [r12,T1]'s
 
715
               ---------------------*/
 
716
              if (num[2]) { /* Let's see if we need to write out something */
 
717
                total_te_count[2] += num[2];
 
718
                iwl_buf_wrt_struct_nocut(&TEOUT[2], tot_data[2], num[2]);
 
719
              }
 
720
 
 
721
              /*-----------------------------------
 
722
                Convert [r12,T2]'s into [r12,T1]'s
 
723
               -----------------------------------*/
 
724
              for(i=0;i<num[3];i++)
 
725
                if (tot_data[3][i].i == tot_data[3][i].k && tot_data[3][i].j == tot_data[3][i].l) {
 
726
                  /* Do NOT write out if ij == kl because (ij|[r12,T1]|ij) = (ij|[r12,T2]|ij) */
 
727
                  tot_data[3][i].val = 0.0;
 
728
                  total_te_count[2]--;
 
729
                  continue;
 
730
                }
 
731
                else {
 
732
                  SWAP(tot_data[3][i].i,tot_data[3][i].k);
 
733
                  SWAP(tot_data[3][i].j,tot_data[3][i].l);
 
734
                }
 
735
              if (num[3]) { /* Let's see if we need to write out something */
 
736
                total_te_count[2] += num[3];
 
737
                iwl_buf_wrt_struct(&TEOUT[2], tot_data[3], num[3], toler);
 
738
              }
 
739
            }
 
740
            
 
741
#if PRINT
 
742
            for(te_type=0;te_type<NUM_TE_TYPES;te_type++) {
 
743
              if (te_type < 3) {
 
744
                for(n=0; n<num[te_type]; n++){
 
745
                  fprintf(teout[te_type], "%5d%5d%5d%5d%20.10lf\n",
 
746
                          abs(tot_data[te_type][n].i)+1, 
 
747
                          tot_data[te_type][n].j+1, 
 
748
                          tot_data[te_type][n].k+1, 
 
749
                          tot_data[te_type][n].l+1, 
 
750
                          tot_data[te_type][n].val);
 
751
                }
 
752
                fflush(teout[te_type]);
 
753
              }
 
754
              else { /* We've swapped ij and kl in [r12,T2] case, remember? */
 
755
                for(n=0; n<num[te_type]; n++){
 
756
                  fprintf(teout[te_type], "%5d%5d%5d%5d%20.10lf\n",
 
757
                          tot_data[te_type][n].k+1, 
 
758
                          tot_data[te_type][n].l+1, 
 
759
                          tot_data[te_type][n].i+1, 
 
760
                          tot_data[te_type][n].j+1, 
 
761
                          tot_data[te_type][n].val);
 
762
                }
 
763
                fflush(teout[te_type]);
 
764
              }
 
765
            }
 
766
#endif
 
767
          } /* end of computing PK-quartets. */
 
768
 
 
769
        } /* end getting unique shell combination */
 
770
 
 
771
  for(te_type=0;te_type<NUM_TE_TYPES-1;te_type++) {
 
772
    iwl_buf_flush(&TEOUT[te_type], 1);
 
773
    iwl_buf_close(&TEOUT[te_type], 1);
 
774
  }
 
775
 
 
776
#if PRINT
 
777
  for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
778
    fclose(teout[te_type]);
 
779
#endif
 
780
 
 
781
 
 
782
  for(te_type=0;te_type<NUM_TE_TYPES-1;te_type++)
 
783
    fprintf(outfile,"    Wrote %d two-electron integrals of %s to IWL file %2d\n",
 
784
            total_te_count[te_type],te_operator[te_type],itapTE[te_type]);
 
785
  fprintf(outfile,"\n");
 
786
 
 
787
  /*---------
 
788
    Clean-up
 
789
   ---------*/
 
790
  free_libr12(&Libr12);
 
791
  free_fjt(&fjt_table);
 
792
  for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
793
    free(tot_data[te_type]);
 
794
  free(sj_arr);
 
795
  free(sk_arr);
 
796
  free(sl_arr);
 
797
  if (Symmetry.nirreps > 1) {
 
798
    for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
799
      free_block(plist_data[te_type]);
 
800
    free(sj_fbf_arr);
 
801
    free(sk_fbf_arr);
 
802
    free(sl_fbf_arr);
 
803
  }
 
804
  if (BasisSet.puream) {
 
805
    if (Symmetry.nirreps > 1)
 
806
      free(puream_data[0]);
 
807
    else
 
808
      for(te_type=0;te_type<NUM_TE_TYPES;te_type++)
 
809
        free(puream_data[te_type]);
 
810
  }
 
811
 
 
812
  return;
 
813
}
 
814
 
 
815
};};