| 1 | //
|
|---|
| 2 | // compute_vxb_a_symm.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 <sstream>
|
|---|
| 30 | #include <stdlib.h>
|
|---|
| 31 | #include <math.h>
|
|---|
| 32 | #include <limits.h>
|
|---|
| 33 |
|
|---|
| 34 | #include <scconfig.h>
|
|---|
| 35 | #include <util/misc/formio.h>
|
|---|
| 36 | #include <util/misc/timer.h>
|
|---|
| 37 | #include <util/class/class.h>
|
|---|
| 38 | #include <util/state/state.h>
|
|---|
| 39 | #include <util/state/state_text.h>
|
|---|
| 40 | #include <util/state/state_bin.h>
|
|---|
| 41 | #include <math/scmat/matrix.h>
|
|---|
| 42 | #include <chemistry/molecule/molecule.h>
|
|---|
| 43 | #include <chemistry/qc/basis/integral.h>
|
|---|
| 44 | #include <chemistry/qc/mbpt/bzerofast.h>
|
|---|
| 45 | #include <chemistry/qc/mbptr12/r12ia.h>
|
|---|
| 46 | #include <chemistry/qc/mbptr12/vxb_eval_info.h>
|
|---|
| 47 | #include <chemistry/qc/mbptr12/pairiter.h>
|
|---|
| 48 | #include <chemistry/qc/mbptr12/r12int_eval.h>
|
|---|
| 49 |
|
|---|
| 50 | using namespace std;
|
|---|
| 51 | using namespace sc;
|
|---|
| 52 |
|
|---|
| 53 | #define PRINT_R12_INTERMED 0
|
|---|
| 54 |
|
|---|
| 55 | void
|
|---|
| 56 | R12IntEval::contrib_to_VXB_a_symm_(const std::string& tform_name)
|
|---|
| 57 | {
|
|---|
| 58 | if (evaluated_)
|
|---|
| 59 | return;
|
|---|
| 60 | LinearR12::ABSMethod abs_method = r12info_->abs_method();
|
|---|
| 61 | Ref<MessageGrp> msg = r12info_->msg();
|
|---|
| 62 | Ref<MemoryGrp> mem = r12info_->mem();
|
|---|
| 63 | Ref<ThreadGrp> thr = r12info_->thr();
|
|---|
| 64 | const int num_te_types = 3;
|
|---|
| 65 | enum te_types {eri=0, r12=1, r12t1=2};
|
|---|
| 66 |
|
|---|
| 67 | tim_enter("mp2-r12a intermeds (symmetric term)");
|
|---|
| 68 |
|
|---|
| 69 | int me = msg->me();
|
|---|
| 70 | int nproc = msg->n();
|
|---|
| 71 |
|
|---|
| 72 | // Do the AO->MO transform
|
|---|
| 73 | Ref<TwoBodyMOIntsTransform> ipjq_tform = get_tform_(tform_name);
|
|---|
| 74 | Ref<R12IntsAcc> ijpq_acc = ipjq_tform->ints_acc();
|
|---|
| 75 | if (ijpq_acc.null() || !ijpq_acc->is_committed())
|
|---|
| 76 | ipjq_tform->compute();
|
|---|
| 77 | if (!ijpq_acc->is_active())
|
|---|
| 78 | ijpq_acc->activate();
|
|---|
| 79 | if (num_te_types != ijpq_acc->num_te_types())
|
|---|
| 80 | throw std::runtime_error("R12IntEval::contrib_to_VXB_a_symm_() -- number of MO integral types is wrong");
|
|---|
| 81 |
|
|---|
| 82 | if (ipjq_tform->space2() != ipjq_tform->space4())
|
|---|
| 83 | throw std::runtime_error("R12IntEval::contrib_to_VXB_a_symm_() -- wrong type of transform is requested (space2 != space4)");
|
|---|
| 84 | Ref<MOIndexSpace> mospace = ipjq_tform->space2();
|
|---|
| 85 |
|
|---|
| 86 | ExEnv::out0() << endl << indent
|
|---|
| 87 | << "Entered \"" << mospace->name()
|
|---|
| 88 | << "\" A (GEBC) intermediates evaluator" << endl;
|
|---|
| 89 | ExEnv::out0() << incindent;
|
|---|
| 90 |
|
|---|
| 91 | int nocc_act = r12info_->nocc_act();
|
|---|
| 92 | int noso = mospace->rank();
|
|---|
| 93 |
|
|---|
| 94 | /*--------------------------------
|
|---|
| 95 | Compute MP2-R12/A intermediates
|
|---|
| 96 | and collect on node0
|
|---|
| 97 | --------------------------------*/
|
|---|
| 98 | ExEnv::out0() << indent << "Begin computation of intermediates" << endl;
|
|---|
| 99 | tim_enter("intermediates");
|
|---|
| 100 | SpatialMOPairIter_eq ij_iter(r12info_->act_occ_space());
|
|---|
| 101 | SpatialMOPairIter_eq kl_iter(r12info_->act_occ_space());
|
|---|
| 102 | int naa = ij_iter.nij_aa(); // Number of alpha-alpha pairs (i > j)
|
|---|
| 103 | int nab = ij_iter.nij_ab(); // Number of alpha-beta pairs
|
|---|
| 104 | if (debug_) {
|
|---|
| 105 | ExEnv::out0() << indent << "naa = " << naa << endl;
|
|---|
| 106 | ExEnv::out0() << indent << "nab = " << nab << endl;
|
|---|
| 107 | }
|
|---|
| 108 |
|
|---|
| 109 | // Compute intermediates
|
|---|
| 110 | if (debug_)
|
|---|
| 111 | ExEnv::out0() << indent << "Ready to compute MP2-R12/A (GEBC) intermediates" << endl;
|
|---|
| 112 |
|
|---|
| 113 | // Compute the number of tasks that have full access to the integrals
|
|---|
| 114 | // and split the work among them
|
|---|
| 115 | vector<int> proc_with_ints;
|
|---|
| 116 | int nproc_with_ints = tasks_with_ints_(ijpq_acc,proc_with_ints);
|
|---|
| 117 |
|
|---|
| 118 |
|
|---|
| 119 | //////////////////////////////////////////////////////////////
|
|---|
| 120 | //
|
|---|
| 121 | // Evaluation of the intermediates proceeds as follows:
|
|---|
| 122 | //
|
|---|
| 123 | // loop over batches of kl, k >= l, 0<=k,l<nocc_act
|
|---|
| 124 | // load (kl|xy), (kl| [T1,r12] |xy), and (lk| [T1,r12] |xy)
|
|---|
| 125 | // (aka kl-sets) into memory
|
|---|
| 126 | //
|
|---|
| 127 | // loop over batches of ij, i>=j, 0<=i,j<nocc_act
|
|---|
| 128 | // load (ij|r12|xy) into memory
|
|---|
| 129 | // (aka ij-sets) into memory
|
|---|
| 130 | // compute V[ij][kl] and T[ij][kl] for all ij and kl in
|
|---|
| 131 | // the "direct product" batch
|
|---|
| 132 | // end ij loop
|
|---|
| 133 | // end kl loop
|
|---|
| 134 | //
|
|---|
| 135 | /////////////////////////////////////////////////////////////////////////////////
|
|---|
| 136 |
|
|---|
| 137 | if (ijpq_acc->has_access(me)) {
|
|---|
| 138 |
|
|---|
| 139 | for(kl_iter.start();int(kl_iter);kl_iter.next()) {
|
|---|
| 140 |
|
|---|
| 141 | const int kl = kl_iter.ij();
|
|---|
| 142 | // Figure out if this task will handle this kl
|
|---|
| 143 | int kl_proc = kl%nproc_with_ints;
|
|---|
| 144 | if (kl_proc != proc_with_ints[me])
|
|---|
| 145 | continue;
|
|---|
| 146 | const int k = kl_iter.i();
|
|---|
| 147 | const int l = kl_iter.j();
|
|---|
| 148 | const int kl_aa = kl_iter.ij_aa();
|
|---|
| 149 | const int kl_ab = kl_iter.ij_ab();
|
|---|
| 150 | const int lk_ab = kl_iter.ij_ba();
|
|---|
| 151 |
|
|---|
| 152 | if (debug_)
|
|---|
| 153 | ExEnv::outn() << indent << "task " << me << ": working on (k,l) = " << k << "," << l << " " << endl;
|
|---|
| 154 |
|
|---|
| 155 | // Get (|1/r12|), (|r12|), and (|[r12,T1]|) integrals
|
|---|
| 156 | tim_enter("MO ints retrieve");
|
|---|
| 157 | double *klxy_buf_eri = ijpq_acc->retrieve_pair_block(k,l,R12IntsAcc::eri);
|
|---|
| 158 | double *klxy_buf_r12 = ijpq_acc->retrieve_pair_block(k,l,R12IntsAcc::r12);
|
|---|
| 159 | double *klxy_buf_r12t1 = ijpq_acc->retrieve_pair_block(k,l,R12IntsAcc::r12t1);
|
|---|
| 160 | double *lkxy_buf_r12t1 = ijpq_acc->retrieve_pair_block(l,k,R12IntsAcc::r12t1);
|
|---|
| 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 *ijxy_buf_r12 = ijpq_acc->retrieve_pair_block(i,j,R12IntsAcc::r12);
|
|---|
| 180 | tim_exit("MO ints retrieve");
|
|---|
| 181 |
|
|---|
| 182 | if (debug_)
|
|---|
| 183 | ExEnv::outn() << indent << "task " << me << ": obtained ij blocks" << endl;
|
|---|
| 184 |
|
|---|
| 185 |
|
|---|
| 186 | tim_enter("MO ints contraction");
|
|---|
| 187 | double Vaa_ijkl, Vab_ijkl, Vab_jikl, Vab_ijlk, Vab_jilk;
|
|---|
| 188 | double Xaa_ijkl, Xab_ijkl, Xab_jikl, Xab_ijlk, Xab_jilk;
|
|---|
| 189 | double Taa_ijkl, Tab_ijkl, Tab_jikl, Tab_ijlk, Tab_jilk;
|
|---|
| 190 | Vaa_ijkl = Vab_ijkl = Vab_jikl = Vab_ijlk = Vab_jilk = 0.0;
|
|---|
| 191 | Xaa_ijkl = Xab_ijkl = Xab_jikl = Xab_ijlk = Xab_jilk = 0.0;
|
|---|
| 192 | Taa_ijkl = Tab_ijkl = Tab_jikl = Tab_ijlk = Tab_jilk = 0.0;
|
|---|
| 193 |
|
|---|
| 194 | const double pfac_xy = 0.5;
|
|---|
| 195 | for(int y=0;y<noso;y++) {
|
|---|
| 196 | for(int x=0;x<noso;x++) {
|
|---|
| 197 | int yx_offset = y*noso+x;
|
|---|
| 198 | int xy_offset = x*noso+y;
|
|---|
| 199 | double ij_r12_xy = ijxy_buf_r12[xy_offset];
|
|---|
| 200 | double ij_r12_yx = ijxy_buf_r12[yx_offset];
|
|---|
| 201 | double kl_eri_xy = klxy_buf_eri[xy_offset];
|
|---|
| 202 | double kl_eri_yx = klxy_buf_eri[yx_offset];
|
|---|
| 203 | Vab_ijkl -= pfac_xy * (ij_r12_xy * kl_eri_xy + ij_r12_yx * kl_eri_yx);
|
|---|
| 204 | if (ij_ab != ji_ab)
|
|---|
| 205 | Vab_jikl -= pfac_xy * (ij_r12_yx * kl_eri_xy + ij_r12_xy * kl_eri_yx);
|
|---|
| 206 | if (kl_ab != lk_ab)
|
|---|
| 207 | Vab_ijlk -= pfac_xy * (ij_r12_xy * kl_eri_yx + ij_r12_yx * kl_eri_xy);
|
|---|
| 208 | if (ij_ab != ji_ab && kl_ab != lk_ab) {
|
|---|
| 209 | Vab_jilk -= pfac_xy * (ij_r12_yx * kl_eri_yx + ij_r12_xy * kl_eri_xy);
|
|---|
| 210 | }
|
|---|
| 211 | if (ij_aa != -1 && kl_aa != -1) {
|
|---|
| 212 | Vaa_ijkl -= pfac_xy * (ij_r12_xy - ij_r12_yx)*(kl_eri_xy - kl_eri_yx);
|
|---|
| 213 | }
|
|---|
| 214 | double kl_r12_xy = klxy_buf_r12[xy_offset];
|
|---|
| 215 | double kl_r12_yx = klxy_buf_r12[yx_offset];
|
|---|
| 216 | Xab_ijkl -= pfac_xy * (ij_r12_xy * kl_r12_xy + ij_r12_yx * kl_r12_yx);
|
|---|
| 217 | if (ij_ab != ji_ab)
|
|---|
| 218 | Xab_jikl -= pfac_xy * (ij_r12_yx * kl_r12_xy + ij_r12_xy * kl_r12_yx);
|
|---|
| 219 | if (kl_ab != lk_ab)
|
|---|
| 220 | Xab_ijlk -= pfac_xy * (ij_r12_xy * kl_r12_yx + ij_r12_yx * kl_r12_xy);
|
|---|
| 221 | if (ij_ab != ji_ab && kl_ab != lk_ab) {
|
|---|
| 222 | Xab_jilk -= pfac_xy * (ij_r12_yx * kl_r12_yx + ij_r12_xy * kl_r12_xy);
|
|---|
| 223 | }
|
|---|
| 224 | if (ij_aa != -1 && kl_aa != -1) {
|
|---|
| 225 | Xaa_ijkl -= pfac_xy * (ij_r12_xy - ij_r12_yx)*(kl_r12_xy - kl_r12_yx);
|
|---|
| 226 | }
|
|---|
| 227 | double kl_r12t1_xy = klxy_buf_r12t1[xy_offset];
|
|---|
| 228 | double kl_r12t1_yx = klxy_buf_r12t1[yx_offset];
|
|---|
| 229 | double lk_r12t1_xy = lkxy_buf_r12t1[xy_offset];
|
|---|
| 230 | double lk_r12t1_yx = lkxy_buf_r12t1[yx_offset];
|
|---|
| 231 | double kl_Tr12_xy = -kl_r12t1_xy-lk_r12t1_yx;
|
|---|
| 232 | double kl_Tr12_yx = -kl_r12t1_yx-lk_r12t1_xy;
|
|---|
| 233 | Tab_ijkl += pfac_xy * (ij_r12_xy * kl_Tr12_xy + ij_r12_yx * kl_Tr12_yx);
|
|---|
| 234 | if (ij_ab != ji_ab)
|
|---|
| 235 | Tab_jikl += pfac_xy * (ij_r12_yx * kl_Tr12_xy + ij_r12_xy * kl_Tr12_yx);
|
|---|
| 236 | if (kl_ab != lk_ab)
|
|---|
| 237 | Tab_ijlk += pfac_xy * (ij_r12_xy * kl_Tr12_yx + ij_r12_yx * kl_Tr12_xy);
|
|---|
| 238 | if (ij_ab != ji_ab && kl_ab != lk_ab) {
|
|---|
| 239 | Tab_jilk += pfac_xy * (ij_r12_yx * kl_Tr12_yx + ij_r12_xy * kl_Tr12_xy);
|
|---|
| 240 | }
|
|---|
| 241 | if (ij_aa != -1 && kl_aa != -1) {
|
|---|
| 242 | Taa_ijkl += pfac_xy * (ij_r12_xy - ij_r12_yx)*(kl_Tr12_xy - kl_Tr12_yx);
|
|---|
| 243 | }
|
|---|
| 244 | }
|
|---|
| 245 | }
|
|---|
| 246 | Vab_.accumulate_element(ij_ab,kl_ab,Vab_ijkl);
|
|---|
| 247 | if (ij_ab != ji_ab)
|
|---|
| 248 | Vab_.accumulate_element(ji_ab,kl_ab,Vab_jikl);
|
|---|
| 249 | if (kl_ab != lk_ab)
|
|---|
| 250 | Vab_.accumulate_element(ij_ab,lk_ab,Vab_ijlk);
|
|---|
| 251 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 252 | Vab_.accumulate_element(ji_ab,lk_ab,Vab_jilk);
|
|---|
| 253 | if (ij_aa != -1 && kl_aa != -1)
|
|---|
| 254 | Vaa_.accumulate_element(ij_aa,kl_aa,Vaa_ijkl);
|
|---|
| 255 | Xab_.accumulate_element(ij_ab,kl_ab,Xab_ijkl);
|
|---|
| 256 | if (ij_ab != ji_ab)
|
|---|
| 257 | Xab_.accumulate_element(ji_ab,kl_ab,Xab_jikl);
|
|---|
| 258 | if (kl_ab != lk_ab)
|
|---|
| 259 | Xab_.accumulate_element(ij_ab,lk_ab,Xab_ijlk);
|
|---|
| 260 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 261 | Xab_.accumulate_element(ji_ab,lk_ab,Xab_jilk);
|
|---|
| 262 | if (ij_aa != -1 && kl_aa != -1)
|
|---|
| 263 | Xaa_.accumulate_element(ij_aa,kl_aa,Xaa_ijkl);
|
|---|
| 264 | Bab_.accumulate_element(ij_ab,kl_ab,Tab_ijkl);
|
|---|
| 265 | if (ij_ab != ji_ab)
|
|---|
| 266 | Bab_.accumulate_element(ji_ab,kl_ab,Tab_jikl);
|
|---|
| 267 | if (kl_ab != lk_ab)
|
|---|
| 268 | Bab_.accumulate_element(ij_ab,lk_ab,Tab_ijlk);
|
|---|
| 269 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 270 | Bab_.accumulate_element(ji_ab,lk_ab,Tab_jilk);
|
|---|
| 271 | if (ij_aa != -1 && kl_aa != -1)
|
|---|
| 272 | Baa_.accumulate_element(ij_aa,kl_aa,Taa_ijkl);
|
|---|
| 273 | tim_exit("MO ints contraction");
|
|---|
| 274 |
|
|---|
| 275 | #if PRINT_R12_INTERMED
|
|---|
| 276 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 277 | printf("Vaa[%d][%d] = %lf\n",ij_aa,kl_aa,Vaa_ij[kl_aa]);
|
|---|
| 278 | printf("Vab[%d][%d] = %lf\n",ij_ab,kl_ab,Vab_ij[kl_ab]);
|
|---|
| 279 | if (ij_ab != ji_ab)
|
|---|
| 280 | printf("Vab[%d][%d] = %lf\n",ji_ab,kl_ab,Vab_ji[kl_ab]);
|
|---|
| 281 | if (kl_ab != lk_ab)
|
|---|
| 282 | printf("Vab[%d][%d] = %lf\n",ij_ab,lk_ab,Vab_ij[lk_ab]);
|
|---|
| 283 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 284 | printf("Vab[%d][%d] = %lf\n",ji_ab,lk_ab,Vab_ji[lk_ab]);
|
|---|
| 285 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 286 | printf("Xaa[%d][%d] = %lf\n",ij_aa,kl_aa,Xaa_ij[kl_aa]);
|
|---|
| 287 | printf("Xab[%d][%d] = %lf\n",ij_ab,kl_ab,Xab_ij[kl_ab]);
|
|---|
| 288 | if (ij_ab != ji_ab)
|
|---|
| 289 | printf("Xab[%d][%d] = %lf\n",ji_ab,kl_ab,Xab_ji[kl_ab]);
|
|---|
| 290 | if (kl_ab != lk_ab)
|
|---|
| 291 | printf("Xab[%d][%d] = %lf\n",ij_ab,lk_ab,Xab_ij[lk_ab]);
|
|---|
| 292 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 293 | printf("Xab[%d][%d] = %lf\n",ji_ab,lk_ab,Xab_ji[lk_ab]);
|
|---|
| 294 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 295 | printf("Taa[%d][%d] = %lf\n",ij_aa,kl_aa,Taa_ij[kl_aa]);
|
|---|
| 296 | printf("Tab[%d][%d] = %lf\n",ij_ab,kl_ab,Tab_ij[kl_ab]);
|
|---|
| 297 | if (ij_ab != ji_ab)
|
|---|
| 298 | printf("Tab[%d][%d] = %lf\n",ji_ab,kl_ab,Tab_ji[kl_ab]);
|
|---|
| 299 | if (kl_ab != lk_ab)
|
|---|
| 300 | printf("Tab[%d][%d] = %lf\n",ij_ab,lk_ab,Tab_ij[lk_ab]);
|
|---|
| 301 | if (ij_ab != ji_ab && kl_ab != lk_ab)
|
|---|
| 302 | printf("Tab[%d][%d] = %lf\n",ji_ab,lk_ab,Tab_ji[lk_ab]);
|
|---|
| 303 | #endif
|
|---|
| 304 | ijpq_acc->release_pair_block(i,j,R12IntsAcc::r12);
|
|---|
| 305 | }
|
|---|
| 306 | ijpq_acc->release_pair_block(k,l,R12IntsAcc::eri);
|
|---|
| 307 | ijpq_acc->release_pair_block(k,l,R12IntsAcc::r12);
|
|---|
| 308 | ijpq_acc->release_pair_block(k,l,R12IntsAcc::r12t1);
|
|---|
| 309 | ijpq_acc->release_pair_block(l,k,R12IntsAcc::r12t1);
|
|---|
| 310 | }
|
|---|
| 311 | }
|
|---|
| 312 | // Tasks that don't do any work here still need to create these timers
|
|---|
| 313 | tim_enter("MO ints retrieve");
|
|---|
| 314 | tim_exit("MO ints retrieve");
|
|---|
| 315 | tim_enter("MO ints contraction");
|
|---|
| 316 | tim_exit("MO ints contraction");
|
|---|
| 317 |
|
|---|
| 318 | tim_exit("intermediates");
|
|---|
| 319 | ExEnv::out0() << indent << "End of computation of intermediates" << endl;
|
|---|
| 320 | ijpq_acc->deactivate();
|
|---|
| 321 |
|
|---|
| 322 | // Symmetrize B intermediate
|
|---|
| 323 | for(int ij=0;ij<naa;ij++)
|
|---|
| 324 | for(int kl=0;kl<=ij;kl++) {
|
|---|
| 325 | double belem = 0.5*(Baa_->get_element(ij,kl) + Baa_->get_element(kl,ij));
|
|---|
| 326 | Baa_->set_element(ij,kl,belem);
|
|---|
| 327 | Baa_->set_element(kl,ij,belem);
|
|---|
| 328 | }
|
|---|
| 329 |
|
|---|
| 330 | for(int ij=0;ij<nab;ij++)
|
|---|
| 331 | for(int kl=0;kl<=ij;kl++) {
|
|---|
| 332 | double belem = 0.5*(Bab_->get_element(ij,kl) + Bab_->get_element(kl,ij));
|
|---|
| 333 | Bab_->set_element(ij,kl,belem);
|
|---|
| 334 | Bab_->set_element(kl,ij,belem);
|
|---|
| 335 | }
|
|---|
| 336 |
|
|---|
| 337 | globally_sum_intermeds_();
|
|---|
| 338 |
|
|---|
| 339 | ExEnv::out0() << decindent;
|
|---|
| 340 | ExEnv::out0() << indent
|
|---|
| 341 | << "Exited \"" << mospace->name()
|
|---|
| 342 | << "\" A (GEBC) intermediates evaluator" << endl;
|
|---|
| 343 |
|
|---|
| 344 | tim_exit("mp2-r12a intermeds (symmetric term)");
|
|---|
| 345 | checkpoint_();
|
|---|
| 346 | }
|
|---|
| 347 |
|
|---|
| 348 | ////////////////////////////////////////////////////////////////////////////
|
|---|
| 349 |
|
|---|
| 350 | // Local Variables:
|
|---|
| 351 | // mode: c++
|
|---|
| 352 | // c-file-style: "CLJ-CONDENSED"
|
|---|
| 353 | // End:
|
|---|
