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