MPQC  3.0.0-alpha
lbgbuild.h
1 //
2 // lbgbuild.h --- definitino of the load-balanced local G matrix builder
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Edward Seidl <seidl@janed.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #ifndef _chemistry_qc_scf_lbgbuild_h
29 #define _chemistry_qc_scf_lbgbuild_h
30 
31 #include <chemistry/qc/scf/gbuild.h>
32 
33 namespace sc {
34 
35 template<class T>
36 class LocalLBGBuild : public GBuild<T> {
37  protected:
38  Ref<MessageGrp> grp_;
39  Ref<TwoBodyInt> tbi_;
40  Ref<Integral> integral_;
42  signed char *pmax;
43 
44  public:
45  LocalLBGBuild(T& t, const Ref<TwoBodyInt>& tbi, const Ref<Integral>& ints,
46  const Ref<GaussianBasisSet>& bs, const Ref<MessageGrp>& g,
47  signed char *pm) :
48  GBuild<T>(t), grp_(g), tbi_(tbi), integral_(ints), gbs_(bs), pmax(pm) {}
49  ~LocalLBGBuild() {}
50 
51  void build_gmat(double accuracy) {
52  Timer tim("ao_gmat");
53  tim.set_default("quartet");
54 
55  double tnint=0;
56  int tol = (int) (log(accuracy)/log(2.0));
57  int me=grp_->me();
58  int nproc = grp_->n();
59 
60  Ref<PetiteList> rpl = integral_->petite_list();
61 
62  // grab references for speed
63  GaussianBasisSet& gbs = *gbs_.pointer();
64  PetiteList& pl = *rpl.pointer();
65  TwoBodyInt& tbi = *tbi_.pointer();
66 
67  tbi.set_redundant(0);
68  const double *intbuf = tbi.buffer();
69 
70  int inds[4];
71 
72  // node zero passes out indices
73  if (me==0) {
74  int i;
75  for (i=0; i < gbs.nshell(); i++) {
76  if (!pl.in_p1(i))
77  continue;
78 
79  inds[0]=i;
80 
81  for (int j=0; j <= i; j++) {
82  int oij = i_offset(i)+j;
83 
84  if (!pl.in_p2(oij))
85  continue;
86 
87  inds[1]=j;
88 
89  tim.enter_default();
90  int from;
91  grp_->recvt(2323, &from, 1);
92  grp_->sendt(from, 3232, inds, 4);
93  tim.exit_default();
94  }
95  }
96 
97  // now clean up
98  inds[0] = inds[1] = inds[2] = inds[3] = -1;
99 
100  for (i=1; i < nproc; i++) {
101  int from;
102  grp_->recvt(2323, &from, 1);
103  grp_->sendt(from, 3232, inds, 4);
104  }
105  }
106 
107  // all other nodes do the work
108  else {
109  do {
110  grp_->sendt(0, 2323, &me, 1);
111  grp_->recvt(3232, inds, 4);
112 
113  int i=inds[0];
114  int j=inds[1];
115 
116  if (i < 0)
117  break;
118 
119  int fi=gbs.shell_to_function(i);
120  int ni=gbs(i).nfunction();
121 
122  int oij = i_offset(i)+j;
123 
124  int fj=gbs.shell_to_function(j);
125  int nj=gbs(j).nfunction();
126  int pmaxij = pmax[oij];
127 
128  for (int k=0; k <= i; k++) {
129  int fk=gbs.shell_to_function(k);
130  int nk=gbs(k).nfunction();
131 
132  int pmaxijk=pmaxij, ptmp;
133  if ((ptmp=pmax[i_offset(i)+k]-2) > pmaxijk) pmaxijk=ptmp;
134  if ((ptmp=pmax[ij_offset(j,k)]-2) > pmaxijk) pmaxijk=ptmp;
135 
136  int okl = i_offset(k);
137  for (int l=0; l <= (k==i?j:k); l++,okl++) {
138 
139  int pmaxijkl = pmaxijk;
140  if ((ptmp=pmax[okl]) > pmaxijkl) pmaxijkl=ptmp;
141  if ((ptmp=pmax[i_offset(i)+l]-2) > pmaxijkl) pmaxijkl=ptmp;
142  if ((ptmp=pmax[ij_offset(j,l)]-2) > pmaxijkl) pmaxijkl=ptmp;
143 
144 
145  if (tbi.log2_shell_bound(i,j,k,l)+pmaxijkl < tol)
146  continue;
147 
148  int qijkl = pl.in_p4(oij,okl,i,j,k,l);
149  if (!qijkl)
150  continue;
151 
152  tim.enter_default();
153  tbi.compute_shell(i,j,k,l);
154  tim.exit_default();
155 
156  int e12 = (i==j);
157  int e34 = (k==l);
158  int e13e24 = (i==k) && (j==l);
159  int e_any = e12||e34||e13e24;
160 
161  int fl=gbs.shell_to_function(l);
162  int nl=gbs(l).nfunction();
163 
164  int ii,jj,kk,ll;
165  int I,J,K,L;
166  int index=0;
167 
168  for (I=0, ii=fi; I < ni; I++, ii++) {
169  for (J=0, jj=fj; J <= (e12 ? I : nj-1); J++, jj++) {
170  for (K=0, kk=fk; K <= (e13e24 ? I : nk-1); K++, kk++) {
171 
172  int lend = (e34 ? ((e13e24)&&(K==I) ? J : K)
173  : ((e13e24)&&(K==I)) ? J : nl-1);
174  for (L=0, ll=fl; L <= lend; L++, ll++, index++) {
175  double pki_int = intbuf[index];
176 
177  if ((pki_int>0?pki_int:-pki_int) < 1.0e-15)
178  continue;
179 
180  if (qijkl > 1)
181  pki_int *= qijkl;
182 
183  if (e_any) {
184  int ij,kl;
185  double val;
186 
187  if (jj == kk) {
188  /*
189  * if i=j=k or j=k=l, then this integral contributes
190  * to J, K1, and K2 of G(ij), so
191  * pkval = (ijkl) - 0.25 * ((ikjl)-(ilkj))
192  * = 0.5 * (ijkl)
193  */
194  if (ii == jj || kk == ll) {
195  ij = i_offset(ii)+jj;
196  kl = i_offset(kk)+ll;
197  val = (ij==kl) ? 0.5*pki_int : pki_int;
198 
199  contribution.cont5(ij,kl,val);
200 
201  } else {
202  /*
203  * if j=k, then this integral contributes
204  * to J and K1 of G(ij)
205  *
206  * pkval = (ijkl) - 0.25 * (ikjl)
207  * = 0.75 * (ijkl)
208  */
209  ij = i_offset(ii)+jj;
210  kl = i_offset(kk)+ll;
211  val = (ij==kl) ? 0.5*pki_int : pki_int;
212 
213  contribution.cont4(ij,kl,val);
214 
215  /*
216  * this integral also contributes to K1 and K2 of
217  * G(il)
218  *
219  * pkval = -0.25 * ((ijkl)+(ikjl))
220  * = -0.5 * (ijkl)
221  */
222  ij = ij_offset(ii,ll);
223  kl = ij_offset(kk,jj);
224  val = (ij==kl) ? 0.5*pki_int : pki_int;
225 
226  contribution.cont3(ij,kl,val);
227  }
228  } else if (ii == kk || jj == ll) {
229  /*
230  * if i=k or j=l, then this integral contributes
231  * to J and K2 of G(ij)
232  *
233  * pkval = (ijkl) - 0.25 * (ilkj)
234  * = 0.75 * (ijkl)
235  */
236  ij = i_offset(ii)+jj;
237  kl = i_offset(kk)+ll;
238  val = (ij==kl) ? 0.5*pki_int : pki_int;
239 
240  contribution.cont4(ij,kl,val);
241 
242  /*
243  * this integral also contributes to K1 and K2 of
244  * G(ik)
245  *
246  * pkval = -0.25 * ((ijkl)+(ilkj))
247  * = -0.5 * (ijkl)
248  */
249  ij = ij_offset(ii,kk);
250  kl = ij_offset(jj,ll);
251  val = (ij==kl) ? 0.5*pki_int : pki_int;
252 
253  contribution.cont3(ij,kl,val);
254 
255  } else {
256  /*
257  * This integral contributes to J of G(ij)
258  *
259  * pkval = (ijkl)
260  */
261  ij = i_offset(ii)+jj;
262  kl = i_offset(kk)+ll;
263  val = (ij==kl) ? 0.5*pki_int : pki_int;
264 
265  contribution.cont1(ij,kl,val);
266 
267  /*
268  * and to K1 of G(ik)
269  *
270  * pkval = -0.25 * (ijkl)
271  */
272  ij = ij_offset(ii,kk);
273  kl = ij_offset(jj,ll);
274  val = (ij==kl) ? 0.5*pki_int : pki_int;
275 
276  contribution.cont2(ij,kl,val);
277 
278  if ((ii != jj) && (kk != ll)) {
279  /*
280  * if i!=j and k!=l, then this integral also
281  * contributes to K2 of G(il)
282  *
283  * pkval = -0.25 * (ijkl)
284  *
285  * note: if we get here, then ik can't equal jl,
286  * so pkval wasn't multiplied by 0.5 above.
287  */
288  ij = ij_offset(ii,ll);
289  kl = ij_offset(kk,jj);
290 
291  contribution.cont2(ij,kl,val);
292  }
293  }
294  } else { // !e_any
295  if (jj == kk) {
296  /*
297  * if j=k, then this integral contributes
298  * to J and K1 of G(ij)
299  *
300  * pkval = (ijkl) - 0.25 * (ikjl)
301  * = 0.75 * (ijkl)
302  */
303  contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
304 
305  /*
306  * this integral also contributes to K1 and K2 of
307  * G(il)
308  *
309  * pkval = -0.25 * ((ijkl)+(ikjl))
310  * = -0.5 * (ijkl)
311  */
312  contribution.cont3(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
313 
314  } else if (ii == kk || jj == ll) {
315  /*
316  * if i=k or j=l, then this integral contributes
317  * to J and K2 of G(ij)
318  *
319  * pkval = (ijkl) - 0.25 * (ilkj)
320  * = 0.75 * (ijkl)
321  */
322  contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
323 
324  /*
325  * this integral also contributes to K1 and K2 of
326  * G(ik)
327  *
328  * pkval = -0.25 * ((ijkl)+(ilkj))
329  * = -0.5 * (ijkl)
330  */
331  contribution.cont3(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
332 
333  } else {
334  /*
335  * This integral contributes to J of G(ij)
336  *
337  * pkval = (ijkl)
338  */
339  contribution.cont1(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
340 
341  /*
342  * and to K1 of G(ik)
343  *
344  * pkval = -0.25 * (ijkl)
345  */
346  contribution.cont2(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
347 
348  /*
349  * and to K2 of G(il)
350  *
351  * pkval = -0.25 * (ijkl)
352  */
353  contribution.cont2(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
354  }
355  }
356  }
357  }
358  }
359  }
360 
361  tnint += (double) ni*nj*nk*nl;
362  }
363  }
364  } while (inds[0] > -1);
365  }
366 
367  grp_->sum(&tnint, 1, 0, 0);
368  ExEnv::out0() << indent << scprintf("%20.0f integrals\n", tnint);
369 
370  tim.exit("ao_gmat");
371  }
372 
373 };
374 
375 }
376 
377 #endif
378 
379 // Local Variables:
380 // mode: c++
381 // c-file-style: "ETS"
382 // End:
sc::TwoBodyInt
This is an abstract base type for classes that compute integrals involving two electrons and 2 functi...
Definition: tbint.h:61
sc::GaussianBasisSet::nshell
unsigned int nshell() const
Return the number of shells.
Definition: gaussbas.h:500
sc::Timer
The Timer class uses RegionTimer to time intervals in an exception safe manner.
Definition: regtime.h:181
sc::Ref
A template class that maintains references counts.
Definition: ref.h:361
sc::Ref::pointer
T * pointer() const
Returns a pointer the reference counted object.
Definition: ref.h:413
sc::GBuild
Definition: gbuild.h:37
sc::TwoBodyInt::compute_shell
virtual void compute_shell(int, int, int, int)=0
Given four shell indices, integrals will be computed and placed in the buffer.
sc::PetiteList
PetiteList is a petite list (see Dupuis & King, IJQC 11,613,(1977) ) that can be used for constructin...
Definition: petite.h:120
sc::LocalLBGBuild
Definition: lbgbuild.h:36
sc::Timer::exit
void exit(const char *region=0)
Exit the current timing region.
sc::TwoBodyInt::set_redundant
virtual void set_redundant(int i)
See redundant().
Definition: tbint.h:168
sc::TwoBodyInt::buffer
virtual const double * buffer(TwoBodyOper::type type=TwoBodyOper::eri) const
The computed shell integrals will be put in the buffer returned by this member.
sc::TwoBodyInt::log2_shell_bound
virtual int log2_shell_bound(int=-1, int=-1, int=-1, int=-1)=0
Return log base 2 of the maximum magnitude of any integral in a shell block obtained from compute_she...
sc::GaussianBasisSet
The GaussianBasisSet class is used describe a basis set composed of atomic gaussian orbitals.
Definition: gaussbas.h:141
sc::ExEnv::out0
static std::ostream & out0()
Return an ostream that writes from node 0.
sc::GaussianBasisSet::shell_to_function
int shell_to_function(int i) const
Return the number of the first function in the given shell.
Definition: gaussbas.h:540
sc::contribution
Definition: petite.h:70
sc::Timer::set_default
void set_default(const char *region)
Default timing regions are provided as a performance optimization.
sc::scprintf
This class allows printf-like output to be sent to an ostream.
Definition: formio.h:97
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14

Generated at Sun Jan 26 2020 23:24:00 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.