← 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_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/Default_Deriv1/te_deriv1_ints.cc
 
1
/*! \file
 
2
    \ingroup CINTS
 
3
    \brief Enter brief description of file here 
 
4
*/
1
5
#include <vector>
2
6
#include <sstream>
3
7
#include <algorithm>
4
 
 
5
 
#include<math.h>
6
 
#include<stdio.h>
7
 
#include<string.h>
 
8
#include<cmath>
 
9
#include<cstring>
 
10
#include<cstdio>
8
11
#include<memory.h>
9
 
#include<stdlib.h>
 
12
#include<cstdlib>
10
13
 
11
 
extern "C" {
12
14
#include<psitypes.h>
13
15
#include<libipv1/ip_lib.h>
14
16
#include<libiwl/iwl.h>
19
21
#include"defines.h"
20
22
#define EXTERN
21
23
#include"global.h"
 
24
#include <stdexcept>
22
25
#include"schwartz.h"
23
26
#include"deriv1_quartet_data.h"
24
27
#include"iwl_tebuf.h"
28
31
#else
29
32
  #include"int_fjt.h"
30
33
#endif
31
 
}
32
34
 
33
35
using std::vector;
34
36
using std::min;
35
37
 
36
 
// For a set of SALCs, we need to be able to select coefficients of all nonzero derivatives
37
 
// of a given shell quartet. It is best done when the relevant subset of the (sparse) SALC coefficient
38
 
// matrix is pre-arranged in an array
39
 
typedef struct {
40
 
  double coef;    // contribution of this cartesian derivative to this SALC
41
 
  int cart_der;   // cartesian derivative (for this quartet) [0,3*num_unique_atoms)
42
 
  int salc;       // SALC index [0,num_atoms*3)
43
 
} cdsalc_elem;
44
 
 
45
 
 
46
 
extern "C" {
47
 
  
48
 
void te_deriv1_ints()
49
 
{
50
 
  const double toler = UserOptions.cutoff;
51
 
  const int num_coords = 3*Molecule.num_atoms;        // total number of coordinates, symmetry-adapted or otherwise
52
 
  // max number of nonzero derivatives for a SO shell quartet = 4 (centers) * 3 (x,y,z) * maximum multiplicity of a nucleus
53
 
  const int max_num_cdsalcs_per_quartet = 12*Symmetry.max_stab_index;
54
 
 
55
 
  /*--- ASCII file to print integrals ---*/
56
 
  FILE *d1eriout;
57
 
 
58
 
  /*--- Various data structures ---*/
59
 
  struct iwlbuf* D1ERIOUT;               /* IWL buffer for target integrals */
60
 
  struct tebuf** tot_data;               /* accum. for contracted integrals */
61
 
  vector<int> num_of_ints_in_totdata;    /* counts how many integrals are in each tot_data */
62
 
  struct shell_pair *sp_ij, *sp_kl;
63
 
  struct unique_shell_pair *usp_ij,*usp_kl;
64
 
  
65
 
  Libderiv_t Libderiv;
66
 
#ifndef USE_TAYLOR_FM
67
 
  double_array_t fjt_table;
68
 
#endif
69
 
  
70
 
  std::vector<int> unique_center;  unique_center.reserve(4);
71
 
  std::vector<int> nuc(4);
72
 
  
73
 
  PSI_INT_LEAST64 total_te_count = 0;
74
 
  vector<PSI_INT_LEAST64> te_count_per_coord(num_coords,0);     // keeps track of how many integrals were written out
75
 
  int ij, kl, ik, jl, ijkl;
76
 
  int ioffset, joffset, koffset, loffset;
77
 
  int count ;
78
 
  int dum;
79
 
  int total_am, am;
80
 
  int orig_am[4];
81
 
  int pkblock_end_index = -1;
82
 
  int i, j, k, l, m, ii, jj, kk, ll;
83
 
  int si, sj, sk, sl ;
84
 
  int sii, sjj, skk, sll , slll;
85
 
  int pi, pj, pk, pl;
86
 
  int max_pj, max_pl;
87
 
  int upk, num_unique_pk;
88
 
  int usi_arr[3], usj_arr[3], usk_arr[3], usl_arr[3];
89
 
  int *sj_arr, *sk_arr, *sl_arr;
90
 
  int *sj_fbf_arr, *sk_fbf_arr, *sl_fbf_arr;
91
 
  int usi,usj,usk,usl;
92
 
  int stab_i,stab_j,stab_k,stab_l,stab_ij,stab_kl;
93
 
  int *R_list, *S_list, *T_list;
94
 
  int R,S,T;
95
 
  int dcr_ij, dcr_kl, dcr_ijkl;
96
 
  int lambda_T = 1;
97
 
  int num_unique_quartets;
98
 
  int plquartet;
99
 
  int max_num_unique_quartets;
100
 
  int max_num_prim_comb;
101
 
 
102
 
  int size, class_size;
103
 
  int max_cart_class_size;
104
 
 
105
 
  int bf_i, bf_j, bf_k, bf_l, so_i, so_j, so_k, so_l, s;
106
 
  int np_i, np_j, np_k, np_l;
107
 
  int ni, nj, nk, nl;
108
 
 
109
 
  int index;
110
 
  int iimax, jjmax, kkmax, llmax;
111
 
  int irrep, npi_ij, npi_kl, npi_ik, npi_jl, ind_offset;
112
 
 
113
 
  int num_prim_comb, p;
114
 
 
115
 
  double AB2, CD2;
116
 
  double pkblock_end_value = 0.0;
117
 
  double temp;
118
 
  
119
 
  vector<double> so_int;
120
 
  {
121
 
    const int max_num_cdsalc_per_quartet = min(max_num_cdsalcs_per_quartet,CDSALCs.nsalcs);
122
 
    so_int.reserve(max_num_cdsalc_per_quartet);
123
 
  }
124
 
  
125
 
  /*---------------
126
 
    Initialization
127
 
   ---------------*/
128
 
#if PRINT
129
 
  eriout = fopen("d1eriout.dat","w");
130
 
#endif
131
 
  D1ERIOUT = new struct iwlbuf[num_coords];
132
 
  for(int c=0; c<num_coords; c++)
133
 
    iwl_buf_init(&D1ERIOUT[c], IOUnits.itapD1ERI_SO+c, toler, 0, 0);
134
 
#ifdef USE_TAYLOR_FM
135
 
  init_Taylor_Fm_Eval(BasisSet.max_am*4-4+DERIV_LVL,UserOptions.cutoff);
136
 
#else
137
 
  init_fjt(BasisSet.max_am*4+DERIV_LVL);
138
 
  init_fjt_table(&fjt_table);
139
 
#endif
140
 
  init_libderiv_base();
141
 
  
142
 
  /*-------------------------
143
 
    Allocate data structures
144
 
   -------------------------*/
145
 
  max_cart_class_size = (ioff[BasisSet.max_am])*
146
 
                        (ioff[BasisSet.max_am])*
147
 
                        (ioff[BasisSet.max_am])*
148
 
                        (ioff[BasisSet.max_am]);
149
 
  max_num_unique_quartets = Symmetry.max_stab_index*
150
 
                            Symmetry.max_stab_index*
151
 
                            Symmetry.max_stab_index;
152
 
  
153
 
  // Final integrals are stored here
154
 
  tot_data = new struct tebuf*[num_coords];
155
 
  for(int c=0; c<num_coords; c++) {
156
 
    tot_data[c] = new struct tebuf[max_num_unique_quartets*max_cart_class_size];
157
 
  }
158
 
  num_of_ints_in_totdata.resize(num_coords);
159
 
  for(int c=0; c<num_coords; c++)
160
 
    num_of_ints_in_totdata[c] = 0;
161
 
  
162
 
  // These arrays are used to hold cartesian AO derivative integrals
163
 
  double* cart_ints[12];
164
 
  for(int i=0; i<12; i++)
165
 
    cart_ints[i] = new double[max_cart_class_size];
166
 
  // plist_ints holds all symmetry-unique AO shell quartets for a given SO shell quartet
167
 
  // plist_ints[i][j] is the pointer to j-th *total* derivative of i-th AO shell quartet which contributes to the current SO quartet
168
 
  // note that only total derivatives are stored, i.e. partial derivatives produced by libderiv are combined
169
 
  // to yield total derivatives
170
 
  double*** plist_ints = new double**[max_num_unique_quartets];
171
 
  for(int i=0; i<max_num_unique_quartets; i++) {
172
 
    plist_ints[i] = new double*[12];
173
 
    for(int j=0; j<12; j++) {
174
 
      plist_ints[i][j] = new double[max_cart_class_size];
175
 
    }
176
 
  }
177
 
  // For every SO shell quartet we need the following information:
178
 
  // 1) irreps of all total derivative operators which give nonzero when applied to the SO quartet
179
 
  // 2) given an irrep, which total derivative operators contribute
180
 
  // For every AO shell quartet in a petit list for a given SO shell quartet we need the following information:
181
 
  // 1) number of nonzero total derivatives (same as number of distinct centers, n, times 3, 1 <= n <= 4)
182
 
  // 2) irreps of all total derivative operators which give nonzero when applied to the AO quartet
183
 
  // 3) given an irrep, which total derivative operators contribute
184
 
 
 
38
 
 
39
namespace psi { namespace CINTS {
185
40
  // For a set of SALCs, we need to be able to select coefficients of all nonzero derivatives
186
41
  // of a given shell quartet. It is best done when the relevant subset of the (sparse) SALC coefficient
187
42
  // matrix is pre-arranged in an array
188
 
  // plist_salcs[s][g][i] is an i-th nonzero contribution to a SALC of irrep g
189
 
  // which results in nonzero when applied to member s of a petite list
190
 
  struct cdsalc_elem_vec {
191
 
    int nelems;
192
 
    cdsalc_elem* elems;
193
 
  };
194
 
 
195
 
  cdsalc_elem_vec** plist_salcs = new cdsalc_elem_vec*[max_num_unique_quartets];
196
 
  for(int i=0; i<max_num_unique_quartets; i++) {
197
 
    plist_salcs[i] = new cdsalc_elem_vec[Symmetry.nirreps];
198
 
    for(int j=0; j<Symmetry.nirreps; j++) {
199
 
      plist_salcs[i][j].nelems=0;
200
 
      plist_salcs[i][j].elems = new cdsalc_elem[max_num_cdsalcs_per_quartet];
201
 
    }
202
 
  }
203
 
  
204
 
  int* salc_all2thisquartet = new int[num_coords];    // maps absolute SALC index to the SALC index for this quartet
205
 
  int* salc_thisquartet2all = new int[max_num_cdsalcs_per_quartet];    // the reverse of the above
206
 
  
207
 
  sj_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
208
 
  sk_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
209
 
  sl_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
210
 
  if (Symmetry.nirreps > 1) {
211
 
    sj_fbf_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
212
 
    sk_fbf_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
213
 
    sl_fbf_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
214
 
  }
215
 
  
216
 
  max_num_prim_comb = (BasisSet.max_num_prims*
217
 
                       BasisSet.max_num_prims)*
218
 
                      (BasisSet.max_num_prims*
219
 
                       BasisSet.max_num_prims);
220
 
 
221
 
  init_libderiv1(&Libderiv, BasisSet.max_am-1, max_num_prim_comb, max_cart_class_size);
222
 
 
223
 
  /*------------------------------------
224
 
    generate all unique shell quartets
225
 
   ------------------------------------*/
226
 
  if(UserOptions.print_lvl >= PRINT_DEBUG) {
227
 
    fprintf(outfile,"  -Electron repulsion integrals:\n\n");
228
 
  }
229
 
  
230
 
  for (int usii=0; usii<Symmetry.num_unique_shells; usii++)
231
 
    for (int usjj=0; usjj<=usii; usjj++)
232
 
      for (int uskk=0; uskk<Symmetry.num_unique_shells; uskk++) {
233
 
        const int usll_max = (usii==uskk ? usjj : uskk);
234
 
        for (int usll=0; usll<=usll_max; usll++){
235
 
          
236
 
          int usi = usii;
237
 
          int usj = usjj;
238
 
          int usk = uskk;
239
 
          int usl = usll;
240
 
          
241
 
          /* place in "ascending" angular mom-
242
 
          my simple way of optimizing PHG recursion (VRR) */
243
 
          /* these first two are good for the HRR */
244
 
          if(BasisSet.shells[Symmetry.us2s[usi]].am < BasisSet.shells[Symmetry.us2s[usj]].am){
245
 
            dum = usi;
246
 
            usi = usj;
247
 
            usj = dum;
248
 
          }
249
 
          if(BasisSet.shells[Symmetry.us2s[usk]].am < BasisSet.shells[Symmetry.us2s[usl]].am){
250
 
            dum = usk;
251
 
            usk = usl;
252
 
            usl = dum;
253
 
          }
254
 
          /* this should be /good/ for the VRR */
255
 
          if(BasisSet.shells[Symmetry.us2s[usi]].am + BasisSet.shells[Symmetry.us2s[usj]].am >
256
 
          BasisSet.shells[Symmetry.us2s[usk]].am + BasisSet.shells[Symmetry.us2s[usl]].am){
257
 
            dum = usi;
258
 
            usi = usk;
259
 
            usk = dum;
260
 
            dum = usj;
261
 
            usj = usl;
262
 
            usl = dum;
263
 
          }
264
 
          
265
 
          si = Symmetry.us2s[usi];
266
 
          sjj = Symmetry.us2s[usj];
267
 
          skk = Symmetry.us2s[usk];
268
 
          sll = Symmetry.us2s[usl];
269
 
          
270
 
          int center_i = BasisSet.shells[si].center-1;
271
 
          int center_j = BasisSet.shells[sjj].center-1;
272
 
          int center_k = BasisSet.shells[skk].center-1;
273
 
          int center_l = BasisSet.shells[sll].center-1;
274
 
          
275
 
          if (Symmetry.nirreps > 1) { /*--- Non-C1 symmetry case ---*/
276
 
            /*--- Generate the petite list of shell quadruplets using DCD approach of Davidson ---*/
277
 
            usp_ij = &(Symmetry.us_pairs[usi][usj]);
278
 
            usp_kl = &(Symmetry.us_pairs[usk][usl]);
279
 
            stab_i = Symmetry.atom_positions[center_i];
280
 
            stab_j = Symmetry.atom_positions[center_j];
281
 
            stab_k = Symmetry.atom_positions[center_k];
282
 
            stab_l = Symmetry.atom_positions[center_l];
283
 
            stab_ij = Symmetry.GnG[stab_i][stab_j];
284
 
            stab_kl = Symmetry.GnG[stab_k][stab_l];
285
 
            R_list = Symmetry.dcr[stab_i][stab_j];
286
 
            S_list = Symmetry.dcr[stab_k][stab_l];
287
 
            T_list = Symmetry.dcr[stab_ij][stab_kl];
288
 
            lambda_T = Symmetry.nirreps/Symmetry.dcr_deg[stab_ij][stab_kl];
289
 
            ni = (BasisSet.puream ? 2*BasisSet.shells[si].am - 1 : ioff[BasisSet.shells[si].am]);
290
 
            nj = (BasisSet.puream ? 2*BasisSet.shells[sjj].am - 1 : ioff[BasisSet.shells[sjj].am]);
291
 
            nk = (BasisSet.puream ? 2*BasisSet.shells[skk].am - 1 : ioff[BasisSet.shells[skk].am]);
292
 
            nl = (BasisSet.puream ? 2*BasisSet.shells[sll].am - 1 : ioff[BasisSet.shells[sll].am]);
293
 
            
294
 
            memset(sj_arr,0,sizeof(int)*max_num_unique_quartets);
295
 
            memset(sk_arr,0,sizeof(int)*max_num_unique_quartets);
296
 
            memset(sl_arr,0,sizeof(int)*max_num_unique_quartets);
297
 
            memset(sj_fbf_arr,0,sizeof(int)*max_num_unique_quartets);
298
 
            memset(sk_fbf_arr,0,sizeof(int)*max_num_unique_quartets);
299
 
            memset(sl_fbf_arr,0,sizeof(int)*max_num_unique_quartets);
300
 
            count = 0;
301
 
            for(dcr_ij=0;dcr_ij<Symmetry.dcr_dim[stab_i][stab_j];dcr_ij++){
302
 
              R = R_list[dcr_ij];
303
 
              sj = BasisSet.shells[sjj].trans_vec[R]-1;
304
 
              for(dcr_ijkl=0;dcr_ijkl<Symmetry.dcr_dim[stab_ij][stab_kl];dcr_ijkl++){
305
 
                T = T_list[dcr_ijkl];
306
 
                sk = BasisSet.shells[skk].trans_vec[T]-1;
307
 
                slll = BasisSet.shells[sll].trans_vec[T]-1;
308
 
                for(dcr_kl=0;dcr_kl<Symmetry.dcr_dim[stab_k][stab_l];dcr_kl++) {
309
 
                  S = S_list[dcr_kl];
310
 
                  sl = BasisSet.shells[slll].trans_vec[S]-1;
311
 
                  
312
 
                  total_am = BasisSet.shells[si].am +
313
 
                  BasisSet.shells[sj].am +
314
 
                  BasisSet.shells[sk].am +
315
 
                  BasisSet.shells[sl].am;
316
 
                  /*-------------------------------------------------------------
317
 
                  Obviously redundant or zero cases should be eliminated here!
318
 
                  Right now only zero case is eliminated. Redundancies arising
319
 
                  in DCD approach when usi == usj etc. may be eliminated too
320
 
                  but lambda_T will have to be replaced by an array (it won't
321
 
                  the same for every shell quartet in petite list anymore).
322
 
                  -------------------------------------------------------------*/
323
 
                  if(!(total_am%2)||
324
 
                    (BasisSet.shells[si].center!=BasisSet.shells[sj].center)||
325
 
                    (BasisSet.shells[sj].center!=BasisSet.shells[sk].center)||
326
 
                    (BasisSet.shells[sk].center!=BasisSet.shells[sl].center)) {
327
 
                    sj_arr[count] = sj;
328
 
                    sk_arr[count] = sk;
329
 
                    sl_arr[count] = sl;
330
 
                    sj_fbf_arr[count] = BasisSet.shells[sj].fbf-1;
331
 
                    sk_fbf_arr[count] = BasisSet.shells[sk].fbf-1;
332
 
                    sl_fbf_arr[count] = BasisSet.shells[sl].fbf-1;
333
 
                    
334
 
                    count++;
335
 
                  }
336
 
                }
337
 
              }
338
 
            } /* petite list is ready to be used */
339
 
            num_unique_quartets = count;
340
 
          }
341
 
          else { /*--- C1 symmetry case ---*/
342
 
            num_unique_quartets = 1;
343
 
            sj_arr[0] = usj;
344
 
            sk_arr[0] = usk;
345
 
            sl_arr[0] = usl;
346
 
            ni = (BasisSet.puream ? 2*BasisSet.shells[usi].am - 1 : ioff[BasisSet.shells[usi].am]);
347
 
            nj = (BasisSet.puream ? 2*BasisSet.shells[usj].am - 1 : ioff[BasisSet.shells[usj].am]);
348
 
            nk = (BasisSet.puream ? 2*BasisSet.shells[usk].am - 1 : ioff[BasisSet.shells[usk].am]);
349
 
            nl = (BasisSet.puream ? 2*BasisSet.shells[usl].am - 1 : ioff[BasisSet.shells[usl].am]);
350
 
            ioffset = BasisSet.shells[usi].fbf - 1;
351
 
            joffset = BasisSet.shells[usj].fbf - 1;
352
 
            koffset = BasisSet.shells[usk].fbf - 1;
353
 
            loffset = BasisSet.shells[usl].fbf - 1;
354
 
          }
355
 
          class_size = ni*nj*nk*nl;
356
 
          
357
 
          /* Determine irreps of all derivative operators whose application to this SO quartet is not identically zero */
358
 
          /* assume Abelian groups only, i.e. at most 8 irreps */
359
 
          int irreps_of_allowed_derivatives = CDSALCs.atom_irreps[center_i] |
360
 
                                              CDSALCs.atom_irreps[center_j] |
361
 
                                              CDSALCs.atom_irreps[center_k] |
362
 
                                              CDSALCs.atom_irreps[center_l];
363
 
          /* Determine the number of all SALC derivative operators which give nonzero when applied to the SO quartet */
364
 
          nuc[0] = center_i;
365
 
          nuc[1] = center_j;
366
 
          nuc[2] = center_k;
367
 
          nuc[3] = center_l;
368
 
          int salc_count = 0;
369
 
          for(int cd=0; cd<num_coords; cd++)
370
 
            salc_all2thisquartet[cd] = -1;
371
 
          for(int c=0; c<4; c++) {
372
 
            int cart_der = 3*nuc[c];
373
 
            for(int xyz=0; xyz<3; xyz++,cart_der++) {
374
 
              const CDSALC_t::cd2salc_map_t& cd2salc_map = CDSALCs.cd2salc_map[cart_der];
375
 
              const int nsalcs = cd2salc_map.nsalcs;
376
 
              for(int s=0; s<nsalcs; s++) {
377
 
                const int salc = cd2salc_map.salcs[s];
378
 
                if(salc_all2thisquartet[salc] == -1) {
379
 
                  salc_all2thisquartet[salc] = salc_count;
380
 
                  salc_thisquartet2all[salc_count] = salc;
381
 
                  ++salc_count;
382
 
                }
383
 
              }
384
 
            }
385
 
          }
386
 
          const int nsalcders_for_this_quartet = salc_count;
387
 
          
388
 
          np_i = BasisSet.shells[si].n_prims;
389
 
          np_j = BasisSet.shells[sjj].n_prims;
390
 
          np_k = BasisSet.shells[skk].n_prims;
391
 
          np_l = BasisSet.shells[sll].n_prims;
392
 
          
393
 
          orig_am[0] = BasisSet.shells[si].am-1;
394
 
          orig_am[1] = BasisSet.shells[sjj].am-1;
395
 
          orig_am[2] = BasisSet.shells[skk].am-1;
396
 
          orig_am[3] = BasisSet.shells[sll].am-1;
397
 
          am = orig_am[0] + orig_am[1] + orig_am[2] + orig_am[3];
398
 
 
399
 
          /*----------------------------------
400
 
            Compute the nonredundant quartets
401
 
           ----------------------------------*/
402
 
          for(plquartet=0;plquartet<num_unique_quartets;plquartet++) {
403
 
            sj = sj_arr[plquartet];
404
 
            sk = sk_arr[plquartet];
405
 
            sl = sl_arr[plquartet];
406
 
            
407
 
            int center_j = BasisSet.shells[sj].center-1;
408
 
            int center_k = BasisSet.shells[sk].center-1;
409
 
            int center_l = BasisSet.shells[sl].center-1;
410
 
            
411
 
            sp_ij = &(BasisSet.shell_pairs[si][sj]);
412
 
            sp_kl = &(BasisSet.shell_pairs[sk][sl]);
413
 
            
414
 
            Libderiv.AB[0] = sp_ij->AB[0];
415
 
            Libderiv.AB[1] = sp_ij->AB[1];
416
 
            Libderiv.AB[2] = sp_ij->AB[2];
417
 
            Libderiv.CD[0] = sp_kl->AB[0];
418
 
            Libderiv.CD[1] = sp_kl->AB[1];
419
 
            Libderiv.CD[2] = sp_kl->AB[2];
420
 
            
421
 
            AB2 = Libderiv.AB[0]*Libderiv.AB[0]+Libderiv.AB[1]*Libderiv.AB[1]+Libderiv.AB[2]*Libderiv.AB[2];
422
 
            CD2 = Libderiv.CD[0]*Libderiv.CD[0]+Libderiv.CD[1]*Libderiv.CD[1]+Libderiv.CD[2]*Libderiv.CD[2];
423
 
            
424
 
            /*--- Compute data for primitive quartets here ---*/
425
 
            num_prim_comb = 0;
426
 
            const double pfac = lambda_T;
 
43
  typedef struct {
 
44
    double coef;    // contribution of this cartesian derivative to this SALC
 
45
    int cart_der;   // cartesian derivative (for this quartet) [0,3*num_unique_atoms)
 
46
    int salc;       // SALC index [0,num_atoms*3)
 
47
  } cdsalc_elem;
 
48
  
 
49
  
 
50
  void te_deriv1_ints()
 
51
  {
 
52
    const double toler = UserOptions.cutoff;
 
53
    const int num_coords = 3*Molecule.num_atoms;        // total number of coordinates, symmetry-adapted or otherwise
 
54
    // max number of nonzero derivatives for a SO shell quartet = 4 (centers) * 3 (x,y,z) * maximum multiplicity of a nucleus
 
55
    const int max_num_cdsalcs_per_quartet = 12*Symmetry.max_stab_index;
 
56
    
 
57
    /*--- ASCII file to print integrals ---*/
 
58
    FILE *d1eriout;
 
59
    
 
60
    /*--- Various data structures ---*/
 
61
    struct iwlbuf* D1ERIOUT;               /* IWL buffer for target integrals */
 
62
    struct tebuf** tot_data;               /* accum. for contracted integrals */
 
63
    vector<int> num_of_ints_in_totdata;    /* counts how many integrals are in each tot_data */
 
64
    struct shell_pair *sp_ij, *sp_kl;
 
65
    struct unique_shell_pair *usp_ij,*usp_kl;
 
66
    
 
67
    Libderiv_t Libderiv;
 
68
#ifndef USE_TAYLOR_FM
 
69
    double_array_t fjt_table;
 
70
#endif
 
71
    
 
72
    std::vector<int> unique_center;  unique_center.reserve(4);
 
73
    std::vector<int> nuc(4);
 
74
    
 
75
    PSI_INT_LEAST64 total_te_count = 0;
 
76
    vector<PSI_INT_LEAST64> te_count_per_coord(num_coords,0);     // keeps track of how many integrals were written out
 
77
    int ij, kl, ik, jl, ijkl;
 
78
    int ioffset, joffset, koffset, loffset;
 
79
    int count ;
 
80
    int dum;
 
81
    int total_am, am;
 
82
    int orig_am[4];
 
83
    int pkblock_end_index = -1;
 
84
    int i, j, k, l, m, ii, jj, kk, ll;
 
85
    int si, sj, sk, sl ;
 
86
    int sii, sjj, skk, sll , slll;
 
87
    int pi, pj, pk, pl;
 
88
    int max_pj, max_pl;
 
89
    int upk, num_unique_pk;
 
90
    int usi_arr[3], usj_arr[3], usk_arr[3], usl_arr[3];
 
91
    int *sj_arr, *sk_arr, *sl_arr;
 
92
    int *sj_fbf_arr, *sk_fbf_arr, *sl_fbf_arr;
 
93
    int usi,usj,usk,usl;
 
94
    int stab_i,stab_j,stab_k,stab_l,stab_ij,stab_kl;
 
95
    int *R_list, *S_list, *T_list;
 
96
    int R,S,T;
 
97
    int dcr_ij, dcr_kl, dcr_ijkl;
 
98
    int lambda_T = 1;
 
99
    int num_unique_quartets;
 
100
    int plquartet;
 
101
    int max_num_unique_quartets;
 
102
    int max_num_prim_comb;
 
103
    
 
104
    int size, class_size;
 
105
    int max_cart_class_size;
 
106
    
 
107
    int bf_i, bf_j, bf_k, bf_l, so_i, so_j, so_k, so_l, s;
 
108
    int np_i, np_j, np_k, np_l;
 
109
    int ni, nj, nk, nl;
 
110
    
 
111
    int index;
 
112
    int iimax, jjmax, kkmax, llmax;
 
113
    int irrep, npi_ij, npi_kl, npi_ik, npi_jl, ind_offset;
 
114
    
 
115
    int num_prim_comb, p;
 
116
    
 
117
    double AB2, CD2;
 
118
    double pkblock_end_value = 0.0;
 
119
    double temp;
 
120
    
 
121
    vector<double> so_int;
 
122
    {
 
123
      const int max_num_cdsalc_per_quartet = min(max_num_cdsalcs_per_quartet,CDSALCs.nsalcs);
 
124
      so_int.reserve(max_num_cdsalc_per_quartet);
 
125
    }
 
126
    
 
127
    /*---------------
 
128
      Initialization
 
129
      ---------------*/
 
130
#if PRINT
 
131
    eriout = fopen("d1eriout.dat","w");
 
132
#endif
 
133
    D1ERIOUT = new struct iwlbuf[num_coords];
 
134
    for(int c=0; c<num_coords; c++)
 
135
      iwl_buf_init(&D1ERIOUT[c], IOUnits.itapD1ERI_SO+c, toler, 0, 0);
 
136
#ifdef USE_TAYLOR_FM
 
137
    init_Taylor_Fm_Eval(BasisSet.max_am*4-4+DERIV_LVL,UserOptions.cutoff);
 
138
#else
 
139
    init_fjt(BasisSet.max_am*4+DERIV_LVL);
 
140
    init_fjt_table(&fjt_table);
 
141
#endif
 
142
    init_libderiv_base();
 
143
    
 
144
    /*-------------------------
 
145
      Allocate data structures
 
146
      -------------------------*/
 
147
    max_cart_class_size = (ioff[BasisSet.max_am])*
 
148
      (ioff[BasisSet.max_am])*
 
149
      (ioff[BasisSet.max_am])*
 
150
      (ioff[BasisSet.max_am]);
 
151
    max_num_unique_quartets = Symmetry.max_stab_index*
 
152
      Symmetry.max_stab_index*
 
153
      Symmetry.max_stab_index;
 
154
    
 
155
    // Final integrals are stored here
 
156
    tot_data = new struct tebuf*[num_coords];
 
157
    for(int c=0; c<num_coords; c++) {
 
158
      tot_data[c] = new struct tebuf[max_num_unique_quartets*max_cart_class_size];
 
159
    }
 
160
    num_of_ints_in_totdata.resize(num_coords);
 
161
    for(int c=0; c<num_coords; c++)
 
162
      num_of_ints_in_totdata[c] = 0;
 
163
    
 
164
    // These arrays are used to hold cartesian AO derivative integrals
 
165
    double* cart_ints[12];
 
166
    for(int i=0; i<12; i++)
 
167
      cart_ints[i] = new double[max_cart_class_size];
 
168
    // plist_ints holds all symmetry-unique AO shell quartets for a given SO shell quartet
 
169
    // plist_ints[i][j] is the pointer to j-th *total* derivative of i-th AO shell quartet which contributes to the current SO quartet
 
170
    // note that only total derivatives are stored, i.e. partial derivatives produced by libderiv are combined
 
171
    // to yield total derivatives
 
172
    double*** plist_ints = new double**[max_num_unique_quartets];
 
173
    for(int i=0; i<max_num_unique_quartets; i++) {
 
174
      plist_ints[i] = new double*[12];
 
175
      for(int j=0; j<12; j++) {
 
176
        plist_ints[i][j] = new double[max_cart_class_size];
 
177
      }
 
178
    }
 
179
    // For every SO shell quartet we need the following information:
 
180
    // 1) irreps of all total derivative operators which give nonzero when applied to the SO quartet
 
181
    // 2) given an irrep, which total derivative operators contribute
 
182
    // For every AO shell quartet in a petit list for a given SO shell quartet we need the following information:
 
183
    // 1) number of nonzero total derivatives (same as number of distinct centers, n, times 3, 1 <= n <= 4)
 
184
    // 2) irreps of all total derivative operators which give nonzero when applied to the AO quartet
 
185
    // 3) given an irrep, which total derivative operators contribute
 
186
    
 
187
    // For a set of SALCs, we need to be able to select coefficients of all nonzero derivatives
 
188
    // of a given shell quartet. It is best done when the relevant subset of the (sparse) SALC coefficient
 
189
    // matrix is pre-arranged in an array
 
190
    // plist_salcs[s][g][i] is an i-th nonzero contribution to a SALC of irrep g
 
191
    // which results in nonzero when applied to member s of a petite list
 
192
    struct cdsalc_elem_vec {
 
193
      int nelems;
 
194
      cdsalc_elem* elems;
 
195
    };
 
196
 
 
197
    cdsalc_elem_vec** plist_salcs = new cdsalc_elem_vec*[max_num_unique_quartets];
 
198
    for(int i=0; i<max_num_unique_quartets; i++) {
 
199
      plist_salcs[i] = new cdsalc_elem_vec[Symmetry.nirreps];
 
200
      for(int j=0; j<Symmetry.nirreps; j++) {
 
201
        plist_salcs[i][j].nelems=0;
 
202
        plist_salcs[i][j].elems = new cdsalc_elem[max_num_cdsalcs_per_quartet];
 
203
      }
 
204
    }
 
205
    
 
206
    int* salc_all2thisquartet = new int[num_coords];    // maps absolute SALC index to the SALC index for this quartet
 
207
    int* salc_thisquartet2all = new int[max_num_cdsalcs_per_quartet];    // the reverse of the above
 
208
    
 
209
    sj_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
210
    sk_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
211
    sl_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
212
    if (Symmetry.nirreps > 1) {
 
213
      sj_fbf_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
214
      sk_fbf_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
215
      sl_fbf_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
 
216
    }
 
217
    
 
218
    max_num_prim_comb = (BasisSet.max_num_prims*
 
219
                         BasisSet.max_num_prims)*
 
220
      (BasisSet.max_num_prims*
 
221
       BasisSet.max_num_prims);
 
222
    
 
223
    init_libderiv1(&Libderiv, BasisSet.max_am-1, max_num_prim_comb, max_cart_class_size);
 
224
    
 
225
    /*------------------------------------
 
226
      generate all unique shell quartets
 
227
      ------------------------------------*/
 
228
    if(UserOptions.print_lvl >= PRINT_DEBUG) {
 
229
      fprintf(outfile,"  -Electron repulsion integrals:\n\n");
 
230
    }
 
231
    
 
232
    for (int usii=0; usii<Symmetry.num_unique_shells; usii++)
 
233
      for (int usjj=0; usjj<=usii; usjj++)
 
234
        for (int uskk=0; uskk<Symmetry.num_unique_shells; uskk++) {
 
235
          const int usll_max = (usii==uskk ? usjj : uskk);
 
236
          for (int usll=0; usll<=usll_max; usll++){
 
237
            
 
238
            int usi = usii;
 
239
            int usj = usjj;
 
240
            int usk = uskk;
 
241
            int usl = usll;
 
242
            
 
243
            /* place in "ascending" angular mom-
 
244
               my simple way of optimizing PHG recursion (VRR) */
 
245
            /* these first two are good for the HRR */
 
246
            if(BasisSet.shells[Symmetry.us2s[usi]].am < BasisSet.shells[Symmetry.us2s[usj]].am){
 
247
              dum = usi;
 
248
              usi = usj;
 
249
              usj = dum;
 
250
            }
 
251
            if(BasisSet.shells[Symmetry.us2s[usk]].am < BasisSet.shells[Symmetry.us2s[usl]].am){
 
252
              dum = usk;
 
253
              usk = usl;
 
254
              usl = dum;
 
255
            }
 
256
            /* this should be /good/ for the VRR */
 
257
            if(BasisSet.shells[Symmetry.us2s[usi]].am + BasisSet.shells[Symmetry.us2s[usj]].am >
 
258
               BasisSet.shells[Symmetry.us2s[usk]].am + BasisSet.shells[Symmetry.us2s[usl]].am){
 
259
              dum = usi;
 
260
              usi = usk;
 
261
              usk = dum;
 
262
              dum = usj;
 
263
              usj = usl;
 
264
              usl = dum;
 
265
            }
 
266
            
 
267
            si = Symmetry.us2s[usi];
 
268
            sjj = Symmetry.us2s[usj];
 
269
            skk = Symmetry.us2s[usk];
 
270
            sll = Symmetry.us2s[usl];
 
271
            
 
272
            int center_i = BasisSet.shells[si].center-1;
 
273
            int center_j = BasisSet.shells[sjj].center-1;
 
274
            int center_k = BasisSet.shells[skk].center-1;
 
275
            int center_l = BasisSet.shells[sll].center-1;
 
276
            
 
277
            if (Symmetry.nirreps > 1) { /*--- Non-C1 symmetry case ---*/
 
278
              /*--- Generate the petite list of shell quadruplets using DCD approach of Davidson ---*/
 
279
              usp_ij = &(Symmetry.us_pairs[usi][usj]);
 
280
              usp_kl = &(Symmetry.us_pairs[usk][usl]);
 
281
              stab_i = Symmetry.atom_positions[center_i];
 
282
              stab_j = Symmetry.atom_positions[center_j];
 
283
              stab_k = Symmetry.atom_positions[center_k];
 
284
              stab_l = Symmetry.atom_positions[center_l];
 
285
              stab_ij = Symmetry.GnG[stab_i][stab_j];
 
286
              stab_kl = Symmetry.GnG[stab_k][stab_l];
 
287
              R_list = Symmetry.dcr[stab_i][stab_j];
 
288
              S_list = Symmetry.dcr[stab_k][stab_l];
 
289
              T_list = Symmetry.dcr[stab_ij][stab_kl];
 
290
              lambda_T = Symmetry.nirreps/Symmetry.dcr_deg[stab_ij][stab_kl];
 
291
              ni = (BasisSet.puream ? 2*BasisSet.shells[si].am - 1 : ioff[BasisSet.shells[si].am]);
 
292
              nj = (BasisSet.puream ? 2*BasisSet.shells[sjj].am - 1 : ioff[BasisSet.shells[sjj].am]);
 
293
              nk = (BasisSet.puream ? 2*BasisSet.shells[skk].am - 1 : ioff[BasisSet.shells[skk].am]);
 
294
              nl = (BasisSet.puream ? 2*BasisSet.shells[sll].am - 1 : ioff[BasisSet.shells[sll].am]);
 
295
              
 
296
              memset(sj_arr,0,sizeof(int)*max_num_unique_quartets);
 
297
              memset(sk_arr,0,sizeof(int)*max_num_unique_quartets);
 
298
              memset(sl_arr,0,sizeof(int)*max_num_unique_quartets);
 
299
              memset(sj_fbf_arr,0,sizeof(int)*max_num_unique_quartets);
 
300
              memset(sk_fbf_arr,0,sizeof(int)*max_num_unique_quartets);
 
301
              memset(sl_fbf_arr,0,sizeof(int)*max_num_unique_quartets);
 
302
              count = 0;
 
303
              for(dcr_ij=0;dcr_ij<Symmetry.dcr_dim[stab_i][stab_j];dcr_ij++){
 
304
                R = R_list[dcr_ij];
 
305
                sj = BasisSet.shells[sjj].trans_vec[R]-1;
 
306
                for(dcr_ijkl=0;dcr_ijkl<Symmetry.dcr_dim[stab_ij][stab_kl];dcr_ijkl++){
 
307
                  T = T_list[dcr_ijkl];
 
308
                  sk = BasisSet.shells[skk].trans_vec[T]-1;
 
309
                  slll = BasisSet.shells[sll].trans_vec[T]-1;
 
310
                  for(dcr_kl=0;dcr_kl<Symmetry.dcr_dim[stab_k][stab_l];dcr_kl++) {
 
311
                    S = S_list[dcr_kl];
 
312
                    sl = BasisSet.shells[slll].trans_vec[S]-1;
 
313
                    
 
314
                    total_am = BasisSet.shells[si].am +
 
315
                      BasisSet.shells[sj].am +
 
316
                      BasisSet.shells[sk].am +
 
317
                      BasisSet.shells[sl].am;
 
318
                    /*-------------------------------------------------------------
 
319
                      Obviously redundant or zero cases should be eliminated here!
 
320
                      Right now only zero case is eliminated. Redundancies arising
 
321
                      in DCD approach when usi == usj etc. may be eliminated too
 
322
                      but lambda_T will have to be replaced by an array (it won't
 
323
                      the same for every shell quartet in petite list anymore).
 
324
                      -------------------------------------------------------------*/
 
325
                    if(!(total_am%2)||
 
326
                       (BasisSet.shells[si].center!=BasisSet.shells[sj].center)||
 
327
                       (BasisSet.shells[sj].center!=BasisSet.shells[sk].center)||
 
328
                       (BasisSet.shells[sk].center!=BasisSet.shells[sl].center)) {
 
329
                      sj_arr[count] = sj;
 
330
                      sk_arr[count] = sk;
 
331
                      sl_arr[count] = sl;
 
332
                      sj_fbf_arr[count] = BasisSet.shells[sj].fbf-1;
 
333
                      sk_fbf_arr[count] = BasisSet.shells[sk].fbf-1;
 
334
                      sl_fbf_arr[count] = BasisSet.shells[sl].fbf-1;
 
335
                      
 
336
                      count++;
 
337
                    }
 
338
                  }
 
339
                }
 
340
              } /* petite list is ready to be used */
 
341
              num_unique_quartets = count;
 
342
            }
 
343
            else { /*--- C1 symmetry case ---*/
 
344
              num_unique_quartets = 1;
 
345
              sj_arr[0] = usj;
 
346
              sk_arr[0] = usk;
 
347
              sl_arr[0] = usl;
 
348
              ni = (BasisSet.puream ? 2*BasisSet.shells[usi].am - 1 : ioff[BasisSet.shells[usi].am]);
 
349
              nj = (BasisSet.puream ? 2*BasisSet.shells[usj].am - 1 : ioff[BasisSet.shells[usj].am]);
 
350
              nk = (BasisSet.puream ? 2*BasisSet.shells[usk].am - 1 : ioff[BasisSet.shells[usk].am]);
 
351
              nl = (BasisSet.puream ? 2*BasisSet.shells[usl].am - 1 : ioff[BasisSet.shells[usl].am]);
 
352
              ioffset = BasisSet.shells[usi].fbf - 1;
 
353
              joffset = BasisSet.shells[usj].fbf - 1;
 
354
              koffset = BasisSet.shells[usk].fbf - 1;
 
355
              loffset = BasisSet.shells[usl].fbf - 1;
 
356
            }
 
357
            class_size = ni*nj*nk*nl;
 
358
            
 
359
            /* Determine irreps of all derivative operators whose application to this SO quartet is not identically zero */
 
360
            /* assume Abelian groups only, i.e. at most 8 irreps */
 
361
            int irreps_of_allowed_derivatives = CDSALCs.atom_irreps[center_i] |
 
362
              CDSALCs.atom_irreps[center_j] |
 
363
              CDSALCs.atom_irreps[center_k] |
 
364
              CDSALCs.atom_irreps[center_l];
 
365
            /* Determine the number of all SALC derivative operators which give nonzero when applied to the SO quartet */
 
366
            nuc[0] = center_i;
 
367
            nuc[1] = center_j;
 
368
            nuc[2] = center_k;
 
369
            nuc[3] = center_l;
 
370
            int salc_count = 0;
 
371
            for(int cd=0; cd<num_coords; cd++)
 
372
              salc_all2thisquartet[cd] = -1;
 
373
            for(int c=0; c<4; c++) {
 
374
              int cart_der = 3*nuc[c];
 
375
              for(int xyz=0; xyz<3; xyz++,cart_der++) {
 
376
                const CDSALC_t::cd2salc_map_t& cd2salc_map = CDSALCs.cd2salc_map[cart_der];
 
377
                const int nsalcs = cd2salc_map.nsalcs;
 
378
                for(int s=0; s<nsalcs; s++) {
 
379
                  const int salc = cd2salc_map.salcs[s];
 
380
                  if(salc_all2thisquartet[salc] == -1) {
 
381
                    salc_all2thisquartet[salc] = salc_count;
 
382
                    salc_thisquartet2all[salc_count] = salc;
 
383
                    ++salc_count;
 
384
                  }
 
385
                }
 
386
              }
 
387
            }
 
388
            const int nsalcders_for_this_quartet = salc_count;
 
389
          
 
390
            np_i = BasisSet.shells[si].n_prims;
 
391
            np_j = BasisSet.shells[sjj].n_prims;
 
392
            np_k = BasisSet.shells[skk].n_prims;
 
393
            np_l = BasisSet.shells[sll].n_prims;
 
394
            
 
395
            orig_am[0] = BasisSet.shells[si].am-1;
 
396
            orig_am[1] = BasisSet.shells[sjj].am-1;
 
397
            orig_am[2] = BasisSet.shells[skk].am-1;
 
398
            orig_am[3] = BasisSet.shells[sll].am-1;
 
399
            am = orig_am[0] + orig_am[1] + orig_am[2] + orig_am[3];
 
400
            
 
401
            /*----------------------------------
 
402
              Compute the nonredundant quartets
 
403
              ----------------------------------*/
 
404
            for(plquartet=0;plquartet<num_unique_quartets;plquartet++) {
 
405
              sj = sj_arr[plquartet];
 
406
              sk = sk_arr[plquartet];
 
407
              sl = sl_arr[plquartet];
 
408
              
 
409
              int center_j = BasisSet.shells[sj].center-1;
 
410
              int center_k = BasisSet.shells[sk].center-1;
 
411
              int center_l = BasisSet.shells[sl].center-1;
 
412
              
 
413
              sp_ij = &(BasisSet.shell_pairs[si][sj]);
 
414
              sp_kl = &(BasisSet.shell_pairs[sk][sl]);
 
415
              
 
416
              Libderiv.AB[0] = sp_ij->AB[0];
 
417
              Libderiv.AB[1] = sp_ij->AB[1];
 
418
              Libderiv.AB[2] = sp_ij->AB[2];
 
419
              Libderiv.CD[0] = sp_kl->AB[0];
 
420
              Libderiv.CD[1] = sp_kl->AB[1];
 
421
              Libderiv.CD[2] = sp_kl->AB[2];
 
422
            
 
423
              AB2 = Libderiv.AB[0]*Libderiv.AB[0]+Libderiv.AB[1]*Libderiv.AB[1]+Libderiv.AB[2]*Libderiv.AB[2];
 
424
              CD2 = Libderiv.CD[0]*Libderiv.CD[0]+Libderiv.CD[1]*Libderiv.CD[1]+Libderiv.CD[2]*Libderiv.CD[2];
 
425
              
 
426
              /*--- Compute data for primitive quartets here ---*/
 
427
              num_prim_comb = 0;
 
428
              const double pfac = lambda_T;
427
429
#define NO_PERM_SYMM 0
428
430
#if NO_PERM_SYMM
429
 
            for (pi = 0; pi < np_i; pi++) {
430
 
              for (pj = 0; pj < np_j; pj++) {
431
 
                for (pk = 0; pk < np_k; pk++) {
432
 
                  for (pl = 0; pl < np_l; pl++){
433
 
                    const int n = 1;
 
431
              for (pi = 0; pi < np_i; pi++) {
 
432
                for (pj = 0; pj < np_j; pj++) {
 
433
                  for (pk = 0; pk < np_k; pk++) {
 
434
                    for (pl = 0; pl < np_l; pl++){
 
435
                      const int n = 1;
434
436
#else
435
 
            for (pi = 0; pi < np_i; pi++) {
436
 
              max_pj = (si == sj) ? pi+1 : np_j;
437
 
              for (pj = 0; pj < max_pj; pj++) {
438
 
                m = (1 + (si == sj && pi != pj));
439
 
                for (pk = 0; pk < np_k; pk++) {
440
 
                  max_pl = (sk == sl) ? pk+1 : np_l;
441
 
                  for (pl = 0; pl < max_pl; pl++){
442
 
                    const int n = m * (1 + (sk == sl && pk != pl));
 
437
                      for (pi = 0; pi < np_i; pi++) {
 
438
                        max_pj = (si == sj) ? pi+1 : np_j;
 
439
                        for (pj = 0; pj < max_pj; pj++) {
 
440
                          m = (1 + (si == sj && pi != pj));
 
441
                          for (pk = 0; pk < np_k; pk++) {
 
442
                            max_pl = (sk == sl) ? pk+1 : np_l;
 
443
                            for (pl = 0; pl < max_pl; pl++){
 
444
                              const int n = m * (1 + (sk == sl && pk != pl));
443
445
#endif
444
446
#ifdef USE_TAYLOR_FM
445
447
                    deriv1_quartet_data(&(Libderiv.PrimQuartet[num_prim_comb++]),
778
780
      free_fjt();
779
781
      #endif
780
782
 
781
 
  return;
782
 
}
783
 
 
784
 
};
785
 
 
 
783
      return;
 
784
            }
 
785
 
 
786
};
 
787
};