← 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_print.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cints/Default_Deriv1/te_deriv1_print.c
1
 
#include <math.h>
2
 
#include <stdio.h>
3
 
#include <string.h>
4
 
#include <memory.h>
5
 
#include <stdlib.h>
6
 
#include <psitypes.h>
7
 
#include <libciomr/libciomr.h>
8
 
#include <libint/libint.h>
9
 
#include <libderiv/libderiv.h>
10
 
#include "defines.h"
11
 
#define EXTERN
12
 
#include "global.h"
13
 
#ifdef USE_TAYLOR_FM
14
 
  #include"taylor_fm_eval.h"
15
 
#else
16
 
  #include"int_fjt.h"
17
 
#endif
18
 
#include "deriv1_quartet_data.h"
19
 
#include "small_fns.h"
20
 
 
21
 
void te_deriv1_print(void)
22
 
{
23
 
  Libderiv_t Libderiv;                    /* Integrals library object */
24
 
#ifndef USE_TAYLOR_FM
25
 
  double_array_t fjt_table;               /* table of auxiliary function F_m(u) for each primitive combination */
26
 
#endif
27
 
 
28
 
  int max_class_size;
29
 
  int max_cart_class_size;
30
 
  int max_num_prim_comb, num_prim_comb;
31
 
  int max_num_unique_quartets;
32
 
  int *sj_arr, *sk_arr, *sl_arr;
33
 
  int usii, usjj, uskk, usll;
34
 
  PSI_INT_LEAST64 quartet_index;
35
 
  int sii, sjj, skk, sll;
36
 
  int num_unique_quartets;
37
 
  int ni, nj, nk, nl;
38
 
  int si, sj, sk, sl;
39
 
  int plquartet, dum;
40
 
  int switch_ij, switch_kl, switch_ijkl;
41
 
  int np_i, np_j, np_k, np_l;
42
 
  int quartet_size;
43
 
  int orig_am[4], am;
44
 
  struct shell_pair *sp_ij, *sp_kl;
45
 
  double AB2, CD2;
46
 
  int pi, pj, pk, pl;
47
 
  int si_fao, sj_fao, sk_fao, sl_fao;
48
 
  int ao_i, ao_j, ao_k, ao_l;
49
 
  int count;
50
 
  int center_i, center_j, center_k, center_l;
51
 
  int iout, jout, kout, lout;
52
 
  int i, k, ij, kl;
53
 
 
54
 
  double *PrintFourInd;
55
 
  double **derivs;
56
 
  FILE **out;
57
 
  char lbl[20];
58
 
 
59
 
  /* some error checking up front */
60
 
  if(Symmetry.nirreps > 1 || BasisSet.puream) {
61
 
    fprintf(outfile, "\nThe derivative integral printing code will almost certainly fail\n");
62
 
    fprintf(outfile, "for non-C1 and/or puream cases.  Exiting.\n");
63
 
    punt("no symmetry or puream cases allowed when printing deriv. ints.");
64
 
  }
65
 
 
66
 
  /*---------------
67
 
    Initialization
68
 
    ---------------*/
69
 
#ifndef USE_TAYLOR_FM
70
 
  init_fjt_table(&fjt_table);
71
 
#endif
72
 
 
73
 
  max_cart_class_size = ioff[BasisSet.max_am]*ioff[BasisSet.max_am]*ioff[BasisSet.max_am]*ioff[BasisSet.max_am];
74
 
  max_class_size = max_cart_class_size;
75
 
  max_num_prim_comb = (BasisSet.max_num_prims*BasisSet.max_num_prims)*
76
 
    (BasisSet.max_num_prims*BasisSet.max_num_prims);
77
 
  init_libderiv1(&Libderiv,BasisSet.max_am-1,max_num_prim_comb,max_class_size);
78
 
  PrintFourInd = init_array(max_cart_class_size);
79
 
  derivs = block_matrix(Molecule.num_atoms*3,max_cart_class_size);
80
 
  out = (FILE **) malloc(Molecule.num_atoms * 3 * sizeof(FILE *));
81
 
 
82
 
  for(i=0; i < Molecule.num_atoms*3; i++) {
83
 
    sprintf(lbl,"eri%d.dat", i);
84
 
    ffile(&out[i], lbl, 0);
85
 
  }
86
 
 
87
 
  max_num_unique_quartets = Symmetry.max_stab_index*
88
 
    Symmetry.max_stab_index*
89
 
    Symmetry.max_stab_index;
90
 
  sj_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
91
 
  sk_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
92
 
  sl_arr = (int *)malloc(sizeof(int)*max_num_unique_quartets);
93
 
 
94
 
  for (usii=0; usii<Symmetry.num_unique_shells; usii++)
95
 
    for (usjj=0; usjj<Symmetry.num_unique_shells; usjj++)
96
 
      for (uskk=0; uskk<Symmetry.num_unique_shells; uskk++)
97
 
        for (usll=0; usll<Symmetry.num_unique_shells; usll++, quartet_index++){
98
 
 
99
 
          sii = Symmetry.us2s[usii];
100
 
          sjj = Symmetry.us2s[usjj];
101
 
          skk = Symmetry.us2s[uskk];
102
 
          sll = Symmetry.us2s[usll];
103
 
 
104
 
          num_unique_quartets = 1;
105
 
          sj_arr[0] = usjj;
106
 
          sk_arr[0] = uskk;
107
 
          sl_arr[0] = usll;
108
 
 
109
 
          /*----------------------------------
110
 
            Compute the nonredundant quartets
111
 
            ----------------------------------*/
112
 
          for(plquartet=0;plquartet<num_unique_quartets;plquartet++) {
113
 
            si = Symmetry.us2s[usii];
114
 
            sj = sj_arr[plquartet];
115
 
            sk = sk_arr[plquartet];
116
 
            sl = sl_arr[plquartet];
117
 
 
118
 
            /*--- Skip this quartet if all four centers are the same ---*/
119
 
            if (BasisSet.shells[si].center == BasisSet.shells[sj].center &&
120
 
                BasisSet.shells[si].center == BasisSet.shells[sk].center &&
121
 
                BasisSet.shells[si].center == BasisSet.shells[sl].center &&
122
 
                DERIV_LVL == 1)
123
 
              continue;
124
 
 
125
 
            switch_ij = 0;
126
 
            switch_kl = 0;
127
 
            switch_ijkl = 0;
128
 
            /* place in "ascending" angular mom-
129
 
               my simple way of optimizing PHG recursion (VRR) */
130
 
            /* these first two are good for the HRR */
131
 
            if(BasisSet.shells[si].am < BasisSet.shells[sj].am){
132
 
              dum = si;
133
 
              si = sj;
134
 
              sj = dum;
135
 
              switch_ij = 1;
136
 
            }
137
 
            if(BasisSet.shells[sk].am < BasisSet.shells[sl].am){
138
 
              dum = sk;
139
 
              sk = sl;
140
 
              sl = dum;
141
 
              switch_kl = 1;
142
 
            }
143
 
            /* this should be /good/ for the VRR */
144
 
            if(BasisSet.shells[si].am + BasisSet.shells[sj].am > BasisSet.shells[sk].am + BasisSet.shells[sl].am){
145
 
              dum = si;
146
 
              si = sk;
147
 
              sk = dum;
148
 
              dum = sj;
149
 
              sj = sl;
150
 
              sl = dum;
151
 
              switch_ijkl = 1;
152
 
            }
153
 
 
154
 
            ni = ioff[BasisSet.shells[si].am];
155
 
            nj = ioff[BasisSet.shells[sj].am];
156
 
            nk = ioff[BasisSet.shells[sk].am];
157
 
            nl = ioff[BasisSet.shells[sl].am];
158
 
            quartet_size = ni*nj*nk*nl;
159
 
 
160
 
            np_i = BasisSet.shells[si].n_prims;
161
 
            np_j = BasisSet.shells[sj].n_prims;
162
 
            np_k = BasisSet.shells[sk].n_prims;
163
 
            np_l = BasisSet.shells[sl].n_prims;
164
 
 
165
 
            orig_am[0] = BasisSet.shells[si].am-1;
166
 
            orig_am[1] = BasisSet.shells[sj].am-1;
167
 
            orig_am[2] = BasisSet.shells[sk].am-1;
168
 
            orig_am[3] = BasisSet.shells[sl].am-1;
169
 
            am = orig_am[0] + orig_am[1] + orig_am[2] + orig_am[3];
170
 
 
171
 
            sp_ij = &(BasisSet.shell_pairs[si][sj]);
172
 
            sp_kl = &(BasisSet.shell_pairs[sk][sl]);
173
 
 
174
 
            Libderiv.AB[0] = sp_ij->AB[0];
175
 
            Libderiv.AB[1] = sp_ij->AB[1];
176
 
            Libderiv.AB[2] = sp_ij->AB[2];
177
 
            Libderiv.CD[0] = sp_kl->AB[0];
178
 
            Libderiv.CD[1] = sp_kl->AB[1];
179
 
            Libderiv.CD[2] = sp_kl->AB[2];
180
 
            
181
 
            AB2 = Libderiv.AB[0]*Libderiv.AB[0]+
182
 
              Libderiv.AB[1]*Libderiv.AB[1]+
183
 
              Libderiv.AB[2]*Libderiv.AB[2];
184
 
            CD2 = Libderiv.CD[0]*Libderiv.CD[0]+
185
 
              Libderiv.CD[1]*Libderiv.CD[1]+
186
 
              Libderiv.CD[2]*Libderiv.CD[2];
187
 
 
188
 
            num_prim_comb = 0;
189
 
            for (pi = 0; pi < np_i; pi++) {
190
 
              for (pj = 0; pj < np_j; pj++) {
191
 
                for (pk = 0; pk < np_k; pk++) {
192
 
                  for (pl = 0; pl < np_l; pl++){
193
 
#ifdef USE_TAYLOR_FM
194
 
                    deriv1_quartet_data(&(Libderiv.PrimQuartet[num_prim_comb++]),
195
 
                                        NULL, AB2, CD2,
196
 
                                        sp_ij, sp_kl, am, pi, pj, pk, pl, 1);
197
 
#else
198
 
                    deriv1_quartet_data(&(Libderiv.PrimQuartet[num_prim_comb++]),
199
 
                                        &fjt_table, AB2, CD2,
200
 
                                        sp_ij, sp_kl, am, pi, pj, pk, pl, 1);
201
 
#endif              
202
 
                  }
203
 
                }
204
 
              }
205
 
            }
206
 
 
207
 
            si_fao = BasisSet.shells[si].fao-1;
208
 
            sj_fao = BasisSet.shells[sj].fao-1;
209
 
            sk_fao = BasisSet.shells[sk].fao-1;
210
 
            sl_fao = BasisSet.shells[sl].fao-1;
211
 
            count = 0;
212
 
            for (ao_i = si_fao; ao_i < si_fao+ni; ao_i++)
213
 
              for (ao_j = sj_fao; ao_j < sj_fao+nj; ao_j++)
214
 
                for (ao_k = sk_fao; ao_k < sk_fao+nk; ao_k++)
215
 
                  for (ao_l = sl_fao; ao_l < sl_fao+nl; ao_l++) {
216
 
                    PrintFourInd[count] = GTOs.bf_norm[orig_am[0]][ao_i-si_fao] *
217
 
                      GTOs.bf_norm[orig_am[1]][ao_j-sj_fao] *
218
 
                      GTOs.bf_norm[orig_am[2]][ao_k-sk_fao] *
219
 
                      GTOs.bf_norm[orig_am[3]][ao_l-sl_fao];
220
 
                    count++;
221
 
                  }
222
 
 
223
 
            build_deriv1_eri[orig_am[0]][orig_am[1]][orig_am[2]][orig_am[3]](&Libderiv,num_prim_comb);
224
 
 
225
 
            center_i = BasisSet.shells[si].center-1;
226
 
            center_j = BasisSet.shells[sj].center-1;
227
 
            center_k = BasisSet.shells[sk].center-1;
228
 
            center_l = BasisSet.shells[sl].center-1;
229
 
 
230
 
            zero_mat(derivs, Molecule.num_atoms*3, max_cart_class_size);
231
 
            for(k=0; k < quartet_size; k++) {
232
 
              derivs[center_i*3+0][k] += Libderiv.ABCD[0][k] * PrintFourInd[k];
233
 
              derivs[center_i*3+1][k] += Libderiv.ABCD[1][k] * PrintFourInd[k];
234
 
              derivs[center_i*3+2][k] += Libderiv.ABCD[2][k] * PrintFourInd[k];
235
 
 
236
 
              derivs[center_j*3+0][k] -= (Libderiv.ABCD[0][k] + Libderiv.ABCD[6][k] + Libderiv.ABCD[9][k]) * PrintFourInd[k];
237
 
              derivs[center_j*3+1][k] -= (Libderiv.ABCD[1][k] + Libderiv.ABCD[7][k] + Libderiv.ABCD[10][k]) * PrintFourInd[k];
238
 
              derivs[center_j*3+2][k] -= (Libderiv.ABCD[2][k] + Libderiv.ABCD[8][k] + Libderiv.ABCD[11][k]) * PrintFourInd[k];
239
 
 
240
 
              derivs[center_k*3+0][k] += Libderiv.ABCD[6][k] * PrintFourInd[k];
241
 
              derivs[center_k*3+1][k] += Libderiv.ABCD[7][k] * PrintFourInd[k];
242
 
              derivs[center_k*3+2][k] += Libderiv.ABCD[8][k] * PrintFourInd[k];
243
 
 
244
 
              derivs[center_l*3+0][k] += Libderiv.ABCD[9][k] * PrintFourInd[k];
245
 
              derivs[center_l*3+1][k] += Libderiv.ABCD[10][k] * PrintFourInd[k];
246
 
              derivs[center_l*3+2][k] += Libderiv.ABCD[11][k] * PrintFourInd[k];
247
 
            }
248
 
 
249
 
            for(k=0; k < Molecule.num_atoms*3; k++) {
250
 
              count = 0;
251
 
              for (ao_i = si_fao; ao_i < si_fao+ni; ao_i++)
252
 
                for (ao_j = sj_fao; ao_j < sj_fao+nj; ao_j++)
253
 
                  for (ao_k = sk_fao; ao_k < sk_fao+nk; ao_k++)
254
 
                    for (ao_l = sl_fao; ao_l < sl_fao+nl; ao_l++) {
255
 
                      if(fabs(derivs[k][count]) > 1e-6) {
256
 
                        iout = ao_i;  jout = ao_j; kout = ao_k; lout = ao_l;
257
 
                        if(switch_ijkl) { dum = iout; iout = kout; kout = dum; dum = jout; jout = lout; lout = dum; }
258
 
                        if(switch_kl) { dum = kout; kout = lout; lout = dum; }
259
 
                        if(switch_ij) { dum = iout; iout = jout; jout = dum; }
260
 
                        ij = INDEX(iout,jout); kl = INDEX(kout,lout);
261
 
                        if(iout >= jout && kout >= lout && ij >= kl)
262
 
                          fprintf(out[k], "%3d %3d %3d %3d %20.14f\n", iout, jout, kout, lout, derivs[k][count]);
263
 
                      }
264
 
 
265
 
                      count++;
266
 
                    }
267
 
            }
268
 
 
269
 
          } /* end of loop over symmetry unique quartets */
270
 
 
271
 
        } /* end of unique shell loops */
272
 
 
273
 
  free(sj_arr);
274
 
  free(sk_arr);
275
 
  free(sl_arr);
276
 
  free(PrintFourInd);
277
 
  for(i=0; i < Molecule.num_atoms*3; i++) fclose(out[i]);
278
 
  free(out);
279
 
  free_block(derivs);
280
 
  free_libderiv(&Libderiv);
281
 
#ifndef USE_TAYLOR_FM
282
 
  free_fjt_table(&fjt_table);
283
 
#endif
284
 
 
285
 
  return;
286
 
}