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-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_potentialint_h_
29 #define _psi_src_lib_libmints_potentialint_h_
30 
32 #include "psi4/libmints/basisset.h"
33 #include "psi4/libmints/matrix.h"
34 #include "psi4/libmints/osrecur.h"
35 
36 namespace psi{
37 
38 class GaussianShell;
39 class OneBodyAOInt;
40 class IntegralFactory;
41 class SphericalTransform;
42 class Vector3;
43 
56 {
57 public:
58  PCMPotentialInt(std::vector<SphericalTransform>&, std::shared_ptr<BasisSet>, std::shared_ptr<BasisSet>, int deriv=0);
60  template<typename PCMPotentialIntFunctor>
61  void compute(PCMPotentialIntFunctor &functor);
62 };
63 
64 template<typename PCMPotentialIntFunctor>
65 void PCMPotentialInt::compute(PCMPotentialIntFunctor &functor)
66 {
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 
107  double ***vi = potential_recur_->vi();
108 
109  double** Zxyzp = Zxyz_->pointer();
110  int ncharge = Zxyz_->rowspi()[0];
111 
112  for (int atom=0; atom<ncharge; ++atom) {
113  memset(buffer_, 0, s1.ncartesian() * s2.ncartesian() * sizeof(double));
114  double PC[3];
115 
116  double Z = Zxyzp[atom][0];
117 
118  double C[3];
119  C[0] = Zxyzp[atom][1];
120  C[1] = Zxyzp[atom][2];
121  C[2] = Zxyzp[atom][3];
122  for (int p1=0; p1<nprim1; ++p1) {
123  double a1 = s1.exp(p1);
124  double c1 = s1.coef(p1);
125  for (int p2=0; p2<nprim2; ++p2) {
126  double a2 = s2.exp(p2);
127  double c2 = s2.coef(p2);
128  double gamma = a1 + a2;
129  double oog = 1.0/gamma;
130 
131  double PA[3], PB[3], P[3];
132  P[0] = (a1*A[0] + a2*B[0])*oog;
133  P[1] = (a1*A[1] + a2*B[1])*oog;
134  P[2] = (a1*A[2] + a2*B[2])*oog;
135  PA[0] = P[0] - A[0];
136  PA[1] = P[1] - A[1];
137  PA[2] = P[2] - A[2];
138  PB[0] = P[0] - B[0];
139  PB[1] = P[1] - B[1];
140  PB[2] = P[2] - B[2];
141  PC[0] = P[0] - C[0];
142  PC[1] = P[1] - C[1];
143  PC[2] = P[2] - C[2];
144 
145  double over_pf = exp(-a1*a2*AB2*oog) * sqrt(M_PI*oog) * M_PI * oog * c1 * c2;
146 
147 
148  // Do recursion
149  potential_recur_->compute(PA, PB, PC, gamma, am1, am2);
150 
151  ao12 = 0;
152  for(int ii = 0; ii <= am1; ii++) {
153  int l1 = am1 - ii;
154  for(int jj = 0; jj <= ii; jj++) {
155  int m1 = ii - jj;
156  int n1 = jj;
157  /*--- create all am components of sj ---*/
158  for(int kk = 0; kk <= am2; kk++) {
159  int l2 = am2 - kk;
160  for(int ll = 0; ll <= kk; ll++) {
161  int m2 = kk - ll;
162  int n2 = ll;
163 
164  // Compute location in the recursion and store the value
165  int iind = l1 * ixm + m1 * iym + n1 * izm;
166  int jind = l2 * jxm + m2 * jym + n2 * jzm;
167  buffer_[ao12++] += -vi[iind][jind][0] * over_pf * Z;
168  }
169  }
170  }
171  }
172  } // End loop over primitives of shell 2
173  } // End loop over primitives of shell 1
174  ao12 = 0;
175  int ao1 = 0;
176  for(int ii = 0; ii <= am1; ii++) {
177  for(int jj = 0; jj <= ii; jj++) {
178  /*--- create all am components of sj ---*/
179  int ao2 = 0;
180  for(int kk = 0; kk <= am2; kk++) {
181  for(int ll = 0; ll <= kk; ll++) {
182  // Compute location in the recursion
183  double val = buffer_[ao12++];
184  // Hand the work off to the functor
185  functor(ao1+bf1_offset, ao2+bf2_offset, atom, val);
186  ao2++;
187  }
188  }
189  ao1++;
190  }
191  }
192  } // End loop over points
193  bf2_offset += nj;
194  } // End loop over shell 2
195  bf1_offset += ni;
196  } // End loop over shell 1
197 }
198 
200 {
201  public:
205  void operator()(int bf1, int bf2, int center, double integral)
206  {
207  outfile->Printf( "bf1: %3d bf2 %3d center (%5d) integral %16.10f\n", bf1, bf2, center, integral);
208  }
209 };
210 
211 
213 {
218  protected:
220  double **pD_;
222  double *charges_;
223  public:
224  ContractOverDensityFunctor(size_t /*ncenters*/, double *charges, SharedMatrix D):
225  pD_(D->pointer()),
226  charges_(charges)
227  {
228  }
229  void operator()(int bf1, int bf2, int center, double integral)
230  {
231  charges_[center] += pD_[bf1][bf2] * integral;
232  }
233 };
234 
235 
237 {
242  protected:
244  double **pF_;
246  const double *charges_;
247  public:
248  ContractOverChargesFunctor(const double* charges, SharedMatrix F):
249  pF_(F->pointer()),
250  charges_(charges)
251  {
252  if(F->rowdim() != F->coldim())
253  throw PSIEXCEPTION("Invalid Fock matrix in ContractOverCharges");
254  int nbf = F->rowdim();
255  ::memset(pF_[0], 0, nbf*nbf*sizeof(double));
256  }
257 
258  void operator()(int bf1, int bf2, int center, double integral)
259  {
260  pF_[bf1][bf2] += integral * charges_[center];
261  }
262 };
263 
264 } //Namespace
265 #endif
Computes potential integrals. Use an IntegralFactory to create this object.
Definition: potential.h:49
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:235
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:258
Definition: pointgrp.h:105
int am() const
The angular momentum of the given contraction.
Definition: gshell.h:269
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:980
#define M_PI
Definition: include/psi4/physconst.h:86
ContractOverChargesFunctor(const double *charges, SharedMatrix F)
Definition: potentialint.h:248
Definition: potentialint.h:212
double * charges_
The array of charges.
Definition: potentialint.h:222
const double * charges_
The array of charges.
Definition: potentialint.h:246
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:229
PCMPotentialInt(std::vector< SphericalTransform > &, std::shared_ptr< BasisSet >, std::shared_ptr< BasisSet >, int deriv=0)
Definition: potentialint.cc:32
std::shared_ptr< BasisSet > bs2_
Definition: onebody.h:59
Definition: potentialint.h:199
Gaussian orbital shell.
Definition: gshell.h:191
double ** pF_
Pointer to the matrix that will contribute to the 2e part of the Fock matrix.
Definition: potentialint.h:244
void compute(PCMPotentialIntFunctor &functor)
Drives the loops over all shell pairs, to compute integrals.
Definition: potentialint.h:65
std::shared_ptr< Matrix > SharedMatrix
Definition: adc.h:50
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:205
std::shared_ptr< BasisSet > bs1_
Definition: onebody.h:58
double exp(int prim) const
Returns the exponent of the given primitive.
Definition: gshell.h:287
std::shared_ptr< PsiOutStream > outfile
Definition: core.cc:99
int deriv() const
What order of derivative was requested?
Definition: onebody.h:130
double ** pD_
Pointer to the density matrix.
Definition: potentialint.h:220
ContractOverDensityFunctor(size_t, double *charges, SharedMatrix D)
Definition: potentialint.h:224
Definition: potentialint.h:236
double * buffer_
Definition: onebody.h:64
Definition: PsiFileImpl.h:38
int nprimitive() const
The number of primitive Gaussians.
Definition: gshell.cc:215
double coef(int pi) const
Return coefficient of pi&#39;th primitive.
Definition: gshell.h:289
double *** vi() const
Returns the potential integral 3D matrix.
Definition: osrecur.h:127
SharedMatrix Zxyz_
Matrix of coordinates/charges of partial charges.
Definition: potential.h:63
#define PSIEXCEPTION(message)
Definition: exception.h:47
int ncartesian() const
Total number of functions if this shell was Cartesian.
Definition: gshell.h:267
Definition: potentialint.h:55