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