Psi4
sointegral_twobody.h
Go to the documentation of this file.
1 /*
2  * @BEGIN LICENSE
3  *
4  * Psi4: an open-source quantum chemistry software package
5  *
6  * Copyright (c) 2007-2018 The Psi4 Developers.
7  *
8  * The copyrights for code used from other parties are included in
9  * the corresponding files.
10  *
11  * This file is part of Psi4.
12  *
13  * Psi4 is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU Lesser General Public License as published by
15  * the Free Software Foundation, version 3.
16  *
17  * Psi4 is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU Lesser General Public License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public License along
23  * with Psi4; if not, write to the Free Software Foundation, Inc.,
24  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25  *
26  * @END LICENSE
27  */
28 
29 #ifndef _psi_src_lib_libmints_sointegral_h_
30 #define _psi_src_lib_libmints_sointegral_h_
31 
32 #include "typedefs.h"
33 
34 #include "onebody.h"
35 #include "twobody.h"
36 #include "integral.h"
37 #include "sobasis.h"
38 #include "gshell.h"
39 #include "petitelist.h"
40 #include "wavefunction.h"
41 #include "cdsalclist.h"
42 #include "dcd.h"
43 #include "basisset.h"
44 
45 #include "psi4/libqt/qt.h"
46 
47 #include <vector>
48 
49 //#define DebugPrint 1
50 
51 #ifdef _OPENMP
52 #include <omp.h>
53 #endif
54 
55 #ifndef DebugPrint
56 #define DebugPrint 0
57 #endif
58 
59 #if DebugPrint
60 #define dprintf(...) outfile->Printf(__VA_ARGS__)
61 #else
62 #define dprintf(...)
63 #endif
64 
65 namespace psi {
66 class BasisSet;
67 
68 // Only include the following function if Doxygen is running to generate appropriate
69 // documentation.
70 #ifdef DOXYGEN
71 class TwoBodySOIntFunctor {
72  public:
73  void operator()(int pirrep, int pso, int qirrep, int qso, int rirrep, int rso, int sirrep, int sso, double value);
74 };
75 #endif
76 
78  protected:
79  std::vector<std::shared_ptr<TwoBodyAOInt> > tb_;
80  std::shared_ptr<IntegralFactory> integral_;
81 
82  std::shared_ptr<SOBasisSet> b1_;
83  std::shared_ptr<SOBasisSet> b2_;
84  std::shared_ptr<SOBasisSet> b3_;
85  std::shared_ptr<SOBasisSet> b4_;
86 
87  size_t size_;
88  std::vector<double *> buffer_;
89  std::vector<double *> temp_;
90  std::vector<double *> temp2_;
91  std::vector<double **> deriv_;
92 
93  int iirrepoff_[8], jirrepoff_[8], kirrepoff_[8], lirrepoff_[8];
94 
95  std::shared_ptr<PetiteList> petite1_;
96  std::shared_ptr<PetiteList> petite2_;
97  std::shared_ptr<PetiteList> petite3_;
98  std::shared_ptr<PetiteList> petite4_;
99 
100  std::shared_ptr<PointGroup> pg_;
101 
102  std::shared_ptr<DCD> dcd_;
103 
105  int nthread_;
106  std::string comm_;
107  int nproc_;
108  int me_;
109 
110  double cutoff_;
111 
113 
114  template <typename TwoBodySOIntFunctor>
115  void provide_IJKL(int, int, int, int, TwoBodySOIntFunctor &body);
116 
117  template <typename TwoBodySOIntFunctor>
118  void provide_IJKL_deriv1(int ish, int jsh, int ksh, int lsh, TwoBodySOIntFunctor &body);
119 
120  void common_init();
121 
122  public:
123  // Constructor, assuming 1 thread
124  TwoBodySOInt(const std::shared_ptr<TwoBodyAOInt> &, const std::shared_ptr<IntegralFactory> &);
125  // Constructor, using vector of AO objects for threading
126  TwoBodySOInt(const std::vector<std::shared_ptr<TwoBodyAOInt> > &tb,
127  const std::shared_ptr<IntegralFactory> &integral);
128  TwoBodySOInt(const std::shared_ptr<TwoBodyAOInt> &aoint, const std::shared_ptr<IntegralFactory> &intfac,
129  const CdSalcList &cdsalcs);
130  TwoBodySOInt(const std::vector<std::shared_ptr<TwoBodyAOInt> > &tb,
131  const std::shared_ptr<IntegralFactory> &integral, const CdSalcList &cdsalcs);
132 
133  virtual ~TwoBodySOInt();
134 
135  bool only_totally_symmetric() const { return only_totally_symmetric_; }
136  void set_only_totally_symmetric(bool ots) { only_totally_symmetric_ = ots; }
137 
138  std::shared_ptr<SOBasisSet> basis() const;
139  std::shared_ptr<SOBasisSet> basis1() const;
140  std::shared_ptr<SOBasisSet> basis2() const;
141  std::shared_ptr<SOBasisSet> basis3() const;
142  std::shared_ptr<SOBasisSet> basis4() const;
143 
144  const double *buffer(int thread = 0) const { return buffer_[thread]; }
145 
146  void set_cutoff(double ints_tolerance) { cutoff_ = ints_tolerance; }
147 
148  // Normal integrals
149  template <typename TwoBodySOIntFunctor>
150  void compute_shell(const SOShellCombinationsIterator &shellIter, TwoBodySOIntFunctor &body) {
151  compute_shell(shellIter.p(), shellIter.q(), shellIter.r(), shellIter.s(), body);
152  }
153 
154  template <typename TwoBodySOIntFunctor>
155  void compute_shell(int, int, int, int, TwoBodySOIntFunctor &body);
156 
157  // User provides an iterator object and this function will walk through it.
158  // Assumes serial run (nthread = 1)
159  template <typename ShellIter, typename TwoBodySOIntFunctor>
160  void compute_quartets(ShellIter &shellIter, TwoBodySOIntFunctor &body) {
161  for (shellIter->first(); shellIter->is_done() == false; shellIter->next()) {
162  this->compute_shell(shellIter->p(), shellIter->q(), shellIter->r(), shellIter->s(), body);
163  }
164  }
165 
166  // Compute integrals in parallel
167  template <typename TwoBodySOIntFunctor>
168  void compute_integrals(TwoBodySOIntFunctor &functor);
169 
170  // Derivative integrals
171  // User provides an iterator object and this function will walk through it.
172  template <typename ShellIter, typename TwoBodySOIntFunctor>
173  void compute_quartets_deriv1(ShellIter &shellIter, TwoBodySOIntFunctor &body) {
174  for (shellIter->first(); shellIter->is_done() == false; shellIter->next()) {
175  compute_shell_deriv1(shellIter->p(), shellIter->q(), shellIter->r(), shellIter->s(), body);
176  }
177  }
178 
179  template <typename TwoBodySOIntFunctor>
180  void compute_shell_deriv1(const SOShellCombinationsIterator &shellIter, TwoBodySOIntFunctor &body) {
181  compute_shell_deriv1(shellIter.p(), shellIter.q(), shellIter.r(), shellIter.s(), body);
182  }
183 
184  template <typename TwoBodySOIntFunctor>
185  void compute_shell_deriv1(int, int, int, int, TwoBodySOIntFunctor &body);
186 
187  // Compute integrals in parallel
188  template <typename TwoBodySOIntFunctor>
189  void compute_integrals_deriv1(TwoBodySOIntFunctor &functor);
190 
191  template <typename TwoBodySOIntFunctor>
192  int compute_pq_pair_deriv1(const int &p, const int &q, const size_t &pair_number, const TwoBodySOIntFunctor &body) {
193  const_cast<TwoBodySOIntFunctor &>(body).load_tpdm(pair_number);
194  auto shellIter = std::make_shared<SO_RS_Iterator>(p, q, b1_, b2_, b3_, b4_);
195 
196  compute_quartets_deriv1(shellIter, const_cast<TwoBodySOIntFunctor &>(body));
197 
198  return 0;
199  }
200 };
201 
202 template <typename TwoBodySOIntFunctor>
203 void TwoBodySOInt::compute_shell(int uish, int ujsh, int uksh, int ulsh, TwoBodySOIntFunctor &body) {
204  dprintf("uish %d, ujsh %d, uksh %d, ulsh %d\n", uish, ujsh, uksh, ulsh);
205 
206  int thread = 0;
207 #ifdef _OPENMP
208  thread = omp_get_thread_num();
209 #endif
210  // Old call WorldComm->thread_id(pthread_self());
211 
212  // timer_on("TwoBodySOInt::compute_shell overall");
213  mints_timer_on("TwoBodySOInt::compute_shell overall");
214  mints_timer_on("TwoBodySOInt::compute_shell setup");
215 
216  const double *aobuff = tb_[thread]->buffer();
217 
218  const SOTransform &t1 = b1_->sotrans(uish);
219  const SOTransform &t2 = b2_->sotrans(ujsh);
220  const SOTransform &t3 = b3_->sotrans(uksh);
221  const SOTransform &t4 = b4_->sotrans(ulsh);
222 
223  const int nso1 = b1_->nfunction(uish);
224  const int nso2 = b2_->nfunction(ujsh);
225  const int nso3 = b3_->nfunction(uksh);
226  const int nso4 = b4_->nfunction(ulsh);
227  const size_t nso = nso1 * nso2 * nso3 * nso4;
228 
229  const int nao2 = b2_->naofunction(ujsh);
230  const int nao3 = b3_->naofunction(uksh);
231  const int nao4 = b4_->naofunction(ulsh);
232 
233  const int iatom = tb_[thread]->basis1()->shell(t1.aoshell[0].aoshell).ncenter();
234  const int jatom = tb_[thread]->basis2()->shell(t2.aoshell[0].aoshell).ncenter();
235  const int katom = tb_[thread]->basis3()->shell(t3.aoshell[0].aoshell).ncenter();
236  const int latom = tb_[thread]->basis4()->shell(t4.aoshell[0].aoshell).ncenter();
237 
238  int nirrep = b1_->nirrep();
239 
240  mints_timer_on("TwoBodySOInt::compute_shell zero buffer");
241 
242  ::memset(buffer_[thread], 0, nso * sizeof(double));
243 
244  mints_timer_off("TwoBodySOInt::compute_shell zero buffer");
245  mints_timer_off("TwoBodySOInt::compute_shell setup");
246 
247  mints_timer_on("TwoBodySOInt::compute_shell full shell transform");
248 
249  // Get the atomic stablizer (the first symmetry operation that maps the atom
250  // onto itself.
251 
252  // These 3 sections are not shell specific so we can just use petite1_
253  const unsigned short istablizer = petite1_->stablizer(iatom);
254  const unsigned short jstablizer = petite1_->stablizer(jatom);
255  const unsigned short kstablizer = petite1_->stablizer(katom);
256  const unsigned short lstablizer = petite1_->stablizer(latom);
257 
258  const int istabdense = dcd_->bits_to_dense_numbering(istablizer);
259  const int jstabdense = dcd_->bits_to_dense_numbering(jstablizer);
260  const int kstabdense = dcd_->bits_to_dense_numbering(kstablizer);
261  const int lstabdense = dcd_->bits_to_dense_numbering(lstablizer);
262 
263  const int ijstablizer = dcd_->intersection(istabdense, jstabdense);
264  const int klstablizer = dcd_->intersection(kstabdense, lstabdense);
265  const int ijklstablizer = dcd_->intersection(ijstablizer, klstablizer);
266 
267  const int *R_list = dcd_->dcr(istabdense, jstabdense);
268  const int *S_list = dcd_->dcr(kstabdense, lstabdense);
269  const int *T_list = dcd_->dcr(ijstablizer, klstablizer);
270 
271  const int R_size = R_list[0];
272  const int S_size = S_list[0];
273  const int T_size = T_list[0];
274 
275  // Check with Andy on this:
276  int lambda_T = petite1_->nirrep() / dcd_->subgroup_dimensions(ijklstablizer);
277 
278  std::vector<int> sj_arr, sk_arr, sl_arr;
279 
280  int si = petite1_->unique_shell_map(uish, 0);
281  const int siatom = tb_[thread]->basis1()->shell(si).ncenter();
282 
283  dprintf("dcd %d", petite1_->group());
284  dprintf("istab %d, jstab %d, kstab %d, lstab %d, ijstab %d, klstab %d\n", istabdense, jstabdense, kstabdense,
285  lstabdense, ijstablizer, klstablizer);
286  dprintf("R_size %d, S_size %d, T_size %d\n", R_size, S_size, T_size);
287 
288  for (int ij = 1; ij <= R_size; ++ij) {
289  int sj = petite2_->unique_shell_map(ujsh, R_list[ij]);
290  const int sjatom = tb_[thread]->basis2()->shell(sj).ncenter();
291 
292  for (int ijkl = 1; ijkl <= T_size; ++ijkl) {
293  int sk = petite3_->unique_shell_map(uksh, T_list[ijkl]);
294  int llsh = petite4_->unique_shell_map(ulsh, T_list[ijkl]);
295  const int skatom = tb_[thread]->basis3()->shell(sk).ncenter();
296 
297  for (int kl = 1; kl <= S_size; ++kl) {
298  int sl = petite4_->shell_map(llsh, S_list[kl]);
299  const int slatom = tb_[thread]->basis4()->shell(sl).ncenter();
300 
301  // Check AM
302  int total_am = tb_[thread]->basis1()->shell(si).am() + tb_[thread]->basis2()->shell(sj).am() +
303  tb_[thread]->basis3()->shell(sk).am() + tb_[thread]->basis4()->shell(sl).am();
304 
305  if (!(total_am % 2) || (siatom != sjatom) || (sjatom != skatom) || (skatom != slatom)) {
306  sj_arr.push_back(sj);
307  sk_arr.push_back(sk);
308  sl_arr.push_back(sl);
309  }
310  }
311  }
312  }
313 
314  // Compute integral using si, sj_arr, sk_arr, sl_arr
315  // Loop over unique quartets
316  const AOTransform &s1 = b1_->aotrans(si);
317 
318  const unsigned short *ifuncpi = s1.nfuncpi;
319 
320  for (size_t n = 0; n < sj_arr.size(); ++n) {
321  int sj = sj_arr[n];
322  int sk = sk_arr[n];
323  int sl = sl_arr[n];
324 
325  const AOTransform &s2 = b2_->aotrans(sj);
326  const AOTransform &s3 = b3_->aotrans(sk);
327  const AOTransform &s4 = b4_->aotrans(sl);
328 
329  const unsigned short *jfuncpi = s2.nfuncpi;
330  const unsigned short *kfuncpi = s3.nfuncpi;
331  const unsigned short *lfuncpi = s4.nfuncpi;
332 
333  // Compute this unique AO shell
334  // timer_on("Computing the AO shell");
335  tb_[thread]->compute_shell(si, sj, sk, sl);
336  // timer_off("Computing the AO shell");
337 
338  mints_timer_on("TwoBodySOInt::compute_shell actual transform");
339 
340  for (int isym = 0; isym < nirrep; ++isym) {
341  unsigned short nifunc = ifuncpi[isym];
342  for (int itr = 0; itr < nifunc; itr++) {
343  mints_timer_on("itr");
344 
345  const AOTransformFunction &ifunc = s1.soshellpi[isym][itr];
346  double icoef = ifunc.coef;
347  int iaofunc = ifunc.aofunc;
348  int isofunc = ifunc.sofunc;
349  int iaooff = iaofunc;
350  int isooff = isofunc;
351 
352  mints_timer_off("itr");
353 
354  for (int jsym = 0; jsym < nirrep; ++jsym) {
355  unsigned short njfunc = jfuncpi[jsym];
356  for (int jtr = 0; jtr < njfunc; jtr++) {
357  mints_timer_on("jtr");
358 
359  const AOTransformFunction &jfunc = s2.soshellpi[jsym][jtr];
360  double jcoef = jfunc.coef * icoef;
361  int jaofunc = jfunc.aofunc;
362  int jsofunc = jfunc.sofunc;
363  int jaooff = iaooff * nao2 + jaofunc;
364  int jsooff = isooff * nso2 + jsofunc;
365 
366  mints_timer_off("jtr");
367 
368  for (int ksym = 0; ksym < nirrep; ++ksym) {
369  unsigned short nkfunc = kfuncpi[ksym];
370  for (int ktr = 0; ktr < nkfunc; ktr++) {
371  mints_timer_on("ktr");
372 
373  const AOTransformFunction &kfunc = s3.soshellpi[ksym][ktr];
374  double kcoef = kfunc.coef * jcoef;
375  int kaofunc = kfunc.aofunc;
376  int ksofunc = kfunc.sofunc;
377  int kaooff = jaooff * nao3 + kaofunc;
378  int ksooff = jsooff * nso3 + ksofunc;
379 
380  mints_timer_off("ktr");
381 
382  int lsym = isym ^ jsym ^ ksym;
383  unsigned short nlfunc = lfuncpi[lsym];
384  for (int ltr = 0; ltr < nlfunc; ltr++) {
385  mints_timer_on("ltr");
386 
387  const AOTransformFunction &lfunc = s4.soshellpi[lsym][ltr];
388  double lcoef = lfunc.coef * kcoef;
389  int laofunc = lfunc.aofunc;
390  int lsofunc = lfunc.sofunc;
391  int laooff = kaooff * nao4 + laofunc;
392  int lsooff = ksooff * nso4 + lsofunc;
393 
394  mints_timer_off("ltr");
395  mints_timer_on("transform");
396 
397  dprintf("lambda_T %d, lcoef %lf, aobuff %lf\n", lambda_T, lcoef, aobuff[laooff]);
398  buffer_[thread][lsooff] += lambda_T * lcoef * aobuff[laooff];
399 
400  mints_timer_off("transform");
401  }
402  }
403  }
404  }
405  }
406  }
407  }
408 
409  mints_timer_off("TwoBodySOInt::compute_shell actual transform");
410  }
411 
412  mints_timer_off("TwoBodySOInt::compute_shell full shell transform");
413 
414  provide_IJKL(uish, ujsh, uksh, ulsh, body);
415 
416  mints_timer_off("TwoBodySOInt::compute_shell overall");
417  // timer_off("TwoBodySOInt::compute_shell overall");
418 }
419 
420 template <typename TwoBodySOIntFunctor>
421 void TwoBodySOInt::provide_IJKL(int ish, int jsh, int ksh, int lsh, TwoBodySOIntFunctor &body) {
422  int thread = 0;
423 // Old call WorldComm->thread_id(pthread_self());
424 #ifdef _OPENMP
425  thread = omp_get_thread_num();
426 #endif
427 
428  mints_timer_on("TwoBodySOInt::provide_IJKL overall");
429  // timer_on("TwoBodySOInt::provide_IJKL overall");
430 
431  int nso2 = b2_->nfunction(jsh);
432  int nso3 = b3_->nfunction(ksh);
433  int nso4 = b4_->nfunction(lsh);
434 
435  int n1 = b1_->nfunction(ish);
436  int n2 = b2_->nfunction(jsh);
437  int n3 = b3_->nfunction(ksh);
438  int n4 = b4_->nfunction(lsh);
439 
440  int itr;
441  int jtr;
442  int ktr;
443  int ltr;
444 
445  for (itr = 0; itr < n1; itr++) {
446  int ifunc = b1_->function(ish) + itr;
447  int isym = b1_->irrep(ifunc);
448  int irel = b1_->function_within_irrep(ifunc);
449  int iabs = iirrepoff_[isym] + irel;
450  int isooff = itr;
451 
452  for (jtr = 0; jtr < n2; jtr++) {
453  int jfunc = b2_->function(jsh) + jtr;
454  int jsym = b2_->irrep(jfunc);
455  int jrel = b2_->function_within_irrep(jfunc);
456  int jabs = jirrepoff_[jsym] + jrel;
457  int jsooff = isooff * nso2 + jtr;
458 
459  for (ktr = 0; ktr < n3; ktr++) {
460  int kfunc = b3_->function(ksh) + ktr;
461  int ksym = b3_->irrep(kfunc);
462  int krel = b3_->function_within_irrep(kfunc);
463  int kabs = kirrepoff_[ksym] + krel;
464  int ksooff = jsooff * nso3 + ktr;
465 
466  for (ltr = 0; ltr < n4; ltr++) {
467  int lfunc = b4_->function(lsh) + ltr;
468  int lsym = b4_->irrep(lfunc);
469  int lrel = b4_->function_within_irrep(lfunc);
470  int labs = lirrepoff_[lsym] + lrel;
471  int lsooff = ksooff * nso4 + ltr;
472 
473  int iiabs = iabs;
474  int jjabs = jabs;
475  int kkabs = kabs;
476  int llabs = labs;
477 
478  int iiirrep = isym;
479  int jjirrep = jsym;
480  int kkirrep = ksym;
481  int llirrep = lsym;
482 
483  int iirel = irel;
484  int jjrel = jrel;
485  int kkrel = krel;
486  int llrel = lrel;
487 
488  if (std::fabs(buffer_[thread][lsooff]) > cutoff_) {
489  if (ish == jsh) {
490  if (iabs < jabs) continue;
491 
492  if (ksh == lsh) {
493  if (kabs < labs) continue;
494  if (INDEX2(iabs, jabs) < INDEX2(kabs, labs)) {
495  if (ish == ksh) // IIII case
496  continue;
497  else { // IIJJ case
498  SWAP_INDEX(ii, kk);
499  SWAP_INDEX(jj, ll);
500  }
501  }
502  } else { // IIJK case
503  if (labs > kabs) {
504  SWAP_INDEX(kk, ll);
505  }
506  if (INDEX2(iabs, jabs) < INDEX2(kabs, labs)) {
507  SWAP_INDEX(ii, kk);
508  SWAP_INDEX(jj, ll);
509  }
510  }
511  } else {
512  if (ksh == lsh) { // IJKK case
513  if (kabs < labs) continue;
514  if (iabs < jabs) {
515  SWAP_INDEX(ii, jj);
516  }
517  if (INDEX2(iabs, jabs) < INDEX2(kabs, labs)) {
518  SWAP_INDEX(ii, kk);
519  SWAP_INDEX(jj, ll);
520  }
521  } else { // IJIJ case
522  if (ish == ksh && jsh == lsh && INDEX2(iabs, jabs) < INDEX2(kabs, labs)) continue;
523  // IJKL case
524  if (iabs < jabs) {
525  SWAP_INDEX(ii, jj);
526  }
527  if (kabs < labs) {
528  SWAP_INDEX(kk, ll);
529  }
530  if (INDEX2(iabs, jabs) < INDEX2(kabs, labs)) {
531  SWAP_INDEX(ii, kk);
532  SWAP_INDEX(jj, ll);
533  }
534  }
535  }
536 
537  mints_timer_on("TwoBodySOInt::provide_IJKL functor");
538 
539  // func off/on
540  body(iiabs, jjabs, kkabs, llabs, iiirrep, iirel, jjirrep, jjrel, kkirrep, kkrel, llirrep, llrel,
541  buffer_[thread][lsooff]);
542 
543  mints_timer_off("TwoBodySOInt::provide_IJKL functor");
544  }
545  }
546  }
547  }
548  }
549  // timer_off("TwoBodySOInt::provide_IJKL overall");
550  mints_timer_off("TwoBodySOInt::provide_IJKL overall");
551 }
552 
553 template <typename TwoBodySOIntFunctor>
554 void TwoBodySOInt::compute_shell_deriv1(int uish, int ujsh, int uksh, int ulsh, TwoBodySOIntFunctor &body) {
555  int thread = 0;
556  // Old call: WorldComm->thread_id(pthread_self());
557 
558  const double *aobuffer = tb_[thread]->buffer();
559 
560  const SOTransform &t1 = b1_->sotrans(uish);
561  const SOTransform &t2 = b2_->sotrans(ujsh);
562  const SOTransform &t3 = b3_->sotrans(uksh);
563  const SOTransform &t4 = b4_->sotrans(ulsh);
564 
565  const int nso1 = b1_->nfunction(uish);
566  const int nso2 = b2_->nfunction(ujsh);
567  const int nso3 = b3_->nfunction(uksh);
568  const int nso4 = b4_->nfunction(ulsh);
569  const size_t nso = nso1 * nso2 * nso3 * nso4;
570 
571  // These are needed, because the buffers for each perturbation are spaced by
572  // quartets expressed in terms of Cartesian basis functions, regardless of puream.
573  const int nao1 = b1_->naofunction(uish);
574  const int nao2 = b2_->naofunction(ujsh);
575  const int nao3 = b3_->naofunction(uksh);
576  const int nao4 = b4_->naofunction(ulsh);
577  const size_t nao = nao1 * nao2 * nao3 * nao4;
578 
579  const int iatom = tb_[thread]->basis1()->shell(t1.aoshell[0].aoshell).ncenter();
580  const int jatom = tb_[thread]->basis2()->shell(t2.aoshell[0].aoshell).ncenter();
581  const int katom = tb_[thread]->basis3()->shell(t3.aoshell[0].aoshell).ncenter();
582  const int latom = tb_[thread]->basis4()->shell(t4.aoshell[0].aoshell).ncenter();
583 
584  // These 3 sections are not shell specific so we can just use petite1_
585  const unsigned short istablizer = petite1_->stablizer(iatom);
586  const unsigned short jstablizer = petite1_->stablizer(jatom);
587  const unsigned short kstablizer = petite1_->stablizer(katom);
588  const unsigned short lstablizer = petite1_->stablizer(latom);
589 
590  const int istabdense = dcd_->bits_to_dense_numbering(istablizer);
591  const int jstabdense = dcd_->bits_to_dense_numbering(jstablizer);
592  const int kstabdense = dcd_->bits_to_dense_numbering(kstablizer);
593  const int lstabdense = dcd_->bits_to_dense_numbering(lstablizer);
594 
595  const int ijstablizer = dcd_->intersection(istabdense, jstabdense);
596  const int klstablizer = dcd_->intersection(kstabdense, lstabdense);
597  const int ijklstablizer = dcd_->intersection(ijstablizer, klstablizer);
598 
599  const int *R_list = dcd_->dcr(istabdense, jstabdense);
600  const int *S_list = dcd_->dcr(kstabdense, lstabdense);
601  const int *T_list = dcd_->dcr(ijstablizer, klstablizer);
602 
603  const int R_size = R_list[0];
604  const int S_size = S_list[0];
605  const int T_size = T_list[0];
606 
607  int lambda_T = petite1_->nirrep() / dcd_->subgroup_dimensions(ijklstablizer);
608 
609  std::vector<int> sj_arr, sk_arr, sl_arr;
610 
611  int si = petite1_->unique_shell_map(uish, 0);
612  const int siatom = tb_[thread]->basis1()->shell(si).ncenter();
613 
614  for (int ij = 1; ij <= R_size; ++ij) {
615  int sj = petite2_->unique_shell_map(ujsh, R_list[ij]);
616  const int sjatom = tb_[thread]->basis2()->shell(sj).ncenter();
617 
618  for (int ijkl = 1; ijkl <= T_size; ++ijkl) {
619  int sk = petite3_->unique_shell_map(uksh, T_list[ijkl]);
620  int llsh = petite4_->unique_shell_map(ulsh, T_list[ijkl]);
621  const int skatom = tb_[thread]->basis3()->shell(sk).ncenter();
622 
623  for (int kl = 1; kl <= S_size; ++kl) {
624  int sl = petite4_->shell_map(llsh, S_list[kl]);
625  const int slatom = tb_[thread]->basis4()->shell(sl).ncenter();
626 
627  // Check AM
628  int total_am = tb_[thread]->basis1()->shell(si).am() + tb_[thread]->basis2()->shell(sj).am() +
629  tb_[thread]->basis3()->shell(sk).am() + tb_[thread]->basis4()->shell(sl).am();
630 
631  // if(!(total_am%2)||
632  // (BasisSet.shells[si].center!=BasisSet.shells[sj].center)||
633  // (BasisSet.shells[sj].center!=BasisSet.shells[sk].center)||
634  // (BasisSet.shells[sk].center!=BasisSet.shells[sl].center)) {
635 
636  // outfile->Printf( "total_am %d siatom %d sjatom %d skatom %d slatom %d\n",
637  // total_am, siatom, sjatom, skatom, slatom);
638  if (!(total_am % 2) || (siatom != sjatom) || (sjatom != skatom) || (skatom != slatom)) {
639  // outfile->Printf( "\tadding\n");
640  sj_arr.push_back(sj);
641  sk_arr.push_back(sk);
642  sl_arr.push_back(sl);
643  }
644  }
645  }
646  }
647 
648  // Obtain SALC transformation objects
649  // This probably won't work. I'll probably need
650  // the SALC for the petite list of shells.
651  // If so, this will move into the for loop below.
652 
653  // Compute integral using si, sj_arr, sk_arr, sl_arr
654  // Loop over unique quartets
655  const AOTransform &s1 = b1_->aotrans(si);
656  const unsigned short *ifuncpi = s1.nfuncpi;
657  const CdSalcWRTAtom &c1 = cdsalcs_->atom_salc(siatom);
658 
659  // Zero out SALC memory
660  for (size_t i = 0; i < cdsalcs_->ncd(); ++i) ::memset(deriv_[thread][i], 0, sizeof(double) * nso);
661 
662  double pfac = 1.0;
663  // if (uish == ujsh)
664  // pfac *= 0.5;
665  // if (uksh == ulsh)
666  // pfac *= 0.5;
667  // if (uish == uksh && ujsh == ulsh || uish == ulsh && ujsh == uksh)
668  // pfac *= 0.5;
669 
670  for (size_t n = 0; n < sj_arr.size(); ++n) {
671  int sj = sj_arr[n];
672  int sk = sk_arr[n];
673  int sl = sl_arr[n];
674 
675  const AOTransform &s2 = b2_->aotrans(sj);
676  const AOTransform &s3 = b3_->aotrans(sk);
677  const AOTransform &s4 = b4_->aotrans(sl);
678 
679  const unsigned short *jfuncpi = s2.nfuncpi;
680  const unsigned short *kfuncpi = s3.nfuncpi;
681  const unsigned short *lfuncpi = s4.nfuncpi;
682 
683  const CdSalcWRTAtom &c2 = cdsalcs_->atom_salc(tb_[thread]->basis2()->shell(sj).ncenter());
684  const CdSalcWRTAtom &c3 = cdsalcs_->atom_salc(tb_[thread]->basis3()->shell(sk).ncenter());
685  const CdSalcWRTAtom &c4 = cdsalcs_->atom_salc(tb_[thread]->basis4()->shell(sl).ncenter());
686 
687  // Compute this unique AO shell
688  tb_[thread]->compute_shell_deriv1(si, sj, sk, sl);
689 
690  int nirrep = b1_->nirrep();
691 
692  for (int isym = 0; isym < nirrep; ++isym) {
693  unsigned short nifunc = ifuncpi[isym];
694  for (int itr = 0; itr < nifunc; itr++) {
695  mints_timer_on("itr");
696 
697  const AOTransformFunction &ifunc = s1.soshellpi[isym][itr];
698  double icoef = ifunc.coef;
699  int iaofunc = ifunc.aofunc;
700  int isofunc = ifunc.sofunc;
701  int iaooff = iaofunc;
702  int isooff = isofunc;
703 
704  mints_timer_off("itr");
705 
706  for (int jsym = 0; jsym < nirrep; ++jsym) {
707  unsigned short njfunc = jfuncpi[jsym];
708  for (int jtr = 0; jtr < njfunc; jtr++) {
709  mints_timer_on("jtr");
710 
711  const AOTransformFunction &jfunc = s2.soshellpi[jsym][jtr];
712  double jcoef = jfunc.coef * icoef;
713  int jaofunc = jfunc.aofunc;
714  int jsofunc = jfunc.sofunc;
715  int jaooff = iaooff * nao2 + jaofunc;
716  int jsooff = isooff * nso2 + jsofunc;
717 
718  mints_timer_off("jtr");
719 
720  for (int ksym = 0; ksym < nirrep; ++ksym) {
721  unsigned short nkfunc = kfuncpi[ksym];
722  for (int ktr = 0; ktr < nkfunc; ktr++) {
723  mints_timer_on("ktr");
724 
725  const AOTransformFunction &kfunc = s3.soshellpi[ksym][ktr];
726  double kcoef = kfunc.coef * jcoef;
727  int kaofunc = kfunc.aofunc;
728  int ksofunc = kfunc.sofunc;
729  int kaooff = jaooff * nao3 + kaofunc;
730  int ksooff = jsooff * nso3 + ksofunc;
731 
732  mints_timer_off("ktr");
733 
734  int lsym = isym ^ jsym ^ ksym;
735  unsigned short nlfunc = lfuncpi[lsym];
736  for (int ltr = 0; ltr < nlfunc; ltr++) {
737  mints_timer_on("ltr");
738 
739  const AOTransformFunction &lfunc = s4.soshellpi[lsym][ltr];
740  double lcoef = lfunc.coef * kcoef;
741  int laofunc = lfunc.aofunc;
742  int lsofunc = lfunc.sofunc;
743  int laooff = kaooff * nao4 + laofunc;
744  int lsooff = ksooff * nso4 + lsofunc;
745 
746  int total_symmetry = ifunc.irrep ^ jfunc.irrep ^ kfunc.irrep ^ lfunc.irrep;
747 
748  // If we're only interested in totally symmetric derivatives, skip all others.
749  if (only_totally_symmetric_ == true && total_symmetry != 0) continue;
750 
751  // OK, the integral at 12 aobuff[laooff] needs
752  // to have lambda_T * lcoef * salcCoef applied to
753  // it and saved to salc index
754  // Special case is center B which must be determined
755  // from centers A, C, and D.
756  double A[3], B[3], C[3], D[3];
757 
758  A[0] = aobuffer[0 * nao + laooff];
759  A[1] = aobuffer[1 * nao + laooff];
760  A[2] = aobuffer[2 * nao + laooff];
761  C[0] = aobuffer[3 * nao + laooff];
762  C[1] = aobuffer[4 * nao + laooff];
763  C[2] = aobuffer[5 * nao + laooff];
764  D[0] = aobuffer[6 * nao + laooff];
765  D[1] = aobuffer[7 * nao + laooff];
766  D[2] = aobuffer[8 * nao + laooff];
767 
768  // Use translational invariance to determine B
769  B[0] = -(A[0] + C[0] + D[0]);
770  B[1] = -(A[1] + C[1] + D[1]);
771  B[2] = -(A[2] + C[2] + D[2]);
772 
773  A[0] *= lambda_T * pfac * lcoef;
774  A[1] *= lambda_T * pfac * lcoef;
775  A[2] *= lambda_T * pfac * lcoef;
776  B[0] *= lambda_T * pfac * lcoef;
777  B[1] *= lambda_T * pfac * lcoef;
778  B[2] *= lambda_T * pfac * lcoef;
779  C[0] *= lambda_T * pfac * lcoef;
780  C[1] *= lambda_T * pfac * lcoef;
781  C[2] *= lambda_T * pfac * lcoef;
782  D[0] *= lambda_T * pfac * lcoef;
783  D[1] *= lambda_T * pfac * lcoef;
784  D[2] *= lambda_T * pfac * lcoef;
785 
786  dprintf(
787  "so' derivatives: A[0] %+lf A[1] %+lf A[2] %+lf\n"
788  " B[0] %+lf B[1] %+lf B[2] %+lf\n"
789  " C[0] %+lf C[1] %+lf C[2] %+lf\n"
790  " D[0] %+lf D[1] %+lf D[2] %+lf\n"
791  "lsooff: %d\n"
792  "iirrep: %d jirrep: %d kirrep: %d lirrep: %d combined: %d\n",
793  A[0], A[1], A[2], B[0], B[1], B[2], C[0], C[1], C[2], D[0], D[1], D[2], lsooff,
794  ifunc.irrep, jfunc.irrep, kfunc.irrep, lfunc.irrep,
795  ifunc.irrep ^ jfunc.irrep ^ kfunc.irrep ^ lfunc.irrep);
796 
797  // For each center apply the so transform and salc coef.
798  // Ax
799  for (int nx = 0; nx < c1.nx(); ++nx) {
800  const CdSalcWRTAtom::Component element = c1.x(nx);
801  double temp = element.coef * A[0];
802  dprintf("Ax SALC#%d pfac %lf, A[0] %lf, contr %lf\n", element.salc,
803  element.coef, A[0], temp);
804  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
805  deriv_[thread][element.salc][lsooff] += temp;
806  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
807  }
808 
809  // Ay
810  for (int ny = 0; ny < c1.ny(); ++ny) {
811  const CdSalcWRTAtom::Component element = c1.y(ny);
812  double temp = element.coef * A[1];
813  dprintf("Ay SALC#%d pfac %lf, A[1] %lf, contr %lf\n", element.salc,
814  element.coef, A[1], temp);
815  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
816  deriv_[thread][element.salc][lsooff] += temp;
817  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
818  }
819 
820  // Az
821  for (int nz = 0; nz < c1.nz(); ++nz) {
822  const CdSalcWRTAtom::Component element = c1.z(nz);
823  double temp = element.coef * A[2];
824  dprintf("Az SALC#%d pfac %lf, A[2] %lf, contr %lf\n", element.salc,
825  element.coef, A[2], temp);
826  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
827  deriv_[thread][element.salc][lsooff] += temp;
828  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
829  }
830 
831  // Bx
832  for (int nx = 0; nx < c2.nx(); ++nx) {
833  const CdSalcWRTAtom::Component element = c2.x(nx);
834  double temp = element.coef * B[0];
835  dprintf("Bx SALC#%d pfac %lf, B[0] %lf, contr %lf\n", element.salc,
836  element.coef, B[0], temp);
837  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
838  deriv_[thread][element.salc][lsooff] += temp;
839  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
840  }
841 
842  // By
843  for (int ny = 0; ny < c2.ny(); ++ny) {
844  const CdSalcWRTAtom::Component element = c2.y(ny);
845  double temp = element.coef * B[1];
846  dprintf("By SALC#%d pfac %lf, B[1] %lf, contr %lf\n", element.salc,
847  element.coef, B[1], temp);
848  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
849  deriv_[thread][element.salc][lsooff] += temp;
850  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
851  }
852 
853  // Bz
854  for (int nz = 0; nz < c2.nz(); ++nz) {
855  const CdSalcWRTAtom::Component element = c2.z(nz);
856  double temp = element.coef * B[2];
857  dprintf("Bz SALC#%d pfac %lf, B[2] %lf, contr %lf\n", element.salc,
858  element.coef, B[2], temp);
859  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
860  deriv_[thread][element.salc][lsooff] += temp;
861  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
862  }
863 
864  // Cx
865  for (int nx = 0; nx < c3.nx(); ++nx) {
866  const CdSalcWRTAtom::Component element = c3.x(nx);
867  double temp = element.coef * C[0];
868  dprintf("Cx SALC#%d pfac %lf, C[0] %lf, contr %lf\n", element.salc,
869  element.coef, C[0], temp);
870  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
871  deriv_[thread][element.salc][lsooff] += temp;
872  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
873  }
874 
875  // Cy
876  for (int ny = 0; ny < c3.ny(); ++ny) {
877  const CdSalcWRTAtom::Component element = c3.y(ny);
878  double temp = element.coef * C[1];
879  dprintf("Cy SALC#%d pfac %lf, C[1] %lf, contr %lf\n", element.salc,
880  element.coef, C[1], temp);
881  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
882  deriv_[thread][element.salc][lsooff] += temp;
883  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
884  }
885 
886  // Cz
887  for (int nz = 0; nz < c3.nz(); ++nz) {
888  const CdSalcWRTAtom::Component element = c3.z(nz);
889  double temp = element.coef * C[2];
890  dprintf("Cz SALC#%d pfac %lf, C[2] %lf, contr %lf\n", element.salc,
891  element.coef, C[2], temp);
892  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
893  deriv_[thread][element.salc][lsooff] += temp;
894  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
895  }
896 
897  // Dx
898  for (int nx = 0; nx < c4.nx(); ++nx) {
899  const CdSalcWRTAtom::Component element = c4.x(nx);
900  double temp = element.coef * D[0];
901  dprintf("Dx SALC#%d pfac %lf, D[0] %lf, contr %lf\n", element.salc,
902  element.coef, D[0], temp);
903  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
904  deriv_[thread][element.salc][lsooff] += temp;
905  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
906  }
907 
908  // Dy
909  for (int ny = 0; ny < c4.ny(); ++ny) {
910  const CdSalcWRTAtom::Component element = c4.y(ny);
911  double temp = element.coef * D[1];
912  dprintf("Dy SALC#%d pfac %lf, D[1] %lf, contr %lf\n", element.salc,
913  element.coef, D[1], temp);
914  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
915  deriv_[thread][element.salc][lsooff] += temp;
916  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
917  }
918 
919  // Dz
920  for (int nz = 0; nz < c4.nz(); ++nz) {
921  const CdSalcWRTAtom::Component element = c4.z(nz);
922  double temp = element.coef * D[2];
923  dprintf("Dz SALC#%d pfac %lf, D[2] %lf, contr %lf\n", element.salc,
924  element.coef, D[2], temp);
925  if (total_symmetry == element.irrep && std::fabs(temp) > cutoff_)
926  deriv_[thread][element.salc][lsooff] += temp;
927  dprintf(" val: %lf\n", deriv_[thread][element.salc][lsooff]);
928  }
929 
930  dprintf("\n");
931  }
932  }
933  }
934  }
935  }
936  }
937  }
938  }
939 
940  provide_IJKL_deriv1(uish, ujsh, uksh, ulsh, body);
941 
942 } // function
943 
944 template <typename TwoBodySOIntFunctor>
945 void TwoBodySOInt::provide_IJKL_deriv1(int ish, int jsh, int ksh, int lsh, TwoBodySOIntFunctor &body) {
946  int thread = 0;
947  // Old call: WorldComm->thread_id(pthread_self());
948 
949  mints_timer_on("TwoBodySOInt::provide_IJKL overall");
950 
951  int nso2 = b2_->nfunction(jsh);
952  int nso3 = b3_->nfunction(ksh);
953  int nso4 = b4_->nfunction(lsh);
954 
955  int n1 = b1_->nfunction(ish);
956  int n2 = b2_->nfunction(jsh);
957  int n3 = b3_->nfunction(ksh);
958  int n4 = b4_->nfunction(lsh);
959 
960  int itr;
961  int jtr;
962  int ktr;
963  int ltr;
964 
965  for (itr = 0; itr < n1; itr++) {
966  int ifunc = b1_->function(ish) + itr;
967  int isym = b1_->irrep(ifunc);
968  int irel = b1_->function_within_irrep(ifunc);
969  int iabs = iirrepoff_[isym] + irel;
970  int isooff = itr;
971 
972  for (jtr = 0; jtr < n2; jtr++) {
973  int jfunc = b2_->function(jsh) + jtr;
974  int jsym = b2_->irrep(jfunc);
975  int jrel = b2_->function_within_irrep(jfunc);
976  int jabs = jirrepoff_[jsym] + jrel;
977  int jsooff = isooff * nso2 + jtr;
978 
979  for (ktr = 0; ktr < n3; ktr++) {
980  int kfunc = b3_->function(ksh) + ktr;
981  int ksym = b3_->irrep(kfunc);
982  int krel = b3_->function_within_irrep(kfunc);
983  int kabs = kirrepoff_[ksym] + krel;
984  int ksooff = jsooff * nso3 + ktr;
985 
986  for (ltr = 0; ltr < n4; ltr++) {
987  int lfunc = b4_->function(lsh) + ltr;
988  int lsym = b4_->irrep(lfunc);
989  int lrel = b4_->function_within_irrep(lfunc);
990  int labs = lirrepoff_[lsym] + lrel;
991  int lsooff = ksooff * nso4 + ltr;
992 
993  // Only totally symmetric pertubations are considered here!
994  if (isym ^ jsym ^ ksym ^ lsym) continue;
995 
996  int iiabs = iabs;
997  int jjabs = jabs;
998  int kkabs = kabs;
999  int llabs = labs;
1000 
1001  int iiirrep = isym;
1002  int jjirrep = jsym;
1003  int kkirrep = ksym;
1004  int llirrep = lsym;
1005 
1006  int iirel = irel;
1007  int jjrel = jrel;
1008  int kkrel = krel;
1009  int llrel = lrel;
1010 
1011  if (ish == jsh) {
1012  if (iabs < jabs) continue;
1013 
1014  if (ksh == lsh) {
1015  if (kabs < labs) continue;
1016  if (INDEX2(iabs, jabs) < INDEX2(kabs, labs)) {
1017  if (ish == ksh) // IIII case
1018  continue;
1019  else { // IIJJ case
1020  SWAP_INDEX(ii, kk);
1021  SWAP_INDEX(jj, ll);
1022  }
1023  }
1024  } else { // IIJK case
1025  if (labs > kabs) {
1026  SWAP_INDEX(kk, ll);
1027  }
1028  if (INDEX2(iabs, jabs) < INDEX2(kabs, labs)) {
1029  SWAP_INDEX(ii, kk);
1030  SWAP_INDEX(jj, ll);
1031  }
1032  }
1033  } else {
1034  if (ksh == lsh) { // IJKK case
1035  if (kabs < labs) continue;
1036  if (iabs < jabs) {
1037  SWAP_INDEX(ii, jj);
1038  }
1039  if (INDEX2(iabs, jabs) < INDEX2(kabs, labs)) {
1040  SWAP_INDEX(ii, kk);
1041  SWAP_INDEX(jj, ll);
1042  }
1043  } else { // IJIJ case
1044  if (ish == ksh && jsh == lsh && INDEX2(iabs, jabs) < INDEX2(kabs, labs)) continue;
1045  // IJKL case
1046  if (iabs < jabs) {
1047  SWAP_INDEX(ii, jj);
1048  }
1049  if (kabs < labs) {
1050  SWAP_INDEX(kk, ll);
1051  }
1052  if (INDEX2(iabs, jabs) < INDEX2(kabs, labs)) {
1053  SWAP_INDEX(ii, kk);
1054  SWAP_INDEX(jj, ll);
1055  }
1056  }
1057  }
1058 
1059  mints_timer_on("TwoBodySOInt::provide_IJKL functor");
1060  for (size_t i = 0; i < cdsalcs_->ncd(); ++i) {
1061  if (std::fabs(deriv_[thread][i][lsooff]) > cutoff_)
1062  body(i, iiabs, jjabs, kkabs, llabs, iiirrep, iirel, jjirrep, jjrel, kkirrep, kkrel, llirrep,
1063  llrel, deriv_[thread][i][lsooff]);
1064  }
1065  body.next_tpdm_element();
1066 
1067  mints_timer_off("TwoBodySOInt::provide_IJKL functor");
1068  }
1069  }
1070  }
1071  }
1072  mints_timer_off("TwoBodySOInt::provide_IJKL overall");
1073 }
1074 
1075 template <typename TwoBodySOIntFunctor>
1076 void TwoBodySOInt::compute_integrals(TwoBodySOIntFunctor &functor) {
1077  if (comm_ == "MADNESS") {
1078  throw PSIEXCEPTION(
1079  "PSI4 was not built with MADNESS. "
1080  "Please rebuild PSI4 with MADNESS, or "
1081  "change your COMMUNICATOR "
1082  "environment variable to MPI or LOCAL.\n");
1083  } else {
1084  std::shared_ptr<SOShellCombinationsIterator> shellIter =
1085  std::make_shared<SOShellCombinationsIterator>(b1_, b2_, b3_, b4_);
1086  this->compute_quartets(shellIter, functor);
1087  }
1088 }
1089 
1090 template <typename TwoBodySOIntFunctor>
1091 void TwoBodySOInt::compute_integrals_deriv1(TwoBodySOIntFunctor &functor) {
1093  throw PSIEXCEPTION(
1094  "The way the TPDM is stored and iterated enables only totally symmetric"
1095  " perturbations to be considered right now!");
1096 
1097  if (comm_ == "MADNESS") {
1098  } else {
1099  auto PQIter = std::make_shared<SO_PQ_Iterator>(b1_);
1100  size_t pair_number = 0;
1101  for (PQIter->first(); PQIter->is_done() == false; PQIter->next()) {
1102  compute_pq_pair_deriv1<TwoBodySOIntFunctor>(PQIter->p(), PQIter->q(), pair_number, functor);
1103  pair_number++;
1104  }
1105  }
1106 }
1107 
1108 typedef std::shared_ptr<OneBodySOInt> SharedOneBodySOInt;
1109 typedef std::shared_ptr<TwoBodySOInt> SharedTwoBodySOInt;
1110 }
1111 
1112 #endif // _psi_src_lib_libmints_sointegral_h_
void compute_shell_deriv1(const SOShellCombinationsIterator &shellIter, TwoBodySOIntFunctor &body)
Definition: sointegral_twobody.h:180
std::shared_ptr< SOBasisSet > basis2() const
Definition: sointegral.cc:459
std::shared_ptr< SOBasisSet > basis4() const
Definition: sointegral.cc:463
SOTransformShell * aoshell
The SOTransformShell object for each AO.
Definition: sobasis.h:98
Definition: sointegral_twobody.h:77
std::shared_ptr< BasisSet > basis
Definition: dx_write.cc:59
std::shared_ptr< PetiteList > petite1_
Definition: sointegral_twobody.h:95
int iirrepoff_[8]
Definition: sointegral_twobody.h:93
Definition: cdsalclist.h:124
int s() const
Definition: integral.h:309
int nz() const
Definition: cdsalclist.h:115
std::shared_ptr< IntegralFactory > integral_
Definition: sointegral_twobody.h:80
int lirrepoff_[8]
Definition: sointegral_twobody.h:93
int compute_pq_pair_deriv1(const int &p, const int &q, const size_t &pair_number, const TwoBodySOIntFunctor &body)
Definition: sointegral_twobody.h:192
int r() const
Definition: integral.h:308
int irrep
Definition: sobasis.h:64
std::shared_ptr< SOBasisSet > basis3() const
Definition: sointegral.cc:461
std::vector< double ** > deriv_
Definition: sointegral_twobody.h:91
Definition: pointgrp.h:106
void compute_integrals(TwoBodySOIntFunctor &functor)
Definition: sointegral_twobody.h:1076
Definition: sobasis.h:59
int kirrepoff_[8]
Definition: sointegral_twobody.h:93
int nao
Definition: dx_write.cc:56
bool only_totally_symmetric_
Definition: sointegral_twobody.h:104
#define mints_timer_off(a)
Definition: typedefs.h:41
int nx() const
Definition: cdsalclist.h:113
std::shared_ptr< OneBodySOInt > SharedOneBodySOInt
Definition: sointegral_twobody.h:1108
void compute_shell(const SOShellCombinationsIterator &shellIter, TwoBodySOIntFunctor &body)
Definition: sointegral_twobody.h:150
int irrep
Definition: cdsalclist.h:87
std::shared_ptr< PetiteList > petite3_
Definition: sointegral_twobody.h:97
#define mints_timer_on(a)
Definition: typedefs.h:40
unsigned short nfuncpi[8]
Definition: sobasis.h:111
void compute_integrals_deriv1(TwoBodySOIntFunctor &functor)
Definition: sointegral_twobody.h:1091
int aoshell
The number of the AO shell from which these functions come.
Definition: sobasis.h:77
Definition: cdsalclist.h:84
const Component & x(int n) const
Definition: cdsalclist.h:117
bool only_totally_symmetric() const
Definition: sointegral_twobody.h:135
#define INDEX2(i, j)
Definition: PK_workers.h:37
const CdSalcWRTAtom & atom_salc(int i) const
Definition: cdsalclist.h:178
std::shared_ptr< SOBasisSet > b4_
Definition: sointegral_twobody.h:85
#define SWAP_INDEX(a, b)
Definition: typedefs.h:66
Definition: sobasis.h:107
#define dprintf(...)
Definition: sointegral_twobody.h:62
std::vector< std::shared_ptr< TwoBodyAOInt > > tb_
Definition: sointegral_twobody.h:79
Definition: cdsalclist.h:82
std::shared_ptr< SOBasisSet > b3_
Definition: sointegral_twobody.h:84
int nproc_
Definition: sointegral_twobody.h:107
std::string comm_
Definition: sointegral_twobody.h:106
std::vector< double * > buffer_
Definition: sointegral_twobody.h:88
int q() const
Definition: integral.h:307
int aofunc
Definition: sobasis.h:62
int me_
Definition: sointegral_twobody.h:108
Header file for the Quantum Trio LibraryDavid Sherrill 1994.
std::shared_ptr< PetiteList > petite4_
Definition: sointegral_twobody.h:98
double cutoff_
Definition: sointegral_twobody.h:110
common_init()
int nso
Definition: dx_write.cc:56
std::shared_ptr< SOBasisSet > b1_
Definition: sointegral_twobody.h:82
int jirrepoff_[8]
Definition: sointegral_twobody.h:93
void compute_quartets(ShellIter &shellIter, TwoBodySOIntFunctor &body)
Definition: sointegral_twobody.h:160
void compute_quartets_deriv1(ShellIter &shellIter, TwoBodySOIntFunctor &body)
Definition: sointegral_twobody.h:173
#define PSI_API
Definition: pragma.h:128
std::shared_ptr< SOBasisSet > b2_
Definition: sointegral_twobody.h:83
int nthread_
Definition: sointegral_twobody.h:105
void set_only_totally_symmetric(bool ots)
Definition: sointegral_twobody.h:136
const Component & y(int n) const
Definition: cdsalclist.h:118
const double * buffer(int thread=0) const
Definition: sointegral_twobody.h:144
std::vector< double * > temp_
Definition: sointegral_twobody.h:89
const CdSalcList * cdsalcs_
Definition: sointegral_twobody.h:112
void set_cutoff(double ints_tolerance)
Definition: sointegral_twobody.h:146
std::shared_ptr< DCD > dcd_
Definition: sointegral_twobody.h:102
std::shared_ptr< PointGroup > pg_
Definition: sointegral_twobody.h:100
int ny() const
Definition: cdsalclist.h:114
size_t ncd() const
Definition: cdsalclist.h:165
int sofunc
Definition: sobasis.h:63
void provide_IJKL(int, int, int, int, TwoBodySOIntFunctor &body)
Definition: sointegral_twobody.h:421
double coef
Definition: sobasis.h:61
int salc
Definition: cdsalclist.h:88
int p() const
Definition: integral.h:306
std::vector< double * > temp2_
Definition: sointegral_twobody.h:90
const Component & z(int n) const
Definition: cdsalclist.h:119
#define PSIEXCEPTION(message)
Definition: exception.h:50
size_t size_
Definition: sointegral_twobody.h:87
double coef
Definition: cdsalclist.h:86
std::shared_ptr< TwoBodySOInt > SharedTwoBodySOInt
Definition: sointegral_twobody.h:1109
std::vector< AOTransformFunction > soshellpi[8]
Definition: sobasis.h:110
void provide_IJKL_deriv1(int ish, int jsh, int ksh, int lsh, TwoBodySOIntFunctor &body)
Definition: sointegral_twobody.h:945
std::shared_ptr< PetiteList > petite2_
Definition: sointegral_twobody.h:96
Definition: sobasis.h:92
Definition: integral.h:271