← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/mpqc/oneiric

« back to all changes in this revision

Viewing changes to src/lib/chemistry/qc/mbptr12/multipole_ints.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2005-11-27 11:41:49 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20051127114149-zgz9r3gk50w8ww2q
Tags: 2.3.0-1
http://bugs.debian.org/339232
* New upstream release.
* debian/rules (SONAME): Activate awk snippet for automatic so-name
  detection again, resulting in a bump to `7' and making a `c2a' for
  the C++ allocator change unnecessary; closes: #339232.
* debian/patches/00list (08_gcc-4.0_fixes): Removed, no longer needed.
* debian/rules (test): Remove workarounds, do not abort build if tests
  fail.
* debian/ref: Removed.
* debian/control.in (libsc): Added Conflict against libsc6c2.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/chemistry/qc/mbptr12/multipole_ints.cc
38
38
using namespace sc;
39
39
 
40
40
void
41
 
R12IntEvalInfo::compute_multipole_ints(RefSymmSCMatrix& MX, RefSymmSCMatrix& MY, RefSymmSCMatrix& MZ,
42
 
                                       RefSymmSCMatrix& MXX, RefSymmSCMatrix& MYY, RefSymmSCMatrix& MZZ)
 
41
R12IntEvalInfo::compute_multipole_ints(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2,
 
42
                                       RefSCMatrix& MX, RefSCMatrix& MY, RefSCMatrix& MZ,
 
43
                                       RefSCMatrix& MXX, RefSCMatrix& MYY, RefSCMatrix& MZZ)
43
44
{
44
 
  Ref<PetiteList> pl = ref_->integral()->petite_list();
45
 
  int nshell = bs_->nshell();
46
 
 
47
 
  // Have to convert scf_vec_ to unblocked form, get_subblock, and convert back to blocked form
48
 
  // because get_subblock doesn't like blocked matrices
49
 
  RefSCDimension aodim = pl->AO_basisdim();
50
 
  RefSCDimension modim;
51
 
  {
52
 
    int nmo = scf_vec_.rowdim().n();
53
 
    int blksize[1];  blksize[0] = nmo;
54
 
    modim = new SCDimension(nmo,1,blksize,"Sorted MO dimension");
55
 
  }
56
 
  RefSCMatrix unblvec(modim,aodim,bs_->matrixkit());
57
 
  unblvec->convert(scf_vec_.pointer());
58
 
  RefSCMatrix subblock = unblvec.get_subblock(nfzc_,nocc_-1,0,bs_->nbasis()-1);
59
 
  RefSCDimension occactdim;
60
 
  {
61
 
    int blksize[1]; blksize[0] = nocc_ - nfzc_;
62
 
    occactdim = new SCDimension(nocc_-nfzc_,1,blksize,"Active occupied MO dimension");
63
 
    occactdim->blocks()->set_subdim(0,new SCDimension(nocc_-nfzc_,"occact MO dim"));
64
 
  }
65
 
  RefSCMatrix OccAct_Vec(occactdim,aodim,bs_->so_matrixkit());
66
 
  OccAct_Vec->convert(subblock.pointer());
67
 
  subblock = 0;
68
 
  unblvec=0;
69
 
 
70
 
  Ref<OneBodyInt> m1_ints = integral_->dipole(0);
71
 
 
72
 
  // form AO dipole moment matrices
73
 
  RefSymmSCMatrix mx(pl->AO_basisdim(), bs_->so_matrixkit());
74
 
  RefSymmSCMatrix my(pl->AO_basisdim(), bs_->so_matrixkit());
75
 
  RefSymmSCMatrix mz(pl->AO_basisdim(), bs_->so_matrixkit());
 
45
  if (!space1->integral()->equiv(space2->integral()))
 
46
    throw ProgrammingError("two MOIndexSpaces use incompatible Integral factories");
 
47
  const Ref<GaussianBasisSet> bs1 = space1->basis();
 
48
  const Ref<GaussianBasisSet> bs2 = space2->basis();
 
49
  const bool bs1_eq_bs2 = (bs1 == bs2);
 
50
  int nshell1 = bs1->nshell();
 
51
  int nshell2 = bs2->nshell();
 
52
 
 
53
  RefSCMatrix vec1t = space1->coefs().t();
 
54
  RefSCMatrix vec2 = space2->coefs();
 
55
 
 
56
  Ref<Integral> localints = space1->integral()->clone();
 
57
  localints->set_basis(bs1,bs2);
 
58
 
 
59
  Ref<OneBodyInt> m1_ints = localints->dipole(0);
 
60
  Ref<OneBodyInt> m2_ints = localints->quadrupole(0);
 
61
 
 
62
  // form AO moment matrices
 
63
  RefSCDimension aodim1 = vec1t.coldim();
 
64
  RefSCDimension aodim2 = vec2.rowdim();
 
65
  Ref<SCMatrixKit> aokit = bs1->so_matrixkit();
 
66
  RefSCMatrix mx(aodim1, aodim2, aokit);
 
67
  RefSCMatrix my(aodim1, aodim2, aokit);
 
68
  RefSCMatrix mz(aodim1, aodim2, aokit);
 
69
  RefSCMatrix mxx(aodim1, aodim2, aokit);
 
70
  RefSCMatrix myy(aodim1, aodim2, aokit);
 
71
  RefSCMatrix mzz(aodim1, aodim2, aokit);
76
72
  mx.assign(0.0);
77
73
  my.assign(0.0);
78
74
  mz.assign(0.0);
 
75
  mxx.assign(0.0);
 
76
  myy.assign(0.0);
 
77
  mzz.assign(0.0);
 
78
  
 
79
  for(int sh1=0; sh1<nshell1; sh1++) {
 
80
    int bf1_offset = bs1->shell_to_function(sh1);
 
81
    int nbf1 = bs1->shell(sh1).nfunction();
 
82
 
 
83
    int sh2max;
 
84
    if (bs1_eq_bs2)
 
85
      sh2max = sh1;
 
86
    else
 
87
      sh2max = nshell2-1;
79
88
    
80
 
  for(int sh1=0; sh1<nshell; sh1++) {
81
 
    int bf1_offset = bs_->shell_to_function(sh1);
82
 
    int nbf1 = bs_->shell(sh1).nfunction();
83
 
    for(int sh2=0; sh2<=sh1; sh2++) {
84
 
      int bf2_offset = bs_->shell_to_function(sh2);
85
 
      int nbf2 = bs_->shell(sh2).nfunction();
 
89
    for(int sh2=0; sh2<=sh2max; sh2++) {
 
90
      int bf2_offset = bs2->shell_to_function(sh2);
 
91
      int nbf2 = bs2->shell(sh2).nfunction();
86
92
      
87
93
      m1_ints->compute_shell(sh1,sh2);
88
94
      const double *m1intsptr = m1_ints->buffer();
 
95
 
 
96
      m2_ints->compute_shell(sh1,sh2);
 
97
      const double *m2intsptr = m2_ints->buffer();
 
98
 
89
99
      int bf1_index = bf1_offset;
90
 
      for(int bf1=0; bf1<nbf1; bf1++, bf1_index++, m1intsptr+=3*nbf2) {
 
100
      for(int bf1=0; bf1<nbf1; bf1++, bf1_index++, m1intsptr+=3*nbf2, m2intsptr+=6*nbf2) {
91
101
        int bf2_index = bf2_offset;
92
102
        const double *ptr1 = m1intsptr;
 
103
        const double *ptr2 = m2intsptr;
93
104
        int bf2max;
94
 
        if (sh1 != sh2)
 
105
        if (bs1_eq_bs2 && sh1 == sh2)
 
106
          bf2max = bf1;
 
107
        else
95
108
          bf2max = nbf2-1;
96
 
        else
97
 
          bf2max = bf1;
98
109
        for(int bf2=0; bf2<=bf2max; bf2++, bf2_index++) {
 
110
 
99
111
          mx.set_element(bf1_index, bf2_index, *(ptr1++));
100
112
          my.set_element(bf1_index, bf2_index, *(ptr1++));
101
113
          mz.set_element(bf1_index, bf2_index, *(ptr1++));
102
 
        }
 
114
 
 
115
          mxx.set_element(bf1_index, bf2_index, *(ptr2++));
 
116
          ptr2 += 2;
 
117
          myy.set_element(bf1_index, bf2_index, *(ptr2++));
 
118
          ptr2++;
 
119
          mzz.set_element(bf1_index, bf2_index, *(ptr2++));
 
120
 
 
121
        }
103
122
      }
104
123
    }
105
124
  }
106
125
 
107
126
  // and clean up a bit
108
127
  m1_ints = 0;
 
128
  m2_ints = 0;
 
129
 
 
130
  // Symmetrize matrices, if necessary
 
131
  if (bs1_eq_bs2) {
 
132
 
 
133
    const int nbasis = bs1->nbasis();
 
134
    
 
135
    for(int bf1=0; bf1<nbasis; bf1++)
 
136
      for(int bf2=0; bf2<=bf1; bf2++) {
 
137
        mx(bf2,bf1) = mx(bf1,bf2);
 
138
        my(bf2,bf1) = my(bf1,bf2);
 
139
        mz(bf2,bf1) = mz(bf1,bf2);
 
140
        mxx(bf2,bf1) = mxx(bf1,bf2);
 
141
        myy(bf2,bf1) = myy(bf1,bf2);
 
142
        mzz(bf2,bf1) = mzz(bf1,bf2);
 
143
      }
 
144
 
 
145
  }
 
146
      
109
147
 
110
148
  // finally, transform
111
 
  MX = bs_->so_matrixkit()->symmmatrix(occactdim);
112
 
  MY = bs_->so_matrixkit()->symmmatrix(occactdim);
113
 
  MZ = bs_->so_matrixkit()->symmmatrix(occactdim);
114
 
  MX.assign(0.0);
115
 
  MY.assign(0.0);
116
 
  MZ.assign(0.0);
117
 
  MX.accumulate_transform(OccAct_Vec,mx);
 
149
  MX = vec1t * mx * vec2;
 
150
  MY = vec1t * my * vec2;
 
151
  MZ = vec1t * mz * vec2;
 
152
  MXX = vec1t * mxx * vec2;
 
153
  MYY = vec1t * myy * vec2;
 
154
  MZZ = vec1t * mzz * vec2;
 
155
  
 
156
  // and clean up a bit
118
157
  mx = 0;
119
 
  MY.accumulate_transform(OccAct_Vec,my);
120
158
  my = 0;
121
 
  MZ.accumulate_transform(OccAct_Vec,mz);
122
159
  mz = 0;
123
 
 
124
 
  // same for quadrupole integrals
125
 
  Ref<OneBodyInt> m2_ints = integral_->quadrupole(0);
126
 
  RefSymmSCMatrix mxx(pl->AO_basisdim(), bs_->so_matrixkit());
127
 
  RefSymmSCMatrix myy(pl->AO_basisdim(), bs_->so_matrixkit());
128
 
  RefSymmSCMatrix mzz(pl->AO_basisdim(), bs_->so_matrixkit());
129
 
  mxx.assign(0.0);
130
 
  myy.assign(0.0);
131
 
  mzz.assign(0.0);
132
 
    
133
 
  for(int sh1=0; sh1<nshell; sh1++) {
134
 
    int bf1_offset = bs_->shell_to_function(sh1);
135
 
    int nbf1 = bs_->shell(sh1).nfunction();
136
 
    for(int sh2=0; sh2<=sh1; sh2++) {
137
 
      int bf2_offset = bs_->shell_to_function(sh2);
138
 
      int nbf2 = bs_->shell(sh2).nfunction();
139
 
      
140
 
      m2_ints->compute_shell(sh1,sh2);
141
 
      const double *m2intsptr = m2_ints->buffer();
142
 
      int bf1_index = bf1_offset;
143
 
      for(int bf1=0; bf1<nbf1; bf1++, bf1_index++, m2intsptr+=6*nbf2) {
144
 
        int bf2_index = bf2_offset;
145
 
        const double *ptr2 = m2intsptr;
146
 
        int bf2max;
147
 
        if (sh1 != sh2)
148
 
          bf2max = nbf2-1;
149
 
        else
150
 
          bf2max = bf1;
151
 
        for(int bf2=0; bf2<=bf2max; bf2++, bf2_index++) {
152
 
          mxx.set_element(bf1_index, bf2_index, *(ptr2++));
153
 
          ptr2 += 2;
154
 
          myy.set_element(bf1_index, bf2_index, *(ptr2++));
155
 
          ptr2++;
156
 
          mzz.set_element(bf1_index, bf2_index, *(ptr2++));
157
 
        }
158
 
      }
159
 
    }
160
 
  }
161
 
  // and clean up a bit
162
 
  m2_ints = 0;
163
 
  pl = 0;
164
 
 
165
 
  // transform
166
 
  MXX = bs_->so_matrixkit()->symmmatrix(occactdim);
167
 
  MYY = bs_->so_matrixkit()->symmmatrix(occactdim);
168
 
  MZZ = bs_->so_matrixkit()->symmmatrix(occactdim);
169
 
  MXX.assign(0.0);
170
 
  MYY.assign(0.0);
171
 
  MZZ.assign(0.0);
172
 
  MXX.accumulate_transform(OccAct_Vec,mxx);
173
160
  mxx = 0;
174
 
  MYY.accumulate_transform(OccAct_Vec,myy);
175
161
  myy = 0;
176
 
  MZZ.accumulate_transform(OccAct_Vec,mzz);
177
162
  mzz = 0;
178
163
 
179
 
  if (debug_ > 1) {
180
 
    MX.print("mu(X) in active occupied MO basis");
181
 
    MY.print("mu(Y) in active occupied MO basis");
182
 
    MZ.print("mu(Z) in active occupied MO basis");
183
 
    MXX.print("mu(XX) in active occupied MO basis");
184
 
    MYY.print("mu(YY) in active occupied MO basis");
185
 
    MZZ.print("mu(ZZ) in active occupied MO basis");
186
 
  }
187
 
 
188
 
  OccAct_Vec = 0;
189
 
 
190
 
  return;
 
164
  //if (debug_ > 1) {
 
165
  //  MX.print("mu(X)");
 
166
  //  MY.print("mu(Y)");
 
167
  //  MZ.print("mu(Z)");
 
168
  //  MXX.print("mu(XX)");
 
169
  //  MYY.print("mu(YY)");
 
170
  //  MZZ.print("mu(ZZ)");
 
171
  //}
 
172
}
 
173
 
 
174
 
 
175
void
 
176
R12IntEvalInfo::compute_overlap_ints(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2,
 
177
                                     RefSCMatrix& S)
 
178
{
 
179
  if (!space1->integral()->equiv(space2->integral()))
 
180
    throw ProgrammingError("two MOIndexSpaces use incompatible Integral factories");
 
181
  const Ref<GaussianBasisSet> bs1 = space1->basis();
 
182
  const Ref<GaussianBasisSet> bs2 = space2->basis();
 
183
  const bool bs1_eq_bs2 = (bs1 == bs2);
 
184
  int nshell1 = bs1->nshell();
 
185
  int nshell2 = bs2->nshell();
 
186
 
 
187
  RefSCMatrix vec1t = space1->coefs().t();
 
188
  RefSCMatrix vec2 = space2->coefs();
 
189
 
 
190
  Ref<Integral> localints = space1->integral()->clone();
 
191
  localints->set_basis(bs1,bs2);
 
192
 
 
193
  Ref<OneBodyInt> ov_ints = localints->overlap();
 
194
 
 
195
  // form AO moment matrices
 
196
  RefSCDimension aodim1 = vec1t.coldim();
 
197
  RefSCDimension aodim2 = vec2.rowdim();
 
198
  Ref<SCMatrixKit> aokit = bs1->so_matrixkit();
 
199
  RefSCMatrix s(aodim1, aodim2, aokit);
 
200
  s.assign(0.0);
 
201
 
 
202
  for(int sh1=0; sh1<nshell1; sh1++) {
 
203
    int bf1_offset = bs1->shell_to_function(sh1);
 
204
    int nbf1 = bs1->shell(sh1).nfunction();
 
205
 
 
206
    int sh2max;
 
207
    if (bs1_eq_bs2)
 
208
      sh2max = sh1;
 
209
    else
 
210
      sh2max = nshell2-1;
 
211
 
 
212
    for(int sh2=0; sh2<=sh2max; sh2++) {
 
213
      int bf2_offset = bs2->shell_to_function(sh2);
 
214
      int nbf2 = bs2->shell(sh2).nfunction();
 
215
 
 
216
      ov_ints->compute_shell(sh1,sh2);
 
217
      const double *ovintsptr = ov_ints->buffer();
 
218
 
 
219
      int bf1_index = bf1_offset;
 
220
      for(int bf1=0; bf1<nbf1; bf1++, bf1_index++, ovintsptr+=nbf2) {
 
221
        int bf2_index = bf2_offset;
 
222
        const double *ptr = ovintsptr;
 
223
        int bf2max;
 
224
        if (bs1_eq_bs2 && sh1 == sh2)
 
225
          bf2max = bf1;
 
226
        else
 
227
          bf2max = nbf2-1;
 
228
        for(int bf2=0; bf2<=bf2max; bf2++, bf2_index++) {
 
229
 
 
230
          s.set_element(bf1_index, bf2_index, *(ptr++));
 
231
 
 
232
        }
 
233
      }
 
234
    }
 
235
  }
 
236
 
 
237
  // and clean up a bit
 
238
  ov_ints = 0;
 
239
 
 
240
  // Symmetrize matrices, if necessary
 
241
  if (bs1_eq_bs2) {
 
242
 
 
243
    const int nbasis = bs1->nbasis();
 
244
 
 
245
    for(int bf1=0; bf1<nbasis; bf1++)
 
246
      for(int bf2=0; bf2<=bf1; bf2++) {
 
247
        s(bf2,bf1) = s(bf1,bf2);
 
248
      }
 
249
 
 
250
  }
 
251
 
 
252
 
 
253
    // finally, transform
 
254
    S = vec1t * s * vec2;
 
255
 
 
256
    // and clean up a bit
 
257
    s = 0;
 
258
 
 
259
    //if (debug_ > 1) {
 
260
    //  S.print("Overlap");
 
261
    //}
191
262
}
192
263
 
193
264
///////////////////////////////////////////////////////////////