← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/mpqc/oneiric

« back to all changes in this revision

Viewing changes to src/lib/chemistry/qc/mbptr12/gbc_contribs.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2005-11-27 11:41:49 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20051127114149-zgz9r3gk50w8ww2q
Tags: 2.3.0-1
http://bugs.debian.org/339232
* New upstream release.
* debian/rules (SONAME): Activate awk snippet for automatic so-name
  detection again, resulting in a bump to `7' and making a `c2a' for
  the C++ allocator change unnecessary; closes: #339232.
* debian/patches/00list (08_gcc-4.0_fixes): Removed, no longer needed.
* debian/rules (test): Remove workarounds, do not abort build if tests
  fail.
* debian/ref: Removed.
* debian/control.in (libsc): Added Conflict against libsc6c2.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/chemistry/qc/mbptr12/gbc_contribs.cc
 
1
//
 
2
// gbc_contribs.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 <string.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/mbptr12/blas.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
void
 
54
R12IntEval::compute_B_gbc_1_()
 
55
{
 
56
  if (abs_method_ == LinearR12::ABS_ABS || abs_method_ == LinearR12::ABS_ABSPlus)
 
57
    throw std::runtime_error("R12IntEval::compute_B_gbc_1_() -- B(GBC1) term can only be computed using a CABS (or CABS+) approach");
 
58
 
 
59
  if (evaluated_)
 
60
    return;
 
61
 
 
62
  Ref<TwoBodyMOIntsTransform> ipjq_tform = get_tform_("(ip|jq)");
 
63
  Ref<R12IntsAcc> ijpq_acc = ipjq_tform->ints_acc();
 
64
  if (!ijpq_acc->is_committed())
 
65
    ipjq_tform->compute();
 
66
  if (!ijpq_acc->is_active())
 
67
    ijpq_acc->activate();
 
68
 
 
69
  tim_enter("B(GBC1) intermediate");
 
70
 
 
71
  const int num_te_types = 2;
 
72
  
 
73
  Ref<MessageGrp> msg = r12info()->msg();
 
74
  int me = msg->me();
 
75
  int nproc = msg->n();
 
76
  ExEnv::out0() << endl << indent
 
77
    << "Entered B(GBC1) intermediate evaluator" << endl;
 
78
  ExEnv::out0() << incindent;
 
79
 
 
80
  RefSCMatrix B_gbc1_aa = Baa_.clone();  B_gbc1_aa.assign(0.0);
 
81
  RefSCMatrix B_gbc1_ab = Bab_.clone();  B_gbc1_ab.assign(0.0);
 
82
 
 
83
  Ref<MOIndexSpace> mo_space = r12info_->obs_space();
 
84
  Ref<MOIndexSpace> occ_space = r12info_->occ_space();
 
85
  Ref<MOIndexSpace> act_occ_space = r12info_->act_occ_space();
 
86
  Ref<MOIndexSpace> vir_space = r12info_->vir_space();
 
87
  Ref<MOIndexSpace> ribs_space = r12info_->ribs_space();
 
88
  form_focc_space_();
 
89
  Ref<MOIndexSpace> focc_space = focc_space_;
 
90
 
 
91
  const int noso = r12info_->mo_space()->rank();
 
92
  const int nocc = r12info_->nocc();
 
93
  const int nvir = noso - nocc;
 
94
  const int nribs = ribs_space->rank();
 
95
 
 
96
  //
 
97
  // Do the AO->MO transform for (act_occ occ|r12|act_occ ribs) and (act_occ focc|r12|act_occ ribs)
 
98
  //
 
99
  Ref<MOIntsTransformFactory> tfactory = r12info_->tfactory();
 
100
 
 
101
  tfactory->set_spaces(act_occ_space,occ_space,
 
102
                       act_occ_space,ribs_space);
 
103
  Ref<TwoBodyMOIntsTransform> imjA_tform = tfactory->twobody_transform_13("(im|jA)");
 
104
  imjA_tform->set_num_te_types(num_te_types);
 
105
  imjA_tform->compute();
 
106
  Ref<R12IntsAcc> ijmA_acc = imjA_tform->ints_acc();
 
107
 
 
108
  tfactory->set_spaces(act_occ_space,focc_space,
 
109
                       act_occ_space,ribs_space);
 
110
  Ref<TwoBodyMOIntsTransform> iMfjA_tform = tfactory->twobody_transform_13("(iMf|jA)");
 
111
  iMfjA_tform->set_num_te_types(num_te_types);
 
112
  iMfjA_tform->compute();
 
113
  Ref<R12IntsAcc> ijMfA_acc = iMfjA_tform->ints_acc();
 
114
  
 
115
  SpatialMOPairIter_eq ij_iter(act_occ_space);
 
116
  SpatialMOPairIter_eq kl_iter(act_occ_space);
 
117
  int naa = ij_iter.nij_aa();          // Number of alpha-alpha pairs (i > j)
 
118
  int nab = ij_iter.nij_ab();          // Number of alpha-beta pairs
 
119
  if (debug_) {
 
120
    ExEnv::out0() << indent << "naa = " << naa << endl;
 
121
    ExEnv::out0() << indent << "nab = " << nab << endl;
 
122
  }
 
123
 
 
124
  // Compute the number of tasks that have full access to the integrals
 
125
  // and split the work among them
 
126
  vector<int> proc_with_ints;
 
127
  int nproc_with_ints = tasks_with_ints_(ijMfA_acc,proc_with_ints);
 
128
 
 
129
  // Compute the first half of the term
 
130
  const int nbraket = nocc*nribs;
 
131
 
 
132
#if 1
 
133
  for(kl_iter.start();int(kl_iter);kl_iter.next()) {
 
134
 
 
135
    const int kl = kl_iter.ij();
 
136
    // Figure out if this task will handle this kl
 
137
    int kl_proc = kl%nproc_with_ints;
 
138
    if (kl_proc != proc_with_ints[me])
 
139
      continue;
 
140
    const int k = kl_iter.i();
 
141
    const int l = kl_iter.j();
 
142
    const int kl_aa = kl_iter.ij_aa();
 
143
    const int kl_ab = kl_iter.ij_ab();
 
144
    const int lk_ab = kl_iter.ij_ba();
 
145
 
 
146
    if (debug_)
 
147
      ExEnv::outn() << indent << "task " << me << ": working on (k,l) = " << k << "," << l << " " << endl;
 
148
 
 
149
    // Get (|r12|) integrals
 
150
    tim_enter("MO ints retrieve");
 
151
    double *klMfA_buf_r12 = ijMfA_acc->retrieve_pair_block(k,l,R12IntsAcc::r12);
 
152
    double *lkMfA_buf_r12 = ijMfA_acc->retrieve_pair_block(l,k,R12IntsAcc::r12);
 
153
    tim_exit("MO ints retrieve");
 
154
 
 
155
    if (debug_)
 
156
      ExEnv::outn() << indent << "task " << me << ": obtained kl blocks" << endl;
 
157
 
 
158
    for(ij_iter.start(kl+1);int(ij_iter);ij_iter.next()) {
 
159
 
 
160
      const int ij = ij_iter.ij();
 
161
      const int i = ij_iter.i();
 
162
      const int j = ij_iter.j();
 
163
      const int ij_aa = ij_iter.ij_aa();
 
164
      const int ij_ab = ij_iter.ij_ab();
 
165
      const int ji_ab = ij_iter.ij_ba();
 
166
 
 
167
      if (debug_)
 
168
        ExEnv::outn() << indent << "task " << me << ": working on (i,j) = " << i << "," << j << " " << endl;
 
169
 
 
170
      // Get (|r12|) integrals
 
171
      tim_enter("MO ints retrieve");
 
172
      double *ijmA_buf_r12 = ijmA_acc->retrieve_pair_block(i,j,R12IntsAcc::r12);
 
173
      double *jimA_buf_r12 = ijmA_acc->retrieve_pair_block(j,i,R12IntsAcc::r12);
 
174
      tim_exit("MO ints retrieve");
 
175
 
 
176
      if (debug_)
 
177
        ExEnv::outn() << indent << "task " << me << ": obtained ij blocks" << endl;
 
178
 
 
179
      double rr_klij = 0.0;
 
180
      double rr_lkji = 0.0;
 
181
      double rr_klji = 0.0;
 
182
      double rr_lkij = 0.0;
 
183
 
 
184
      const int unit_stride = 1;
 
185
      rr_klij = F77_DDOT(&nbraket,klMfA_buf_r12,&unit_stride,ijmA_buf_r12,&unit_stride);
 
186
      if (kl_ab != lk_ab && ij_ab != ji_ab) {
 
187
        rr_lkji = F77_DDOT(&nbraket,lkMfA_buf_r12,&unit_stride,jimA_buf_r12,&unit_stride);
 
188
      }
 
189
      else
 
190
        rr_lkji = rr_klij;
 
191
      B_gbc1_ab.set_element(kl_ab,ij_ab,-(rr_klij+rr_lkji));
 
192
      B_gbc1_ab.set_element(lk_ab,ji_ab,-(rr_klij+rr_lkji));
 
193
      
 
194
      if (kl_ab != lk_ab)
 
195
        rr_lkij = F77_DDOT(&nbraket,lkMfA_buf_r12,&unit_stride,ijmA_buf_r12,&unit_stride);
 
196
      else
 
197
        rr_lkij = rr_klij;
 
198
      if (ij_ab != ji_ab)
 
199
        rr_klji += F77_DDOT(&nbraket,klMfA_buf_r12,&unit_stride,jimA_buf_r12,&unit_stride);
 
200
      else
 
201
        rr_klji = rr_klij;
 
202
      B_gbc1_ab.set_element(kl_ab,ji_ab,-(rr_klji+rr_lkij));
 
203
      B_gbc1_ab.set_element(lk_ab,ij_ab,-(rr_klji+rr_lkij));
 
204
 
 
205
      if (ij_aa != -1 && kl_aa != -1)
 
206
        B_gbc1_aa.set_element(kl_aa,ij_aa,-(rr_klij+rr_lkji-rr_klji-rr_lkij));
 
207
        
 
208
      ijmA_acc->release_pair_block(i,j,R12IntsAcc::r12);
 
209
      ijmA_acc->release_pair_block(j,i,R12IntsAcc::r12);
 
210
    }
 
211
 
 
212
    ijMfA_acc->release_pair_block(k,l,R12IntsAcc::r12);
 
213
    ijMfA_acc->release_pair_block(l,k,R12IntsAcc::r12);
 
214
  }
 
215
#endif
 
216
 
 
217
#if 1
 
218
  //
 
219
  // Do the AO->MO transform for (act_occ focc|r12|act_occ vir)
 
220
  //
 
221
  tfactory->set_spaces(act_occ_space,focc_space,
 
222
                       act_occ_space,vir_space);
 
223
  Ref<TwoBodyMOIntsTransform> iMfja_tform = tfactory->twobody_transform_13("(iMf|ja)");
 
224
  iMfja_tform->set_num_te_types(num_te_types);
 
225
  iMfja_tform->compute();
 
226
  Ref<R12IntsAcc> ijMfa_acc = iMfja_tform->ints_acc();
 
227
 
 
228
  nproc_with_ints = tasks_with_ints_(ijMfa_acc,proc_with_ints);
 
229
 
 
230
  // Compute the second half of the term
 
231
  for(kl_iter.start();int(kl_iter);kl_iter.next()) {
 
232
 
 
233
    const int kl = kl_iter.ij();
 
234
    // Figure out if this task will handle this kl
 
235
    int kl_proc = kl%nproc_with_ints;
 
236
    if (kl_proc != proc_with_ints[me])
 
237
      continue;
 
238
    const int k = kl_iter.i();
 
239
    const int l = kl_iter.j();
 
240
    const int kl_aa = kl_iter.ij_aa();
 
241
    const int kl_ab = kl_iter.ij_ab();
 
242
    const int lk_ab = kl_iter.ij_ba();
 
243
 
 
244
    if (debug_)
 
245
      ExEnv::outn() << indent << "task " << me << ": working on (k,l) = " << k << "," << l << " " << endl;
 
246
 
 
247
    // Get (|r12|) integrals
 
248
    tim_enter("MO ints retrieve");
 
249
    double *klMfa_buf_r12 = ijMfa_acc->retrieve_pair_block(k,l,R12IntsAcc::r12);
 
250
    double *lkMfa_buf_r12 = ijMfa_acc->retrieve_pair_block(l,k,R12IntsAcc::r12);
 
251
    tim_exit("MO ints retrieve");
 
252
 
 
253
    if (debug_)
 
254
      ExEnv::outn() << indent << "task " << me << ": obtained kl blocks" << endl;
 
255
 
 
256
    for(ij_iter.start(kl+1);int(ij_iter);ij_iter.next()) {
 
257
 
 
258
      const int ij = ij_iter.ij();
 
259
      const int i = ij_iter.i();
 
260
      const int j = ij_iter.j();
 
261
      const int ij_aa = ij_iter.ij_aa();
 
262
      const int ij_ab = ij_iter.ij_ab();
 
263
      const int ji_ab = ij_iter.ij_ba();
 
264
 
 
265
      if (debug_)
 
266
        ExEnv::outn() << indent << "task " << me << ": working on (i,j) = " << i << "," << j << " " << endl;
 
267
 
 
268
      // Get (|r12|) integrals
 
269
      tim_enter("MO ints retrieve");
 
270
      double *ijpq_buf_r12 = ijpq_acc->retrieve_pair_block(i,j,R12IntsAcc::r12);
 
271
      tim_exit("MO ints retrieve");
 
272
 
 
273
      if (debug_)
 
274
        ExEnv::outn() << indent << "task " << me << ": obtained ij blocks" << endl;
 
275
 
 
276
      int ma = 0;
 
277
      for(int m=0; m<nocc; m++)
 
278
        for(int a=0; a<nvir; a++, ma++) {
 
279
          
 
280
          const int aa = a+nocc;
 
281
          const int ma_offset = m*noso+aa;
 
282
          const int am_offset = aa*noso+m;
 
283
 
 
284
          double rr_klij = 0.0;
 
285
          double rr_lkij = 0.0;
 
286
          double rr_klji = 0.0;
 
287
          double rr_lkji = 0.0;
 
288
          
 
289
          rr_klij = klMfa_buf_r12[ma]*ijpq_buf_r12[ma_offset];
 
290
          rr_lkji = lkMfa_buf_r12[ma]*ijpq_buf_r12[am_offset];
 
291
          B_gbc1_ab.accumulate_element(kl_ab,ij_ab,-(rr_klij+rr_lkji));
 
292
          if (kl_ab != lk_ab && ij_ab != ji_ab) {
 
293
            B_gbc1_ab.accumulate_element(lk_ab,ji_ab,-(rr_lkji+rr_klij));
 
294
          }
 
295
          
 
296
          rr_lkij = lkMfa_buf_r12[ma]*ijpq_buf_r12[ma_offset];
 
297
          rr_klji = klMfa_buf_r12[ma]*ijpq_buf_r12[am_offset];
 
298
          if (kl_ab != lk_ab)
 
299
            B_gbc1_ab.accumulate_element(lk_ab,ij_ab,-(rr_lkij+rr_klji));
 
300
          if (ij_ab != ji_ab)
 
301
            B_gbc1_ab.accumulate_element(kl_ab,ji_ab,-(rr_klji+rr_lkij));
 
302
          
 
303
          if (ij_aa != -1 && kl_aa != -1)
 
304
            B_gbc1_aa.accumulate_element(kl_aa,ij_aa,-(rr_klij-rr_lkij-rr_klji+rr_lkji));
 
305
        }
 
306
      
 
307
      ijpq_acc->release_pair_block(i,j,R12IntsAcc::r12);
 
308
    }
 
309
 
 
310
    ijMfa_acc->release_pair_block(k,l,R12IntsAcc::r12);
 
311
    ijMfa_acc->release_pair_block(l,k,R12IntsAcc::r12);
 
312
  }
 
313
#endif
 
314
 
 
315
  if (debug_ > 1) {
 
316
    B_gbc1_aa.print("Alpha-alpha B(GBC1) contribution");
 
317
    B_gbc1_ab.print("Alpha-beta B(GBC1) contribution");
 
318
  }
 
319
  // Symmetrize the B contribution
 
320
  B_gbc1_aa.scale(0.5);
 
321
  B_gbc1_ab.scale(0.5);
 
322
  RefSCMatrix B_gbc1_aa_t = B_gbc1_aa.t();
 
323
  Baa_.accumulate(B_gbc1_aa); Baa_.accumulate(B_gbc1_aa_t);
 
324
  RefSCMatrix B_gbc1_ab_t = B_gbc1_ab.t();
 
325
  Bab_.accumulate(B_gbc1_ab); Bab_.accumulate(B_gbc1_ab_t);
 
326
 
 
327
  globally_sum_intermeds_();
 
328
 
 
329
  ExEnv::out0() << decindent;
 
330
  ExEnv::out0() << indent << "Exited B(GBC1) intermediate evaluator" << endl;
 
331
 
 
332
  tim_exit("B(GBC1) intermediate");
 
333
}
 
334
 
 
335
 
 
336
void
 
337
R12IntEval::compute_B_gbc_2_()
 
338
{
 
339
  if (abs_method_ == LinearR12::ABS_ABS || abs_method_ == LinearR12::ABS_ABSPlus)
 
340
    throw std::runtime_error("R12IntEval::compute_B_gbc_2_() -- B(GBC2) term can only be computed using a CABS (or CABS+) approach");
 
341
 
 
342
  if (evaluated_)
 
343
    return;
 
344
 
 
345
  Ref<TwoBodyMOIntsTransform> ipjq_tform = get_tform_("(ip|jq)");
 
346
  Ref<R12IntsAcc> ijpq_acc = ipjq_tform->ints_acc();
 
347
  if (!ijpq_acc->is_committed())
 
348
    ipjq_tform->compute();
 
349
  if (!ijpq_acc->is_active())
 
350
    ijpq_acc->activate();
 
351
 
 
352
  tim_enter("B(GBC2) intermediate");
 
353
 
 
354
  const int num_te_types = 2; // only integrals of r_{12} are needed
 
355
  Ref<MessageGrp> msg = r12info_->msg();
 
356
  int me = msg->me();
 
357
  int nproc = msg->n();
 
358
  ExEnv::out0() << endl << indent
 
359
    << "Entered B(GBC2) intermediate evaluator" << endl;
 
360
  ExEnv::out0() << incindent;
 
361
 
 
362
  RefSCMatrix X_ijklF_ab = Bab_.clone();
 
363
  RefSCMatrix B_gbc2_aa = Baa_.clone();
 
364
  RefSCMatrix B_gbc2_ab = Bab_.clone();
 
365
  X_ijklF_ab.assign(0.0);
 
366
  B_gbc2_aa.assign(0.0);
 
367
  B_gbc2_ab.assign(0.0);
 
368
 
 
369
  Ref<MOIndexSpace> obs_space = r12info_->obs_space();
 
370
  Ref<MOIndexSpace> occ_space = r12info_->occ_space();
 
371
  Ref<MOIndexSpace> act_occ_space = r12info_->act_occ_space();
 
372
  Ref<MOIndexSpace> ribs_space = r12info_->ribs_space();
 
373
  form_factocc_space_();
 
374
  Ref<MOIndexSpace> factocc_space = factocc_space_;
 
375
 
 
376
  const int nocc = r12info_->nocc();
 
377
  const int nribs = ribs_space->rank();
 
378
 
 
379
  // compute r_{12}^2 operator in act.occ.pair/act.occ.-focc. basis
 
380
  RefSCMatrix R2 = compute_r2_(act_occ_space,factocc_space);
 
381
#if 1
 
382
  // Compute contribution X += (r^2)_{ij}^{k l_f}
 
383
  if (me == 0)
 
384
    X_ijklF_ab.accumulate(R2);
 
385
#endif
 
386
 
 
387
  //
 
388
  // Compute contribution X -= r_{ij}^{\alpha'm} r_{m\alpha'}^{k l_f}
 
389
  //                         + r_{ji}^{\alpha'm} r_{\alpha'm}^{k l_f}
 
390
  //
 
391
  Ref<MOIntsTransformFactory> tfactory = r12info_->tfactory();
 
392
 
 
393
  tfactory->set_spaces(act_occ_space,occ_space,
 
394
                       act_occ_space,ribs_space);
 
395
  Ref<TwoBodyMOIntsTransform> imjA_tform = tfactory->twobody_transform_13("(im|jA)");
 
396
  imjA_tform->set_num_te_types(num_te_types);
 
397
  imjA_tform->compute();
 
398
  Ref<R12IntsAcc> ijmA_acc = imjA_tform->ints_acc();
 
399
 
 
400
  tfactory->set_spaces(act_occ_space,occ_space,
 
401
                       factocc_space,ribs_space);
 
402
  Ref<TwoBodyMOIntsTransform> kmlfA_tform = tfactory->twobody_transform_13("(km|lfA)");
 
403
  kmlfA_tform->set_num_te_types(num_te_types);
 
404
  kmlfA_tform->compute();
 
405
  Ref<R12IntsAcc> klfmA_acc = kmlfA_tform->ints_acc();
 
406
 
 
407
  tfactory->set_spaces(factocc_space,occ_space,
 
408
                       act_occ_space,ribs_space);
 
409
  Ref<TwoBodyMOIntsTransform> lfmkA_tform = tfactory->twobody_transform_13("(lfm|kA)");
 
410
  lfmkA_tform->set_num_te_types(num_te_types);
 
411
  lfmkA_tform->compute();
 
412
  Ref<R12IntsAcc> lfkmA_acc = lfmkA_tform->ints_acc();
 
413
 
 
414
  // Compute the number of tasks that have full access to the integrals
 
415
  // and split the work among them
 
416
  vector<int> proc_with_ints;
 
417
  int nproc_with_ints = tasks_with_ints_(ijmA_acc,proc_with_ints);
 
418
 
 
419
  SpatialMOPairIter_eq ij_iter(act_occ_space);
 
420
  SpatialMOPairIter_eq kl_iter(act_occ_space);
 
421
 
 
422
  int nbraket = nocc*nribs;
 
423
#if 1
 
424
  for(kl_iter.start();int(kl_iter);kl_iter.next()) {
 
425
 
 
426
    const int kl = kl_iter.ij();
 
427
    // Figure out if this task will handle this kl
 
428
    int kl_proc = kl%nproc_with_ints;
 
429
    if (kl_proc != proc_with_ints[me])
 
430
      continue;
 
431
    const int k = kl_iter.i();
 
432
    const int l = kl_iter.j();
 
433
    const int kl_ab = kl_iter.ij_ab();
 
434
    const int lk_ab = kl_iter.ij_ba();
 
435
 
 
436
    if (debug_)
 
437
      ExEnv::outn() << indent << "task " << me << ": working on (k,l) = " << k << "," << l << " " << endl;
 
438
 
 
439
    // Get (|r12|) integrals
 
440
    tim_enter("MO ints retrieve");
 
441
    double *klfmA_buf_r12 = klfmA_acc->retrieve_pair_block(k,l,R12IntsAcc::r12);
 
442
    double *lfkmA_buf_r12 = lfkmA_acc->retrieve_pair_block(l,k,R12IntsAcc::r12);
 
443
    double *lkfmA_buf_r12 = klfmA_acc->retrieve_pair_block(l,k,R12IntsAcc::r12);
 
444
    double *kflmA_buf_r12 = lfkmA_acc->retrieve_pair_block(k,l,R12IntsAcc::r12);
 
445
    tim_exit("MO ints retrieve");
 
446
 
 
447
    if (debug_)
 
448
      ExEnv::outn() << indent << "task " << me << ": obtained kl blocks" << endl;
 
449
 
 
450
    for(ij_iter.start();int(ij_iter);ij_iter.next()) {
 
451
 
 
452
      const int ij = ij_iter.ij();
 
453
      const int i = ij_iter.i();
 
454
      const int j = ij_iter.j();
 
455
      const int ij_ab = ij_iter.ij_ab();
 
456
      const int ji_ab = ij_iter.ij_ba();
 
457
 
 
458
      if (debug_)
 
459
        ExEnv::outn() << indent << "task " << me << ": working on (i,j) = " << i << "," << j << " " << endl;
 
460
 
 
461
      // Get (|r12|) integrals
 
462
      tim_enter("MO ints retrieve");
 
463
      double *ijmA_buf_r12 = ijmA_acc->retrieve_pair_block(i,j,R12IntsAcc::r12);
 
464
      double *jimA_buf_r12 = ijmA_acc->retrieve_pair_block(j,i,R12IntsAcc::r12);
 
465
      tim_exit("MO ints retrieve");
 
466
 
 
467
      if (debug_)
 
468
        ExEnv::outn() << indent << "task " << me << ": obtained ij blocks" << endl;
 
469
 
 
470
      double X_ijklf = 0.0;
 
471
      double X_jiklf = 0.0;
 
472
      double X_ijlkf = 0.0;
 
473
      double X_jilkf = 0.0;
 
474
 
 
475
      const int unit_stride = 1;
 
476
      X_ijklf += F77_DDOT(&nbraket,lfkmA_buf_r12,&unit_stride,jimA_buf_r12,&unit_stride);
 
477
      X_ijklf += F77_DDOT(&nbraket,klfmA_buf_r12,&unit_stride,ijmA_buf_r12,&unit_stride);
 
478
      X_ijklF_ab.accumulate_element(ij_ab,kl_ab,-X_ijklf);
 
479
      if (kl_ab != lk_ab) {
 
480
        X_ijlkf += F77_DDOT(&nbraket,kflmA_buf_r12,&unit_stride,jimA_buf_r12,&unit_stride);
 
481
        X_ijlkf += F77_DDOT(&nbraket,lkfmA_buf_r12,&unit_stride,ijmA_buf_r12,&unit_stride);
 
482
        X_ijklF_ab.accumulate_element(ij_ab,lk_ab,-X_ijlkf);
 
483
      }
 
484
      if (ij_ab != ji_ab) {
 
485
        X_jiklf += F77_DDOT(&nbraket,lfkmA_buf_r12,&unit_stride,ijmA_buf_r12,&unit_stride);
 
486
        X_jiklf += F77_DDOT(&nbraket,klfmA_buf_r12,&unit_stride,jimA_buf_r12,&unit_stride);
 
487
        X_ijklF_ab.accumulate_element(ji_ab,kl_ab,-X_jiklf);
 
488
        if (kl_ab != lk_ab) {
 
489
          X_jilkf += F77_DDOT(&nbraket,kflmA_buf_r12,&unit_stride,ijmA_buf_r12,&unit_stride);
 
490
          X_jilkf += F77_DDOT(&nbraket,lkfmA_buf_r12,&unit_stride,jimA_buf_r12,&unit_stride);
 
491
          X_ijklF_ab.accumulate_element(ji_ab,lk_ab,-X_jilkf);
 
492
        }
 
493
      }
 
494
      
 
495
      ijmA_acc->release_pair_block(i,j,R12IntsAcc::r12);
 
496
      ijmA_acc->release_pair_block(j,i,R12IntsAcc::r12);
 
497
    }
 
498
 
 
499
    klfmA_acc->release_pair_block(k,l,R12IntsAcc::r12);
 
500
    lfkmA_acc->release_pair_block(l,k,R12IntsAcc::r12);
 
501
    klfmA_acc->release_pair_block(l,k,R12IntsAcc::r12);
 
502
    lfkmA_acc->release_pair_block(k,l,R12IntsAcc::r12);
 
503
  }
 
504
#endif
 
505
 
 
506
#if 1
 
507
  //
 
508
  // Compute contribution X -= r_{ij}^{pq} r_{pq}^{k l_f}
 
509
  //
 
510
  tfactory->set_spaces(act_occ_space,obs_space,
 
511
                       factocc_space,obs_space);
 
512
  Ref<TwoBodyMOIntsTransform> kplfq_tform = tfactory->twobody_transform_13("(kp|lfq)");
 
513
  kplfq_tform->set_num_te_types(num_te_types);
 
514
  kplfq_tform->compute();
 
515
  Ref<R12IntsAcc> klfpq_acc = kplfq_tform->ints_acc();
 
516
 
 
517
  nbraket = obs_space->rank() * obs_space->rank();
 
518
  for(kl_iter.start();int(kl_iter);kl_iter.next()) {
 
519
 
 
520
    const int kl = kl_iter.ij();
 
521
    // Figure out if this task will handle this kl
 
522
    int kl_proc = kl%nproc_with_ints;
 
523
    if (kl_proc != proc_with_ints[me])
 
524
      continue;
 
525
    const int k = kl_iter.i();
 
526
    const int l = kl_iter.j();
 
527
    const int kl_ab = kl_iter.ij_ab();
 
528
    const int lk_ab = kl_iter.ij_ba();
 
529
 
 
530
    if (debug_)
 
531
      ExEnv::outn() << indent << "task " << me << ": working on (k,l) = " << k << "," << l << " " << endl;
 
532
 
 
533
    // Get (|r12|) integrals
 
534
    tim_enter("MO ints retrieve");
 
535
    double *klfpq_buf_r12 = klfpq_acc->retrieve_pair_block(k,l,R12IntsAcc::r12);
 
536
    double *lkfpq_buf_r12 = klfpq_acc->retrieve_pair_block(l,k,R12IntsAcc::r12);
 
537
    tim_exit("MO ints retrieve");
 
538
 
 
539
    if (debug_)
 
540
      ExEnv::outn() << indent << "task " << me << ": obtained kl blocks" << endl;
 
541
 
 
542
    for(ij_iter.start();int(ij_iter);ij_iter.next()) {
 
543
 
 
544
      const int ij = ij_iter.ij();
 
545
      const int i = ij_iter.i();
 
546
      const int j = ij_iter.j();
 
547
      const int ij_ab = ij_iter.ij_ab();
 
548
      const int ji_ab = ij_iter.ij_ba();
 
549
 
 
550
      if (debug_)
 
551
        ExEnv::outn() << indent << "task " << me << ": working on (i,j) = " << i << "," << j << " " << endl;
 
552
 
 
553
      // Get (|r12|) integrals
 
554
      tim_enter("MO ints retrieve");
 
555
      double *ijpq_buf_r12 = ijpq_acc->retrieve_pair_block(i,j,R12IntsAcc::r12);
 
556
      double *jipq_buf_r12 = ijpq_acc->retrieve_pair_block(j,i,R12IntsAcc::r12);
 
557
      tim_exit("MO ints retrieve");
 
558
 
 
559
      if (debug_)
 
560
        ExEnv::outn() << indent << "task " << me << ": obtained ij blocks" << endl;
 
561
 
 
562
      double X_ijklf = 0.0;
 
563
      double X_jiklf = 0.0;
 
564
      double X_ijlkf = 0.0;
 
565
      double X_jilkf = 0.0;
 
566
 
 
567
      const int unit_stride = 1;
 
568
      X_ijklf += F77_DDOT(&nbraket,klfpq_buf_r12,&unit_stride,ijpq_buf_r12,&unit_stride);
 
569
      X_ijklF_ab.accumulate_element(ij_ab,kl_ab,-X_ijklf);
 
570
      if (kl_ab != lk_ab) {
 
571
        X_ijlkf += F77_DDOT(&nbraket,lkfpq_buf_r12,&unit_stride,ijpq_buf_r12,&unit_stride);
 
572
        X_ijklF_ab.accumulate_element(ij_ab,lk_ab,-X_ijlkf);
 
573
      }
 
574
      if (ij_ab != ji_ab) {
 
575
        X_jiklf += F77_DDOT(&nbraket,klfpq_buf_r12,&unit_stride,jipq_buf_r12,&unit_stride);
 
576
        X_ijklF_ab.accumulate_element(ji_ab,kl_ab,-X_jiklf);
 
577
        if (kl_ab != lk_ab) {
 
578
          X_jilkf += F77_DDOT(&nbraket,lkfpq_buf_r12,&unit_stride,jipq_buf_r12,&unit_stride);
 
579
          X_ijklF_ab.accumulate_element(ji_ab,lk_ab,-X_jilkf);
 
580
        }
 
581
      }
 
582
      
 
583
      ijpq_acc->release_pair_block(i,j,R12IntsAcc::r12);
 
584
      ijpq_acc->release_pair_block(j,i,R12IntsAcc::r12);
 
585
    }
 
586
 
 
587
    klfpq_acc->release_pair_block(k,l,R12IntsAcc::r12);
 
588
    klfpq_acc->release_pair_block(l,k,R12IntsAcc::r12);
 
589
  }
 
590
  globally_sum_scmatrix_(X_ijklF_ab);
 
591
#endif
 
592
 
 
593
  //
 
594
  // Compute B_gbc2 = X_ijklF + X_jilkF :
 
595
  // B_gbc2_ab_ijkl = X_ijklF_ab + X_jilkF_ab
 
596
  // B_gbc2_aa_ijkl = X_ijklF_aa + X_jilkF_aa = X_ijklF_ab - X_jiklF_ab + X_jilkF_ab - X_ijlkF_ab
 
597
  //
 
598
 
 
599
  for(kl_iter.start();int(kl_iter);kl_iter.next()) {
 
600
 
 
601
    const int kl_aa = kl_iter.ij_aa();
 
602
    const int kl_ab = kl_iter.ij_ab();
 
603
    const int lk_ab = kl_iter.ij_ba();
 
604
 
 
605
    for(ij_iter.start();int(ij_iter);ij_iter.next()) {
 
606
 
 
607
      const int ij_aa = ij_iter.ij_aa();
 
608
      const int ij_ab = ij_iter.ij_ab();
 
609
      const int ji_ab = ij_iter.ij_ba();
 
610
 
 
611
      const double B_ab_ijkl = X_ijklF_ab.get_element(ij_ab,kl_ab) + X_ijklF_ab.get_element(ji_ab,lk_ab);
 
612
      const double B_ab_ijlk = X_ijklF_ab.get_element(ij_ab,lk_ab) + X_ijklF_ab.get_element(ji_ab,kl_ab);
 
613
      const double B_ab_jikl = B_ab_ijlk;
 
614
      const double B_ab_jilk = B_ab_ijkl;
 
615
 
 
616
      B_gbc2_ab.set_element( ij_ab, kl_ab, B_ab_ijkl);
 
617
      if (kl_ab != lk_ab)
 
618
        B_gbc2_ab.set_element( ij_ab, lk_ab, B_ab_ijlk);
 
619
      if (ij_ab != ji_ab) {
 
620
        B_gbc2_ab.set_element( ji_ab, kl_ab, B_ab_jikl);
 
621
        if (kl_ab != lk_ab)
 
622
          B_gbc2_ab.set_element( ji_ab, lk_ab, B_ab_jilk);
 
623
      }
 
624
 
 
625
      if (ij_aa != -1 && kl_aa != -1) {
 
626
        B_gbc2_aa.set_element( ij_aa, kl_aa, B_ab_ijkl - B_ab_jikl);
 
627
      }
 
628
 
 
629
    }
 
630
  }
 
631
 
 
632
  if (debug_ > 1) {
 
633
    B_gbc2_aa.print("Alpha-alpha B(GBC2) contribution");
 
634
    B_gbc2_ab.print("Alpha-beta B(GBC2) contribution");
 
635
  }
 
636
  // Symmetrize the B contribution
 
637
  B_gbc2_aa.scale(0.5);
 
638
  B_gbc2_ab.scale(0.5);
 
639
  RefSCMatrix B_gbc2_aa_t = B_gbc2_aa.t();
 
640
  Baa_.accumulate(B_gbc2_aa); Baa_.accumulate(B_gbc2_aa_t);
 
641
  RefSCMatrix B_gbc2_ab_t = B_gbc2_ab.t();
 
642
  Bab_.accumulate(B_gbc2_ab); Bab_.accumulate(B_gbc2_ab_t);
 
643
 
 
644
  globally_sum_intermeds_();
 
645
 
 
646
  ExEnv::out0() << decindent;
 
647
  ExEnv::out0() << indent << "Exited B(GBC2) intermediate evaluator" << endl;
 
648
 
 
649
  tim_exit("B(GBC2) intermediate");
 
650
}
 
651
 
 
652
////////////////////////////////////////////////////////////////////////////
 
653
 
 
654
// Local Variables:
 
655
// mode: c++
 
656
// c-file-style: "CLJ-CONDENSED"
 
657
// End: