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 
33 #include "psi4/libmints/basisset.h"
34 #include "psi4/libmints/matrix.h"
35 #include "psi4/libmints/osrecur.h"
36 
37 namespace psi{
38 
39 class GaussianShell;
40 class SphericalTransform;
41 
54 {
55 public:
56  PCMPotentialInt(std::vector<SphericalTransform>&, std::shared_ptr<BasisSet>, std::shared_ptr<BasisSet>, int deriv=0);
58  template<typename PCMPotentialIntFunctor>
59  void compute(PCMPotentialIntFunctor &functor);
60 };
61 
62 template<typename PCMPotentialIntFunctor>
63 void PCMPotentialInt::compute(PCMPotentialIntFunctor &functor)
64 {
65  // Do not worry about zeroing out result
66  int ns1 = bs1_->nshell();
67  int ns2 = bs2_->nshell();
68  int bf1_offset = 0;
69  for (int i=0; i<ns1; ++i) {
70  const GaussianShell& s1 = bs1_->shell(i);
71  int ni = s1.ncartesian();
72  int bf2_offset = 0;
73  for (int j=0; j<ns2; ++j) {
74  const GaussianShell& s2 = bs2_->shell(j);
75  int nj = s2.ncartesian();
76  // Compute the shell
77 
78  int ao12;
79  int am1 = s1.am();
80  int am2 = s2.am();
81  int nprim1 = s1.nprimitive();
82  int nprim2 = s2.nprimitive();
83  double A[3], B[3];
84  A[0] = s1.center()[0];
85  A[1] = s1.center()[1];
86  A[2] = s1.center()[2];
87  B[0] = s2.center()[0];
88  B[1] = s2.center()[1];
89  B[2] = s2.center()[2];
90 
91  int izm = 1;
92  int iym = am1 + 1;
93  int ixm = iym * iym;
94  int jzm = 1;
95  int jym = am2 + 1;
96  int jxm = jym * jym;
97 
98  // compute intermediates
99  double AB2 = 0.0;
100  AB2 += (A[0] - B[0]) * (A[0] - B[0]);
101  AB2 += (A[1] - B[1]) * (A[1] - B[1]);
102  AB2 += (A[2] - B[2]) * (A[2] - B[2]);
103 
104 
105  double ***vi = potential_recur_->vi();
106 
107  double** Zxyzp = Zxyz_->pointer();
108  int ncharge = Zxyz_->rowspi()[0];
109 
110  for (int atom=0; atom<ncharge; ++atom) {
111  memset(buffer_, 0, s1.ncartesian() * s2.ncartesian() * sizeof(double));
112  double PC[3];
113 
114  double Z = Zxyzp[atom][0];
115 
116  double C[3];
117  C[0] = Zxyzp[atom][1];
118  C[1] = Zxyzp[atom][2];
119  C[2] = Zxyzp[atom][3];
120  for (int p1=0; p1<nprim1; ++p1) {
121  double a1 = s1.exp(p1);
122  double c1 = s1.coef(p1);
123  for (int p2=0; p2<nprim2; ++p2) {
124  double a2 = s2.exp(p2);
125  double c2 = s2.coef(p2);
126  double gamma = a1 + a2;
127  double oog = 1.0/gamma;
128 
129  double PA[3], PB[3], P[3];
130  P[0] = (a1*A[0] + a2*B[0])*oog;
131  P[1] = (a1*A[1] + a2*B[1])*oog;
132  P[2] = (a1*A[2] + a2*B[2])*oog;
133  PA[0] = P[0] - A[0];
134  PA[1] = P[1] - A[1];
135  PA[2] = P[2] - A[2];
136  PB[0] = P[0] - B[0];
137  PB[1] = P[1] - B[1];
138  PB[2] = P[2] - B[2];
139  PC[0] = P[0] - C[0];
140  PC[1] = P[1] - C[1];
141  PC[2] = P[2] - C[2];
142 
143  double over_pf = exp(-a1*a2*AB2*oog) * sqrt(M_PI*oog) * M_PI * oog * c1 * c2;
144 
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 {
199  public:
203  void operator()(int bf1, int bf2, int center, double integral)
204  {
205  outfile->Printf( "bf1: %3d bf2 %3d center (%5d) integral %16.10f\n", bf1, bf2, center, integral);
206  }
207 };
208 
209 
211 {
216  protected:
218  double **pD_;
220  double *charges_;
221  public:
222  ContractOverDensityFunctor(size_t /*ncenters*/, double *charges, SharedMatrix D):
223  pD_(D->pointer()),
224  charges_(charges)
225  {
226  }
227  void operator()(int bf1, int bf2, int center, double integral)
228  {
229  charges_[center] += pD_[bf1][bf2] * integral;
230  }
231 };
232 
233 
235 {
240  protected:
242  double **pF_;
244  const double *charges_;
245  public:
246  ContractOverChargesFunctor(const double* charges, SharedMatrix F):
247  pF_(F->pointer()),
248  charges_(charges)
249  {
250  if(F->rowdim() != F->coldim())
251  throw PSIEXCEPTION("Invalid Fock matrix in ContractOverCharges");
252  int nbf = F->rowdim();
253  ::memset(pF_[0], 0, nbf*nbf*sizeof(double));
254  }
255 
256  void operator()(int bf1, int bf2, int center, double integral)
257  {
258  pF_[bf1][bf2] += integral * charges_[center];
259  }
260 };
261 
262 } //Namespace
263 #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:289
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:256
Definition: pointgrp.h:106
int am() const
The angular momentum of the given contraction.
Definition: gshell.h:329
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:246
Definition: potentialint.h:210
double * charges_
The array of charges.
Definition: potentialint.h:220
const double * charges_
The array of charges.
Definition: potentialint.h:244
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:227
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:197
Gaussian orbital shell.
Definition: gshell.h:222
double ** pF_
Pointer to the matrix that will contribute to the 2e part of the Fock matrix.
Definition: potentialint.h:242
void compute(PCMPotentialIntFunctor &functor)
Drives the loops over all shell pairs, to compute integrals.
Definition: potentialint.h:63
std::shared_ptr< Matrix > SharedMatrix
Definition: adc.h:49
void operator()(int bf1, int bf2, int center, double integral)
Definition: potentialint.h:203
std::shared_ptr< BasisSet > bs1_
Definition: onebody.h:58
double exp(int prim) const
Returns the exponent of the given primitive.
Definition: gshell.h:347
std::shared_ptr< PsiOutStream > outfile
Definition: core.cc:100
int deriv() const
What order of derivative was requested?
Definition: onebody.h:130
double ** pD_
Pointer to the density matrix.
Definition: potentialint.h:218
ContractOverDensityFunctor(size_t, double *charges, SharedMatrix D)
Definition: potentialint.h:222
Definition: potentialint.h:234
double * buffer_
Definition: onebody.h:64
Definition: PsiFileImpl.h:39
int nprimitive() const
The number of primitive Gaussians.
Definition: gshell.cc:252
double coef(int pi) const
Return coefficient of pi&#39;th primitive.
Definition: gshell.h:349
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:327
Definition: potentialint.h:53