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