| [0b990d] | 1 | //
|
|---|
| 2 | // compute_a_gebc.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 |
|
|---|
| 49 | using namespace std;
|
|---|
| 50 | using namespace sc;
|
|---|
| 51 |
|
|---|
| 52 | #define PRINT4Q_MP2 0
|
|---|
| 53 | #define PRINT_R12_INTERMED 0
|
|---|
| 54 |
|
|---|
| 55 | #define COMPUTE_AB_BLOCK_ONLY 0
|
|---|
| 56 | #define COMPUTE_MA_BLOCK_ONLY 0
|
|---|
| 57 |
|
|---|
| 58 | void
|
|---|
| 59 | R12IntEval::obs_contrib_to_VXB_gebc_vbseqobs_()
|
|---|
| 60 | {
|
|---|
| 61 | if (evaluated_)
|
|---|
| 62 | return;
|
|---|
| 63 | LinearR12::ABSMethod abs_method = r12info_->abs_method();
|
|---|
| 64 | Ref<MessageGrp> msg = r12info_->msg();
|
|---|
| 65 | Ref<MemoryGrp> mem = r12info_->mem();
|
|---|
| 66 | Ref<ThreadGrp> thr = r12info_->thr();
|
|---|
| 67 | const int num_te_types = 3;
|
|---|
| 68 | enum te_types {eri=0, r12=1, r12t1=2};
|
|---|
| 69 |
|
|---|
| 70 | tim_enter("mp2-r12a intermeds");
|
|---|
| 71 |
|
|---|
| 72 | int me = msg->me();
|
|---|
| 73 | int nproc = msg->n();
|
|---|
| 74 |
|
|---|
| 75 | ExEnv::out0() << endl << indent
|
|---|
| 76 | << "Entered OBS A (GEBC) intermediates evaluator" << endl;
|
|---|
| 77 | ExEnv::out0() << incindent;
|
|---|
| 78 |
|
|---|
| 79 | // Do the AO->MO transform
|
|---|
| 80 | Ref<TwoBodyMOIntsTransform> ipjq_tform = get_tform_("(ip|jq)");
|
|---|
| 81 | Ref<R12IntsAcc> ipjq_acc = ipjq_tform->ints_acc();
|
|---|
| 82 | if (!ipjq_acc->is_committed()) {
|
|---|
| 83 | ipjq_tform->set_num_te_types(num_te_types);
|
|---|
| 84 | ipjq_tform->compute();
|
|---|
| 85 | }
|
|---|
| 86 | if (num_te_types != ipjq_acc->num_te_types())
|
|---|
| 87 | throw std::runtime_error("R12IntEval::obs_contrib_to_VXB_gebc() -- number of MO integral types is wrong");
|
|---|
| 88 |
|
|---|
| 89 | int nocc = r12info_->nocc();
|
|---|
| 90 | int nocc_act = r12info_->nocc_act();
|
|---|
| 91 | int nfzc = r12info_->nfzc();
|
|---|
| 92 | int nfzv = r12info_->nfzv();
|
|---|
| 93 | int noso = r12info_->mo_space()->rank();
|
|---|
| 94 | int nvir = noso - nocc;
|
|---|
| 95 |
|
|---|
| 96 | /*--------------------------------
|
|---|
| 97 | Compute MP2-R12/A intermediates
|
|---|
| 98 | and collect on node0
|
|---|
| 99 | --------------------------------*/
|
|---|
| 100 | ExEnv::out0() << indent << "Begin computation of intermediates" << endl;
|
|---|
| 101 | tim_enter("intermediates");
|
|---|
| 102 | SpatialMOPairIter_eq ij_iter(r12info_->act_occ_space());
|
|---|
| 103 | SpatialMOPairIter_eq kl_iter(r12info_->act_occ_space());
|
|---|
| 104 | int naa = ij_iter.nij_aa(); // Number of alpha-alpha pairs (i > j)
|
|---|
| 105 | int nab = ij_iter.nij_ab(); // Number of alpha-beta pairs
|
|---|
| 106 | if (debug_) {
|
|---|
| 107 | ExEnv::out0() << indent << "naa = " << naa << endl;
|
|---|
| 108 | ExEnv::out0() << indent << "nab = " << nab << endl;
|
|---|
| 109 | }
|
|---|
| 110 |
|
|---|
| 111 | // Compute intermediates
|
|---|
| 112 | if (debug_)
|
|---|
| 113 | ExEnv::out0() << indent << "Ready to compute MP2-R12/A (GEBC) intermediates" << endl;
|
|---|
| 114 |
|
|---|
| 115 | // Compute the number of tasks that have full access to the integrals
|
|---|
| 116 | // and split the work among them
|
|---|
| 117 | vector<int> proc_with_ints;
|
|---|
| 118 | int nproc_with_ints = tasks_with_ints_(ipjq_acc,proc_with_ints);
|
|---|
| 119 |
|
|---|
| 120 |
|
|---|
| 121 | //////////////////////////////////////////////////////////////
|
|---|
| 122 | //
|
|---|
| 123 | // Evaluation of the intermediates proceeds as follows:
|
|---|
| 124 | //
|
|---|
| 125 | // loop over batches of kl, k >= l, 0<=k,l<nocc_act
|
|---|
| 126 | // load (kl|xy), (kl| [T1,r12] |xy), and (lk| [T1,r12] |xy)
|
|---|
| 127 | // (aka kl-sets) into memory
|
|---|
| 128 | //
|
|---|
| 129 | // loop over batches of ij, i>=j, 0<=i,j<nocc_act
|
|---|
| 130 | // load (ij|r12|xy) into memory
|
|---|
| 131 | // (aka ij-sets) into memory
|
|---|
| 132 | // compute V[ij][kl] and T[ij][kl] for all ij and kl in
|
|---|
| 133 | // the "direct product" batch
|
|---|
| 134 | // end ij loop
|
|---|
| 135 | // end kl loop
|
|---|
| 136 | //
|
|---|
| 137 | /////////////////////////////////////////////////////////////////////////////////
|
|---|
| 138 |
|
|---|
| 139 | bool two_basis_form = (r12info_->basis() != r12info_->basis_ri());
|
|---|
| 140 | double pfac_xy_1, pfac_xy_2;
|
|---|
| 141 | if (two_basis_form &&
|
|---|
| 142 | ( abs_method == LinearR12::ABS_ABS ||
|
|---|
| 143 | abs_method == LinearR12::ABS_ABSPlus ) ) {
|
|---|
| 144 | pfac_xy_1 = 0.5;
|
|---|
| 145 | pfac_xy_2 = -0.5;
|
|---|
| 146 | }
|
|---|
| 147 | else {
|
|---|
| 148 | pfac_xy_1 = 0.5;
|
|---|
| 149 | pfac_xy_2 = 0.5;
|
|---|
| 150 | }
|
|---|
| 151 |
|
|---|
| 152 | if (ipjq_acc->has_access(me)) {
|
|---|
| 153 |
|
|---|
| 154 | for(kl_iter.start();int(kl_iter);kl_iter.next()) {
|
|---|
| 155 |
|
|---|
| 156 | const int kl = kl_iter.ij();
|
|---|
| 157 | // Figure out if this task will handle this kl
|
|---|
| 158 | int kl_proc = kl%nproc_with_ints;
|
|---|
| 159 | if (kl_proc != proc_with_ints[me])
|
|---|
| 160 | continue;
|
|---|
| 161 | const int k = kl_iter.i();
|
|---|
| 162 | const int l = kl_iter.j();
|
|---|
| 163 | const int kl_aa = kl_iter.ij_aa();
|
|---|
| 164 | const int kl_ab = kl_iter.ij_ab();
|
|---|
| 165 | const int lk_ab = kl_iter.ij_ba();
|
|---|
| 166 |
|
|---|
| 167 | if (debug_)
|
|---|
| 168 | ExEnv::outn() << indent << "task " << me << ": working on (k,l) = " << k << "," << l << " " << endl;
|
|---|
| 169 |
|
|---|
| 170 | // Get (|1/r12|), (|r12|), and (|[r12,T1]|) integrals
|
|---|
| 171 | tim_enter("MO ints retrieve");
|
|---|
| 172 | double *klxy_buf_eri = ipjq_acc->retrieve_pair_block(k,l,R12IntsAcc::eri);
|
|---|
| 173 | double *klxy_buf_r12 = ipjq_acc->retrieve_pair_block(k,l,R12IntsAcc::r12);
|
|---|
| 174 | double *klxy_buf_r12t1 = ipjq_acc->retrieve_pair_block(k,l,R12IntsAcc::r12t1);
|
|---|
| 175 | double *lkxy_buf_r12t1 = ipjq_acc->retrieve_pair_block(l,k,R12IntsAcc::r12t1);
|
|---|
| 176 | tim_exit("MO ints retrieve");
|
|---|
| 177 |
|
|---|
| 178 | if (debug_)
|
|---|
| 179 | ExEnv::outn() << indent << "task " << me << ": obtained kl blocks" << endl;
|
|---|
| 180 |
|
|---|
| 181 | // Compute MP2 energies
|
|---|
| 182 | RefDiagSCMatrix act_occ_evals = r12info_->act_occ_space()->evals();
|
|---|
| 183 | RefDiagSCMatrix all_evals = r12info_->obs_space()->evals();
|
|---|
| 184 | double emp2_aa = 0.0;
|
|---|
| 185 | double emp2_ab = 0.0;
|
|---|
| 186 | for(int a=nocc; a<noso; a++) {
|
|---|
| 187 | for(int b=nocc; b<noso; b++) {
|
|---|
| 188 | const int ab_offset = a*noso+b;
|
|---|
| 189 | const int ba_offset = b*noso+a;
|
|---|
| 190 | const double oo_delta_ijab = 1.0/(act_occ_evals(k)+act_occ_evals(l)-all_evals(a)-all_evals(b));
|
|---|
| 191 | const double eri_kalb = klxy_buf_eri[ab_offset];
|
|---|
| 192 | const double eri_kbla = klxy_buf_eri[ba_offset];
|
|---|
| 193 | emp2_ab += 0.5*(eri_kalb * eri_kalb + eri_kbla * eri_kbla) * oo_delta_ijab;
|
|---|
| 194 | if (kl_aa != -1) {
|
|---|
| 195 | emp2_aa += (eri_kalb - eri_kbla) * (eri_kalb - eri_kbla) * oo_delta_ijab;
|
|---|
| 196 | }
|
|---|
| 197 | }
|
|---|
| 198 | }
|
|---|
| 199 | emp2pair_ab_.accumulate_element(kl_ab,emp2_ab);
|
|---|
| 200 | if (kl_ab != lk_ab)
|
|---|
| 201 | emp2pair_ab_.accumulate_element(lk_ab,emp2_ab);
|
|---|
| 202 | if (kl_aa != -1)
|
|---|
| 203 | emp2pair_aa_.accumulate_element(kl_aa,emp2_aa);
|
|---|
| 204 |
|
|---|
| 205 | // 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
|
|---|
| 206 | for(ij_iter.start(kl+1);int(ij_iter);ij_iter.next()) {
|
|---|
| 207 |
|
|---|
| 208 | const int i = ij_iter.i();
|
|---|
| 209 | const int j = ij_iter.j();
|
|---|
| 210 | const int ij_aa = ij_iter.ij_aa();
|
|---|
| 211 | const int ij_ab = ij_iter.ij_ab();
|
|---|
| 212 | const int ji_ab = ij_iter.ij_ba();
|
|---|
| 213 |
|
|---|
| 214 | if (debug_)
|
|---|
| 215 | ExEnv::outn() << indent << "task " << me << ": (k,l) = " << k << "," << l << ": (i,j) = " << i << "," << j << endl;
|
|---|
| 216 |
|
|---|
| 217 | tim_enter("MO ints retrieve");
|
|---|
| 218 | double *ijxy_buf_r12 = ipjq_acc->retrieve_pair_block(i,j,R12IntsAcc::r12);
|
|---|
| 219 | tim_exit("MO ints retrieve");
|
|---|
| 220 |
|
|---|
| 221 | if (debug_)
|
|---|
| 222 | ExEnv::outn() << indent << "task " << me << ": obtained ij blocks" << endl;
|
|---|
| 223 |
|
|---|
| 224 |
|
|---|
| 225 | tim_enter("MO ints contraction");
|
|---|
| 226 | double Vaa_ijkl, Vab_ijkl, Vab_jikl, Vab_ijlk, Vab_jilk;
|
|---|
| 227 | double Xaa_ijkl, Xab_ijkl, Xab_jikl, Xab_ijlk, Xab_jilk;
|
|---|
| 228 | double Taa_ijkl, Tab_ijkl, Tab_jikl, Tab_ijlk, Tab_jilk;
|
|---|
| 229 | Vaa_ijkl = Vab_ijkl = Vab_jikl = Vab_ijlk = Vab_jilk = 0.0;
|
|---|
| 230 | Xaa_ijkl = Xab_ijkl = Xab_jikl = Xab_ijlk = Xab_jilk = 0.0;
|
|---|
| 231 | Taa_ijkl = Tab_ijkl = Tab_jikl = Tab_ijlk = Tab_jilk = 0.0;
|
|---|
| 232 |
|
|---|
| 233 | for(int y=0;y<noso;y++) {
|
|---|
| 234 | double pfac_xy;
|
|---|
| 235 | if (y >= nocc)
|
|---|
| 236 | pfac_xy = pfac_xy_1;
|
|---|
| 237 | else
|
|---|
| 238 | pfac_xy = pfac_xy_2;
|
|---|
| 239 | for(int x=0;x<noso;x++) {
|
|---|
| 240 |
|
|---|
| 241 | #if COMPUTE_AB_BLOCK_ONLY
|
|---|
| 242 | if (y < nocc || x < nocc) {
|
|---|
| 243 | pfac_xy = 0.0;
|
|---|
| 244 | }
|
|---|
| 245 | else {
|
|---|
| 246 | if (y >= nocc)
|
|---|
| 247 | pfac_xy = pfac_xy_1;
|
|---|
| 248 | else
|
|---|
| 249 | pfac_xy = pfac_xy_2;
|
|---|
| 250 | }
|
|---|
| 251 | #endif
|
|---|
| 252 | #if COMPUTE_MA_BLOCK_ONLY
|
|---|
| 253 | if ((y < nocc && x < nocc) ||
|
|---|
| 254 | (y >= nocc && x >= nocc)) {
|
|---|
| 255 | pfac_xy = 0.0;
|
|---|
| 256 | }
|
|---|
| 257 | else {
|
|---|
| 258 | pfac_xy = 0.5;
|
|---|
| 259 | }
|
|---|
| 260 | #endif
|
|---|
| 261 |
|
|---|
| 262 | int yx_offset = y*noso+x;
|
|---|
| 263 | int xy_offset = x*noso+y;
|
|---|
| 264 | double ij_r12_xy = ijxy_buf_r12[xy_offset];
|
|---|
| 265 | double ij_r12_yx = ijxy_buf_r12[yx_offset];
|
|---|
| 266 | double kl_eri_xy = klxy_buf_eri[xy_offset];
|
|---|
| 267 | double kl_eri_yx = klxy_buf_eri[yx_offset];
|
|---|
| 268 | Vab_ijkl -= pfac_xy * (ij_r12_xy * kl_eri_xy + ij_r12_yx * kl_eri_yx);
|
|---|
| 269 | if (ij_ab != ji_ab)
|
|---|
| 270 | Vab_jikl -= pfac_xy * (ij_r12_yx * kl_eri_xy + ij_r12_xy * kl_eri_yx);
|
|---|
| 271 | if (kl_ab != lk_ab)
|
|---|
| 272 | Vab_ijlk -= pfac_xy * (ij_r12_xy * kl_eri_yx + ij_r12_yx * kl_eri_xy);
|
|---|
| 273 | if (ij_ab != ji_ab && kl_ab != lk_ab) {
|
|---|
| 274 | Vab_jilk -= pfac_xy * (ij_r12_yx * kl_eri_yx + ij_r12_xy * kl_eri_xy);
|
|---|
| 275 | }
|
|---|
| 276 | if (ij_aa != -1 && kl_aa != -1) {
|
|---|
| 277 | Vaa_ijkl -= pfac_xy * (ij_r12_xy - ij_r12_yx)*(kl_eri_xy - kl_eri_yx);
|
|---|
| 278 | }
|
|---|
| 279 | double kl_r12_xy = klxy_buf_r12[xy_offset];
|
|---|
| 280 | double kl_r12_yx = klxy_buf_r12[yx_offset];
|
|---|
| 281 | Xab_ijkl -= pfac_xy * (ij_r12_xy * kl_r12_xy + ij_r12_yx * kl_r12_yx);
|
|---|
| 282 | if (ij_ab != ji_ab)
|
|---|
| 283 | Xab_jikl -= pfac_xy * (ij_r12_yx * kl_r12_xy + ij_r12_xy * kl_r12_yx);
|
|---|
| 284 | if (kl_ab != lk_ab)
|
|---|
| 285 | Xab_ijlk -= pfac_xy * (ij_r12_xy * kl_r12_yx + ij_r12_yx * kl_r12_xy);
|
|---|
| 286 | if (ij_ab != ji_ab && kl_ab != lk_ab) {
|
|---|
| 287 | Xab_jilk -= pfac_xy * (ij_r12_yx * kl_r12_yx + ij_r12_xy * kl_r12_xy);
|
|---|
| 288 | }
|
|---|
| 289 | if (ij_aa != -1 && kl_aa != -1) {
|
|---|
| 290 | Xaa_ijkl -= pfac_xy * (ij_r12_xy - ij_r12_yx)*(kl_r12_xy - kl_r12_yx);
|
|---|
| 291 | }
|
|---|
| 292 | double kl_r12t1_xy = klxy_buf_r12t1[xy_offset];
|
|---|
| 293 | double kl_r12t1_yx = klxy_buf_r12t1[yx_offset];
|
|---|
| 294 | double lk_r12t1_xy = lkxy_buf_r12t1[xy_offset];
|
|---|
| 295 | double lk_r12t1_yx = lkxy_buf_r12t1[yx_offset];
|
|---|
| 296 | double kl_Tr12_xy = -kl_r12t1_xy-lk_r12t1_yx;
|
|---|
| 297 | double kl_Tr12_yx = -kl_r12t1_yx-lk_r12t1_xy;
|
|---|
| 298 | Tab_ijkl += pfac_xy * (ij_r12_xy * kl_Tr12_xy + ij_r12_yx * kl_Tr12_yx);
|
|---|
| 299 | if (ij_ab != ji_ab)
|
|---|
| 300 | Tab_jikl += pfac_xy * (ij_r12_yx * kl_Tr12_xy + ij_r12_xy * kl_Tr12_yx);
|
|---|
| 301 | if (kl_ab != lk_ab)
|
|---|
| 302 | Tab_ijlk += pfac_xy * (ij_r12_xy * kl_Tr12_yx + ij_r12_yx * kl_Tr12_xy);
|
|---|
| 303 | if (ij_ab != ji_ab && kl_ab != lk_ab) {
|
|---|
| 304 | Tab_jilk += pfac_xy * (ij_r12_yx * kl_Tr12_yx + ij_r12_xy * kl_Tr12_xy);
|
|---|
| 305 | }
|
|---|
| 306 | if (ij_aa != -1 && kl_aa != -1) {
|
|---|
| 307 | Taa_ijkl += pfac_xy * (ij_r12_xy - ij_r12_yx)*(kl_Tr12_xy - kl_Tr12_yx);
|
|---|
| 308 | }
|
|---|
| 309 | }
|
|---|
| 310 | }
|
|---|
| 311 | Vab_.accumulate_element(ij_ab,kl_ab,Vab_ijkl);
|
|---|
| 312 | if (ij_ab != ji_ab)
|
|---|
| 313 | Vab_.accumulate_element(ji_ab,kl_ab,Vab_jikl);
|
|---|
| 314 | if (kl_ab != lk_ab)
|
|---|
| 315 | Vab_.accumulate_element(ij_ab,lk_ab,Vab_ijlk);
|
|---|
| 316 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 317 | Vab_.accumulate_element(ji_ab,lk_ab,Vab_jilk);
|
|---|
| 318 | if (ij_aa != -1 && kl_aa != -1)
|
|---|
| 319 | Vaa_.accumulate_element(ij_aa,kl_aa,Vaa_ijkl);
|
|---|
| 320 | Xab_.accumulate_element(ij_ab,kl_ab,Xab_ijkl);
|
|---|
| 321 | if (ij_ab != ji_ab)
|
|---|
| 322 | Xab_.accumulate_element(ji_ab,kl_ab,Xab_jikl);
|
|---|
| 323 | if (kl_ab != lk_ab)
|
|---|
| 324 | Xab_.accumulate_element(ij_ab,lk_ab,Xab_ijlk);
|
|---|
| 325 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 326 | Xab_.accumulate_element(ji_ab,lk_ab,Xab_jilk);
|
|---|
| 327 | if (ij_aa != -1 && kl_aa != -1)
|
|---|
| 328 | Xaa_.accumulate_element(ij_aa,kl_aa,Xaa_ijkl);
|
|---|
| 329 | Bab_.accumulate_element(ij_ab,kl_ab,Tab_ijkl);
|
|---|
| 330 | if (ij_ab != ji_ab)
|
|---|
| 331 | Bab_.accumulate_element(ji_ab,kl_ab,Tab_jikl);
|
|---|
| 332 | if (kl_ab != lk_ab)
|
|---|
| 333 | Bab_.accumulate_element(ij_ab,lk_ab,Tab_ijlk);
|
|---|
| 334 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 335 | Bab_.accumulate_element(ji_ab,lk_ab,Tab_jilk);
|
|---|
| 336 | if (ij_aa != -1 && kl_aa != -1)
|
|---|
| 337 | Baa_.accumulate_element(ij_aa,kl_aa,Taa_ijkl);
|
|---|
| 338 | tim_exit("MO ints contraction");
|
|---|
| 339 |
|
|---|
| 340 | #if PRINT_R12_INTERMED
|
|---|
| 341 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 342 | printf("Vaa[%d][%d] = %lf\n",ij_aa,kl_aa,Vaa_ij[kl_aa]);
|
|---|
| 343 | printf("Vab[%d][%d] = %lf\n",ij_ab,kl_ab,Vab_ij[kl_ab]);
|
|---|
| 344 | if (ij_ab != ji_ab)
|
|---|
| 345 | printf("Vab[%d][%d] = %lf\n",ji_ab,kl_ab,Vab_ji[kl_ab]);
|
|---|
| 346 | if (kl_ab != lk_ab)
|
|---|
| 347 | printf("Vab[%d][%d] = %lf\n",ij_ab,lk_ab,Vab_ij[lk_ab]);
|
|---|
| 348 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 349 | printf("Vab[%d][%d] = %lf\n",ji_ab,lk_ab,Vab_ji[lk_ab]);
|
|---|
| 350 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 351 | printf("Xaa[%d][%d] = %lf\n",ij_aa,kl_aa,Xaa_ij[kl_aa]);
|
|---|
| 352 | printf("Xab[%d][%d] = %lf\n",ij_ab,kl_ab,Xab_ij[kl_ab]);
|
|---|
| 353 | if (ij_ab != ji_ab)
|
|---|
| 354 | printf("Xab[%d][%d] = %lf\n",ji_ab,kl_ab,Xab_ji[kl_ab]);
|
|---|
| 355 | if (kl_ab != lk_ab)
|
|---|
| 356 | printf("Xab[%d][%d] = %lf\n",ij_ab,lk_ab,Xab_ij[lk_ab]);
|
|---|
| 357 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 358 | printf("Xab[%d][%d] = %lf\n",ji_ab,lk_ab,Xab_ji[lk_ab]);
|
|---|
| 359 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 360 | printf("Taa[%d][%d] = %lf\n",ij_aa,kl_aa,Taa_ij[kl_aa]);
|
|---|
| 361 | printf("Tab[%d][%d] = %lf\n",ij_ab,kl_ab,Tab_ij[kl_ab]);
|
|---|
| 362 | if (ij_ab != ji_ab)
|
|---|
| 363 | printf("Tab[%d][%d] = %lf\n",ji_ab,kl_ab,Tab_ji[kl_ab]);
|
|---|
| 364 | if (kl_ab != lk_ab)
|
|---|
| 365 | printf("Tab[%d][%d] = %lf\n",ij_ab,lk_ab,Tab_ij[lk_ab]);
|
|---|
| 366 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 367 | printf("Tab[%d][%d] = %lf\n",ji_ab,lk_ab,Tab_ji[lk_ab]);
|
|---|
| 368 | #endif
|
|---|
| 369 | ipjq_acc->release_pair_block(i,j,R12IntsAcc::r12);
|
|---|
| 370 | }
|
|---|
| 371 | ipjq_acc->release_pair_block(k,l,R12IntsAcc::eri);
|
|---|
| 372 | ipjq_acc->release_pair_block(k,l,R12IntsAcc::r12);
|
|---|
| 373 | ipjq_acc->release_pair_block(k,l,R12IntsAcc::r12t1);
|
|---|
| 374 | ipjq_acc->release_pair_block(l,k,R12IntsAcc::r12t1);
|
|---|
| 375 | }
|
|---|
| 376 | }
|
|---|
| 377 | // Tasks that don't do any work here still need to create these timers
|
|---|
| 378 | tim_enter("MO ints retrieve");
|
|---|
| 379 | tim_exit("MO ints retrieve");
|
|---|
| 380 | tim_enter("MO ints contraction");
|
|---|
| 381 | tim_exit("MO ints contraction");
|
|---|
| 382 |
|
|---|
| 383 | tim_exit("intermediates");
|
|---|
| 384 | ExEnv::out0() << indent << "End of computation of intermediates" << endl;
|
|---|
| 385 | ipjq_acc->deactivate();
|
|---|
| 386 |
|
|---|
| 387 | // Symmetrize B intermediate
|
|---|
| 388 | for(int ij=0;ij<naa;ij++)
|
|---|
| 389 | for(int kl=0;kl<=ij;kl++) {
|
|---|
| 390 | double belem = 0.5*(Baa_->get_element(ij,kl) + Baa_->get_element(kl,ij));
|
|---|
| 391 | Baa_->set_element(ij,kl,belem);
|
|---|
| 392 | Baa_->set_element(kl,ij,belem);
|
|---|
| 393 | }
|
|---|
| 394 |
|
|---|
| 395 | for(int ij=0;ij<nab;ij++)
|
|---|
| 396 | for(int kl=0;kl<=ij;kl++) {
|
|---|
| 397 | double belem = 0.5*(Bab_->get_element(ij,kl) + Bab_->get_element(kl,ij));
|
|---|
| 398 | Bab_->set_element(ij,kl,belem);
|
|---|
| 399 | Bab_->set_element(kl,ij,belem);
|
|---|
| 400 | }
|
|---|
| 401 |
|
|---|
| 402 | globally_sum_intermeds_();
|
|---|
| 403 |
|
|---|
| 404 | ExEnv::out0() << decindent;
|
|---|
| 405 | ExEnv::out0() << indent << "Exited OBS A (GEBC) intermediates evaluator" << endl;
|
|---|
| 406 |
|
|---|
| 407 | tim_exit("mp2-r12a intermeds");
|
|---|
| 408 | checkpoint_();
|
|---|
| 409 | }
|
|---|
| 410 |
|
|---|
| 411 | ////////////////////////////////////////////////////////////////////////////
|
|---|
| 412 |
|
|---|
| 413 | // Local Variables:
|
|---|
| 414 | // mode: c++
|
|---|
| 415 | // c-file-style: "CLJ-CONDENSED"
|
|---|
| 416 | // End:
|
|---|
