← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/Tools/compute_eri.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/Tools/compute_eri.cc
1
 
 
2
 
extern "C" {
3
 
#include <stdio.h>
4
 
#include <string.h>
5
 
#include <stdlib.h>
 
1
/*! \file compute_eri.cc
 
2
    \ingroup CINTS
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
 
 
6
#include <cstdio>
 
7
#include <cstring>
 
8
#include <cstdlib>
 
9
 
6
10
#include <libint/libint.h>
7
11
#include <libqt/qt.h>
8
12
#include "defines.h"
14
18
  #include"taylor_fm_eval.h"
15
19
#else
16
20
  #include"int_fjt.h"
 
21
  #include"fjt.h"
17
22
#endif
18
23
#include "quartet_data.h"
19
24
#include "symmetrize.h"
20
25
#include "small_fns.h"
21
 
}
22
 
 
23
 
static void inline switch_ij(int& a, int& b) {int dum = a; a = b; b = dum;};
24
 
 
25
 
extern "C" int compute_eri(double* target,
26
 
                           Libint_t* Libint, int& si, int& sj, int& sk, int& sl,
27
 
                           int& inc1, int& inc2, int& inc3, int& inc4, const bool do_not_permute)
28
 
{
29
 
 
 
26
 
 
27
namespace psi { 
 
28
  namespace CINTS {
 
29
    
 
30
    static void inline switch_ij(int& a, int& b) {int dum = a; a = b; b = dum;};
 
31
    
 
32
    int compute_eri(double* target,
 
33
                    Libint_t* Libint, int& si, int& sj, int& sk, int& sl,
 
34
                    int& inc1, int& inc2, int& inc3, int& inc4, const bool do_not_permute)
 
35
    {
 
36
      
30
37
#ifndef USE_TAYLOR_FM
31
 
  static double_array_t fjt_table;
32
 
  init_fjt_table(&fjt_table);
 
38
      static double_array_t fjt_table;
 
39
      init_fjt_table(&fjt_table);
33
40
#endif
34
 
 
35
 
  /* place in "ascending" angular mom-
36
 
     my simple way of optimizing PHG recursion (VRR) */
37
 
 
38
 
  bool switch_bra_ket = false;
39
 
  /* this should be /good/ for the VRR */
40
 
  if ( BasisSet.shells[si].am + BasisSet.shells[sj].am + inc1 + inc2 >
41
 
       BasisSet.shells[sk].am + BasisSet.shells[sl].am + inc3 + inc4 ) {
42
 
    switch_bra_ket = true;
43
 
    switch_ij(si,sk);
44
 
    switch_ij(sj,sl);
45
 
    switch_ij(inc1, inc3);
46
 
    switch_ij(inc2, inc4);
47
 
  }
48
 
 
49
 
  /* these two are good for the HRR */
50
 
  bool switch_bra = false;
51
 
  if(BasisSet.shells[si].am + inc1 < BasisSet.shells[sj].am + inc2){
52
 
    switch_bra = true;
53
 
    switch_ij(si,sj);
54
 
    switch_ij(inc1,inc2);
55
 
  }
56
 
  bool switch_ket = false;
57
 
  if(BasisSet.shells[sk].am + inc3 < BasisSet.shells[sl].am + inc4){
58
 
    switch_ket = true;
59
 
    switch_ij(sk,sl);
60
 
    switch_ij(inc3,inc4);
61
 
  }
62
 
 
63
 
  int am1 = BasisSet.shells[si].am - 1 + inc1;
64
 
  int am2 = BasisSet.shells[sj].am - 1 + inc2;
65
 
  int am3 = BasisSet.shells[sk].am - 1 + inc3;
66
 
  int am4 = BasisSet.shells[sl].am - 1 + inc4;
67
 
  int am = am1 + am2 + am3 + am4;
68
 
  
69
 
  int c1 = BasisSet.shells[si].center - 1;
70
 
  int c2 = BasisSet.shells[sj].center - 1;
71
 
  int c3 = BasisSet.shells[sk].center - 1;
72
 
  int c4 = BasisSet.shells[sl].center - 1;
73
 
 
74
 
  // If all centers are the same and angular momentum of functions is odd -- the ERI will be 0 by AM selection
75
 
  if ( (c1 == c2) && (c1 == c3) && (c1 == c4) && (am%2 == 1)) {
76
 
    return 0;
77
 
  }
78
 
 
79
 
  struct shell_pair* sp_ij = &(BasisSet.shell_pairs[si][sj]);
80
 
  struct shell_pair* sp_kl = &(BasisSet.shell_pairs[sk][sl]);
81
 
 
82
 
  Libint->AB[0] = sp_ij->AB[0];
83
 
  Libint->AB[1] = sp_ij->AB[1];
84
 
  Libint->AB[2] = sp_ij->AB[2];
85
 
  Libint->CD[0] = sp_kl->AB[0];
86
 
  Libint->CD[1] = sp_kl->AB[1];
87
 
  Libint->CD[2] = sp_kl->AB[2];
88
 
            
89
 
  double AB2 = Libint->AB[0]*Libint->AB[0]+
90
 
               Libint->AB[1]*Libint->AB[1]+
91
 
               Libint->AB[2]*Libint->AB[2];
92
 
  double CD2 = Libint->CD[0]*Libint->CD[0]+
93
 
               Libint->CD[1]*Libint->CD[1]+
94
 
               Libint->CD[2]*Libint->CD[2];
95
 
 
96
 
  /*--- Compute data for primitive quartets here ---*/
97
 
  int np1 = BasisSet.shells[si].n_prims;
98
 
  int np2 = BasisSet.shells[sj].n_prims;
99
 
  int np3 = BasisSet.shells[sk].n_prims;
100
 
  int np4 = BasisSet.shells[sl].n_prims;
101
 
 
102
 
  /*--- Compute data for primitive quartets here ---*/
103
 
  int num_prim_comb = 0;
104
 
//  bool si_eq_sj = (si == sj && am1 == am2);
105
 
//  bool sk_eq_sl = (sk == sl && am3 == am4);
106
 
  bool si_eq_sj = false;
107
 
  bool sk_eq_sl = false;
108
 
  for (int p1 = 0; p1 < np1; p1++) {
109
 
    int max_p2 = si_eq_sj ? p1+1 : np2;
110
 
    for (int p2 = 0; p2 < max_p2; p2++) {
111
 
      int m = (1 + (si_eq_sj && p1 != p2));
112
 
      for (int p3 = 0; p3 < np3; p3++) {
113
 
        int max_p4 = sk_eq_sl ? p3+1 : np4;
114
 
        for (int p4 = 0; p4 < max_p4; p4++){
115
 
          int n = m * (1 + (sk_eq_sl && p3 != p4));
 
41
      
 
42
      /* place in "ascending" angular mom-
 
43
         my simple way of optimizing PHG recursion (VRR) */
 
44
      
 
45
      bool switch_bra_ket = false;
 
46
      /* this should be /good/ for the VRR */
 
47
      if ( BasisSet.shells[si].am + BasisSet.shells[sj].am + inc1 + inc2 >
 
48
           BasisSet.shells[sk].am + BasisSet.shells[sl].am + inc3 + inc4 ) {
 
49
        switch_bra_ket = true;
 
50
        switch_ij(si,sk);
 
51
        switch_ij(sj,sl);
 
52
        switch_ij(inc1, inc3);
 
53
        switch_ij(inc2, inc4);
 
54
      }
 
55
      
 
56
      /* these two are good for the HRR */
 
57
      bool switch_bra = false;
 
58
      if(BasisSet.shells[si].am + inc1 < BasisSet.shells[sj].am + inc2){
 
59
        switch_bra = true;
 
60
        switch_ij(si,sj);
 
61
        switch_ij(inc1,inc2);
 
62
      }
 
63
      bool switch_ket = false;
 
64
      if(BasisSet.shells[sk].am + inc3 < BasisSet.shells[sl].am + inc4){
 
65
        switch_ket = true;
 
66
        switch_ij(sk,sl);
 
67
        switch_ij(inc3,inc4);
 
68
      }
 
69
      
 
70
      int am1 = BasisSet.shells[si].am - 1 + inc1;
 
71
      int am2 = BasisSet.shells[sj].am - 1 + inc2;
 
72
      int am3 = BasisSet.shells[sk].am - 1 + inc3;
 
73
      int am4 = BasisSet.shells[sl].am - 1 + inc4;
 
74
      int am = am1 + am2 + am3 + am4;
 
75
      
 
76
      int c1 = BasisSet.shells[si].center - 1;
 
77
      int c2 = BasisSet.shells[sj].center - 1;
 
78
      int c3 = BasisSet.shells[sk].center - 1;
 
79
      int c4 = BasisSet.shells[sl].center - 1;
 
80
      
 
81
      // If all centers are the same and angular momentum of functions is odd -- the ERI will be 0 by AM selection
 
82
      if ( (c1 == c2) && (c1 == c3) && (c1 == c4) && (am%2 == 1)) {
 
83
        return 0;
 
84
      }
 
85
      
 
86
      struct shell_pair* sp_ij = &(BasisSet.shell_pairs[si][sj]);
 
87
      struct shell_pair* sp_kl = &(BasisSet.shell_pairs[sk][sl]);
 
88
      
 
89
      Libint->AB[0] = sp_ij->AB[0];
 
90
      Libint->AB[1] = sp_ij->AB[1];
 
91
      Libint->AB[2] = sp_ij->AB[2];
 
92
      Libint->CD[0] = sp_kl->AB[0];
 
93
      Libint->CD[1] = sp_kl->AB[1];
 
94
      Libint->CD[2] = sp_kl->AB[2];
 
95
      
 
96
      double AB2 = 
 
97
        Libint->AB[0]*Libint->AB[0]+
 
98
        Libint->AB[1]*Libint->AB[1]+
 
99
        Libint->AB[2]*Libint->AB[2];
 
100
      double CD2 = 
 
101
        Libint->CD[0]*Libint->CD[0]+
 
102
        Libint->CD[1]*Libint->CD[1]+
 
103
        Libint->CD[2]*Libint->CD[2];
 
104
      
 
105
      /*--- Compute data for primitive quartets here ---*/
 
106
      int np1 = BasisSet.shells[si].n_prims;
 
107
      int np2 = BasisSet.shells[sj].n_prims;
 
108
      int np3 = BasisSet.shells[sk].n_prims;
 
109
      int np4 = BasisSet.shells[sl].n_prims;
 
110
      
 
111
      /*--- Compute data for primitive quartets here ---*/
 
112
      int num_prim_comb = 0;
 
113
      //  bool si_eq_sj = (si == sj && am1 == am2);
 
114
      //  bool sk_eq_sl = (sk == sl && am3 == am4);
 
115
      bool si_eq_sj = false;
 
116
      bool sk_eq_sl = false;
 
117
      for (int p1 = 0; p1 < np1; p1++) {
 
118
        int max_p2 = si_eq_sj ? p1+1 : np2;
 
119
        for (int p2 = 0; p2 < max_p2; p2++) {
 
120
          int m = (1 + (si_eq_sj && p1 != p2));
 
121
          for (int p3 = 0; p3 < np3; p3++) {
 
122
            int max_p4 = sk_eq_sl ? p3+1 : np4;
 
123
            for (int p4 = 0; p4 < max_p4; p4++){
 
124
              int n = m * (1 + (sk_eq_sl && p3 != p4));
116
125
#ifdef USE_TAYLOR_FM
117
 
          quartet_data(&(Libint->PrimQuartet[num_prim_comb++]), NULL, AB2, CD2,
118
 
                       sp_ij, sp_kl, am, p1, p2, p3, p4, n);
119
 
#else
120
 
          quartet_data(&(Libint->PrimQuartet[num_prim_comb++]), &fjt_table, AB2, CD2,
121
 
                       sp_ij, sp_kl, am, p1, p2, p3, p4, n);
122
 
#endif
123
 
        }
124
 
      }
125
 
    }
126
 
  }
127
 
 
128
 
  int nao[4];
129
 
  nao[0] = ioff[am1+1];
130
 
  nao[1] = ioff[am2+1];
131
 
  nao[2] = ioff[am3+1];
132
 
  nao[3] = ioff[am4+1];
133
 
  int nbra = nao[0] * nao[1];
134
 
  int nket = nao[2] * nao[3];
135
 
  int size = nbra * nket;
136
 
  
137
 
  // In general we need a couple extra buffers to be able to copy and permute integrals
138
 
  static int max_size = ioff[BasisSet.max_am]*ioff[BasisSet.max_am]*ioff[BasisSet.max_am]*ioff[BasisSet.max_am];
139
 
  static double* perm_buf = new double[max_size];
140
 
  static double* copy_buf = new double[max_size];
141
 
  if (size > max_size) {
142
 
    delete[] perm_buf;
143
 
    delete[] copy_buf;
144
 
    perm_buf = new double[size];
145
 
    copy_buf = new double[size];
146
 
    max_size = size;
147
 
  }
148
 
#ifdef NONDOUBLE_INTS
149
 
  double* raw_data = copy_buf;
150
 
#endif
151
 
 
152
 
  if (am) {
153
 
#ifdef NONDOUBLE_INTS
154
 
    REALTYPE* target_ptr = build_eri[am1][am2][am3][am4](Libint,num_prim_comb);
155
 
    for(i=0;i<size;i++)
156
 
      raw_data[i] = (double) target_ptr[i];
157
 
#else
158
 
    double* raw_data = build_eri[am1][am2][am3][am4](Libint,num_prim_comb);
159
 
#endif
160
 
 
161
 
    //
162
 
    // Permute shells back, if necessary
163
 
    if (!do_not_permute) {
164
 
      int num_perm = (int) switch_bra_ket + (int) switch_bra + (int) switch_ket;
165
 
    
166
 
      double *sbuf1, *sbuf2, *sbuf3;
167
 
      double *tbuf1, *tbuf2, *tbuf3;
168
 
      // in general, permutations will require careful orchestration of
169
 
      // the utilization of the buffers. If no permutations are necessary -- just copy the data
170
 
      if (num_perm == 0) {
171
 
        for(int i=0; i<size; i++)
172
 
          target[i] = raw_data[i];
173
 
      }
174
 
      // otherwise -- do the dance
175
 
      else if (num_perm == 1) {
176
 
        sbuf1 = sbuf2 = sbuf3 = raw_data;
177
 
        tbuf1 = tbuf2 = tbuf3 = target;
178
 
      }
179
 
      else if (num_perm == 2) {
180
 
        if (!switch_bra_ket) {
181
 
          sbuf1 = raw_data;
182
 
          tbuf1 = sbuf2 = perm_buf;
183
 
          tbuf2 = target;
184
 
        }
185
 
        else if (!switch_ket) {
186
 
          sbuf1 = raw_data;
187
 
          tbuf1 = sbuf3 = perm_buf;
188
 
          tbuf3 = target;
189
 
        }
190
 
        else if (!switch_bra) {
191
 
          sbuf2 = raw_data;
192
 
          tbuf2 = sbuf3 = perm_buf;
193
 
          tbuf3 = target;
194
 
        }
195
 
      }
196
 
      else if (num_perm == 3) {
197
 
        sbuf1 = raw_data;
198
 
        tbuf1 = sbuf2 = perm_buf;
199
 
        tbuf2 = sbuf3 = raw_data;
200
 
        tbuf3 = target;
201
 
      }
202
 
 
203
 
      // Note that the order of permutations is opposite to the order of shell index permutations
204
 
      // at the beginning of this function
205
 
      if (switch_bra) {
206
 
        ijkl_to_jikl(sbuf1, tbuf1, nao[0], nao[1], nket);
207
 
        switch_ij(si,sj);
208
 
        switch_ij(inc1,inc2);
209
 
      }
210
 
      if (switch_ket) {
211
 
        ijkl_to_ijlk(sbuf2, tbuf2, nbra, nao[2], nao[3]);
212
 
        switch_ij(sk,sl);
213
 
        switch_ij(inc3,inc4);
214
 
      }
215
 
      if (switch_bra_ket) {
216
 
        ijkl_to_klij(sbuf3, tbuf3, nbra, nket);
217
 
        switch_ij(si,sk);
218
 
        switch_ij(sj,sl);
219
 
        switch_ij(inc1, inc3);
220
 
        switch_ij(inc2, inc4);
221
 
      }
222
 
    }
223
 
  }
224
 
  else {
225
 
    REALTYPE temp = 0.0;
226
 
    prim_data* prim_quartet = Libint->PrimQuartet;
227
 
    for(int p=0;p<num_prim_comb;p++,prim_quartet++) {
228
 
      temp += prim_quartet->F[0];
229
 
    }
230
 
    target[0] = (double) temp;
231
 
  }
232
 
 
233
 
  return size;
 
126
              quartet_data(&(Libint->PrimQuartet[num_prim_comb++]), NULL, AB2, CD2,
 
127
                           sp_ij, sp_kl, am, p1, p2, p3, p4, n);
 
128
#else
 
129
              quartet_data(&(Libint->PrimQuartet[num_prim_comb++]), &fjt_table, AB2, CD2,
 
130
                           sp_ij, sp_kl, am, p1, p2, p3, p4, n);
 
131
#endif
 
132
            }
 
133
          }
 
134
        }
 
135
      }
 
136
      
 
137
      int nao[4];
 
138
      nao[0] = ioff[am1+1];
 
139
      nao[1] = ioff[am2+1];
 
140
      nao[2] = ioff[am3+1];
 
141
      nao[3] = ioff[am4+1];
 
142
      int nbra = nao[0] * nao[1];
 
143
      int nket = nao[2] * nao[3];
 
144
      int size = nbra * nket;
 
145
      
 
146
      // In general we need a couple extra buffers to be able to copy and permute integrals
 
147
      static int max_size = ioff[BasisSet.max_am]*ioff[BasisSet.max_am]*ioff[BasisSet.max_am]*ioff[BasisSet.max_am];
 
148
      static double* perm_buf = new double[max_size];
 
149
      static double* copy_buf = new double[max_size];
 
150
      if (size > max_size) {
 
151
        delete[] perm_buf;
 
152
        delete[] copy_buf;
 
153
        perm_buf = new double[size];
 
154
        copy_buf = new double[size];
 
155
        max_size = size;
 
156
      }
 
157
#ifdef NONDOUBLE_INTS
 
158
      double* raw_data = copy_buf;
 
159
#endif
 
160
      
 
161
      if (am) {
 
162
#ifdef NONDOUBLE_INTS
 
163
        REALTYPE* target_ptr = build_eri[am1][am2][am3][am4](Libint,num_prim_comb);
 
164
        for(i=0;i<size;i++)
 
165
          raw_data[i] = (double) target_ptr[i];
 
166
#else
 
167
        double* raw_data = build_eri[am1][am2][am3][am4](Libint,num_prim_comb);
 
168
#endif
 
169
        
 
170
        //
 
171
        // Permute shells back, if necessary
 
172
        if (!do_not_permute) {
 
173
          int num_perm = (int) switch_bra_ket + (int) switch_bra + (int) switch_ket;
 
174
          
 
175
          double *sbuf1, *sbuf2, *sbuf3;
 
176
          double *tbuf1, *tbuf2, *tbuf3;
 
177
          // in general, permutations will require careful orchestration of
 
178
          // the utilization of the buffers. If no permutations are necessary -- just copy the data
 
179
          if (num_perm == 0) {
 
180
            for(int i=0; i<size; i++)
 
181
              target[i] = raw_data[i];
 
182
          }
 
183
          // otherwise -- do the dance
 
184
          else if (num_perm == 1) {
 
185
            sbuf1 = sbuf2 = sbuf3 = raw_data;
 
186
            tbuf1 = tbuf2 = tbuf3 = target;
 
187
          }
 
188
          else if (num_perm == 2) {
 
189
            if (!switch_bra_ket) {
 
190
              sbuf1 = raw_data;
 
191
              tbuf1 = sbuf2 = perm_buf;
 
192
              tbuf2 = target;
 
193
            }
 
194
            else if (!switch_ket) {
 
195
              sbuf1 = raw_data;
 
196
              tbuf1 = sbuf3 = perm_buf;
 
197
              tbuf3 = target;
 
198
            }
 
199
            else if (!switch_bra) {
 
200
              sbuf2 = raw_data;
 
201
              tbuf2 = sbuf3 = perm_buf;
 
202
              tbuf3 = target;
 
203
            }
 
204
          }
 
205
          else if (num_perm == 3) {
 
206
            sbuf1 = raw_data;
 
207
            tbuf1 = sbuf2 = perm_buf;
 
208
            tbuf2 = sbuf3 = raw_data;
 
209
            tbuf3 = target;
 
210
          }
 
211
          
 
212
          // Note that the order of permutations is opposite to the order of shell index permutations
 
213
          // at the beginning of this function
 
214
          if (switch_bra) {
 
215
            ijkl_to_jikl(sbuf1, tbuf1, nao[0], nao[1], nket);
 
216
            switch_ij(si,sj);
 
217
            switch_ij(inc1,inc2);
 
218
          }
 
219
          if (switch_ket) {
 
220
            ijkl_to_ijlk(sbuf2, tbuf2, nbra, nao[2], nao[3]);
 
221
            switch_ij(sk,sl);
 
222
            switch_ij(inc3,inc4);
 
223
          }
 
224
          if (switch_bra_ket) {
 
225
            ijkl_to_klij(sbuf3, tbuf3, nbra, nket);
 
226
            switch_ij(si,sk);
 
227
            switch_ij(sj,sl);
 
228
            switch_ij(inc1, inc3);
 
229
            switch_ij(inc2, inc4);
 
230
          }
 
231
        }
 
232
      }
 
233
      else {
 
234
        REALTYPE temp = 0.0;
 
235
        prim_data* prim_quartet = Libint->PrimQuartet;
 
236
        for(int p=0;p<num_prim_comb;p++,prim_quartet++) {
 
237
          temp += prim_quartet->F[0];
 
238
        }
 
239
        target[0] = (double) temp;
 
240
      }
 
241
      
 
242
      return size;
 
243
    }
 
244
  }
234
245
}