Psi4
potentialint.h
Go to the documentation of this file.
1 /*
2  * @BEGIN LICENSE
3  *
4  * Psi4: an open-source quantum chemistry software package
5  *
6  * Copyright (c) 2007-2018 The Psi4 Developers.
7  *
8  * The copyrights for code used from other parties are included in
9  * the corresponding files.
10  *
11  * This file is part of Psi4.
12  *
13  * Psi4 is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU Lesser General Public License as published by
15  * the Free Software Foundation, version 3.
16  *
17  * Psi4 is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU Lesser General Public License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public License along
23  * with Psi4; if not, write to the Free Software Foundation, Inc.,
24  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25  *
26  * @END LICENSE
27  */
28 
29 #ifndef _psi_src_lib_libmints_potentialint_h_
30 #define _psi_src_lib_libmints_potentialint_h_
31 
32 #include "psi4/libmints/gshell.h"
33 #include "psi4/libmints/basisset.h"
35 #include "psi4/libmints/matrix.h"
36 #include "psi4/libmints/osrecur.h"
38 
39 namespace psi {
40 
41 class GaussianShell;
42 class SphericalTransform;
43 class BasisSet;
44 
56 class PCMPotentialInt : public PotentialInt {
57  public:
58  PCMPotentialInt(std::vector<SphericalTransform> &, std::shared_ptr<BasisSet>, std::shared_ptr<BasisSet>,
59  int deriv = 0);
61  template <typename PCMPotentialIntFunctor>
62  void compute(PCMPotentialIntFunctor &functor);
63 };
64 
65 template <typename PCMPotentialIntFunctor>
66 void PCMPotentialInt::compute(PCMPotentialIntFunctor &functor) {
67  // Do not worry about zeroing out result
68  int ns1 = bs1_->nshell();
69  int ns2 = bs2_->nshell();
70  int bf1_offset = 0;
71  for (int i = 0; i < ns1; ++i) {
72  const GaussianShell &s1 = bs1_->shell(i);
73  int ni = s1.ncartesian();
74  int bf2_offset = 0;
75  for (int j = 0; j < ns2; ++j) {
76  const GaussianShell &s2 = bs2_->shell(j);
77  int nj = s2.ncartesian();
78  // Compute the shell
79 
80  int ao12;
81  int am1 = s1.am();
82  int am2 = s2.am();
83  int nprim1 = s1.nprimitive();
84  int nprim2 = s2.nprimitive();
85  double A[3], B[3];
86  A[0] = s1.center()[0];
87  A[1] = s1.center()[1];
88  A[2] = s1.center()[2];
89  B[0] = s2.center()[0];
90  B[1] = s2.center()[1];
91  B[2] = s2.center()[2];
92 
93  int izm = 1;
94  int iym = am1 + 1;
95  int ixm = iym * iym;
96  int jzm = 1;
97  int jym = am2 + 1;
98  int jxm = jym * jym;
99 
100  // compute intermediates
101  double AB2 = 0.0;
102  AB2 += (A[0] - B[0]) * (A[0] - B[0]);
103  AB2 += (A[1] - B[1]) * (A[1] - B[1]);
104  AB2 += (A[2] - B[2]) * (A[2] - B[2]);
105 
106  double ***vi = potential_recur_->vi();
107 
108  double **Zxyzp = Zxyz_->pointer();
109  int ncharge = Zxyz_->rowspi()[0];
110 
111  for (int atom = 0; atom < ncharge; ++atom) {
112  memset(buffer_, 0, s1.ncartesian() * s2.ncartesian() * sizeof(double));
113  double PC[3];
114 
115  double Z = Zxyzp[atom][0];
116 
117  double C[3];
118  C[0] = Zxyzp[atom][1];
119  C[1] = Zxyzp[atom][2];
120  C[2] = Zxyzp[atom][3];
121  for (int p1 = 0; p1 < nprim1; ++p1) {
122  double a1 = s1.exp(p1);
123  double c1 = s1.coef(p1);
124  for (int p2 = 0; p2 < nprim2; ++p2) {
125  double a2 = s2.exp(p2);
126  double c2 = s2.coef(p2);
127  double gamma = a1 + a2;
128  double oog = 1.0 / gamma;
129 
130  double PA[3], PB[3], P[3];
131  P[0] = (a1 * A[0] + a2 * B[0]) * oog;
132  P[1] = (a1 * A[1] + a2 * B[1]) * oog;
133  P[2] = (a1 * A[2] + a2 * B[2]) * oog;
134  PA[0] = P[0] - A[0];
135  PA[1] = P[1] - A[1];
136  PA[2] = P[2] - A[2];
137  PB[0] = P[0] - B[0];
138  PB[1] = P[1] - B[1];
139  PB[2] = P[2] - B[2];
140  PC[0] = P[0] - C[0];
141  PC[1] = P[1] - C[1];
142  PC[2] = P[2] - C[2];
143 
144  double over_pf = exp(-a1 * a2 * AB2 * oog) * sqrt(M_PI * oog) * M_PI * oog * c1 * c2;
145 
146  // Do recursion
147  potential_recur_->compute(PA, PB, PC, gamma, am1, am2);
148 
149  ao12 = 0;
150  for (int ii = 0; ii <= am1; ii++) {
151  int l1 = am1 - ii;
152  for (int jj = 0; jj <= ii; jj++) {
153  int m1 = ii - jj;
154  int n1 = jj;
155  /*--- create all am components of sj ---*/
156  for (int kk = 0; kk <= am2; kk++) {
157  int l2 = am2 - kk;
158  for (int ll = 0; ll <= kk; ll++) {
159  int m2 = kk - ll;
160  int n2 = ll;
161 
162  // Compute location in the recursion and store the value
163  int iind = l1 * ixm + m1 * iym + n1 * izm;
164  int jind = l2 * jxm + m2 * jym + n2 * jzm;
165  buffer_[ao12++] += -vi[iind][jind][0] * over_pf * Z;
166  }
167  }
168  }
169  }
170  } // End loop over primitives of shell 2
171  } // End loop over primitives of shell 1
172  ao12 = 0;
173  int ao1 = 0;
174  for (int ii = 0; ii <= am1; ii++) {
175  for (int jj = 0; jj <= ii; jj++) {
176  /*--- create all am components of sj ---*/
177  int ao2 = 0;
178  for (int kk = 0; kk <= am2; kk++) {
179  for (int ll = 0; ll <= kk; ll++) {
180  // Compute location in the recursion
181  double val = buffer_[ao12++];
182  // Hand the work off to the functor
183  functor(ao1 + bf1_offset, ao2 + bf2_offset, atom, val);
184  ao2++;
185  }
186  }
187  ao1++;
188  }
189  }
190  } // End loop over points
191  bf2_offset += nj;
192  } // End loop over shell 2
193  bf1_offset += ni;
194  } // End loop over shell 1
195 }
196 
198  public:
202  void operator()(int bf1, int bf2, int center, double integral) {
203  outfile->Printf("bf1: %3d bf2 %3d center (%5d) integral %16.10f\n", bf1, bf2, center, integral);
204  }
205 };
206 
212  protected:
214  double **pD_;
216  double *charges_;
217 
218  public:
219  ContractOverDensityFunctor(size_t /*ncenters*/, double *charges, SharedMatrix D)
220  : pD_(D->pointer()), charges_(charges) {}
221  void operator()(int bf1, int bf2, int center, double integral) { charges_[center] += pD_[bf1][bf2] * integral; }
222 };
223 
229  protected:
231  double **pF_;
233  const double *charges_;
234 
235  public:
236  ContractOverChargesFunctor(const double *charges, SharedMatrix F) : pF_(F->pointer()), charges_(charges) {
237  if (F->rowdim() != F->coldim()) throw PSIEXCEPTION("Invalid Fock matrix in ContractOverCharges");
238  int nbf = F->rowdim();
239  ::memset(pF_[0], 0, nbf * nbf * sizeof(double));
240  }
241 
242  void operator()(int bf1, int bf2, int center, double integral) { pF_[bf1][bf2] += integral * charges_[center]; }
243 };
244 
245 } // namespace psi
246 #endif
Computes potential integrals. Use an IntegralFactory to create this object.
Definition: potential.h:50
ObaraSaikaTwoCenterVIRecursion * potential_recur_
Recursion object that does the heavy lifting.
Definition: potential.h:60
const double * center() const
Returns the center of the Molecule this shell is on.
Definition: gshell.cc:236
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:242
Definition: pointgrp.h:104
int am() const
The angular momentum of the given contraction.
Definition: gshell.h:261
virtual void compute(double PA[3], double PB[3], double PC[3], double zeta, int am1, int am2)
Computes the potential integral 3D matrix using the data provided.
Definition: osrecur.cc:1387
#define M_PI
Definition: include/psi4/physconst.h:428
ContractOverChargesFunctor(const double *charges, SharedMatrix F)
Definition: potentialint.h:236
std::shared_ptr< PsiOutStream > outfile
Definition: core.cc:102
Definition: potentialint.h:207
double * charges_
The array of charges.
Definition: potentialint.h:216
const double * charges_
The array of charges.
Definition: potentialint.h:233
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:221
PCMPotentialInt(std::vector< SphericalTransform > &, std::shared_ptr< BasisSet >, std::shared_ptr< BasisSet >, int deriv=0)
Definition: potentialint.cc:33
std::shared_ptr< BasisSet > bs2_
Definition: onebody.h:58
Definition: potentialint.h:197
Gaussian orbital shell.
Definition: gshell.h:168
double ** pF_
Pointer to the matrix that will contribute to the 2e part of the Fock matrix.
Definition: potentialint.h:231
void compute(PCMPotentialIntFunctor &functor)
Drives the loops over all shell pairs, to compute integrals.
Definition: potentialint.h:66
std::shared_ptr< Matrix > SharedMatrix
Definition: adc.h:49
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:202
std::shared_ptr< BasisSet > bs1_
Definition: onebody.h:57
double exp(int prim) const
Returns the exponent of the given primitive.
Definition: gshell.h:279
int deriv() const
What order of derivative was requested?
Definition: onebody.h:130
double ** pD_
Pointer to the density matrix.
Definition: potentialint.h:214
ContractOverDensityFunctor(size_t, double *charges, SharedMatrix D)
Definition: potentialint.h:219
Definition: potentialint.h:224
double * buffer_
Definition: onebody.h:63
int nprimitive() const
The number of primitive Gaussians.
Definition: gshell.cc:206
double coef(int pi) const
Return coefficient of pi&#39;th primitive.
Definition: gshell.h:281
double *** vi() const
Returns the potential integral 3D matrix.
Definition: osrecur.h:125
SharedMatrix Zxyz_
Matrix of coordinates/charges of partial charges.
Definition: potential.h:63
#define PSIEXCEPTION(message)
Definition: exception.h:50
int ncartesian() const
Total number of functions if this shell was Cartesian.
Definition: gshell.h:259
Definition: potentialint.h:56