lgbuild.h
//
// lgbuild.h --- definition of the local G matrix builder
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Edward Seidl <seidl@janed.com>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB. If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//
#ifndef _chemistry_qc_scf_lgbuild_h
#define _chemistry_qc_scf_lgbuild_h
#ifdef __GNUC__
#pragma interface
#endif
#undef SCF_CHECK_INTS
#undef SCF_CHECK_BOUNDS
#undef SCF_DONT_USE_BOUNDS
#include <scconfig.h>
#include <chemistry/qc/scf/gbuild.h>
namespace sc {
template<class T>
class LocalGBuild : public GBuild<T> {
public:
double tnint;
protected:
MessageGrp *grp_;
TwoBodyInt *tbi_;
GaussianBasisSet *gbs_;
PetiteList *rpl_;
signed char * restrictxx pmax;
int threadno_;
int nthread_;
double accuracy_;
public:
LocalGBuild(T& t, const Ref<TwoBodyInt>& tbi, const Ref<PetiteList>& rpl,
const Ref<GaussianBasisSet>& bs, const Ref<MessageGrp>& g,
signed char *pm, double acc, int nt=1, int tn=0) :
GBuild<T>(t),
pmax(pm), nthread_(nt), threadno_(tn), accuracy_(acc)
{
grp_ = g.pointer();
tbi_ = tbi.pointer();
rpl_ = rpl.pointer();
gbs_ = bs.pointer();
}
~LocalGBuild() {}
void run() {
int tol = (int) (log(accuracy_)/log(2.0));
int me=grp_->me();
int nproc = grp_->n();
// grab references for speed
GaussianBasisSet& gbs = *gbs_;
PetiteList& pl = *rpl_;
TwoBodyInt& tbi = *tbi_;
tbi.set_redundant(0);
const double *intbuf = tbi.buffer();
tnint=0;
sc_int_least64_t threadind=0;
sc_int_least64_t ijklind=0;
for (int i=0; i < gbs.nshell(); i++) {
if (!pl.in_p1(i))
continue;
int fi=gbs.shell_to_function(i);
int ni=gbs(i).nfunction();
for (int j=0; j <= i; j++) {
int oij = i_offset(i)+j;
if (!pl.in_p2(oij))
continue;
int fj=gbs.shell_to_function(j);
int nj=gbs(j).nfunction();
int pmaxij = pmax[oij];
for (int k=0; k <= i; k++, ijklind++) {
if (ijklind%nproc != me)
continue;
threadind++;
if (threadind % nthread_ != threadno_)
continue;
int fk=gbs.shell_to_function(k);
int nk=gbs(k).nfunction();
int pmaxijk=pmaxij, ptmp;
if ((ptmp=pmax[i_offset(i)+k]-1) > pmaxijk) pmaxijk=ptmp;
if ((ptmp=pmax[ij_offset(j,k)]-1) > pmaxijk) pmaxijk=ptmp;
int okl = i_offset(k);
for (int l=0; l <= (k==i?j:k); l++,okl++) {
int pmaxijkl = pmaxijk;
if ((ptmp=pmax[okl]) > pmaxijkl) pmaxijkl=ptmp;
if ((ptmp=pmax[i_offset(i)+l]-1) > pmaxijkl) pmaxijkl=ptmp;
if ((ptmp=pmax[ij_offset(j,l)]-1) > pmaxijkl) pmaxijkl=ptmp;
int qijkl = pl.in_p4(oij,okl,i,j,k,l);
if (!qijkl)
continue;
#ifdef SCF_CHECK_BOUNDS
double intbound = pow(2.0,double(tbi.log2_shell_bound(i,j,k,l)));
double pbound = pow(2.0,double(pmaxijkl));
intbound *= qijkl;
GBuild<T>::contribution.set_bound(intbound, pbound);
#else
# ifndef SCF_DONT_USE_BOUNDS
if (tbi.log2_shell_bound(i,j,k,l)+pmaxijkl < tol)
continue;
# endif
#endif
//tim_enter_default();
tbi.compute_shell(i,j,k,l);
//tim_exit_default();
int e12 = (i==j);
int e34 = (k==l);
int e13e24 = (i==k) && (j==l);
int e_any = e12||e34||e13e24;
int fl=gbs.shell_to_function(l);
int nl=gbs(l).nfunction();
int ii,jj,kk,ll;
int I,J,K,L;
int index=0;
for (I=0, ii=fi; I < ni; I++, ii++) {
for (J=0, jj=fj; J <= (e12 ? I : nj-1); J++, jj++) {
for (K=0, kk=fk; K <= (e13e24 ? I : nk-1); K++, kk++) {
int lend = (e34 ? ((e13e24)&&(K==I) ? J : K)
: ((e13e24)&&(K==I)) ? J : nl-1);
for (L=0, ll=fl; L <= lend; L++, ll++, index++) {
double pki_int = intbuf[index];
if ((pki_int>0?pki_int:-pki_int) < 1.0e-15)
continue;
#ifdef SCF_CHECK_INTS
#ifdef HAVE_ISNAN
if (isnan(pki_int))
abort();
#endif
#endif
if (qijkl > 1)
pki_int *= qijkl;
if (e_any) {
int ij,kl;
double val;
if (jj == kk) {
/*
* if i=j=k or j=k=l, then this integral contributes
* to J, K1, and K2 of G(ij), so
* pkval = (ijkl) - 0.25 * ((ikjl)-(ilkj))
* = 0.5 * (ijkl)
*/
if (ii == jj || kk == ll) {
ij = i_offset(ii)+jj;
kl = i_offset(kk)+ll;
val = (ij==kl) ? 0.5*pki_int : pki_int;
GBuild<T>::contribution.cont5(ij,kl,val);
} else {
/*
* if j=k, then this integral contributes
* to J and K1 of G(ij)
*
* pkval = (ijkl) - 0.25 * (ikjl)
* = 0.75 * (ijkl)
*/
ij = i_offset(ii)+jj;
kl = i_offset(kk)+ll;
val = (ij==kl) ? 0.5*pki_int : pki_int;
GBuild<T>::contribution.cont4(ij,kl,val);
/*
* this integral also contributes to K1 and K2 of
* G(il)
*
* pkval = -0.25 * ((ijkl)+(ikjl))
* = -0.5 * (ijkl)
*/
ij = ij_offset(ii,ll);
kl = ij_offset(kk,jj);
val = (ij==kl) ? 0.5*pki_int : pki_int;
GBuild<T>::contribution.cont3(ij,kl,val);
}
} else if (ii == kk || jj == ll) {
/*
* if i=k or j=l, then this integral contributes
* to J and K2 of G(ij)
*
* pkval = (ijkl) - 0.25 * (ilkj)
* = 0.75 * (ijkl)
*/
ij = i_offset(ii)+jj;
kl = i_offset(kk)+ll;
val = (ij==kl) ? 0.5*pki_int : pki_int;
GBuild<T>::contribution.cont4(ij,kl,val);
/*
* this integral also contributes to K1 and K2 of
* G(ik)
*
* pkval = -0.25 * ((ijkl)+(ilkj))
* = -0.5 * (ijkl)
*/
ij = ij_offset(ii,kk);
kl = ij_offset(jj,ll);
val = (ij==kl) ? 0.5*pki_int : pki_int;
GBuild<T>::contribution.cont3(ij,kl,val);
} else {
/*
* This integral contributes to J of G(ij)
*
* pkval = (ijkl)
*/
ij = i_offset(ii)+jj;
kl = i_offset(kk)+ll;
val = (ij==kl) ? 0.5*pki_int : pki_int;
GBuild<T>::contribution.cont1(ij,kl,val);
/*
* and to K1 of G(ik)
*
* pkval = -0.25 * (ijkl)
*/
ij = ij_offset(ii,kk);
kl = ij_offset(jj,ll);
val = (ij==kl) ? 0.5*pki_int : pki_int;
GBuild<T>::contribution.cont2(ij,kl,val);
if ((ii != jj) && (kk != ll)) {
/*
* if i!=j and k!=l, then this integral also
* contributes to K2 of G(il)
*
* pkval = -0.25 * (ijkl)
*
* note: if we get here, then ik can't equal jl,
* so pkval wasn't multiplied by 0.5 above.
*/
ij = ij_offset(ii,ll);
kl = ij_offset(kk,jj);
GBuild<T>::contribution.cont2(ij,kl,val);
}
}
} else { // !e_any
if (jj == kk) {
/*
* if j=k, then this integral contributes
* to J and K1 of G(ij)
*
* pkval = (ijkl) - 0.25 * (ikjl)
* = 0.75 * (ijkl)
*/
GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
/*
* this integral also contributes to K1 and K2 of
* G(il)
*
* pkval = -0.25 * ((ijkl)+(ikjl))
* = -0.5 * (ijkl)
*/
GBuild<T>::contribution.cont3(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
} else if (ii == kk || jj == ll) {
/*
* if i=k or j=l, then this integral contributes
* to J and K2 of G(ij)
*
* pkval = (ijkl) - 0.25 * (ilkj)
* = 0.75 * (ijkl)
*/
GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
/*
* this integral also contributes to K1 and K2 of
* G(ik)
*
* pkval = -0.25 * ((ijkl)+(ilkj))
* = -0.5 * (ijkl)
*/
GBuild<T>::contribution.cont3(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
} else {
/*
* This integral contributes to J of G(ij)
*
* pkval = (ijkl)
*/
GBuild<T>::contribution.cont1(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
/*
* and to K1 of G(ik)
*
* pkval = -0.25 * (ijkl)
*/
GBuild<T>::contribution.cont2(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
/*
* and to K2 of G(il)
*
* pkval = -0.25 * (ijkl)
*/
GBuild<T>::contribution.cont2(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
}
}
}
}
}
}
tnint += (double) ni*nj*nk*nl;
}
}
}
}
}
};
}
#endif
// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End: