source: ThirdParty/mpqc_open/src/lib/chemistry/qc/scf/lgbuild.h@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 10.6 KB
Line 
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
42namespace sc {
43
44template<class T>
45class 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 }
72 ~LocalGBuild() {}
73
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:
Note: See TracBrowser for help on using the repository browser.