Main Page Class Hierarchy Compound List File List Compound 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 <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: