source: ThirdParty/mpqc_open/src/lib/chemistry/qc/mbptr12/compute_a_gebc_abs1.cc

Candidate_v1.6.1
Last change on this file was 860145, checked in by Frederik Heber <heber@…>, 8 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'

  • Property mode set to 100644
File size: 14.4 KB
Line 
1//
2// compute_a_gebc_abs1.cc
3//
4// Copyright (C) 2004 Edward Valeev
5//
6// Author: Edward Valeev <edward.valeev@chemistry.gatech.edu>
7// Maintainer: EV
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#include <stdexcept>
29#include <stdlib.h>
30#include <math.h>
31#include <limits.h>
32
33#include <scconfig.h>
34#include <util/misc/formio.h>
35#include <util/misc/timer.h>
36#include <util/class/class.h>
37#include <util/state/state.h>
38#include <util/state/state_text.h>
39#include <util/state/state_bin.h>
40#include <math/scmat/matrix.h>
41#include <chemistry/molecule/molecule.h>
42#include <chemistry/qc/basis/integral.h>
43#include <chemistry/qc/mbpt/bzerofast.h>
44#include <chemistry/qc/mbptr12/r12ia.h>
45#include <chemistry/qc/mbptr12/vxb_eval_info.h>
46#include <chemistry/qc/mbptr12/pairiter.h>
47#include <chemistry/qc/mbptr12/r12int_eval.h>
48
49using namespace std;
50using namespace sc;
51
52#define PRINT_R12_INTERMED 0
53
54void
55R12IntEval::abs1_contrib_to_VXB_gebc_()
56{
57 if (evaluated_)
58 return;
59 LinearR12::ABSMethod abs_method = r12info_->abs_method();
60 Ref<MessageGrp> msg = r12info_->msg();
61 Ref<MemoryGrp> mem = r12info_->mem();
62 Ref<ThreadGrp> thr = r12info_->thr();
63 const int num_te_types = 4;
64 enum te_types {eri=0, r12=1, r12t1=2, r12t2=3};
65
66 tim_enter("mp2-r12a intermeds");
67
68 int me = msg->me();
69 int nproc = msg->n();
70
71 ExEnv::out0() << endl << indent
72 << "Entered ABS A (GEBC) intermediates evaluator" << endl;
73 ExEnv::out0() << incindent;
74
75 // Do the AO->MO transform
76 Ref<MOIntsTransformFactory> tfactory = r12info_->tfactory();
77 tfactory->set_spaces(r12info_->act_occ_space(),r12info_->occ_space(),
78 r12info_->act_occ_space(),r12info_->ribs_space());
79 Ref<TwoBodyMOIntsTransform> ikjy_tform = tfactory->twobody_transform_13("(ik|jy)");
80 ikjy_tform->set_num_te_types(num_te_types);
81 ikjy_tform->compute();
82 Ref<R12IntsAcc> ijky_acc = ikjy_tform->ints_acc();
83 if (num_te_types != ijky_acc->num_te_types())
84 throw std::runtime_error("R12IntEval::obs_contrib_to_VXB_gebc() -- number of MO integral types is wrong");
85
86 const int nocc = r12info_->nocc();
87 const int noso_ri = r12info_->ribs_space()->rank();
88
89 /*--------------------------------
90 Compute MP2-R12/A intermediates
91 and collect on node0
92 --------------------------------*/
93 ExEnv::out0() << indent << "Begin computation of intermediates" << endl;
94 tim_enter("intermediates");
95 SpatialMOPairIter_eq ij_iter(r12info_->act_occ_space());
96 SpatialMOPairIter_eq kl_iter(r12info_->act_occ_space());
97 int naa = ij_iter.nij_aa(); // Number of alpha-alpha pairs (i > j)
98 int nab = ij_iter.nij_ab(); // Number of alpha-beta pairs
99 if (debug_) {
100 ExEnv::out0() << indent << "naa = " << naa << endl;
101 ExEnv::out0() << indent << "nab = " << nab << endl;
102 }
103
104 // Compute intermediates
105 if (debug_)
106 ExEnv::out0() << indent << "Ready to compute MP2-R12/A (GEBC) intermediates" << endl;
107
108 // Compute the number of tasks that have full access to the integrals
109 // and split the work among them
110 vector<int> proc_with_ints;
111 int nproc_with_ints = tasks_with_ints_(ijky_acc,proc_with_ints);
112
113
114 //////////////////////////////////////////////////////////////
115 //
116 // Evaluation of the intermediates proceeds as follows:
117 //
118 // loop over batches of kl, k >= l, 0<=k,l<nocc_act
119 // load (kl|xy), (kl| [T1,r12] |xy), and (lk| [T1,r12] |xy)
120 // (aka kl-sets) into memory
121 //
122 // loop over batches of ij, i>=j, 0<=i,j<nocc_act
123 // load (ij|r12|xy) into memory
124 // (aka ij-sets) into memory
125 // compute V[ij][kl] and T[ij][kl] for all ij and kl in
126 // the "direct product" batch
127 // end ij loop
128 // end kl loop
129 //
130 /////////////////////////////////////////////////////////////////////////////////
131
132 if (ijky_acc->has_access(me)) {
133
134 for(kl_iter.start();int(kl_iter);kl_iter.next()) {
135
136 const int kl = kl_iter.ij();
137 // Figure out if this task will handle this kl
138 int kl_proc = kl%nproc_with_ints;
139 if (kl_proc != proc_with_ints[me])
140 continue;
141 const int k = kl_iter.i();
142 const int l = kl_iter.j();
143 const int kl_aa = kl_iter.ij_aa();
144 const int kl_ab = kl_iter.ij_ab();
145 const int lk_ab = kl_iter.ij_ba();
146
147 if (debug_)
148 ExEnv::outn() << indent << "task " << me << ": working on (k,l) = " << k << "," << l << " " << endl;
149
150 // Get (|1/r12|), (|r12|), (|[r12,T1]|), and (|[r12,T2]|) integrals
151 tim_enter("MO ints retrieve");
152 double *klox_buf_eri = ijky_acc->retrieve_pair_block(k,l,R12IntsAcc::eri);
153 double *klox_buf_r12 = ijky_acc->retrieve_pair_block(k,l,R12IntsAcc::r12);
154 double *klox_buf_r12t1 = ijky_acc->retrieve_pair_block(k,l,R12IntsAcc::r12t1);
155 double *klox_buf_r12t2 = ijky_acc->retrieve_pair_block(k,l,R12IntsAcc::r12t2);
156
157 double *lkox_buf_eri = ijky_acc->retrieve_pair_block(l,k,R12IntsAcc::eri);
158 double *lkox_buf_r12 = ijky_acc->retrieve_pair_block(l,k,R12IntsAcc::r12);
159 double *lkox_buf_r12t1 = ijky_acc->retrieve_pair_block(l,k,R12IntsAcc::r12t1);
160 double *lkox_buf_r12t2 = ijky_acc->retrieve_pair_block(l,k,R12IntsAcc::r12t2);
161 tim_exit("MO ints retrieve");
162
163 if (debug_)
164 ExEnv::outn() << indent << "task " << me << ": obtained kl blocks" << endl;
165
166 // to avoid every task hitting same ij at the same time, stagger ij-accesses, i.e. each kl task will start with ij=kl+1
167 for(ij_iter.start(kl+1);int(ij_iter);ij_iter.next()) {
168
169 const int i = ij_iter.i();
170 const int j = ij_iter.j();
171 const int ij_aa = ij_iter.ij_aa();
172 const int ij_ab = ij_iter.ij_ab();
173 const int ji_ab = ij_iter.ij_ba();
174
175 if (debug_)
176 ExEnv::outn() << indent << "task " << me << ": (k,l) = " << k << "," << l << ": (i,j) = " << i << "," << j << endl;
177
178 tim_enter("MO ints retrieve");
179 double *ijox_buf_r12 = ijky_acc->retrieve_pair_block(i,j,R12IntsAcc::r12);
180 double *jiox_buf_r12 = ijky_acc->retrieve_pair_block(j,i,R12IntsAcc::r12);
181 tim_exit("MO ints retrieve");
182
183 if (debug_)
184 ExEnv::outn() << indent << "task " << me << ": obtained ij blocks" << endl;
185
186
187 tim_enter("MO ints contraction");
188 double Vaa_ijkl, Vab_ijkl, Vab_jikl, Vab_ijlk, Vab_jilk;
189 double Xaa_ijkl, Xab_ijkl, Xab_jikl, Xab_ijlk, Xab_jilk;
190 double Taa_ijkl, Tab_ijkl, Tab_jikl, Tab_ijlk, Tab_jilk;
191 Vaa_ijkl = Vab_ijkl = Vab_jikl = Vab_ijlk = Vab_jilk = 0.0;
192 Xaa_ijkl = Xab_ijkl = Xab_jikl = Xab_ijlk = Xab_jilk = 0.0;
193 Taa_ijkl = Tab_ijkl = Tab_jikl = Tab_ijlk = Tab_jilk = 0.0;
194 for(int o=0;o<nocc;o++) {
195 const double pfac_xy = 1.0;
196 for(int x=0;x<noso_ri;x++) {
197 int ox_offset = o*noso_ri + x;
198 double ij_r12_ox = ijox_buf_r12[ox_offset];
199 double ji_r12_ox = jiox_buf_r12[ox_offset];
200 double kl_eri_ox = klox_buf_eri[ox_offset];
201 double lk_eri_ox = lkox_buf_eri[ox_offset];
202 Vab_ijkl -= pfac_xy * (ij_r12_ox * kl_eri_ox + ji_r12_ox * lk_eri_ox);
203 if (ij_ab != ji_ab)
204 Vab_jikl -= pfac_xy * (ji_r12_ox * kl_eri_ox + ij_r12_ox * lk_eri_ox);
205 if (kl_ab != lk_ab)
206 Vab_ijlk -= pfac_xy * (ij_r12_ox * lk_eri_ox + ji_r12_ox * kl_eri_ox);
207 if (ij_ab != ji_ab && kl_ab != lk_ab) {
208 Vab_jilk -= pfac_xy * (ij_r12_ox * kl_eri_ox + ji_r12_ox * lk_eri_ox);
209 }
210 if (ij_aa != -1 && kl_aa != -1) {
211 Vaa_ijkl -= pfac_xy * (ij_r12_ox - ji_r12_ox)*(kl_eri_ox - lk_eri_ox);
212 }
213 double kl_r12_ox = klox_buf_r12[ox_offset];
214 double lk_r12_ox = lkox_buf_r12[ox_offset];
215 Xab_ijkl -= pfac_xy * (ij_r12_ox * kl_r12_ox + ji_r12_ox * lk_r12_ox);
216 if (ij_ab != ji_ab)
217 Xab_jikl -= pfac_xy * (ji_r12_ox * kl_r12_ox + ij_r12_ox * lk_r12_ox);
218 if (kl_ab != lk_ab)
219 Xab_ijlk -= pfac_xy * (ij_r12_ox * lk_r12_ox + ji_r12_ox * kl_r12_ox);
220 if (ij_ab != ji_ab && kl_ab != lk_ab) {
221 Xab_jilk -= pfac_xy * (ij_r12_ox * kl_r12_ox + ji_r12_ox * lk_r12_ox);
222 }
223 if (ij_aa != -1 && kl_aa != -1) {
224 Xaa_ijkl -= pfac_xy * (ij_r12_ox - ji_r12_ox)*(kl_r12_ox - lk_r12_ox);
225 }
226 double kl_r12t1_ox = klox_buf_r12t1[ox_offset];
227 double kl_r12t2_ox = klox_buf_r12t2[ox_offset];
228 double lk_r12t1_ox = lkox_buf_r12t1[ox_offset];
229 double lk_r12t2_ox = lkox_buf_r12t2[ox_offset];
230 double kl_Tr12_ox = -kl_r12t1_ox-kl_r12t2_ox;
231 double lk_Tr12_ox = -lk_r12t1_ox-lk_r12t2_ox;
232 Tab_ijkl += pfac_xy * (ij_r12_ox * kl_Tr12_ox + ji_r12_ox * lk_Tr12_ox);
233 if (ij_ab != ji_ab)
234 Tab_jikl += pfac_xy * (ji_r12_ox * kl_Tr12_ox + ij_r12_ox * lk_Tr12_ox);
235 if (kl_ab != lk_ab)
236 Tab_ijlk += pfac_xy * (ij_r12_ox * lk_Tr12_ox + ji_r12_ox * kl_Tr12_ox);
237 if (ij_ab != ji_ab && kl_ab != lk_ab) {
238 Tab_jilk += pfac_xy * (ij_r12_ox * kl_Tr12_ox + ji_r12_ox * lk_Tr12_ox);
239 }
240 if (ij_aa != -1 && kl_aa != -1) {
241 Taa_ijkl += pfac_xy * (ij_r12_ox - ji_r12_ox)*(kl_Tr12_ox - lk_Tr12_ox);
242 }
243 }
244 }
245 Vab_.accumulate_element(ij_ab,kl_ab,Vab_ijkl);
246 if (ij_ab != ji_ab)
247 Vab_.accumulate_element(ji_ab,kl_ab,Vab_jikl);
248 if (kl_ab != lk_ab)
249 Vab_.accumulate_element(ij_ab,lk_ab,Vab_ijlk);
250 if (ij_ab != ji_ab && kl_ab != lk_ab)
251 Vab_.accumulate_element(ji_ab,lk_ab,Vab_jilk);
252 if (ij_aa != -1 && kl_aa != -1)
253 Vaa_.accumulate_element(ij_aa,kl_aa,Vaa_ijkl);
254 Xab_.accumulate_element(ij_ab,kl_ab,Xab_ijkl);
255 if (ij_ab != ji_ab)
256 Xab_.accumulate_element(ji_ab,kl_ab,Xab_jikl);
257 if (kl_ab != lk_ab)
258 Xab_.accumulate_element(ij_ab,lk_ab,Xab_ijlk);
259 if (ij_ab != ji_ab && kl_ab != lk_ab)
260 Xab_.accumulate_element(ji_ab,lk_ab,Xab_jilk);
261 if (ij_aa != -1 && kl_aa != -1)
262 Xaa_.accumulate_element(ij_aa,kl_aa,Xaa_ijkl);
263 Bab_.accumulate_element(ij_ab,kl_ab,Tab_ijkl);
264 if (ij_ab != ji_ab)
265 Bab_.accumulate_element(ji_ab,kl_ab,Tab_jikl);
266 if (kl_ab != lk_ab)
267 Bab_.accumulate_element(ij_ab,lk_ab,Tab_ijlk);
268 if (ij_ab != ji_ab && kl_ab != lk_ab)
269 Bab_.accumulate_element(ji_ab,lk_ab,Tab_jilk);
270 if (ij_aa != -1 && kl_aa != -1)
271 Baa_.accumulate_element(ij_aa,kl_aa,Taa_ijkl);
272 tim_exit("MO ints contraction");
273
274#if PRINT_R12_INTERMED
275 if (ij_ab != ji_ab && kl_ab != lk_ab)
276 printf("Vaa[%d][%d] = %lf\n",ij_aa,kl_aa,Vaa_ij[kl_aa]);
277 printf("Vab[%d][%d] = %lf\n",ij_ab,kl_ab,Vab_ij[kl_ab]);
278 if (ij_ab != ji_ab)
279 printf("Vab[%d][%d] = %lf\n",ji_ab,kl_ab,Vab_ji[kl_ab]);
280 if (kl_ab != lk_ab)
281 printf("Vab[%d][%d] = %lf\n",ij_ab,lk_ab,Vab_ij[lk_ab]);
282 if (ij_ab != ji_ab && kl_ab != lk_ab)
283 printf("Vab[%d][%d] = %lf\n",ji_ab,lk_ab,Vab_ji[lk_ab]);
284 if (ij_ab != ji_ab && kl_ab != lk_ab)
285 printf("Xaa[%d][%d] = %lf\n",ij_aa,kl_aa,Xaa_ij[kl_aa]);
286 printf("Xab[%d][%d] = %lf\n",ij_ab,kl_ab,Xab_ij[kl_ab]);
287 if (ij_ab != ji_ab)
288 printf("Xab[%d][%d] = %lf\n",ji_ab,kl_ab,Xab_ji[kl_ab]);
289 if (kl_ab != lk_ab)
290 printf("Xab[%d][%d] = %lf\n",ij_ab,lk_ab,Xab_ij[lk_ab]);
291 if (ij_ab != ji_ab && kl_ab != lk_ab)
292 printf("Xab[%d][%d] = %lf\n",ji_ab,lk_ab,Xab_ji[lk_ab]);
293 if (ij_ab != ji_ab && kl_ab != lk_ab)
294 printf("Taa[%d][%d] = %lf\n",ij_aa,kl_aa,Taa_ij[kl_aa]);
295 printf("Tab[%d][%d] = %lf\n",ij_ab,kl_ab,Tab_ij[kl_ab]);
296 if (ij_ab != ji_ab)
297 printf("Tab[%d][%d] = %lf\n",ji_ab,kl_ab,Tab_ji[kl_ab]);
298 if (kl_ab != lk_ab)
299 printf("Tab[%d][%d] = %lf\n",ij_ab,lk_ab,Tab_ij[lk_ab]);
300 if (ij_ab != ji_ab && kl_ab != lk_ab)
301 printf("Tab[%d][%d] = %lf\n",ji_ab,lk_ab,Tab_ji[lk_ab]);
302#endif
303 ijky_acc->release_pair_block(i,j,R12IntsAcc::r12);
304 ijky_acc->release_pair_block(j,i,R12IntsAcc::r12);
305 }
306 ijky_acc->release_pair_block(k,l,R12IntsAcc::eri);
307 ijky_acc->release_pair_block(k,l,R12IntsAcc::r12);
308 ijky_acc->release_pair_block(k,l,R12IntsAcc::r12t1);
309 ijky_acc->release_pair_block(k,l,R12IntsAcc::r12t2);
310 ijky_acc->release_pair_block(l,k,R12IntsAcc::eri);
311 ijky_acc->release_pair_block(l,k,R12IntsAcc::r12);
312 ijky_acc->release_pair_block(l,k,R12IntsAcc::r12t1);
313 ijky_acc->release_pair_block(l,k,R12IntsAcc::r12t2);
314 }
315 }
316 // Tasks that don't do any work here still need to create these timers
317 tim_enter("MO ints retrieve");
318 tim_exit("MO ints retrieve");
319 tim_enter("MO ints contraction");
320 tim_exit("MO ints contraction");
321
322 tim_exit("intermediates");
323 ExEnv::out0() << indent << "End of computation of intermediates" << endl;
324 ijky_acc->deactivate();
325
326 // Symmetrize B intermediate
327 for(int ij=0;ij<naa;ij++)
328 for(int kl=0;kl<=ij;kl++) {
329 double belem = 0.5*(Baa_->get_element(ij,kl) + Baa_->get_element(kl,ij));
330 Baa_->set_element(ij,kl,belem);
331 Baa_->set_element(kl,ij,belem);
332 }
333
334 for(int ij=0;ij<nab;ij++)
335 for(int kl=0;kl<=ij;kl++) {
336 double belem = 0.5*(Bab_->get_element(ij,kl) + Bab_->get_element(kl,ij));
337 Bab_->set_element(ij,kl,belem);
338 Bab_->set_element(kl,ij,belem);
339 }
340
341 globally_sum_intermeds_();
342
343 ExEnv::out0() << decindent;
344 ExEnv::out0() << indent << "Exited ABS A (GEBC) intermediates evaluator" << endl;
345
346 tim_exit("mp2-r12a intermeds");
347 checkpoint_();
348}
349
350////////////////////////////////////////////////////////////////////////////
351
352// Local Variables:
353// mode: c++
354// c-file-style: "CLJ-CONDENSED"
355// End:
Note: See TracBrowser for help on using the repository browser.