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