- Main Page
- Related Pages
- Classes
- Files

# 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 #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: