Main Page | Class Hierarchy | Class List | File List | Class Members | Related Pages

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 <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 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 * if i=j=k or j=k=l, then this integral contributes 00190 * to J, K1, and K2 of G(ij), so 00191 * pkval = (ijkl) - 0.25 * ((ikjl)-(ilkj)) 00192 * = 0.5 * (ijkl) 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 * if j=k, then this integral contributes 00204 * to J and K1 of G(ij) 00205 * 00206 * pkval = (ijkl) - 0.25 * (ikjl) 00207 * = 0.75 * (ijkl) 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 * this integral also contributes to K1 and K2 of 00217 * G(il) 00218 * 00219 * pkval = -0.25 * ((ijkl)+(ikjl)) 00220 * = -0.5 * (ijkl) 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 * if i=k or j=l, then this integral contributes 00231 * to J and K2 of G(ij) 00232 * 00233 * pkval = (ijkl) - 0.25 * (ilkj) 00234 * = 0.75 * (ijkl) 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 * this integral also contributes to K1 and K2 of 00244 * G(ik) 00245 * 00246 * pkval = -0.25 * ((ijkl)+(ilkj)) 00247 * = -0.5 * (ijkl) 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 * This integral contributes to J of G(ij) 00258 * 00259 * pkval = (ijkl) 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 * and to K1 of G(ik) 00269 * 00270 * pkval = -0.25 * (ijkl) 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 * if i!=j and k!=l, then this integral also 00281 * contributes to K2 of G(il) 00282 * 00283 * pkval = -0.25 * (ijkl) 00284 * 00285 * note: if we get here, then ik can't equal jl, 00286 * so pkval wasn't multiplied by 0.5 above. 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 { // !e_any 00295 if (jj == kk) { 00296 /* 00297 * if j=k, then this integral contributes 00298 * to J and K1 of G(ij) 00299 * 00300 * pkval = (ijkl) - 0.25 * (ikjl) 00301 * = 0.75 * (ijkl) 00302 */ 00303 GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int); 00304 00305 /* 00306 * this integral also contributes to K1 and K2 of 00307 * G(il) 00308 * 00309 * pkval = -0.25 * ((ijkl)+(ikjl)) 00310 * = -0.5 * (ijkl) 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 * if i=k or j=l, then this integral contributes 00317 * to J and K2 of G(ij) 00318 * 00319 * pkval = (ijkl) - 0.25 * (ilkj) 00320 * = 0.75 * (ijkl) 00321 */ 00322 GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int); 00323 00324 /* 00325 * this integral also contributes to K1 and K2 of 00326 * G(ik) 00327 * 00328 * pkval = -0.25 * ((ijkl)+(ilkj)) 00329 * = -0.5 * (ijkl) 00330 */ 00331 GBuild<T>::contribution.cont3(ij_offset(ii,kk),ij_offset(jj,ll),pki_int); 00332 00333 } else { 00334 /* 00335 * This integral contributes to J of G(ij) 00336 * 00337 * pkval = (ijkl) 00338 */ 00339 GBuild<T>::contribution.cont1(i_offset(ii)+jj,i_offset(kk)+ll,pki_int); 00340 00341 /* 00342 * and to K1 of G(ik) 00343 * 00344 * pkval = -0.25 * (ijkl) 00345 */ 00346 GBuild<T>::contribution.cont2(ij_offset(ii,kk),ij_offset(jj,ll),pki_int); 00347 00348 /* 00349 * and to K2 of G(il) 00350 * 00351 * pkval = -0.25 * (ijkl) 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 // Local Variables: 00374 // mode: c++ 00375 // c-file-style: "ETS" 00376 // End:

Generated at Sat Dec 18 15:14:21 2004 for MPQC 2.2.3 using the documentation package Doxygen 1.3.7-20040617.