← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

Viewing changes to src/bin/input/overlap_float.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_float.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
#include "float.h"
 
17
#include "linalg.h"
 
18
 
 
19
namespace psi { namespace input {
 
20
 
 
21
/*--- External function declaration ---*/
 
22
extern void OI_OSrecurs_float(FLOAT** OIX, FLOAT** OIY, FLOAT** OIZ, FLOAT PA[3], FLOAT PB[3],
 
23
                              FLOAT gamma, int lmaxi, int lmaxj);
 
24
 
 
25
/*------------------------------------------------------
 
26
  This function computes the overlap matrix S12 between
 
27
  the new and old bases (in SO basis)
 
28
 ------------------------------------------------------*/
 
29
FLOAT **overlap_new_old_float()
 
30
{
 
31
  struct coordinates A, B;
 
32
  FLOAT P[3], PA[3], PB[3];
 
33
 
 
34
  int ii, jj, kk, ll;
 
35
  int maxsz1, maxsz2;
 
36
  int x1, x2, y1, y2, z1, z2;
 
37
  int s1, s2, ao1, ao2, l1, l2, n1, n2, p1, p2,
 
38
      prim1, prim2, fprim1, fprim2, nprim1, nprim2, atom1, atom2, bf1, bf2;
 
39
  FLOAT exp1, exp2, cc1, cc2, norm1;
 
40
  FLOAT ab2, pfac, gam, oog;
 
41
  double *ptr1, *ptr2;
 
42
  FLOAT **stemp, **S12_AO, **S12, **OIX, **OIY, **OIZ;
 
43
  FLOAT **tmpmat;
 
44
  FLOAT **usotao_FLOAT, **oldusotao_FLOAT_transp;
 
45
 
 
46
  S12_AO = create_matrix(num_ao, Oldcalc.num_ao);
 
47
  maxsz1 = ioff[max_angmom+1];
 
48
  maxsz2 = ioff[Oldcalc.max_angmom+1];
 
49
  stemp = create_matrix(maxsz1,maxsz2);
 
50
 
 
51
  OIX = create_matrix(max_angmom+1,Oldcalc.max_angmom+1);
 
52
  OIY = create_matrix(max_angmom+1,Oldcalc.max_angmom+1);
 
53
  OIZ = create_matrix(max_angmom+1,Oldcalc.max_angmom+1);
 
54
  
 
55
  /*--- Loop over shells ---*/
 
56
  for(s1=0;s1<num_shells;s1++) {
 
57
      ao1 = first_ao_shell[s1];
 
58
      nprim1 = nprim_in_shell[s1];
 
59
      fprim1 = first_prim_shell[s1];
 
60
      bf1 = first_ao_shell[s1];
 
61
      l1 = shell_ang_mom[s1];
 
62
      n1 = ioff[l1+1];
 
63
      atom1 = shell_nucleus[s1];
 
64
      A.x = geometry[atom1][0];
 
65
      A.y = geometry[atom1][1];
 
66
      A.z = geometry[atom1][2];
 
67
      for(s2=0;s2<Oldcalc.num_shells;s2++) {
 
68
          ao2 = Oldcalc.first_ao_shell[s2];
 
69
          nprim2 = Oldcalc.nprim_in_shell[s2];
 
70
          fprim2 = Oldcalc.first_prim_shell[s2];
 
71
          bf2 = Oldcalc.first_ao_shell[s2];
 
72
          l2 = Oldcalc.shell_ang_mom[s2];
 
73
          n2 = ioff[l2+1];
 
74
          atom2 = Oldcalc.shell_nucleus[s2];
 
75
          B.x = Oldcalc.geometry[atom2][0];
 
76
          B.y = Oldcalc.geometry[atom2][1];
 
77
          B.z = Oldcalc.geometry[atom2][2];
 
78
          
 
79
          ab2  = (A.x-B.x)*(A.x-B.x);
 
80
          ab2 += (A.y-B.y)*(A.y-B.y);
 
81
          ab2 += (A.z-B.z)*(A.z-B.z);
 
82
 
 
83
          /*--- zero the temporary storage for accumulating contractions ---*/
 
84
          memset(stemp[0],0,sizeof(FLOAT)*maxsz1*maxsz2);
 
85
          
 
86
          /*--- Loop over primitives here ---*/
 
87
          for(p1=0,prim1=fprim1;p1<nprim1;p1++,prim1++) {
 
88
              exp1 = exponents[prim1];
 
89
              cc1 = contr_coeff[prim1];
 
90
              for(p2=0,prim2=fprim2;p2<nprim2;p2++,prim2++) {
 
91
                  exp2 = Oldcalc.exponents[prim2];
 
92
                  cc2 = Oldcalc.contr_coeff[prim2][l2];
 
93
 
 
94
                  gam = exp1 + exp2;
 
95
                  oog = 1.0/gam;
 
96
                  P[0] = (exp1*A.x + exp2*B.x)*oog;
 
97
                  P[1] = (exp1*A.y + exp2*B.y)*oog;
 
98
                  P[2] = (exp1*A.z + exp2*B.z)*oog;
 
99
                  PA[0] = P[0] - A.x;
 
100
                  PA[1] = P[1] - A.y;
 
101
                  PA[2] = P[2] - A.z;
 
102
                  PB[0] = P[0] - B.x;
 
103
                  PB[1] = P[1] - B.y;
 
104
                  PB[2] = P[2] - B.z;
 
105
 
 
106
                  pfac = EXP(-exp1*exp2*ab2*oog)*SQRT(M_PI*oog)*M_PI*oog*cc1*cc2;
 
107
 
 
108
                  OI_OSrecurs_float(OIX,OIY,OIZ,PA,PB,gam,l1,l2);
 
109
 
 
110
                  /*--- create all am components of si ---*/
 
111
                  ao1 = 0;
 
112
                  for(ii = 0; ii <= l1; ii++){
 
113
                      x1 = l1 - ii;
 
114
                      for(jj = 0; jj <= ii; jj++){
 
115
                          y1 = ii - jj;
 
116
                          z1 = jj;
 
117
                          /*--- create all am components of sj ---*/
 
118
                          ao2 = 0;
 
119
                          for(kk = 0; kk <= l2; kk++){
 
120
                              x2 = l2 - kk;
 
121
                              for(ll = 0; ll <= kk; ll++){
 
122
                                  y2 = kk - ll;
 
123
                                  z2 = ll;
 
124
                                  
 
125
                                  stemp[ao1][ao2] += pfac*OIX[x1][x2]*OIY[y1][y2]*OIZ[z1][z2];
 
126
                                  
 
127
                                  ao2++;
 
128
                              }
 
129
                          }
 
130
                          ao1++;
 
131
                      }
 
132
                  } /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
 
133
              }
 
134
          } /*--- end primitive contraction ---*/
 
135
 
 
136
          /*--- Normalize the contracted integrals and put them into S12 ---*/
 
137
          ptr1 = GTOs.bf_norm[l1];
 
138
          ptr2 = GTOs.bf_norm[l2];
 
139
          for(ii=0; ii<n1; ii++) {
 
140
              norm1 = ptr1[ii];
 
141
              for(jj=0; jj<n2; jj++) {
 
142
                  S12_AO[bf1+ii][bf2+jj] = stemp[ii][jj]*norm1*ptr2[jj];
 
143
              }
 
144
          }
 
145
          
 
146
      }
 
147
  }   /*--- This shell pair is done ---*/
 
148
 
 
149
  delete_matrix(OIX);
 
150
  delete_matrix(OIY);
 
151
  delete_matrix(OIZ);
 
152
  delete_matrix(stemp);
 
153
 
 
154
  /*--- transform to SO basis ---*/
 
155
  tmpmat = create_matrix(num_so,Oldcalc.num_ao);
 
156
  usotao_FLOAT = convert_matrix(usotao,num_so,num_ao,0);
 
157
  oldusotao_FLOAT_transp = convert_matrix(Oldcalc.usotao,Oldcalc.num_so,Oldcalc.num_ao,1);
 
158
  matrix_mult(usotao_FLOAT,num_so,num_ao,S12_AO,num_ao,Oldcalc.num_ao,tmpmat);
 
159
  delete_matrix(S12_AO);
 
160
  S12 = create_matrix(num_so,Oldcalc.num_so);
 
161
  matrix_mult(tmpmat,num_so,Oldcalc.num_ao,oldusotao_FLOAT_transp,Oldcalc.num_ao,Oldcalc.num_so,S12);
 
162
  delete_matrix(tmpmat);
 
163
 
 
164
  return S12;
 
165
}
 
166
 
 
167
}} // namespace psi::input