← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/input/overlap.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/input/overlap.cc
 
1
/*! \file
 
2
    \ingroup INPUT
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <cstring>
 
8
#include <cmath>
 
9
#include <libciomr/libciomr.h>
 
10
 
 
11
#include "input.h"
 
12
#include "defines.h"
 
13
#define EXTERN
 
14
#include "global.h"
 
15
 
 
16
namespace psi { namespace input {
 
17
 
 
18
/*--- External function declaration ---*/
 
19
extern void OI_OSrecurs(double **OIX, double **OIY, double **OIZ, struct coordinates PA, struct coordinates PB,
 
20
                        double gamma, int lmaxi, int lmaxj);
 
21
 
 
22
/*------------------------------------------------------
 
23
  This function computes the overlap matrix S12 between
 
24
  the new and old bases (in SO basis)
 
25
 ------------------------------------------------------*/
 
26
double **overlap_new_old()
 
27
{
 
28
  struct coordinates A, B, P, PA, PB;
 
29
 
 
30
  int ii, jj, kk, ll;
 
31
  int maxsz1, maxsz2;
 
32
  int x1, x2, y1, y2, z1, z2;
 
33
  int s1, s2, ao1, ao2, l1, l2, n1, n2, p1, p2,
 
34
    prim1, prim2, fprim1, fprim2, nprim1, nprim2, atom1, atom2, bf1, bf2;
 
35
  double exp1, exp2, cc1, cc2, norm1;
 
36
  double ab2, pfac, gam, oog;
 
37
  double *ptr1, *ptr2;
 
38
  double **stemp, **S12_AO, **S12, **OIX, **OIY, **OIZ;
 
39
  double **tmpmat;
 
40
 
 
41
  S12_AO = block_matrix(num_ao, Oldcalc.num_ao);
 
42
  maxsz1 = ioff[max_angmom+1];
 
43
  maxsz2 = ioff[Oldcalc.max_angmom+1];
 
44
  stemp = block_matrix(maxsz1,maxsz2);
 
45
 
 
46
  OIX = block_matrix(max_angmom+1,Oldcalc.max_angmom+1);
 
47
  OIY = block_matrix(max_angmom+1,Oldcalc.max_angmom+1);
 
48
  OIZ = block_matrix(max_angmom+1,Oldcalc.max_angmom+1);
 
49
  
 
50
  /*--- Loop over shells ---*/
 
51
  for(s1=0;s1<num_shells;s1++) {
 
52
    ao1 = first_ao_shell[s1];
 
53
    nprim1 = nprim_in_shell[s1];
 
54
    fprim1 = first_prim_shell[s1];
 
55
    bf1 = first_ao_shell[s1];
 
56
    l1 = shell_ang_mom[s1];
 
57
    n1 = ioff[l1+1];
 
58
    atom1 = shell_nucleus[s1];
 
59
    A.x = geometry[atom1][0];
 
60
    A.y = geometry[atom1][1];
 
61
    A.z = geometry[atom1][2];
 
62
    for(s2=0;s2<Oldcalc.num_shells;s2++) {
 
63
      ao2 = Oldcalc.first_ao_shell[s2];
 
64
      nprim2 = Oldcalc.nprim_in_shell[s2];
 
65
      fprim2 = Oldcalc.first_prim_shell[s2];
 
66
      bf2 = Oldcalc.first_ao_shell[s2];
 
67
      l2 = Oldcalc.shell_ang_mom[s2];
 
68
      n2 = ioff[l2+1];
 
69
      atom2 = Oldcalc.shell_nucleus[s2];
 
70
      B.x = Oldcalc.geometry[atom2][0];
 
71
      B.y = Oldcalc.geometry[atom2][1];
 
72
      B.z = Oldcalc.geometry[atom2][2];
 
73
          
 
74
      ab2  = (A.x-B.x)*(A.x-B.x);
 
75
      ab2 += (A.y-B.y)*(A.y-B.y);
 
76
      ab2 += (A.z-B.z)*(A.z-B.z);
 
77
 
 
78
      /*--- zero the temporary storage for accumulating contractions ---*/
 
79
      memset(stemp[0],0,sizeof(double)*maxsz1*maxsz2);
 
80
          
 
81
      /*--- Loop over primitives here ---*/
 
82
      for(p1=0,prim1=fprim1;p1<nprim1;p1++,prim1++) {
 
83
        exp1 = exponents[prim1];
 
84
        cc1 = contr_coeff[prim1];
 
85
        for(p2=0,prim2=fprim2;p2<nprim2;p2++,prim2++) {
 
86
          exp2 = Oldcalc.exponents[prim2];
 
87
          cc2 = Oldcalc.contr_coeff[prim2][l2];
 
88
 
 
89
          gam = exp1 + exp2;
 
90
          oog = 1.0/gam;
 
91
          P.x = (exp1*A.x + exp2*B.x)*oog;
 
92
          P.y = (exp1*A.y + exp2*B.y)*oog;
 
93
          P.z = (exp1*A.z + exp2*B.z)*oog;
 
94
          PA.x = P.x - A.x;
 
95
          PA.y = P.y - A.y;
 
96
          PA.z = P.z - A.z;
 
97
          PB.x = P.x - B.x;
 
98
          PB.y = P.y - B.y;
 
99
          PB.z = P.z - B.z;
 
100
 
 
101
          pfac = exp(-exp1*exp2*ab2*oog)*sqrt(M_PI*oog)*M_PI*oog*cc1*cc2;
 
102
 
 
103
          OI_OSrecurs(OIX,OIY,OIZ,PA,PB,gam,l1,l2);
 
104
 
 
105
          /*--- create all am components of si ---*/
 
106
          ao1 = 0;
 
107
          for(ii = 0; ii <= l1; ii++){
 
108
            x1 = l1 - ii;
 
109
            for(jj = 0; jj <= ii; jj++){
 
110
              y1 = ii - jj;
 
111
              z1 = jj;
 
112
              /*--- create all am components of sj ---*/
 
113
              ao2 = 0;
 
114
              for(kk = 0; kk <= l2; kk++){
 
115
                x2 = l2 - kk;
 
116
                for(ll = 0; ll <= kk; ll++){
 
117
                  y2 = kk - ll;
 
118
                  z2 = ll;
 
119
                                  
 
120
                  stemp[ao1][ao2] += pfac*OIX[x1][x2]*OIY[y1][y2]*OIZ[z1][z2];
 
121
                                  
 
122
                  ao2++;
 
123
                }
 
124
              }
 
125
              ao1++;
 
126
            }
 
127
          } /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
 
128
        }
 
129
      } /*--- end primitive contraction ---*/
 
130
 
 
131
      /*--- Normalize the contracted integrals and put them into S12 ---*/
 
132
      ptr1 = GTOs.bf_norm[l1];
 
133
      ptr2 = GTOs.bf_norm[l2];
 
134
      for(ii=0; ii<n1; ii++) {
 
135
        norm1 = ptr1[ii];
 
136
        for(jj=0; jj<n2; jj++) {
 
137
          S12_AO[bf1+ii][bf2+jj] = stemp[ii][jj]*norm1*ptr2[jj];
 
138
        }
 
139
      }
 
140
          
 
141
    }
 
142
  }   /*--- This shell pair is done ---*/
 
143
 
 
144
  free_block(OIX);
 
145
  free_block(OIY);
 
146
  free_block(OIZ);
 
147
  free_block(stemp);
 
148
 
 
149
  /*--- transform to SO basis ---*/
 
150
  tmpmat = block_matrix(num_so,Oldcalc.num_ao);
 
151
  mmult(usotao,0,S12_AO,0,tmpmat,0,num_so,num_ao,Oldcalc.num_ao,0);
 
152
  free_block(S12_AO);
 
153
  S12 = block_matrix(num_so,Oldcalc.num_so);
 
154
  mmult(tmpmat,0,Oldcalc.usotao,1,S12,0,num_so,Oldcalc.num_ao,Oldcalc.num_so,0);
 
155
  free_block(tmpmat);
 
156
 
 
157
  if (print_lvl >= DEBUGPRINT) {
 
158
    fprintf(outfile,"  -Overlap matrix between the new and old basis (in SO basis):\n");
 
159
    print_mat(S12,num_so,Oldcalc.num_so,outfile);
 
160
  }
 
161
  
 
162
  return S12;
 
163
}
 
164
 
 
165
 
 
166
/*--------------------------------------------------
 
167
  This function computes the overlap matrix S11 for
 
168
  the new basis set (in SO basis)
 
169
 --------------------------------------------------*/
 
170
double **overlap()
 
171
{
 
172
  struct coordinates A, B, P, PA, PB;
 
173
 
 
174
  int ii, jj, kk, ll;
 
175
  int maxsz1, maxsz2;
 
176
  int x1, x2, y1, y2, z1, z2;
 
177
  int s1, s2, ao1, ao2, l1, l2, n1, n2, p1, p2,
 
178
      prim1, prim2, fprim1, fprim2, nprim1, nprim2, atom1, atom2, bf1, bf2;
 
179
  double exp1, exp2, cc1, cc2, norm1;
 
180
  double ab2, pfac, gam, oog;
 
181
  double *ptr1, *ptr2;
 
182
  double **stemp, **S11_AO, **S11, **OIX, **OIY, **OIZ;
 
183
  double **tmpmat;
 
184
 
 
185
  S11_AO = block_matrix(num_ao, num_ao);
 
186
  maxsz1 = ioff[max_angmom+1];
 
187
  maxsz2 = ioff[max_angmom+1];
 
188
  stemp = block_matrix(maxsz1,maxsz2);
 
189
 
 
190
  OIX = block_matrix(max_angmom+1,max_angmom+1);
 
191
  OIY = block_matrix(max_angmom+1,max_angmom+1);
 
192
  OIZ = block_matrix(max_angmom+1,max_angmom+1);
 
193
  
 
194
  /*--- Loop over shells ---*/
 
195
  for(s1=0;s1<num_shells;s1++) {
 
196
      ao1 = first_ao_shell[s1];
 
197
      nprim1 = nprim_in_shell[s1];
 
198
      fprim1 = first_prim_shell[s1];
 
199
      bf1 = first_ao_shell[s1];
 
200
      l1 = shell_ang_mom[s1];
 
201
      n1 = ioff[l1+1];
 
202
      atom1 = shell_nucleus[s1];
 
203
      A.x = geometry[atom1][0];
 
204
      A.y = geometry[atom1][1];
 
205
      A.z = geometry[atom1][2];
 
206
      for(s2=0;s2<num_shells;s2++) {
 
207
          ao2 = first_ao_shell[s2];
 
208
          nprim2 = nprim_in_shell[s2];
 
209
          fprim2 = first_prim_shell[s2];
 
210
          bf2 = first_ao_shell[s2];
 
211
          l2 = shell_ang_mom[s2];
 
212
          n2 = ioff[l2+1];
 
213
          atom2 = shell_nucleus[s2];
 
214
          B.x = geometry[atom2][0];
 
215
          B.y = geometry[atom2][1];
 
216
          B.z = geometry[atom2][2];
 
217
          
 
218
          ab2  = (A.x-B.x)*(A.x-B.x);
 
219
          ab2 += (A.y-B.y)*(A.y-B.y);
 
220
          ab2 += (A.z-B.z)*(A.z-B.z);
 
221
 
 
222
          /*--- zero the temporary storage for accumulating contractions ---*/
 
223
          memset(stemp[0],0,sizeof(double)*maxsz1*maxsz2);
 
224
 
 
225
          /*--- Loop over primitives here ---*/
 
226
          for(p1=0,prim1=fprim1;p1<nprim1;p1++,prim1++) {
 
227
              exp1 = exponents[prim1];
 
228
              cc1 = contr_coeff[prim1];
 
229
              for(p2=0,prim2=fprim2;p2<nprim2;p2++,prim2++) {
 
230
                  exp2 = exponents[prim2];
 
231
                  cc2 = contr_coeff[prim2];
 
232
 
 
233
                  gam = exp1 + exp2;
 
234
                  oog = 1.0/gam;
 
235
                  P.x = (exp1*A.x + exp2*B.x)*oog;
 
236
                  P.y = (exp1*A.y + exp2*B.y)*oog;
 
237
                  P.z = (exp1*A.z + exp2*B.z)*oog;
 
238
                  PA.x = P.x - A.x;
 
239
                  PA.y = P.y - A.y;
 
240
                  PA.z = P.z - A.z;
 
241
                  PB.x = P.x - B.x;
 
242
                  PB.y = P.y - B.y;
 
243
                  PB.z = P.z - B.z;
 
244
 
 
245
                  pfac = exp(-exp1*exp2*ab2*oog)*sqrt(M_PI*oog)*M_PI*oog*cc1*cc2;
 
246
 
 
247
                  OI_OSrecurs(OIX,OIY,OIZ,PA,PB,gam,l1,l2);
 
248
 
 
249
                  /*--- create all am components of si ---*/
 
250
                  ao1 = 0;
 
251
                  for(ii = 0; ii <= l1; ii++){
 
252
                      x1 = l1 - ii;
 
253
                      for(jj = 0; jj <= ii; jj++){
 
254
                          y1 = ii - jj;
 
255
                          z1 = jj;
 
256
                          /*--- create all am components of sj ---*/
 
257
                          ao2 = 0;
 
258
                          for(kk = 0; kk <= l2; kk++){
 
259
                              x2 = l2 - kk;
 
260
                              for(ll = 0; ll <= kk; ll++){
 
261
                                  y2 = kk - ll;
 
262
                                  z2 = ll;
 
263
                                  
 
264
                                  stemp[ao1][ao2] += pfac*OIX[x1][x2]*OIY[y1][y2]*OIZ[z1][z2];
 
265
                                  
 
266
                                  ao2++;
 
267
                              }
 
268
                          }
 
269
                          ao1++;
 
270
                      }
 
271
                  } /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
 
272
              }
 
273
          } /*--- end primitive contraction ---*/
 
274
 
 
275
          /*--- Normalize the contracted integrals and put them into S11 ---*/
 
276
          ptr1 = GTOs.bf_norm[l1];
 
277
          ptr2 = GTOs.bf_norm[l2];
 
278
          for(ii=0; ii<n1; ii++) {
 
279
              norm1 = ptr1[ii];
 
280
              for(jj=0; jj<n2; jj++) {
 
281
                  S11_AO[bf1+ii][bf2+jj] = stemp[ii][jj]*norm1*ptr2[jj];
 
282
              }
 
283
          }
 
284
      }
 
285
  }   /*--- This shell pair is done ---*/
 
286
 
 
287
  free_block(OIX);
 
288
  free_block(OIY);
 
289
  free_block(OIZ);
 
290
  free_block(stemp);
 
291
 
 
292
  /*--- transform to SO basis ---*/
 
293
  tmpmat = block_matrix(num_so,num_ao);
 
294
  mmult(usotao,0,S11_AO,0,tmpmat,0,num_so,num_ao,num_ao,0);
 
295
  free_block(S11_AO);
 
296
  S11 = block_matrix(num_so,num_so);
 
297
  mmult(tmpmat,0,usotao,1,S11,0,num_so,num_ao,num_so,0);
 
298
  free_block(tmpmat);
 
299
 
 
300
  if (print_lvl >= DEBUGPRINT) {
 
301
      fprintf(outfile,"  -Overlap matrix in the new basis (in SO basis):\n");
 
302
      print_mat(S11,num_so,num_so,outfile);
 
303
  }
 
304
  
 
305
  return S11;
 
306
}
 
307
 
 
308
}} // namespace psi::input