lbgbuild.h
//
// lbgbuild.h --- definitino of the load-balanced 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_lbgbuild_h
#define _chemistry_qc_scf_lbgbuild_h
#ifdef __GNUC__
#pragma interface
#endif
#include <chemistry/qc/scf/gbuild.h>
namespace sc {
template<class T>
class LocalLBGBuild : public GBuild<T> {
protected:
Ref<MessageGrp> grp_;
Ref<TwoBodyInt> tbi_;
Ref<Integral> integral_;
Ref<GaussianBasisSet> gbs_;
signed char *pmax;
public:
LocalLBGBuild(T& t, const Ref<TwoBodyInt>& tbi, const Ref<Integral>& ints,
const Ref<GaussianBasisSet>& bs, const Ref<MessageGrp>& g,
signed char *pm) :
GBuild<T>(t), grp_(g), tbi_(tbi), integral_(ints), gbs_(bs), pmax(pm) {}
~LocalLBGBuild() {}
void build_gmat(double accuracy) {
tim_enter("ao_gmat");
tim_set_default("quartet");
double tnint=0;
int tol = (int) (log(accuracy)/log(2.0));
int me=grp_->me();
int nproc = grp_->n();
Ref<PetiteList> rpl = integral_->petite_list();
// grab references for speed
GaussianBasisSet& gbs = *gbs_.pointer();
PetiteList& pl = *rpl.pointer();
TwoBodyInt& tbi = *tbi_.pointer();
tbi.set_redundant(0);
const double *intbuf = tbi.buffer();
int inds[4];
// node zero passes out indices
if (me==0) {
int i;
for (i=0; i < gbs.nshell(); i++) {
if (!pl.in_p1(i))
continue;
inds[0]=i;
for (int j=0; j <= i; j++) {
int oij = i_offset(i)+j;
if (!pl.in_p2(oij))
continue;
inds[1]=j;
tim_enter_default();
int from;
grp_->recvt(2323, &from, 1);
grp_->sendt(from, 3232, inds, 4);
tim_exit_default();
}
}
// now clean up
inds[0] = inds[1] = inds[2] = inds[3] = -1;
for (i=1; i < nproc; i++) {
int from;
grp_->recvt(2323, &from, 1);
grp_->sendt(from, 3232, inds, 4);
}
}
// all other nodes do the work
else {
do {
grp_->sendt(0, 2323, &me, 1);
grp_->recvt(3232, inds, 4);
int i=inds[0];
int j=inds[1];
if (i < 0)
break;
int fi=gbs.shell_to_function(i);
int ni=gbs(i).nfunction();
int oij = i_offset(i)+j;
int fj=gbs.shell_to_function(j);
int nj=gbs(j).nfunction();
int pmaxij = pmax[oij];
for (int k=0; k <= i; k++) {
int fk=gbs.shell_to_function(k);
int nk=gbs(k).nfunction();
int pmaxijk=pmaxij, ptmp;
if ((ptmp=pmax[i_offset(i)+k]-2) > pmaxijk) pmaxijk=ptmp;
if ((ptmp=pmax[ij_offset(j,k)]-2) > 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]-2) > pmaxijkl) pmaxijkl=ptmp;
if ((ptmp=pmax[ij_offset(j,l)]-2) > pmaxijkl) pmaxijkl=ptmp;
if (tbi.log2_shell_bound(i,j,k,l)+pmaxijkl < tol)
continue;
int qijkl = pl.in_p4(oij,okl,i,j,k,l);
if (!qijkl)
continue;
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;
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;
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;
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;
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;
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;
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;
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;
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);
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)
*/
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)
*/
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)
*/
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)
*/
contribution.cont3(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
} else {
/*
* This integral contributes to J of G(ij)
*
* pkval = (ijkl)
*/
contribution.cont1(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
/*
* and to K1 of G(ik)
*
* pkval = -0.25 * (ijkl)
*/
contribution.cont2(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
/*
* and to K2 of G(il)
*
* pkval = -0.25 * (ijkl)
*/
contribution.cont2(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
}
}
}
}
}
}
tnint += (double) ni*nj*nk*nl;
}
}
} while (inds[0] > -1);
}
grp_->sum(&tnint, 1, 0, 0);
ExEnv::out0() << indent << scprintf("%20.0f integrals\n", tnint);
tim_exit("ao_gmat");
}
};
}
#endif
// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End: