MPQC
2.3.1
|
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 <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 // grab references for speed 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 //tim_enter_default(); 00147 tbi.compute_shell(i,j,k,l); 00148 //tim_exit_default(); 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 #ifdef HAVE_ISNAN 00177 if (isnan(pki_int)) 00178 abort(); 00179 #endif 00180 #endif 00181 00182 if (qijkl > 1) 00183 pki_int *= qijkl; 00184 00185 if (e_any) { 00186 int ij,kl; 00187 double val; 00188 00189 if (jj == kk) { 00190 /* 00191 * if i=j=k or j=k=l, then this integral contributes 00192 * to J, K1, and K2 of G(ij), so 00193 * pkval = (ijkl) - 0.25 * ((ikjl)-(ilkj)) 00194 * = 0.5 * (ijkl) 00195 */ 00196 if (ii == jj || kk == ll) { 00197 ij = i_offset(ii)+jj; 00198 kl = i_offset(kk)+ll; 00199 val = (ij==kl) ? 0.5*pki_int : pki_int; 00200 00201 GBuild<T>::contribution.cont5(ij,kl,val); 00202 00203 } else { 00204 /* 00205 * if j=k, then this integral contributes 00206 * to J and K1 of G(ij) 00207 * 00208 * pkval = (ijkl) - 0.25 * (ikjl) 00209 * = 0.75 * (ijkl) 00210 */ 00211 ij = i_offset(ii)+jj; 00212 kl = i_offset(kk)+ll; 00213 val = (ij==kl) ? 0.5*pki_int : pki_int; 00214 00215 GBuild<T>::contribution.cont4(ij,kl,val); 00216 00217 /* 00218 * this integral also contributes to K1 and K2 of 00219 * G(il) 00220 * 00221 * pkval = -0.25 * ((ijkl)+(ikjl)) 00222 * = -0.5 * (ijkl) 00223 */ 00224 ij = ij_offset(ii,ll); 00225 kl = ij_offset(kk,jj); 00226 val = (ij==kl) ? 0.5*pki_int : pki_int; 00227 00228 GBuild<T>::contribution.cont3(ij,kl,val); 00229 } 00230 } else if (ii == kk || jj == ll) { 00231 /* 00232 * if i=k or j=l, then this integral contributes 00233 * to J and K2 of G(ij) 00234 * 00235 * pkval = (ijkl) - 0.25 * (ilkj) 00236 * = 0.75 * (ijkl) 00237 */ 00238 ij = i_offset(ii)+jj; 00239 kl = i_offset(kk)+ll; 00240 val = (ij==kl) ? 0.5*pki_int : pki_int; 00241 00242 GBuild<T>::contribution.cont4(ij,kl,val); 00243 00244 /* 00245 * this integral also contributes to K1 and K2 of 00246 * G(ik) 00247 * 00248 * pkval = -0.25 * ((ijkl)+(ilkj)) 00249 * = -0.5 * (ijkl) 00250 */ 00251 ij = ij_offset(ii,kk); 00252 kl = ij_offset(jj,ll); 00253 val = (ij==kl) ? 0.5*pki_int : pki_int; 00254 00255 GBuild<T>::contribution.cont3(ij,kl,val); 00256 00257 } else { 00258 /* 00259 * This integral contributes to J of G(ij) 00260 * 00261 * pkval = (ijkl) 00262 */ 00263 ij = i_offset(ii)+jj; 00264 kl = i_offset(kk)+ll; 00265 val = (ij==kl) ? 0.5*pki_int : pki_int; 00266 00267 GBuild<T>::contribution.cont1(ij,kl,val); 00268 00269 /* 00270 * and to K1 of G(ik) 00271 * 00272 * pkval = -0.25 * (ijkl) 00273 */ 00274 ij = ij_offset(ii,kk); 00275 kl = ij_offset(jj,ll); 00276 val = (ij==kl) ? 0.5*pki_int : pki_int; 00277 00278 GBuild<T>::contribution.cont2(ij,kl,val); 00279 00280 if ((ii != jj) && (kk != ll)) { 00281 /* 00282 * if i!=j and k!=l, then this integral also 00283 * contributes to K2 of G(il) 00284 * 00285 * pkval = -0.25 * (ijkl) 00286 * 00287 * note: if we get here, then ik can't equal jl, 00288 * so pkval wasn't multiplied by 0.5 above. 00289 */ 00290 ij = ij_offset(ii,ll); 00291 kl = ij_offset(kk,jj); 00292 00293 GBuild<T>::contribution.cont2(ij,kl,val); 00294 } 00295 } 00296 } else { // !e_any 00297 if (jj == kk) { 00298 /* 00299 * if j=k, then this integral contributes 00300 * to J and K1 of G(ij) 00301 * 00302 * pkval = (ijkl) - 0.25 * (ikjl) 00303 * = 0.75 * (ijkl) 00304 */ 00305 GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int); 00306 00307 /* 00308 * this integral also contributes to K1 and K2 of 00309 * G(il) 00310 * 00311 * pkval = -0.25 * ((ijkl)+(ikjl)) 00312 * = -0.5 * (ijkl) 00313 */ 00314 GBuild<T>::contribution.cont3(ij_offset(ii,ll),ij_offset(kk,jj),pki_int); 00315 00316 } else if (ii == kk || jj == ll) { 00317 /* 00318 * if i=k or j=l, then this integral contributes 00319 * to J and K2 of G(ij) 00320 * 00321 * pkval = (ijkl) - 0.25 * (ilkj) 00322 * = 0.75 * (ijkl) 00323 */ 00324 GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int); 00325 00326 /* 00327 * this integral also contributes to K1 and K2 of 00328 * G(ik) 00329 * 00330 * pkval = -0.25 * ((ijkl)+(ilkj)) 00331 * = -0.5 * (ijkl) 00332 */ 00333 GBuild<T>::contribution.cont3(ij_offset(ii,kk),ij_offset(jj,ll),pki_int); 00334 00335 } else { 00336 /* 00337 * This integral contributes to J of G(ij) 00338 * 00339 * pkval = (ijkl) 00340 */ 00341 GBuild<T>::contribution.cont1(i_offset(ii)+jj,i_offset(kk)+ll,pki_int); 00342 00343 /* 00344 * and to K1 of G(ik) 00345 * 00346 * pkval = -0.25 * (ijkl) 00347 */ 00348 GBuild<T>::contribution.cont2(ij_offset(ii,kk),ij_offset(jj,ll),pki_int); 00349 00350 /* 00351 * and to K2 of G(il) 00352 * 00353 * pkval = -0.25 * (ijkl) 00354 */ 00355 GBuild<T>::contribution.cont2(ij_offset(ii,ll),ij_offset(kk,jj),pki_int); 00356 } 00357 } 00358 } 00359 } 00360 } 00361 } 00362 00363 tnint += (double) ni*nj*nk*nl; 00364 } 00365 } 00366 } 00367 } 00368 } 00369 }; 00370 00371 } 00372 00373 #endif 00374 00375 // Local Variables: 00376 // mode: c++ 00377 // c-file-style: "ETS" 00378 // End: