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 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 
57 {
58 public:
59  PCMPotentialInt(std::vector<SphericalTransform>&, std::shared_ptr<BasisSet>, std::shared_ptr<BasisSet>, 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 {
68  // Do not worry about zeroing out result
69  int ns1 = bs1_->nshell();
70  int ns2 = bs2_->nshell();
71  int bf1_offset = 0;
72  for (int i=0; i<ns1; ++i) {
73  const GaussianShell& s1 = bs1_->shell(i);
74  int ni = s1.ncartesian();
75  int bf2_offset = 0;
76  for (int j=0; j<ns2; ++j) {
77  const GaussianShell& s2 = bs2_->shell(j);
78  int nj = s2.ncartesian();
79  // Compute the shell
80 
81  int ao12;
82  int am1 = s1.am();
83  int am2 = s2.am();
84  int nprim1 = s1.nprimitive();
85  int nprim2 = s2.nprimitive();
86  double A[3], B[3];
87  A[0] = s1.center()[0];
88  A[1] = s1.center()[1];
89  A[2] = s1.center()[2];
90  B[0] = s2.center()[0];
91  B[1] = s2.center()[1];
92  B[2] = s2.center()[2];
93 
94  int izm = 1;
95  int iym = am1 + 1;
96  int ixm = iym * iym;
97  int jzm = 1;
98  int jym = am2 + 1;
99  int jxm = jym * jym;
100 
101  // compute intermediates
102  double AB2 = 0.0;
103  AB2 += (A[0] - B[0]) * (A[0] - B[0]);
104  AB2 += (A[1] - B[1]) * (A[1] - B[1]);
105  AB2 += (A[2] - B[2]) * (A[2] - B[2]);
106 
107 
108  double ***vi = potential_recur_->vi();
109 
110  double** Zxyzp = Zxyz_->pointer();
111  int ncharge = Zxyz_->rowspi()[0];
112 
113  for (int atom=0; atom<ncharge; ++atom) {
114  memset(buffer_, 0, s1.ncartesian() * s2.ncartesian() * sizeof(double));
115  double PC[3];
116 
117  double Z = Zxyzp[atom][0];
118 
119  double C[3];
120  C[0] = Zxyzp[atom][1];
121  C[1] = Zxyzp[atom][2];
122  C[2] = Zxyzp[atom][3];
123  for (int p1=0; p1<nprim1; ++p1) {
124  double a1 = s1.exp(p1);
125  double c1 = s1.coef(p1);
126  for (int p2=0; p2<nprim2; ++p2) {
127  double a2 = s2.exp(p2);
128  double c2 = s2.coef(p2);
129  double gamma = a1 + a2;
130  double oog = 1.0/gamma;
131 
132  double PA[3], PB[3], P[3];
133  P[0] = (a1*A[0] + a2*B[0])*oog;
134  P[1] = (a1*A[1] + a2*B[1])*oog;
135  P[2] = (a1*A[2] + a2*B[2])*oog;
136  PA[0] = P[0] - A[0];
137  PA[1] = P[1] - A[1];
138  PA[2] = P[2] - A[2];
139  PB[0] = P[0] - B[0];
140  PB[1] = P[1] - B[1];
141  PB[2] = P[2] - B[2];
142  PC[0] = P[0] - C[0];
143  PC[1] = P[1] - C[1];
144  PC[2] = P[2] - C[2];
145 
146  double over_pf = exp(-a1*a2*AB2*oog) * sqrt(M_PI*oog) * M_PI * oog * c1 * c2;
147 
148 
149  // Do recursion
150  potential_recur_->compute(PA, PB, PC, gamma, am1, am2);
151 
152  ao12 = 0;
153  for(int ii = 0; ii <= am1; ii++) {
154  int l1 = am1 - ii;
155  for(int jj = 0; jj <= ii; jj++) {
156  int m1 = ii - jj;
157  int n1 = jj;
158  /*--- create all am components of sj ---*/
159  for(int kk = 0; kk <= am2; kk++) {
160  int l2 = am2 - kk;
161  for(int ll = 0; ll <= kk; ll++) {
162  int m2 = kk - ll;
163  int n2 = ll;
164 
165  // Compute location in the recursion and store the value
166  int iind = l1 * ixm + m1 * iym + n1 * izm;
167  int jind = l2 * jxm + m2 * jym + n2 * jzm;
168  buffer_[ao12++] += -vi[iind][jind][0] * over_pf * Z;
169  }
170  }
171  }
172  }
173  } // End loop over primitives of shell 2
174  } // End loop over primitives of shell 1
175  ao12 = 0;
176  int ao1 = 0;
177  for(int ii = 0; ii <= am1; ii++) {
178  for(int jj = 0; jj <= ii; jj++) {
179  /*--- create all am components of sj ---*/
180  int ao2 = 0;
181  for(int kk = 0; kk <= am2; kk++) {
182  for(int ll = 0; ll <= kk; ll++) {
183  // Compute location in the recursion
184  double val = buffer_[ao12++];
185  // Hand the work off to the functor
186  functor(ao1+bf1_offset, ao2+bf2_offset, atom, val);
187  ao2++;
188  }
189  }
190  ao1++;
191  }
192  }
193  } // End loop over points
194  bf2_offset += nj;
195  } // End loop over shell 2
196  bf1_offset += ni;
197  } // End loop over shell 1
198 }
199 
201 {
202  public:
206  void operator()(int bf1, int bf2, int center, double integral)
207  {
208  outfile->Printf( "bf1: %3d bf2 %3d center (%5d) integral %16.10f\n", bf1, bf2, center, integral);
209  }
210 };
211 
212 
214 {
219  protected:
221  double **pD_;
223  double *charges_;
224  public:
225  ContractOverDensityFunctor(size_t /*ncenters*/, double *charges, SharedMatrix D):
226  pD_(D->pointer()),
227  charges_(charges)
228  {
229  }
230  void operator()(int bf1, int bf2, int center, double integral)
231  {
232  charges_[center] += pD_[bf1][bf2] * integral;
233  }
234 };
235 
236 
238 {
243  protected:
245  double **pF_;
247  const double *charges_;
248  public:
249  ContractOverChargesFunctor(const double* charges, SharedMatrix F):
250  pF_(F->pointer()),
251  charges_(charges)
252  {
253  if(F->rowdim() != F->coldim())
254  throw PSIEXCEPTION("Invalid Fock matrix in ContractOverCharges");
255  int nbf = F->rowdim();
256  ::memset(pF_[0], 0, nbf*nbf*sizeof(double));
257  }
258 
259  void operator()(int bf1, int bf2, int center, double integral)
260  {
261  pF_[bf1][bf2] += integral * charges_[center];
262  }
263 };
264 
265 } //Namespace
266 #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:61
const double * center() const
Returns the center of the Molecule this shell is on.
Definition: gshell.cc:263
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:259
Definition: pointgrp.h:106
int am() const
The angular momentum of the given contraction.
Definition: gshell.h:294
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:981
#define M_PI
Definition: include/psi4/physconst.h:87
ContractOverChargesFunctor(const double *charges, SharedMatrix F)
Definition: potentialint.h:249
Definition: potentialint.h:213
double * charges_
The array of charges.
Definition: potentialint.h:223
const double * charges_
The array of charges.
Definition: potentialint.h:247
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:230
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:59
Definition: potentialint.h:200
Gaussian orbital shell.
Definition: gshell.h:189
double ** pF_
Pointer to the matrix that will contribute to the 2e part of the Fock matrix.
Definition: potentialint.h:245
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:206
std::shared_ptr< BasisSet > bs1_
Definition: onebody.h:58
double exp(int prim) const
Returns the exponent of the given primitive.
Definition: gshell.h:312
std::shared_ptr< PsiOutStream > outfile
Definition: core.cc:102
int deriv() const
What order of derivative was requested?
Definition: onebody.h:130
double ** pD_
Pointer to the density matrix.
Definition: potentialint.h:221
ContractOverDensityFunctor(size_t, double *charges, SharedMatrix D)
Definition: potentialint.h:225
Definition: potentialint.h:237
double * buffer_
Definition: onebody.h:64
int nprimitive() const
The number of primitive Gaussians.
Definition: gshell.cc:226
double coef(int pi) const
Return coefficient of pi&#39;th primitive.
Definition: gshell.h:314
double *** vi() const
Returns the potential integral 3D matrix.
Definition: osrecur.h:128
SharedMatrix Zxyz_
Matrix of coordinates/charges of partial charges.
Definition: potential.h:64
#define PSIEXCEPTION(message)
Definition: exception.h:48
int ncartesian() const
Total number of functions if this shell was Cartesian.
Definition: gshell.h:292
Definition: potentialint.h:56