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_lbgbuild_h
00029 #define _chemistry_qc_scf_lbgbuild_h
00030
00031 #ifdef __GNUC__
00032 #pragma interface
00033 #endif
00034
00035 #include <chemistry/qc/scf/gbuild.h>
00036
00037 namespace sc {
00038
00039 template<class T>
00040 class LocalLBGBuild : public GBuild<T> {
00041 protected:
00042 Ref<MessageGrp> grp_;
00043 Ref<TwoBodyInt> tbi_;
00044 Ref<Integral> integral_;
00045 Ref<GaussianBasisSet> gbs_;
00046 signed char *pmax;
00047
00048 public:
00049 LocalLBGBuild(T& t, const Ref<TwoBodyInt>& tbi, const Ref<Integral>& ints,
00050 const Ref<GaussianBasisSet>& bs, const Ref<MessageGrp>& g,
00051 signed char *pm) :
00052 GBuild<T>(t), grp_(g), tbi_(tbi), integral_(ints), gbs_(bs), pmax(pm) {}
00053 ~LocalLBGBuild() {}
00054
00055 void build_gmat(double accuracy) {
00056 tim_enter("ao_gmat");
00057 tim_set_default("quartet");
00058
00059 double tnint=0;
00060 int tol = (int) (log(accuracy)/log(2.0));
00061 int me=grp_->me();
00062 int nproc = grp_->n();
00063
00064 Ref<PetiteList> rpl = integral_->petite_list();
00065
00066
00067 GaussianBasisSet& gbs = *gbs_.pointer();
00068 PetiteList& pl = *rpl.pointer();
00069 TwoBodyInt& tbi = *tbi_.pointer();
00070
00071 tbi.set_redundant(0);
00072 const double *intbuf = tbi.buffer();
00073
00074 int inds[4];
00075
00076
00077 if (me==0) {
00078 int i;
00079 for (i=0; i < gbs.nshell(); i++) {
00080 if (!pl.in_p1(i))
00081 continue;
00082
00083 inds[0]=i;
00084
00085 for (int j=0; j <= i; j++) {
00086 int oij = i_offset(i)+j;
00087
00088 if (!pl.in_p2(oij))
00089 continue;
00090
00091 inds[1]=j;
00092
00093 tim_enter_default();
00094 int from;
00095 grp_->recvt(2323, &from, 1);
00096 grp_->sendt(from, 3232, inds, 4);
00097 tim_exit_default();
00098 }
00099 }
00100
00101
00102 inds[0] = inds[1] = inds[2] = inds[3] = -1;
00103
00104 for (i=1; i < nproc; i++) {
00105 int from;
00106 grp_->recvt(2323, &from, 1);
00107 grp_->sendt(from, 3232, inds, 4);
00108 }
00109 }
00110
00111
00112 else {
00113 do {
00114 grp_->sendt(0, 2323, &me, 1);
00115 grp_->recvt(3232, inds, 4);
00116
00117 int i=inds[0];
00118 int j=inds[1];
00119
00120 if (i < 0)
00121 break;
00122
00123 int fi=gbs.shell_to_function(i);
00124 int ni=gbs(i).nfunction();
00125
00126 int oij = i_offset(i)+j;
00127
00128 int fj=gbs.shell_to_function(j);
00129 int nj=gbs(j).nfunction();
00130 int pmaxij = pmax[oij];
00131
00132 for (int k=0; k <= i; k++) {
00133 int fk=gbs.shell_to_function(k);
00134 int nk=gbs(k).nfunction();
00135
00136 int pmaxijk=pmaxij, ptmp;
00137 if ((ptmp=pmax[i_offset(i)+k]-2) > pmaxijk) pmaxijk=ptmp;
00138 if ((ptmp=pmax[ij_offset(j,k)]-2) > pmaxijk) pmaxijk=ptmp;
00139
00140 int okl = i_offset(k);
00141 for (int l=0; l <= (k==i?j:k); l++,okl++) {
00142
00143 int pmaxijkl = pmaxijk;
00144 if ((ptmp=pmax[okl]) > pmaxijkl) pmaxijkl=ptmp;
00145 if ((ptmp=pmax[i_offset(i)+l]-2) > pmaxijkl) pmaxijkl=ptmp;
00146 if ((ptmp=pmax[ij_offset(j,l)]-2) > pmaxijkl) pmaxijkl=ptmp;
00147
00148
00149 if (tbi.log2_shell_bound(i,j,k,l)+pmaxijkl < tol)
00150 continue;
00151
00152 int qijkl = pl.in_p4(oij,okl,i,j,k,l);
00153 if (!qijkl)
00154 continue;
00155
00156 tim_enter_default();
00157 tbi.compute_shell(i,j,k,l);
00158 tim_exit_default();
00159
00160 int e12 = (i==j);
00161 int e34 = (k==l);
00162 int e13e24 = (i==k) && (j==l);
00163 int e_any = e12||e34||e13e24;
00164
00165 int fl=gbs.shell_to_function(l);
00166 int nl=gbs(l).nfunction();
00167
00168 int ii,jj,kk,ll;
00169 int I,J,K,L;
00170 int index=0;
00171
00172 for (I=0, ii=fi; I < ni; I++, ii++) {
00173 for (J=0, jj=fj; J <= (e12 ? I : nj-1); J++, jj++) {
00174 for (K=0, kk=fk; K <= (e13e24 ? I : nk-1); K++, kk++) {
00175
00176 int lend = (e34 ? ((e13e24)&&(K==I) ? J : K)
00177 : ((e13e24)&&(K==I)) ? J : nl-1);
00178 for (L=0, ll=fl; L <= lend; L++, ll++, index++) {
00179 double pki_int = intbuf[index];
00180
00181 if ((pki_int>0?pki_int:-pki_int) < 1.0e-15)
00182 continue;
00183
00184 if (qijkl > 1)
00185 pki_int *= qijkl;
00186
00187 if (e_any) {
00188 int ij,kl;
00189 double val;
00190
00191 if (jj == kk) {
00192
00193
00194
00195
00196
00197
00198 if (ii == jj || kk == ll) {
00199 ij = i_offset(ii)+jj;
00200 kl = i_offset(kk)+ll;
00201 val = (ij==kl) ? 0.5*pki_int : pki_int;
00202
00203 contribution.cont5(ij,kl,val);
00204
00205 } else {
00206
00207
00208
00209
00210
00211
00212
00213 ij = i_offset(ii)+jj;
00214 kl = i_offset(kk)+ll;
00215 val = (ij==kl) ? 0.5*pki_int : pki_int;
00216
00217 contribution.cont4(ij,kl,val);
00218
00219
00220
00221
00222
00223
00224
00225
00226 ij = ij_offset(ii,ll);
00227 kl = ij_offset(kk,jj);
00228 val = (ij==kl) ? 0.5*pki_int : pki_int;
00229
00230 contribution.cont3(ij,kl,val);
00231 }
00232 } else if (ii == kk || jj == ll) {
00233
00234
00235
00236
00237
00238
00239
00240 ij = i_offset(ii)+jj;
00241 kl = i_offset(kk)+ll;
00242 val = (ij==kl) ? 0.5*pki_int : pki_int;
00243
00244 contribution.cont4(ij,kl,val);
00245
00246
00247
00248
00249
00250
00251
00252
00253 ij = ij_offset(ii,kk);
00254 kl = ij_offset(jj,ll);
00255 val = (ij==kl) ? 0.5*pki_int : pki_int;
00256
00257 contribution.cont3(ij,kl,val);
00258
00259 } else {
00260
00261
00262
00263
00264
00265 ij = i_offset(ii)+jj;
00266 kl = i_offset(kk)+ll;
00267 val = (ij==kl) ? 0.5*pki_int : pki_int;
00268
00269 contribution.cont1(ij,kl,val);
00270
00271
00272
00273
00274
00275
00276 ij = ij_offset(ii,kk);
00277 kl = ij_offset(jj,ll);
00278 val = (ij==kl) ? 0.5*pki_int : pki_int;
00279
00280 contribution.cont2(ij,kl,val);
00281
00282 if ((ii != jj) && (kk != ll)) {
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292 ij = ij_offset(ii,ll);
00293 kl = ij_offset(kk,jj);
00294
00295 contribution.cont2(ij,kl,val);
00296 }
00297 }
00298 } else {
00299 if (jj == kk) {
00300
00301
00302
00303
00304
00305
00306
00307 contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
00308
00309
00310
00311
00312
00313
00314
00315
00316 contribution.cont3(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
00317
00318 } else if (ii == kk || jj == ll) {
00319
00320
00321
00322
00323
00324
00325
00326 contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
00327
00328
00329
00330
00331
00332
00333
00334
00335 contribution.cont3(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
00336
00337 } else {
00338
00339
00340
00341
00342
00343 contribution.cont1(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
00344
00345
00346
00347
00348
00349
00350 contribution.cont2(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
00351
00352
00353
00354
00355
00356
00357 contribution.cont2(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
00358 }
00359 }
00360 }
00361 }
00362 }
00363 }
00364
00365 tnint += (double) ni*nj*nk*nl;
00366 }
00367 }
00368 } while (inds[0] > -1);
00369 }
00370
00371 grp_->sum(&tnint, 1, 0, 0);
00372 ExEnv::out0() << indent << scprintf("%20.0f integrals\n", tnint);
00373
00374 tim_exit("ao_gmat");
00375 }
00376
00377 };
00378
00379 }
00380
00381 #endif
00382
00383
00384
00385
00386