lbgbuild.h
00001 //
00002 // lbgbuild.h --- definitino of the load-balanced local G matrix builder
00003 //
00004 // Copyright (C) 1996 Limit Point Systems, Inc.
00005 //
00006 // Author: Edward Seidl <seidl@janed.com>
00007 // Maintainer: LPS
00008 //
00009 // This file is part of the SC Toolkit.
00010 //
00011 // The SC Toolkit is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Library General Public License as published by
00013 // the Free Software Foundation; either version 2, or (at your option)
00014 // any later version.
00015 //
00016 // The SC Toolkit is distributed in the hope that it will be useful,
00017 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00019 // GNU Library General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Library General Public License
00022 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
00023 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
00024 //
00025 // The U.S. Government is granted a limited license as per AL 91-7.
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 // grab references for speed
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 // node zero passes out indices
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 // now clean up
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 // all other nodes do the work
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 * if i=j=k or j=k=l, then this integral contributes
00194 * to J, K1, and K2 of G(ij), so
00195 * pkval = (ijkl) - 0.25 * ((ikjl)-(ilkj))
00196 * = 0.5 * (ijkl)
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 * if j=k, then this integral contributes
00208 * to J and K1 of G(ij)
00209 *
00210 * pkval = (ijkl) - 0.25 * (ikjl)
00211 * = 0.75 * (ijkl)
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 * this integral also contributes to K1 and K2 of
00221 * G(il)
00222 *
00223 * pkval = -0.25 * ((ijkl)+(ikjl))
00224 * = -0.5 * (ijkl)
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 * if i=k or j=l, then this integral contributes
00235 * to J and K2 of G(ij)
00236 *
00237 * pkval = (ijkl) - 0.25 * (ilkj)
00238 * = 0.75 * (ijkl)
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 * this integral also contributes to K1 and K2 of
00248 * G(ik)
00249 *
00250 * pkval = -0.25 * ((ijkl)+(ilkj))
00251 * = -0.5 * (ijkl)
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 * This integral contributes to J of G(ij)
00262 *
00263 * pkval = (ijkl)
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 * and to K1 of G(ik)
00273 *
00274 * pkval = -0.25 * (ijkl)
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 * if i!=j and k!=l, then this integral also
00285 * contributes to K2 of G(il)
00286 *
00287 * pkval = -0.25 * (ijkl)
00288 *
00289 * note: if we get here, then ik can't equal jl,
00290 * so pkval wasn't multiplied by 0.5 above.
00291 */
00292 ij = ij_offset(ii,ll);
00293 kl = ij_offset(kk,jj);
00294
00295 contribution.cont2(ij,kl,val);
00296 }
00297 }
00298 } else { // !e_any
00299 if (jj == kk) {
00300 /*
00301 * if j=k, then this integral contributes
00302 * to J and K1 of G(ij)
00303 *
00304 * pkval = (ijkl) - 0.25 * (ikjl)
00305 * = 0.75 * (ijkl)
00306 */
00307 contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
00308
00309 /*
00310 * this integral also contributes to K1 and K2 of
00311 * G(il)
00312 *
00313 * pkval = -0.25 * ((ijkl)+(ikjl))
00314 * = -0.5 * (ijkl)
00315 */
00316 contribution.cont3(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
00317
00318 } else if (ii == kk || jj == ll) {
00319 /*
00320 * if i=k or j=l, then this integral contributes
00321 * to J and K2 of G(ij)
00322 *
00323 * pkval = (ijkl) - 0.25 * (ilkj)
00324 * = 0.75 * (ijkl)
00325 */
00326 contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
00327
00328 /*
00329 * this integral also contributes to K1 and K2 of
00330 * G(ik)
00331 *
00332 * pkval = -0.25 * ((ijkl)+(ilkj))
00333 * = -0.5 * (ijkl)
00334 */
00335 contribution.cont3(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
00336
00337 } else {
00338 /*
00339 * This integral contributes to J of G(ij)
00340 *
00341 * pkval = (ijkl)
00342 */
00343 contribution.cont1(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
00344
00345 /*
00346 * and to K1 of G(ik)
00347 *
00348 * pkval = -0.25 * (ijkl)
00349 */
00350 contribution.cont2(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
00351
00352 /*
00353 * and to K2 of G(il)
00354 *
00355 * pkval = -0.25 * (ijkl)
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 // Local Variables:
00384 // mode: c++
00385 // c-file-style: "ETS"
00386 // End:
Generated at Fri Jan 10 08:14:09 2003 for MPQC 2.1.3 using the documentation package Doxygen 1.2.14.