lgbuild.h
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifndef _chemistry_qc_scf_lgbuild_h
00029 #define _chemistry_qc_scf_lgbuild_h
00030
00031 #ifdef __GNUC__
00032 #pragma interface
00033 #endif
00034
00035 #undef SCF_CHECK_INTS
00036 #undef SCF_CHECK_BOUNDS
00037 #undef SCF_DONT_USE_BOUNDS
00038
00039 #include <scconfig.h>
00040 #include <chemistry/qc/scf/gbuild.h>
00041
00042 namespace sc {
00043
00044 template<class T>
00045 class LocalGBuild : public GBuild<T> {
00046 public:
00047 double tnint;
00048
00049 protected:
00050 MessageGrp *grp_;
00051 TwoBodyInt *tbi_;
00052 GaussianBasisSet *gbs_;
00053 PetiteList *rpl_;
00054
00055 signed char * restrictxx pmax;
00056 int threadno_;
00057 int nthread_;
00058 double accuracy_;
00059
00060 public:
00061 LocalGBuild(T& t, const Ref<TwoBodyInt>& tbi, const Ref<PetiteList>& rpl,
00062 const Ref<GaussianBasisSet>& bs, const Ref<MessageGrp>& g,
00063 signed char *pm, double acc, int nt=1, int tn=0) :
00064 GBuild<T>(t),
00065 pmax(pm), nthread_(nt), threadno_(tn), accuracy_(acc)
00066 {
00067 grp_ = g.pointer();
00068 tbi_ = tbi.pointer();
00069 rpl_ = rpl.pointer();
00070 gbs_ = bs.pointer();
00071 }
00072 ~LocalGBuild() {}
00073
00074 void run() {
00075 int tol = (int) (log(accuracy_)/log(2.0));
00076 int me=grp_->me();
00077 int nproc = grp_->n();
00078
00079
00080 GaussianBasisSet& gbs = *gbs_;
00081 PetiteList& pl = *rpl_;
00082 TwoBodyInt& tbi = *tbi_;
00083
00084 tbi.set_redundant(0);
00085 const double *intbuf = tbi.buffer();
00086
00087 tnint=0;
00088 sc_int_least64_t threadind=0;
00089 sc_int_least64_t ijklind=0;
00090
00091 for (int i=0; i < gbs.nshell(); i++) {
00092 if (!pl.in_p1(i))
00093 continue;
00094
00095 int fi=gbs.shell_to_function(i);
00096 int ni=gbs(i).nfunction();
00097
00098 for (int j=0; j <= i; j++) {
00099 int oij = i_offset(i)+j;
00100
00101 if (!pl.in_p2(oij))
00102 continue;
00103
00104 int fj=gbs.shell_to_function(j);
00105 int nj=gbs(j).nfunction();
00106 int pmaxij = pmax[oij];
00107
00108 for (int k=0; k <= i; k++, ijklind++) {
00109 if (ijklind%nproc != me)
00110 continue;
00111
00112 threadind++;
00113 if (threadind % nthread_ != threadno_)
00114 continue;
00115
00116 int fk=gbs.shell_to_function(k);
00117 int nk=gbs(k).nfunction();
00118
00119 int pmaxijk=pmaxij, ptmp;
00120 if ((ptmp=pmax[i_offset(i)+k]-1) > pmaxijk) pmaxijk=ptmp;
00121 if ((ptmp=pmax[ij_offset(j,k)]-1) > pmaxijk) pmaxijk=ptmp;
00122
00123 int okl = i_offset(k);
00124 for (int l=0; l <= (k==i?j:k); l++,okl++) {
00125 int pmaxijkl = pmaxijk;
00126 if ((ptmp=pmax[okl]) > pmaxijkl) pmaxijkl=ptmp;
00127 if ((ptmp=pmax[i_offset(i)+l]-1) > pmaxijkl) pmaxijkl=ptmp;
00128 if ((ptmp=pmax[ij_offset(j,l)]-1) > pmaxijkl) pmaxijkl=ptmp;
00129
00130 int qijkl = pl.in_p4(oij,okl,i,j,k,l);
00131 if (!qijkl)
00132 continue;
00133
00134 #ifdef SCF_CHECK_BOUNDS
00135 double intbound = pow(2.0,double(tbi.log2_shell_bound(i,j,k,l)));
00136 double pbound = pow(2.0,double(pmaxijkl));
00137 intbound *= qijkl;
00138 GBuild<T>::contribution.set_bound(intbound, pbound);
00139 #else
00140 # ifndef SCF_DONT_USE_BOUNDS
00141 if (tbi.log2_shell_bound(i,j,k,l)+pmaxijkl < tol)
00142 continue;
00143 # endif
00144 #endif
00145
00146
00147 tbi.compute_shell(i,j,k,l);
00148
00149
00150 int e12 = (i==j);
00151 int e34 = (k==l);
00152 int e13e24 = (i==k) && (j==l);
00153 int e_any = e12||e34||e13e24;
00154
00155 int fl=gbs.shell_to_function(l);
00156 int nl=gbs(l).nfunction();
00157
00158 int ii,jj,kk,ll;
00159 int I,J,K,L;
00160 int index=0;
00161
00162 for (I=0, ii=fi; I < ni; I++, ii++) {
00163 for (J=0, jj=fj; J <= (e12 ? I : nj-1); J++, jj++) {
00164 for (K=0, kk=fk; K <= (e13e24 ? I : nk-1); K++, kk++) {
00165 int lend = (e34 ? ((e13e24)&&(K==I) ? J : K)
00166 : ((e13e24)&&(K==I)) ? J : nl-1);
00167
00168 for (L=0, ll=fl; L <= lend; L++, ll++, index++) {
00169
00170 double pki_int = intbuf[index];
00171
00172 if ((pki_int>0?pki_int:-pki_int) < 1.0e-15)
00173 continue;
00174
00175 #ifdef SCF_CHECK_INTS
00176 if (isnan(pki_int))
00177 abort();
00178 #endif
00179
00180 if (qijkl > 1)
00181 pki_int *= qijkl;
00182
00183 if (e_any) {
00184 int ij,kl;
00185 double val;
00186
00187 if (jj == kk) {
00188
00189
00190
00191
00192
00193
00194 if (ii == jj || kk == ll) {
00195 ij = i_offset(ii)+jj;
00196 kl = i_offset(kk)+ll;
00197 val = (ij==kl) ? 0.5*pki_int : pki_int;
00198
00199 GBuild<T>::contribution.cont5(ij,kl,val);
00200
00201 } else {
00202
00203
00204
00205
00206
00207
00208
00209 ij = i_offset(ii)+jj;
00210 kl = i_offset(kk)+ll;
00211 val = (ij==kl) ? 0.5*pki_int : pki_int;
00212
00213 GBuild<T>::contribution.cont4(ij,kl,val);
00214
00215
00216
00217
00218
00219
00220
00221
00222 ij = ij_offset(ii,ll);
00223 kl = ij_offset(kk,jj);
00224 val = (ij==kl) ? 0.5*pki_int : pki_int;
00225
00226 GBuild<T>::contribution.cont3(ij,kl,val);
00227 }
00228 } else if (ii == kk || jj == ll) {
00229
00230
00231
00232
00233
00234
00235
00236 ij = i_offset(ii)+jj;
00237 kl = i_offset(kk)+ll;
00238 val = (ij==kl) ? 0.5*pki_int : pki_int;
00239
00240 GBuild<T>::contribution.cont4(ij,kl,val);
00241
00242
00243
00244
00245
00246
00247
00248
00249 ij = ij_offset(ii,kk);
00250 kl = ij_offset(jj,ll);
00251 val = (ij==kl) ? 0.5*pki_int : pki_int;
00252
00253 GBuild<T>::contribution.cont3(ij,kl,val);
00254
00255 } else {
00256
00257
00258
00259
00260
00261 ij = i_offset(ii)+jj;
00262 kl = i_offset(kk)+ll;
00263 val = (ij==kl) ? 0.5*pki_int : pki_int;
00264
00265 GBuild<T>::contribution.cont1(ij,kl,val);
00266
00267
00268
00269
00270
00271
00272 ij = ij_offset(ii,kk);
00273 kl = ij_offset(jj,ll);
00274 val = (ij==kl) ? 0.5*pki_int : pki_int;
00275
00276 GBuild<T>::contribution.cont2(ij,kl,val);
00277
00278 if ((ii != jj) && (kk != ll)) {
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288 ij = ij_offset(ii,ll);
00289 kl = ij_offset(kk,jj);
00290
00291 GBuild<T>::contribution.cont2(ij,kl,val);
00292 }
00293 }
00294 } else {
00295 if (jj == kk) {
00296
00297
00298
00299
00300
00301
00302
00303 GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
00304
00305
00306
00307
00308
00309
00310
00311
00312 GBuild<T>::contribution.cont3(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
00313
00314 } else if (ii == kk || jj == ll) {
00315
00316
00317
00318
00319
00320
00321
00322 GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
00323
00324
00325
00326
00327
00328
00329
00330
00331 GBuild<T>::contribution.cont3(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
00332
00333 } else {
00334
00335
00336
00337
00338
00339 GBuild<T>::contribution.cont1(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
00340
00341
00342
00343
00344
00345
00346 GBuild<T>::contribution.cont2(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
00347
00348
00349
00350
00351
00352
00353 GBuild<T>::contribution.cont2(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
00354 }
00355 }
00356 }
00357 }
00358 }
00359 }
00360
00361 tnint += (double) ni*nj*nk*nl;
00362 }
00363 }
00364 }
00365 }
00366 }
00367 };
00368
00369 }
00370
00371 #endif
00372
00373
00374
00375
00376
Generated at Sat Dec 18 15:14:21 2004 for MPQC
2.2.3 using the documentation package Doxygen
1.3.7-20040617.