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