MPQC  3.0.0-alpha
contract_tbint_tensors_to_obtensor.h
1 //
2 // contract_tbint_tensors_to_obtensor.h
3 //
4 // Copyright (C) 2008 Martin Torheyden
5 //
6 // Author: Martin Torheyden <mtorhey@vt.edu>
7 //
8 // This file is part of the SC Toolkit.
9 //
10 // The SC Toolkit is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Library General Public License as published by
12 // the Free Software Foundation; either version 2, or (at your option)
13 // any later version.
14 //
15 // The SC Toolkit is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Library General Public License for more details.
19 //
20 // You should have received a copy of the GNU Library General Public License
21 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
22 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
23 //
24 // The U.S. Government is granted a limited license as per AL 91-7.
25 //
26 
27 
28 #ifndef _chemistry_qc_mbptr12_contract_tbint_tensors_to_obtensor_h
29 #define _chemistry_qc_mbptr12_contract_tbint_tensors_to_obtensor_h
30 
31 #include <cmath>
32 #include <math/mmisc/pairiter.h>
33 #include <chemistry/qc/lcao/utils.h>
34 #include <chemistry/qc/lcao/utils.impl.h>
35 #include <chemistry/qc/mbptr12/r12int_eval.h>
36 #include <util/misc/print.h>
37 
38 
39 namespace sc {
40 
41  template <typename DataProcess_Bra,
42  typename DataProcess_Ket,
43  bool CorrFactorInBra,
44  bool CorrFactorInKet,
45  bool CorrFactorInInt>
46  void R12IntEval::contract_tbint_tensors_to_obtensor(
47  RefSCMatrix& T,
48  SpinCase2 pairspin,
49  TwoBodyTensorInfo tbtensor_type_bra,
50  TwoBodyTensorInfo tbtensor_type_ket,
51  const Ref<OrbitalSpace>& space1_bra,
52  const Ref<OrbitalSpace>& space1_intb,
53  const Ref<OrbitalSpace>& space2_intb,
54  const Ref<OrbitalSpace>& space3_intb,
55  const Ref<OrbitalSpace>& space1_ket,
56  const Ref<OrbitalSpace>& space1_intk,
57  const Ref<OrbitalSpace>& space2_intk,
58  const Ref<OrbitalSpace>& space3_intk,
59  const Ref<mbptr12::TwoParticleContraction>& tpcontract,
60  const std::vector<std::string>& tformkeys_bra,
61  const std::vector<std::string>& tformkeys_ket) {
62 
63  bool r12coeffs_bra = (tbtensor_type_bra.twobodytensor_type()==TwoBodyTensorInfo::geminalcoeff) ? true : false;
64  bool r12coeffs_ket = (tbtensor_type_ket.twobodytensor_type()==TwoBodyTensorInfo::geminalcoeff) ? true : false;
65 
66  TwoBodyOper::type tbint_type_bra = (r12coeffs_bra) ? corrfactor()->tbint_type_eri() // some dummy integral type
67  : tbtensor_type_bra.twobodyint_type();
68  TwoBodyOper::type tbint_type_ket = (r12coeffs_ket) ? corrfactor()->tbint_type_eri() // some dummy integral type
69  : tbtensor_type_ket.twobodyint_type();
70 
71  // are external spaces of particles 1 and 2 equivalent?
72  const bool part1_strong_equiv_part2 = (space1_bra==space2_intb &&
73  space1_ket==space2_intk);
74 
75  // can external spaces of particles 1 and 2 be equivalent?
76  const bool part1_weak_equiv_part2 = false;
77 
78  bool correct_semantics = ( (space1_intb->rank() == space1_intk->rank()) &&
79  (space2_intb->rank() == space2_intk->rank()) &&
80  (space3_intb->rank() == space3_intk->rank()) );
81 
82  // also dimensions of tpcontract must match those of space2_intb and space3_intb
83  correct_semantics = ( correct_semantics &&
84  (tpcontract->nrow() == space2_intb->rank()) &&
85  (tpcontract->ncol() == space3_intb->rank()) );
86 
87  if (!correct_semantics)
88  throw ProgrammingError("R12IntEval::contract_tbint_tensor_() -- incorrect call semantics",
89  __FILE__,__LINE__);
90 
91  // are internal spaces of particles 1 and 2 equivalent?
92  const bool part1_intequiv_part2 = (space2_intb==space3_intb &&
93  space2_intk==space3_intk);
94 
95  const bool alphabeta = !(part1_strong_equiv_part2 &&
96  part1_intequiv_part2);
97  //const bool alphabeta = (pairspin==AlphaBeta) ? true : false;
98  // NOTE! Even if computing in AlphaBeta, internal sums can be over AlphaAlpha!!!
99 
100  // Using spinorbital iterators means I don't take into account perm symmetry
101  //const SpinCase2 S = (alphabeta ? AlphaBeta : AlphaAlpha);
102  const SpinCase2 S = pairspin;
103 
104  const bool antisymmetrize = !alphabeta;
105 
106  const bool CorrFactorInBraInt = CorrFactorInBra && CorrFactorInInt;
107  const bool CorrFactorInKetInt = CorrFactorInKet && CorrFactorInInt;
108 
109  const unsigned int nbrasets = (CorrFactorInBra ? corrfactor()->nfunctions() : 1);
110  const unsigned int nketsets = (CorrFactorInKet ? corrfactor()->nfunctions() : 1);
111  const unsigned int nintsets = (CorrFactorInInt ? corrfactor()->nfunctions() : 1);
112  const unsigned int nbraintsets = (CorrFactorInBraInt ? UINT_MAX : nbrasets*nintsets);
113  const unsigned int nketintsets = (CorrFactorInKetInt ? UINT_MAX : nketsets*nintsets);
114 
115 
116  //
117  // create transforms, if needed
118  //
119  typedef std::vector<Ref<TwoBodyMOIntsTransform> > tformvec;
120 
121  // bra transforms
122  const size_t num_tforms_bra = tformkeys_bra.size();
123  tformvec transforms_bra(num_tforms_bra);
124  for (unsigned int t = 0; t < num_tforms_bra; ++t) {
125  transforms_bra[t] = moints_runtime4()->get(tformkeys_bra[t]);
126  }
127 
128  // ket transforms
129  const size_t num_tforms_ket = tformkeys_ket.size();
130  tformvec transforms_ket(num_tforms_ket);
131  for (unsigned int t = 0; t < num_tforms_ket; ++t) {
132  transforms_ket[t] = moints_runtime4()->get(tformkeys_ket[t]);
133  }
134 
135  //
136  // Generate contract label
137  //
138  Timer timer("Generic tensor contract");
139  std::string label;
140  {
141  std::ostringstream oss_bra;
142  oss_bra << "<" << space1_bra->id() << " " << space1_intb->id() << (antisymmetrize ? "||" : "|")
143  << space2_intb->id() << " " << space3_intb->id() << ">";
144  const std::string label_bra = oss_bra.str();
145  std::ostringstream oss_ket;
146  oss_ket << "<" << space1_ket->id() << " " << space1_intk->id() << (antisymmetrize ? "||" : "|")
147  << space2_intk->id() << " " << space3_intk->id() << ">";
148  const std::string label_ket = oss_ket.str();
149  std::ostringstream oss;
150  oss << "<" << space1_bra->id() << (antisymmetrize ? "||" : "|")
151  << space1_ket->id() << "> = "
152  << label_bra << " . " << label_ket << "^T";
153  label = oss.str();
154  }
155  ExEnv::out0() << std::endl << indent
156  << "Entered generic contraction (" << label << ")" << std::endl;
157  ExEnv::out0() << incindent;
158 
159  //
160  // Construct maps
161  //
162  // WARNING: Assuming all transforms are over same spaces!!!
163  //
164  Ref<OrbitalSpace> tspace1_bra = transforms_bra[0]->space1();
165  Ref<OrbitalSpace> tspace1_intb = transforms_bra[0]->space3();
166  Ref<OrbitalSpace> tspace2_intb = transforms_bra[0]->space2();
167  Ref<OrbitalSpace> tspace3_intb = transforms_bra[0]->space4();
168  Ref<OrbitalSpace> tspace1_ket = transforms_ket[0]->space1();
169  Ref<OrbitalSpace> tspace1_intk = transforms_ket[0]->space3();
170  Ref<OrbitalSpace> tspace2_intk = transforms_ket[0]->space2();
171  Ref<OrbitalSpace> tspace3_intk = transforms_ket[0]->space4();
172  // maps spaceX to spaceX of the transform
173  std::vector<unsigned int> map1_bra, map1_ket,
174  map1_intb, map2_intb, map3_intb, map1_intk, map2_intk, map3_intk;
175  // maps space3_intb to space2_intb of transform
176  std::vector<unsigned int> map23_intb;
177  // maps space2_intb to space3_intb of transform
178  std::vector<unsigned int> map32_intb;
179  // maps space3_intk to space2_intk of transform
180  std::vector<unsigned int> map23_intk;
181  // maps space2_intk to space3_intk of transform
182  std::vector<unsigned int> map32_intk;
183 
184  { // bra maps
185  map1_bra = *tspace1_bra<<*space1_bra;
186  map1_intb = *tspace1_intb<<*space1_intb;
187  map2_intb = *tspace2_intb<<*space2_intb;
188  map3_intb = *tspace3_intb<<*space3_intb;
189  // Will antisymmetrize the integrals? Then need ijkl AND ijlk
190  if(S!=AlphaBeta) {
191  if (tspace2_intb == tspace3_intb) {
192  map23_intb = map2_intb;
193  map32_intb = map3_intb;
194  }
195  else {
196  map23_intb = *tspace2_intb<<*space3_intb;
197  map32_intb = *tspace3_intb<<*space2_intb;
198  }
199  }
200  }
201 
202  { // ket maps
203  map1_ket = *tspace1_ket<<*space1_ket;
204  map1_intk = *tspace1_intk<<*space1_intk;
205  map2_intk = *tspace2_intk<<*space2_intk;
206  map3_intk = *tspace3_intk<<*space3_intk;
207  // Will antisymmetrize the integrals? Then need ijkl AND ijlk
208  if(S!=AlphaBeta) {
209  if (tspace2_intk == tspace3_intk) {
210  map23_intk = map2_intk;
211  map32_intk = map3_intk;
212  }
213  else {
214  map23_intk = *tspace2_intk<<*space3_intk;
215  map32_intk = *tspace3_intk<<*space2_intk;
216  }
217  }
218  }
219 
220  const unsigned int bratform_block_ncols = tspace3_intb->rank();
221  const unsigned int kettform_block_ncols = tspace3_intk->rank();
222  const RefDiagSCMatrix evals1_bra = space1_bra->evals();
223  const RefDiagSCMatrix evals1_ket = space1_ket->evals();
224  const RefDiagSCMatrix evals1_intb = space1_intb->evals();
225  const RefDiagSCMatrix evals2_intb = space2_intb->evals();
226  const RefDiagSCMatrix evals3_intb = space3_intb->evals();
227  const RefDiagSCMatrix evals1_intk = space1_intk->evals();
228  const RefDiagSCMatrix evals2_intk = space2_intk->evals();
229  const RefDiagSCMatrix evals3_intk = space3_intk->evals();
230 
231  SpinMOPairIter iterint(space2_intb->rank(),(alphabeta ? space3_intb : space2_intb)->rank(),S);
232  // size of one block of |space1_int space2_int>
233  const unsigned int nint = iterint.nij();
234  RefSCDimension space1_bra_dim = space1_bra->dim();
235  RefSCDimension space1_ket_dim = space1_ket->dim();
236  const unsigned int nbra = space1_bra_dim->n();
237  const unsigned int nket = space1_ket_dim->n();
238 
239  // size of integral blocks to contract
240  const unsigned int blksize_int = space2_intb->rank() * space3_intb->rank();
241  double* T_alphap = new double[blksize_int];
242  double* T_qp = new double[blksize_int];
243 
244 
245  int nmo_space1_bra = space1_bra->rank();
246  int nmo_space1_intb = space1_intb->rank();
247  int nmo_space1_ket = space1_ket->rank();
248  int nmo_space1_intk = space1_intk->rank();
249 
250  // outermost loop over contraction blocks to minimize number of activate()/deactivate() calls
251  for(unsigned int fint=0; fint<nintsets; ++fint) {
252  unsigned int fbraoffset = 0;
253  for(unsigned int fbra=0; fbra<nbrasets; ++fbra,fbraoffset+=nbra) {
254  const unsigned int fbraint = fbra*nintsets+fint;
255  Ref<TwoBodyMOIntsTransform> tformb = transforms_bra[fbraint];
256  const Ref<TwoBodyIntDescr>& intdescrb = tformb->intdescr();
257  const unsigned int intsetidx_bra = intdescrb->intset(tbint_type_bra);
258 
259  Ref<DistArray4> accumb = tformb->ints_distarray4();
260  // if transforms have not been computed yet, compute
261  if (accumb.null()) {
262  tformb->compute();
263  accumb = tformb->ints_distarray4();
264  }
265  accumb->activate();
266 
267  unsigned int fketoffset = 0;
268  for(unsigned int fket=0; fket<nketsets; ++fket,fketoffset+=nket) {
269  const unsigned int fketint = fket*nintsets+fint;
270  Ref<TwoBodyMOIntsTransform> tformk = transforms_ket[fketint];
271  const Ref<TwoBodyIntDescr>& intdescrk = tformk->intdescr();
272  const unsigned int intsetidx_ket = intdescrk->intset(tbint_type_ket);
273 
274  Ref<DistArray4> accumk = tformk->ints_distarray4();
275  if (accumk.null()) {
276  tformk->compute();
277  accumk = tformk->ints_distarray4();
278  }
279  accumk->activate();
280 
281  if (debug_ >= DefaultPrintThresholds::diagnostics) {
282  ExEnv::out0() << indent << "Using transforms "
283  << tformb->name() << " and "
284  << tformk->name() << std::endl;
285  }
286 
287  // split work over tasks which have access to integrals
288  // WARNING: assuming same accessibility for both bra and ket transforms
289  std::vector<int> proc_with_ints;
290  const int nproc_with_ints = accumb->tasks_with_access(proc_with_ints);
291  const int me = r12world()->world()->msg()->me();
292  const bool nket_ge_nevals = (nket >= nproc_with_ints);
293 
294  Timer tim_mo_ints_retrieve;
295  tim_mo_ints_retrieve.set_default("MO ints retrieve");
296  if (accumb->has_access(me)) {
297  unsigned int alphapq = 0;
298  for(unsigned int alpha=0; alpha<nmo_space1_bra; alpha++) { // external bra index
299  for(unsigned int p=0; p<nmo_space1_intb; p++) { // internal bra index
300  unsigned int alphap = alpha*nmo_space1_intb + p;
301  bool fetch_alphap_block = false;
302  if (nket_ge_nevals) {
303  fetch_alphap_block = true;
304  }
305  else {
306  const int first_alphap_task = alphapq%nproc_with_ints;
307  const int last_alphap_task = (alphapq+nket-1)%nproc_with_ints;
308  // figure out if for this alphap there is alphapq to be handled by me
309  if ( !(first_alphap_task > proc_with_ints[me] && proc_with_ints[me] > last_alphap_task) )
310  fetch_alphap_block = true;
311  }
312  if(!fetch_alphap_block)
313  continue;
314 
315  const unsigned int alphaalpha = map1_bra[alpha];
316  const unsigned int pp = map1_intb[p];
317 
318  const double *alphap_buf;
319  if(!r12coeffs_bra) {
320  tim_mo_ints_retrieve.enter_default();
321  alphap_buf = accumb->retrieve_pair_block(alphaalpha,pp,intsetidx_bra);
322  tim_mo_ints_retrieve.exit_default();
323  }
324 
325  if (debug_ >= DefaultPrintThresholds::mostO4)
326  ExEnv::outn() << indent << "task " << me << ": obtained alphap blocks" << std::endl;
327 
328  for(unsigned int q=0; q<nmo_space1_ket; q++) { // external ket index
329  unsigned int pq = p*nmo_space1_ket + q; // pq is a rectangular index pair
330  const int alphapq_proc = alphapq%nproc_with_ints;
331  if(pairspin==AlphaBeta) {
332  if (alphapq_proc != proc_with_ints[me])
333  continue;
334  const unsigned int qq = map1_ket[q];
335 
336  if (debug_ >= DefaultPrintThresholds::mostO4)
337  ExEnv::outn() << indent << "task " << me << ": working on <alpha,p | q,p> = <"
338  << alpha << "," << p << " | " << q << "," << p << ">" << std::endl;
339 
340  const double *qp_buf;
341  if(!r12coeffs_ket) {
342  tim_mo_ints_retrieve.enter_default();
343  qp_buf = accumk->retrieve_pair_block(qq,pp,intsetidx_ket);
344  tim_mo_ints_retrieve.exit_default();
345  }
346 
347  if (debug_ >= DefaultPrintThresholds::mostO4)
348  ExEnv::outn() << indent << "task " << me << ": obtained qp blocks" << std::endl;
349 
350  // zero out intblocks
351  memset(T_alphap, 0, blksize_int*sizeof(double));
352  memset(T_qp, 0, blksize_int*sizeof(double));
353 
354  if (debug_ >= DefaultPrintThresholds::mostO4) {
355  ExEnv::out0() << indent << "alpha = " << alpha << " p = " << p << " q = " << q
356  << incindent << std::endl;
357  }
358 
359  for(iterint.start(); iterint; iterint.next()) { // internal pair index. Add \f$ \bar{r}_{p_{\gamma f\gamma} \alpha}^{rs} t_{rs}^{pq} \f$. rs is a rectangular index pair.
360  unsigned int r = iterint.i();
361  unsigned int s = iterint.j();
362  unsigned int rs = iterint.ij();
363 
364  int RS_bra,RS_ket;
365  {
366  const unsigned int rr = map2_intb[r];
367  const unsigned int ss = map3_intb[s];
368  RS_bra = rr*bratform_block_ncols + ss;
369  }
370  {
371  const unsigned int rr = map2_intk[r];
372  const unsigned int ss = map3_intk[s];
373  RS_ket = rr*kettform_block_ncols + ss;
374  }
375  double I_alphaprs;
376  if(!r12coeffs_bra) {
377  I_alphaprs = alphap_buf[RS_bra];
378  }
379  else {
380  I_alphaprs = C_CuspConsistent(alpha,p,r,s,S);
381  }
382  double I_qprs;
383  if(!r12coeffs_ket) {
384  I_qprs = qp_buf[RS_ket];
385  }
386  else {
387  I_qprs = C_CuspConsistent(q,p,r,s,S);
388  }
389 
390  if (debug_ >= DefaultPrintThresholds::mostO6) {
391  ExEnv::out0() << indent << " r = " << r << " s = " << s << std::endl;
392  }
393 
394  if (debug_ >= DefaultPrintThresholds::mostO6) {
395  ExEnv::out0() << indent << " <alphap|rs> = " << I_alphaprs << std::endl
396  << indent << " <qp|rs> = " << I_qprs << std::endl;
397  }
398 
399  double T_alphaprs;
400  if(!r12coeffs_bra) {
401  T_alphaprs = DataProcess_Bra::I2T(I_alphaprs,alpha,p,r,s,
402  evals1_bra,evals2_intb,evals1_intb,evals3_intb);
403  }
404  else {
405  T_alphaprs = I_alphaprs;
406  }
407 
408  double T_qprs;
409  if(!r12coeffs_ket) {
410  T_qprs = DataProcess_Ket::I2T(I_qprs,q,p,r,s,
411  evals1_ket,evals2_intk,evals1_intk,evals3_intk);
412  }
413  else {
414  T_qprs = I_qprs;
415  }
416 
417  if (debug_ >= DefaultPrintThresholds::mostO6) {
418  ExEnv::out0() << indent << " <alphap|T|rs> = " << T_alphaprs << std::endl
419  << indent << " <qp|T|rs> = " << T_qprs << std::endl;
420  }
421 
422  T_alphap[rs] = T_alphaprs;
423  T_qp[rs] = T_qprs;
424 
425  } // loop over iterint, i.e. rs
426 
427  // contract matrices
428  double T_alphapqp = tpcontract->contract(T_alphap,T_qp);
429  T_alphapqp *= 2.0;
430 
431  if (debug_ >= DefaultPrintThresholds::mostO4) {
432  ExEnv::out0() << decindent << indent
433  << " <alphap|qp> = " << T_alphapqp << std::endl;
434  }
435 
436  T.accumulate_element(alpha,q,T_alphapqp);
437 
438  if(!r12coeffs_ket) {
439  accumk->release_pair_block(qq,pp,intsetidx_ket);
440  }
441 
442  alphapq += 1;
443  } // if(pairspin==AlphaBeta)
444  else { // if(pairspin!=AlphaBeta)
445  if(q==p) continue; // q==p is not allowed for pairspin!=AlphaBeta -> skip.
446  if (alphapq_proc != proc_with_ints[me])
447  continue;
448  const unsigned int qq = map1_ket[q];
449 
450  if (debug_ >= DefaultPrintThresholds::mostO4)
451  ExEnv::outn() << indent << "task " << me << ": working on <alpha,p | q,p> = <"
452  << alpha << "," << p << " | " << q << "," << p << ">" << std::endl;
453 
454  const double *qp_buf;
455  const double *pq_buf;
456  if(!r12coeffs_ket) {
457  tim_mo_ints_retrieve.enter_default();
458  if(q>p){
459  qp_buf = accumk->retrieve_pair_block(qq,pp,intsetidx_ket);
460  }
461  else {
462  pq_buf = accumk->retrieve_pair_block(pp,qq,intsetidx_ket);
463  }
464  tim_mo_ints_retrieve.exit_default();
465  }
466 
467  if (debug_ >= DefaultPrintThresholds::mostO4)
468  ExEnv::outn() << indent << "task " << me << ": obtained qp blocks" << std::endl;
469 
470  // zero out intblocks
471  memset(T_alphap, 0, blksize_int*sizeof(double));
472  memset(T_qp, 0, blksize_int*sizeof(double));
473 
474  if (debug_ >= DefaultPrintThresholds::mostO4) {
475  ExEnv::out0() << indent << "alpha = " << alpha << " p = " << p << " q = " << q
476  << incindent << std::endl;
477  }
478 
479  for(iterint.start(); iterint; iterint.next()) { // internal pair index
480  unsigned int r = iterint.i();
481  unsigned int s = iterint.j();
482  unsigned int rs = iterint.ij();
483 
484  int RS_bra,RS_ket;
485  {
486  const unsigned int rr = map2_intb[r];
487  const unsigned int ss = map3_intb[s];
488  RS_bra = rr*bratform_block_ncols + ss;
489  }
490  {
491  const unsigned int rr = map2_intk[r];
492  const unsigned int ss = map3_intk[s];
493  RS_ket = rr*kettform_block_ncols + ss;
494  }
495  double I_alphaprs;
496  if(!r12coeffs_bra) {
497  I_alphaprs = alphap_buf[RS_bra];
498  }
499  else {
500  I_alphaprs = C_CuspConsistent(alpha,p,r,s,S);
501  }
502  // getting integrals with exchanged indices
503  int SR_bra, SR_ket;
504  {
505  const unsigned int rr = map32_intb[r];
506  const unsigned int ss = map23_intb[s];
507  SR_bra = ss*bratform_block_ncols + rr;
508  }
509  {
510  const unsigned int rr = map32_intk[r];
511  const unsigned int ss = map23_intk[s];
512  SR_ket = ss*kettform_block_ncols + rr;
513  }
514 
515  double I_alphapsr;
516  if(!r12coeffs_bra) {
517  I_alphapsr = alphap_buf[SR_bra];
518  }
519  else {
520  I_alphapsr = C_CuspConsistent(alpha,p,s,r,S);
521  }
522 
523  if (debug_ >= DefaultPrintThresholds::mostO6) {
524  ExEnv::out0() << indent << " r = " << r << " s = " << s << std::endl;
525  }
526 
527  if(q>p) {
528  double I_qprs;
529  if(!r12coeffs_ket) {
530  I_qprs = qp_buf[RS_ket];
531  }
532  else {
533  I_qprs = C_CuspConsistent(q,p,r,s,S);
534  }
535 
536  if (debug_ >= DefaultPrintThresholds::mostO6) {
537  ExEnv::out0() << indent << " <alphap|rs> = " << I_alphaprs << std::endl
538  << indent << " <qp|rs> = " << I_qprs << std::endl;
539  }
540 
541  double I_qpsr;
542  if(!r12coeffs_ket) {
543  I_qpsr = qp_buf[SR_ket];
544  }
545  else {
546  I_qpsr = C_CuspConsistent(q,p,s,r,S);
547  }
548 
549  if (debug_ >= DefaultPrintThresholds::mostO6) {
550  ExEnv::out0() << " <alphap|rs> = " << I_alphaprs << std::endl
551  << " <alphap|sr> = " << I_alphapsr << std::endl
552  << " <qp|rs> = " << I_qprs << std::endl
553  << " <qp|sr> = " << I_qpsr << std::endl;
554  }
555 
556  double T_alphaprs;
557  if(!r12coeffs_bra) {
558  T_alphaprs = DataProcess_Bra::I2T(I_alphaprs-I_alphapsr,alpha,p,r,s,
559  evals1_bra,evals2_intb,evals1_intb,evals3_intb);
560  }
561  else {
562  T_alphaprs = I_alphaprs;
563  }
564 
565  double T_qprs;
566  if(!r12coeffs_ket) {
567  T_qprs = DataProcess_Ket::I2T(I_qprs-I_qpsr,q,p,r,s,
568  evals1_ket,evals2_intk,evals1_intk,evals3_intk);
569  }
570  else {
571  T_qprs = I_qprs;
572  }
573 
574  T_alphap[rs] = T_alphaprs;
575  T_qp[rs] = T_qprs;
576 
577  } // if(q>p)
578  else { // if(p>q). Add \f$ \bar{r}_{p_{\gamma f\gamma} \alpha}^{rs} t_{rs}^{pq} \f$. pq is a triangular index.
579  double I_pqrs;
580  if(!r12coeffs_ket) {
581  I_pqrs = pq_buf[RS_ket];
582  }
583  else {
584  I_pqrs = C_CuspConsistent(p,q,r,s,S);
585  }
586 
587  double I_pqsr;
588  if(!r12coeffs_ket) {
589  I_pqsr = pq_buf[SR_ket];
590  }
591  else {
592  I_pqsr = C_CuspConsistent(p,q,s,r,S);
593  }
594 
595  if (debug_ >= DefaultPrintThresholds::mostO6) {
596  ExEnv::out0() << " <alphap|rs> = " << I_alphaprs << std::endl
597  << " <alphap|rs> = " << I_alphapsr << std::endl
598  << " <pq|rs> = " << I_pqrs << std::endl
599  << " <pq|sr> = " << I_pqsr << std::endl;
600  }
601 
602  double T_alphaprs;
603  if(!r12coeffs_bra) {
604  T_alphaprs = DataProcess_Bra::I2T(I_alphaprs-I_alphapsr,alpha,p,r,s,
605  evals1_bra,evals2_intb,evals1_intb,evals3_intb);
606  }
607  else {
608  T_alphaprs = I_alphaprs;
609  }
610 
611  double T_qprs;
612  if(!r12coeffs_ket) {
613  T_qprs = DataProcess_Ket::I2T(I_pqsr-I_pqrs,q,p,r,s,
614  evals1_ket,evals2_intk,evals1_intk,evals3_intk);
615  }
616  else {
617  T_qprs = -I_pqrs;
618  }
619 
620  T_alphap[rs] = T_alphaprs;
621  T_qp[rs] = T_qprs;
622 
623  } // if(p>q)
624  } // loop over iterint, i.e. rs
625 
626  // contract matrices
627  double T_alphapqp = tpcontract->contract(T_alphap,T_qp);
628 
629  if (debug_ >= DefaultPrintThresholds::mostO4) {
630  ExEnv::out0() << decindent << indent
631  << " <alphap|qp> = " << T_alphapqp << std::endl;
632  }
633 
634  T.accumulate_element(alpha,q,T_alphapqp);
635 
636  if(!r12coeffs_ket) {
637  if(q>p) {
638  accumk->release_pair_block(qq,pp,intsetidx_ket);
639  }
640  else {
641  accumk->release_pair_block(pp,qq,intsetidx_ket);
642  }
643  }
644 
645  alphapq += 1;
646  } // if(pairspin!=AlphaBeta)
647  } // loop over q
648 
649  if(fetch_alphap_block && (!r12coeffs_bra)) {
650  accumb->release_pair_block(alphaalpha,pp,intsetidx_bra);
651  }
652 
653  } // loop over p
654  } // loop over alpha
655 
656  }
657 
658  if (accumb != accumk) {
659  if (accumk->data_persistent()) accumk->deactivate();
660  }
661 
662  } // fket loop
663  if (accumb->data_persistent()) accumb->deactivate();
664  } // fbra loop
665  } // fint loop
666 
667  delete[] T_alphap;
668  delete[] T_qp;
669 
670  ExEnv::out0() << decindent;
671  ExEnv::out0() << indent << "Exited generic contraction (" << label << ")" << std::endl;
672 
673  timer.exit();
674  }
675 
676 }
677 
678 #endif /* _chemistry_qc_mbptr12_contract_tbint_tensors_to_obtensor_h */
sc::Timer
The Timer class uses RegionTimer to time intervals in an exception safe manner.
Definition: regtime.h:181
sc::RefSCMatrix
The RefSCMatrix class is a smart pointer to an SCMatrix specialization.
Definition: matrix.h:135
sc::Ref
A template class that maintains references counts.
Definition: ref.h:361
sc::RefDiagSCMatrix
The RefDiagSCMatrix class is a smart pointer to an DiagSCMatrix specialization.
Definition: matrix.h:389
sc::RefSCDimension
The RefSCDimension class is a smart pointer to an SCDimension specialization.
Definition: dim.h:152
sc::Ref::null
bool null() const
Return true if this is a reference to a null object.
Definition: ref.h:423
sc::DefaultPrintThresholds::diagnostics
static const unsigned int diagnostics
Additional diagnostics, e.g., transform progress, etc.
Definition: print.h:39
sc::Timer::exit
void exit(const char *region=0)
Exit the current timing region.
sc::R12IntEval::r12world
const Ref< R12WavefunctionWorld > & r12world() const
the R12World in which this object lives
Definition: r12int_eval.h:642
sc::SpinMOPairIter
SpinMOPairIter iterates over pairs of spinorbitals.
Definition: pairiter.h:276
sc::MOPairIter::nij
int nij() const
Returns the number of pair combinations for this iterator.
Definition: pairiter.h:77
sc::TwoBodyTensorInfo
Provides information about the type of a two body tensor.
Definition: twobodytensorinfo.h:37
sc::ProgrammingError
This is thrown when a situations arises that should be impossible.
Definition: scexception.h:92
sc::TwoBodyOper::type
type
types of known two-body operators
Definition: operator.h:318
sc::R12IntEval::C_CuspConsistent
double C_CuspConsistent(int i, int j, int k, int l, SpinCase2 pairspin)
Returns the cusp consistent coefficient .
sc::DefaultPrintThresholds::mostO4
static const unsigned int mostO4
Print most o^4 quantities.
Definition: print.h:57
sc::ExEnv::out0
static std::ostream & out0()
Return an ostream that writes from node 0.
sc::DefaultPrintThresholds::mostO6
static const unsigned int mostO6
Print most o^6 quantities.
Definition: print.h:75
sc::ExEnv::outn
static std::ostream & outn()
Return an ostream that writes from all nodes.
Definition: exenv.h:81
sc::antisymmetrize
void antisymmetrize(RefSCMatrix &Aanti, const RefSCMatrix &A, const Ref< OrbitalSpace > &bra, const Ref< OrbitalSpace > &ket, bool accumulate=false)
Antisymmetrizes 4-index quantity <ij|A|kl> -> <ij|A|kl> - <ij|A|lk> and saves to Aanti.
sc::Timer::set_default
void set_default(const char *region)
Default timing regions are provided as a performance optimization.
sc::R12IntEval::contract_tbint_tensors_to_obtensor
void contract_tbint_tensors_to_obtensor(RefSCMatrix &T, SpinCase2 pairspin, TwoBodyTensorInfo tbtensor_type_bra, TwoBodyTensorInfo tbtensor_type_ket, const Ref< OrbitalSpace > &space1_bra, const Ref< OrbitalSpace > &space1_intb, const Ref< OrbitalSpace > &space2_intb, const Ref< OrbitalSpace > &space3_intb, const Ref< OrbitalSpace > &space1_ket, const Ref< OrbitalSpace > &space1_intk, const Ref< OrbitalSpace > &space2_intk, const Ref< OrbitalSpace > &space3_intk, const Ref< mbptr12::TwoParticleContraction > &tpcontract, const std::vector< std::string > &tformkeys_bra, const std::vector< std::string > &tformkeys_ket)
<space1bra space1_intb |Tbra| space2_intb space3_intb> * <space2_intk space3_intk |Tket| space1ket sp...
Definition: contract_tbint_tensors_to_obtensor.h:46
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
Generated at Sun Jan 26 2020 23:23:59 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.