MPQC  3.0.0-alpha
lgbuild.h
1 //
2 // lgbuild.h --- definition of the 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_lgbuild_h
29 #define _chemistry_qc_scf_lgbuild_h
30 
31 #undef SCF_CHECK_INTS
32 #undef SCF_CHECK_BOUNDS
33 #undef SCF_DONT_USE_BOUNDS
34 
35 #include <mpqc_config.h>
36 #include <chemistry/qc/scf/gbuild.h>
37 
38 namespace sc {
39 
40 template<class T>
41 class LocalGBuild : public GBuild<T> {
42  public:
43  double tnint;
44  double tinttime;
45 
46  protected:
47  MessageGrp *grp_;
48  TwoBodyInt *tbi_;
49  GaussianBasisSet *gbs_;
50  PetiteList *rpl_;
51 
52  signed char * RESTRICT pmax;
53  int threadno_;
54  int nthread_;
55  double accuracy_;
56 
57  public:
58  LocalGBuild(T& t, const Ref<TwoBodyInt>& tbi, const Ref<PetiteList>& rpl,
59  const Ref<GaussianBasisSet>& bs, const Ref<MessageGrp>& g,
60  signed char *pm, double acc, int nt=1, int tn=0) :
61  GBuild<T>(t),
62  pmax(pm), threadno_(tn), nthread_(nt), accuracy_(acc)
63  {
64  grp_ = g.pointer();
65  tbi_ = tbi.pointer();
66  rpl_ = rpl.pointer();
67  gbs_ = bs.pointer();
68  }
69  ~LocalGBuild() {}
70 
71  void run() {
72  int tol = (int) (log(accuracy_)/log(2.0));
73  int me=grp_->me();
74  int nproc = grp_->n();
75 
76  // grab references for speed
77  GaussianBasisSet& gbs = *gbs_;
78  PetiteList& pl = *rpl_;
79  TwoBodyInt& tbi = *tbi_;
80 
81  tbi.set_redundant(0);
82  const double *intbuf = tbi.buffer();
83 
84  tnint=0;
85  tinttime=0.0;
86  sc_int_least64_t threadind=0;
87  sc_int_least64_t ijklind=0;
88 
89  for (int i=0; i < gbs.nshell(); i++) {
90  if (!pl.in_p1(i))
91  continue;
92 
93  int fi=gbs.shell_to_function(i);
94  int ni=gbs(i).nfunction();
95 
96  for (int j=0; j <= i; j++) {
97  int oij = i_offset(i)+j;
98 
99  if (!pl.in_p2(oij))
100  continue;
101 
102  int fj=gbs.shell_to_function(j);
103  int nj=gbs(j).nfunction();
104  int pmaxij = pmax[oij];
105 
106  for (int k=0; k <= i; k++, ijklind++) {
107  if (ijklind%nproc != me)
108  continue;
109 
110  threadind++;
111  if (threadind % nthread_ != threadno_)
112  continue;
113 
114  int fk=gbs.shell_to_function(k);
115  int nk=gbs(k).nfunction();
116 
117  int pmaxijk=pmaxij, ptmp;
118  if ((ptmp=pmax[i_offset(i)+k]-1) > pmaxijk) pmaxijk=ptmp;
119  if ((ptmp=pmax[ij_offset(j,k)]-1) > pmaxijk) pmaxijk=ptmp;
120 
121  int okl = i_offset(k);
122  for (int l=0; l <= (k==i?j:k); l++,okl++) {
123  int pmaxijkl = pmaxijk;
124  if ((ptmp=pmax[okl]) > pmaxijkl) pmaxijkl=ptmp;
125  if ((ptmp=pmax[i_offset(i)+l]-1) > pmaxijkl) pmaxijkl=ptmp;
126  if ((ptmp=pmax[ij_offset(j,l)]-1) > pmaxijkl) pmaxijkl=ptmp;
127 
128  int qijkl = pl.in_p4(oij,okl,i,j,k,l);
129  if (!qijkl)
130  continue;
131 
132 #ifdef SCF_CHECK_BOUNDS
133  double intbound = pow(2.0,double(tbi.log2_shell_bound(i,j,k,l)));
134  double pbound = pow(2.0,double(pmaxijkl));
135  intbound *= qijkl;
136  GBuild<T>::contribution.set_bound(intbound, pbound);
137 #else
138 # ifndef SCF_DONT_USE_BOUNDS
139  if (tbi.log2_shell_bound(i,j,k,l)+pmaxijkl < tol)
140  continue;
141 # endif
142 #endif
143 
144 #ifdef MPQC_NEW_FEATURES
145  const auto tstart = sc::RegionTimer::get_time_point();
146 #endif
147  tbi.compute_shell(i,j,k,l);
148 #ifdef MPQC_NEW_FEATURES
149  const auto tstop = RegionTimer::get_time_point();
150  const std::chrono::duration<double> time_elapsed = tstop - tstart;
151  tinttime += time_elapsed.count();
152 #endif
153 
154  int e12 = (i==j);
155  int e34 = (k==l);
156  int e13e24 = (i==k) && (j==l);
157  int e_any = e12||e34||e13e24;
158 
159  int fl=gbs.shell_to_function(l);
160  int nl=gbs(l).nfunction();
161 
162  int ii,jj,kk,ll;
163  int I,J,K,L;
164  int index=0;
165 
166  for (I=0, ii=fi; I < ni; I++, ii++) {
167  for (J=0, jj=fj; J <= (e12 ? I : nj-1); J++, jj++) {
168  for (K=0, kk=fk; K <= (e13e24 ? I : nk-1); K++, kk++) {
169  int lend = (e34 ? ((e13e24)&&(K==I) ? J : K)
170  : ((e13e24)&&(K==I)) ? J : nl-1);
171 
172  for (L=0, ll=fl; L <= lend; L++, ll++, index++) {
173 
174  double pki_int = intbuf[index];
175 
176  if ((pki_int>0?pki_int:-pki_int) < 1.0e-15)
177  continue;
178 
179 #ifdef SCF_CHECK_INTS
180 #ifdef HAVE_ISNAN
181  if (isnan(pki_int))
182  abort();
183 #endif
184 #endif
185 
186  if (qijkl > 1)
187  pki_int *= qijkl;
188 
189  if (e_any) {
190  int ij,kl;
191  double val;
192 
193  if (jj == kk) {
194  /*
195  * if i=j=k or j=k=l, then this integral contributes
196  * to J, K1, and K2 of G(ij), so
197  * pkval = (ijkl) - 0.25 * ((ikjl)-(ilkj))
198  * = 0.5 * (ijkl)
199  */
200  if (ii == jj || kk == ll) {
201  ij = i_offset(ii)+jj;
202  kl = i_offset(kk)+ll;
203  val = (ij==kl) ? 0.5*pki_int : pki_int;
204 
205  GBuild<T>::contribution.cont5(ij,kl,val);
206 
207  } else {
208  /*
209  * if j=k, then this integral contributes
210  * to J and K1 of G(ij)
211  *
212  * pkval = (ijkl) - 0.25 * (ikjl)
213  * = 0.75 * (ijkl)
214  */
215  ij = i_offset(ii)+jj;
216  kl = i_offset(kk)+ll;
217  val = (ij==kl) ? 0.5*pki_int : pki_int;
218 
219  GBuild<T>::contribution.cont4(ij,kl,val);
220 
221  /*
222  * this integral also contributes to K1 and K2 of
223  * G(il)
224  *
225  * pkval = -0.25 * ((ijkl)+(ikjl))
226  * = -0.5 * (ijkl)
227  */
228  ij = ij_offset(ii,ll);
229  kl = ij_offset(kk,jj);
230  val = (ij==kl) ? 0.5*pki_int : pki_int;
231 
232  GBuild<T>::contribution.cont3(ij,kl,val);
233  }
234  } else if (ii == kk || jj == ll) {
235  /*
236  * if i=k or j=l, then this integral contributes
237  * to J and K2 of G(ij)
238  *
239  * pkval = (ijkl) - 0.25 * (ilkj)
240  * = 0.75 * (ijkl)
241  */
242  ij = i_offset(ii)+jj;
243  kl = i_offset(kk)+ll;
244  val = (ij==kl) ? 0.5*pki_int : pki_int;
245 
246  GBuild<T>::contribution.cont4(ij,kl,val);
247 
248  /*
249  * this integral also contributes to K1 and K2 of
250  * G(ik)
251  *
252  * pkval = -0.25 * ((ijkl)+(ilkj))
253  * = -0.5 * (ijkl)
254  */
255  ij = ij_offset(ii,kk);
256  kl = ij_offset(jj,ll);
257  val = (ij==kl) ? 0.5*pki_int : pki_int;
258 
259  GBuild<T>::contribution.cont3(ij,kl,val);
260 
261  } else {
262  /*
263  * This integral contributes to J of G(ij)
264  *
265  * pkval = (ijkl)
266  */
267  ij = i_offset(ii)+jj;
268  kl = i_offset(kk)+ll;
269  val = (ij==kl) ? 0.5*pki_int : pki_int;
270 
271  GBuild<T>::contribution.cont1(ij,kl,val);
272 
273  /*
274  * and to K1 of G(ik)
275  *
276  * pkval = -0.25 * (ijkl)
277  */
278  ij = ij_offset(ii,kk);
279  kl = ij_offset(jj,ll);
280  val = (ij==kl) ? 0.5*pki_int : pki_int;
281 
282  GBuild<T>::contribution.cont2(ij,kl,val);
283 
284  if ((ii != jj) && (kk != ll)) {
285  /*
286  * if i!=j and k!=l, then this integral also
287  * contributes to K2 of G(il)
288  *
289  * pkval = -0.25 * (ijkl)
290  *
291  * note: if we get here, then ik can't equal jl,
292  * so pkval wasn't multiplied by 0.5 above.
293  */
294  ij = ij_offset(ii,ll);
295  kl = ij_offset(kk,jj);
296 
297  GBuild<T>::contribution.cont2(ij,kl,val);
298  }
299  }
300  } else { // !e_any
301  if (jj == kk) {
302  /*
303  * if j=k, then this integral contributes
304  * to J and K1 of G(ij)
305  *
306  * pkval = (ijkl) - 0.25 * (ikjl)
307  * = 0.75 * (ijkl)
308  */
309  GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
310 
311  /*
312  * this integral also contributes to K1 and K2 of
313  * G(il)
314  *
315  * pkval = -0.25 * ((ijkl)+(ikjl))
316  * = -0.5 * (ijkl)
317  */
318  GBuild<T>::contribution.cont3(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
319 
320  } else if (ii == kk || jj == ll) {
321  /*
322  * if i=k or j=l, then this integral contributes
323  * to J and K2 of G(ij)
324  *
325  * pkval = (ijkl) - 0.25 * (ilkj)
326  * = 0.75 * (ijkl)
327  */
328  GBuild<T>::contribution.cont4(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
329 
330  /*
331  * this integral also contributes to K1 and K2 of
332  * G(ik)
333  *
334  * pkval = -0.25 * ((ijkl)+(ilkj))
335  * = -0.5 * (ijkl)
336  */
337  GBuild<T>::contribution.cont3(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
338 
339  } else {
340  /*
341  * This integral contributes to J of G(ij)
342  *
343  * pkval = (ijkl)
344  */
345  GBuild<T>::contribution.cont1(i_offset(ii)+jj,i_offset(kk)+ll,pki_int);
346 
347  /*
348  * and to K1 of G(ik)
349  *
350  * pkval = -0.25 * (ijkl)
351  */
352  GBuild<T>::contribution.cont2(ij_offset(ii,kk),ij_offset(jj,ll),pki_int);
353 
354  /*
355  * and to K2 of G(il)
356  *
357  * pkval = -0.25 * (ijkl)
358  */
359  GBuild<T>::contribution.cont2(ij_offset(ii,ll),ij_offset(kk,jj),pki_int);
360  }
361  }
362  }
363  }
364  }
365  }
366 
367  tnint += (double) ni*nj*nk*nl;
368  }
369  }
370  }
371  }
372  }
373 };
374 
375 }
376 
377 #endif
378 
379 // Local Variables:
380 // mode: c++
381 // c-file-style: "ETS"
382 // End:
sc::LocalGBuild::run
void run()
This is called with the Thread is run from a ThreadGrp.
Definition: lgbuild.h:71
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::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::MessageGrp::me
int me()
Returns my processor number. In the range [0,n()).
Definition: message.h:194
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::LocalGBuild::tinttime
double tinttime
total user time spent computing integrals
Definition: lgbuild.h:44
sc::LocalGBuild::tnint
double tnint
total # of integrals computed in this thread
Definition: lgbuild.h:43
sc::MessageGrp::n
int n()
Returns the number of processors.
Definition: message.h:192
sc::LocalGBuild
Definition: lgbuild.h:41
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::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::MessageGrp
The MessageGrp abstract class provides a mechanism for moving data and objects between nodes in a par...
Definition: message.h:120
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.