← Back to branch summary

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

« back to all changes in this revision

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