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