Psi4
ecpint.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) 2016-2017 Robert A. Shaw.
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 /* class ECPInt calculates one-body ECP integrals
30  class AngularIntegral calculates and stores the angular integrals needed for the ECP integration
31  class RadialIntegral abstracts the calculation of the radial integrals needed for the ECP integration,
32  such that if a different approach was desired later, this could be done with minimal alterations to ECPInt.
33 
34  Robert A. Shaw 2016
35 
36  REFERENCES:
37  (Flores06) R. Flores-Moreno et al., J. Comput. Chem. 27 (2006), 1009-1019
38  (MM81) L. E. McMurchie and E. R. Davidson, J. Comp. Phys. 44 (1981), 289 - 301
39 */
40 
41 #ifndef ECPINT_HEAD
42 #define ECPINT_HEAD
43 
44 #include <vector>
45 #include "psi4/libmints/multiarr.h"
47 #include "psi4/libmints/typedefs.h"
48 #include "psi4/libmints/onebody.h"
50 #include "psi4/libmints/bessel.h"
51 
52 namespace psi {
53 
54 class BasisSet;
55 class GaussianShell;
56 class IntegralFactory;
57 class SphericalTransform;
58 
66 static TwoIndex<double> realSphericalHarmonics(int lmax, double x, double phi);
67 
73 struct ShellPairData {
74  int LA, LB, maxLBasis;
75  int ncartA, ncartB;
76  double A[3], B[3];
77  double A2, Am, B2, Bm, RAB2, RABm;
78 };
79 
89  private:
91  int LB, LE;
93  int wDim, maxL;
94 
99 
101  double calcG(int l, int m) const;
102  double calcH1(int i, int j, int l, int m) const;
103  double calcH2(int i, int j, int k, int m) const;
104 
111  ThreeIndex<double> uklm(int lam, int mu) const;
117  ThreeIndex<double> Pijk(int maxI) const;
118 
128  void makeW(FiveIndex<double> &U);
134 
135  public:
137  AngularIntegral();
139  AngularIntegral(int LB, int LE);
145  void init(int LB, int LE);
149  void compute();
150 
152  void clear();
153 
163  double getIntegral(int k, int l, int m, int lam, int mu) const;
175  double getIntegral(int k, int l, int m, int lam, int mu, int rho, int sigma) const;
176 
178  bool isZero(int k, int l, int m, int lam, int mu, double tolerance) const;
180  bool isZero(int k, int l, int m, int lam, int mu, int rho, int sigma, double tolerance) const;
181 };
182 
192  private:
199 
202 
204  double tolerance;
205 
207  static double integrand(double r, double *p, int ix);
208 
217  void buildBessel(std::vector<double> &r, int nr, int maxL, TwoIndex<double> &values, double weight = 1.0);
218 
219  double calcKij(double Na, double Nb, double zeta_a, double zeta_b, double R2) const;
220 
229  void buildU(const GaussianShell &U, int l, int N, GCQuadrature &grid, double *Utab);
230 
242  void buildF(const GaussianShell &shell, double A, int lstart, int lend, std::vector<double> &r, int nr, int start,
243  int end, TwoIndex<double> &F);
244 
255  int integrate(int maxL, int gridSize, TwoIndex<double> &intValues, GCQuadrature &grid, std::vector<double> &values,
256  int offset = 0, int skip = 1);
257 
258  public:
260  RadialIntegral();
261 
271  void init(int maxL, double tol = 1e-15, int small = 256, int large = 1024);
272 
280  void buildParameters(const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data);
281 
294  void type1(int maxL, int N, int offset, const GaussianShell &U, const GaussianShell &shellA,
295  const GaussianShell &shellB, ShellPairData &data, TwoIndex<double> &values);
296 
312  void type2(int lam, int l1start, int l1end, int l2start, int l2end, int N, const GaussianShell &U,
313  const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data, TwoIndex<double> &values);
314 };
315 
324 class ECPInt : public OneBodyAOInt {
325  private:
330 
332  double calcC(int a, int m, double A) const;
333  void makeC(FiveIndex<double> &C, int L, double *A);
334 
336  void type1(const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data,
339  void type2(int l, const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB,
341 
343  void compute_pair(const GaussianShell &shellA, const GaussianShell &shellB);
344 
346  void compute_shell_pair(const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB,
347  TwoIndex<double> &values, int shiftA = 0, int shiftB = 0);
348 
349  public:
355  ECPInt(std::vector<SphericalTransform> &, std::shared_ptr<BasisSet>, std::shared_ptr<BasisSet>, int deriv = 0);
356  virtual ~ECPInt();
357 };
358 
359 class ECPSOInt : public OneBodySOInt {
360  int natom_;
361 
362  public:
363  ECPSOInt(const std::shared_ptr<OneBodyAOInt> &, const std::shared_ptr<IntegralFactory> &);
364  ECPSOInt(const std::shared_ptr<OneBodyAOInt> &, const IntegralFactory *);
365 };
366 
367 } // namespace psi
368 #endif
Definition: onebody.h:55
static double integrand(double r, double *p, int ix)
This integrand simply returns the pretabulated integrand values stored in p given an index ix...
Definition: ecpint.cc:369
void makeW(FiveIndex< double > &U)
Definition: ecpint.cc:219
int LE
Definition: ecpint.h:91
RadialIntegral radInts
The interface to the radial integral calculation.
Definition: ecpint.h:327
int * U
Definition: stringlist.cc:66
AngularIntegral()
Default constructor creates empty object.
Definition: ecpint.cc:306
FiveIndex< double > makeU()
Definition: ecpint.cc:200
double RAB2
Definition: ecpint.h:77
double A2
Definition: ecpint.h:77
double calcC(int a, int m, double A) const
Worker functions for calculating binomial expansion coefficients.
Definition: ecpint.cc:658
static TwoIndex< double > realSphericalHarmonics(int lmax, double x, double phi)
Definition: ecpint.cc:53
Performs adaptive Gauss-Chebyshev quadrature for any given function.
Definition: gaussquad.h:62
Stores the (shifted) angular momenta, number of cartesians in a shell pair, and shifted centers...
Definition: ecpint.h:73
Definition: pointgrp.h:104
FiveIndex< double > W
Stores the type 1 angular integrals.
Definition: ecpint.h:96
double calcKij(double Na, double Nb, double zeta_a, double zeta_b, double R2) const
Definition: ecpint.cc:363
int ncartB
Definition: ecpint.h:75
void init(int maxL, double tol=1e-15, int small=256, int large=1024)
Definition: ecpint.cc:345
ThreeIndex< double > uklm(int lam, int mu) const
Definition: ecpint.cc:141
void type2(int lam, int l1start, int l1end, int l2start, int l2end, int N, const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data, TwoIndex< double > &values)
Definition: ecpint.cc:534
void compute()
Definition: ecpint.cc:315
Computes a modified spherical Bessel function of the first kind.
Definition: bessel.h:57
double getIntegral(int k, int l, int m, int lam, int mu) const
Definition: ecpint.cc:323
double calcH1(int i, int j, int l, int m) const
Definition: ecpint.cc:121
TwoIndex< double > K
Definition: ecpint.h:201
void type1(int maxL, int N, int offset, const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data, TwoIndex< double > &values)
Definition: ecpint.cc:426
GCQuadrature smallGrid
The smaller integration grid, default for the type 2 integrals.
Definition: ecpint.h:196
int maxL
Definition: ecpint.h:93
int ncartA
Definition: ecpint.h:75
RadialIntegral()
Default constructor creates an empty object.
Definition: ecpint.cc:343
double A[3]
Definition: ecpint.h:76
void compute_pair(const GaussianShell &shellA, const GaussianShell &shellB)
Overridden shell-pair integral calculation over all ECP centers.
Definition: ecpint.cc:915
int maxLBasis
Definition: ecpint.h:74
ECPInt(std::vector< SphericalTransform > &, std::shared_ptr< BasisSet >, std::shared_ptr< BasisSet >, int deriv=0)
Definition: ecpint.cc:639
SevenIndex< double > omega
Stores the type 2 angular integrals.
Definition: ecpint.h:98
Definition: sointegral_onebody.h:40
Definition: ecpint.h:324
double RABm
Definition: ecpint.h:77
ECPSOInt(const std::shared_ptr< OneBodyAOInt > &, const std::shared_ptr< IntegralFactory > &)
Definition: ecpint.cc:932
ThreeIndex< double > Pijk(int maxI) const
Definition: ecpint.cc:181
double Bm
Definition: ecpint.h:77
int integrate(int maxL, int gridSize, TwoIndex< double > &intValues, GCQuadrature &grid, std::vector< double > &values, int offset=0, int skip=1)
Definition: ecpint.cc:409
double tolerance
Tolerance for change below which an integral is considered converged.
Definition: ecpint.h:204
void makeOmega(FiveIndex< double > &U)
Definition: ecpint.cc:265
Definition: ecpint.h:359
int LB
Definition: ecpint.h:74
int wDim
Limits for the w-integral indices, and angular momentum indices.
Definition: ecpint.h:93
void type1(const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data, FiveIndex< double > &CA, FiveIndex< double > &CB, TwoIndex< double > &values)
Calculates the type 1 integrals for the given ECP center over the given shell pair.
Definition: ecpint.cc:686
Gaussian orbital shell.
Definition: gshell.h:168
Abstracts the calculation of radial integrals for ECP integration.
Definition: ecpint.h:191
double Am
Definition: ecpint.h:77
void init(int LB, int LE)
Definition: ecpint.cc:308
bool isZero(int k, int l, int m, int lam, int mu, double tolerance) const
is W(k, l, m, lam, mu) zero to within a given tolerance?
Definition: ecpint.cc:328
Calculates and stores the angular integrals needed for ECP integration.
Definition: ecpint.h:88
void makeC(FiveIndex< double > &C, int L, double *A)
Definition: ecpint.cc:665
int LB
Maximum angular momentum of orbital basis and ECP basis, respectively.
Definition: ecpint.h:91
int natom_
Definition: ecpint.h:360
Definition: integral.h:374
TwoIndex< double > P2
Definition: ecpint.h:201
virtual ~ECPInt()
Definition: ecpint.cc:656
void buildParameters(const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data)
Definition: ecpint.cc:371
AngularIntegral angInts
The angular integrals, which can be reused over all ECP centers.
Definition: ecpint.h:329
void compute_shell_pair(const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB, TwoIndex< double > &values, int shiftA=0, int shiftB=0)
Computes the overall ECP integrals over the given ECP center and shell pair.
Definition: ecpint.cc:860
TwoIndex< double > p
Matrices of parameters needed in both type 1 and 2 integrations.
Definition: ecpint.h:201
void type2(int l, const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data, FiveIndex< double > &CA, FiveIndex< double > &CB, ThreeIndex< double > &values)
Calculates the type 2 integrals for the given ECP center over the given shell pair.
Definition: ecpint.cc:769
int deriv() const
What order of derivative was requested?
Definition: onebody.h:130
TwoIndex< double > P
Definition: ecpint.h:201
int LA
Definition: ecpint.h:74
GCQuadrature bigGrid
The larger integration grid for type 1 integrals, and for when the smaller grid fails for type 2 inte...
Definition: ecpint.h:194
double calcH2(int i, int j, int k, int m) const
Definition: ecpint.cc:130
BesselFunction bessie
Modified spherical Bessel function of the first kind.
Definition: ecpint.h:198
void buildBessel(std::vector< double > &r, int nr, int maxL, TwoIndex< double > &values, double weight=1.0)
Definition: ecpint.cc:355
double B[3]
Definition: ecpint.h:76
void clear()
TODO: Clears the W and omega arrays.
Definition: ecpint.cc:321
void buildF(const GaussianShell &shell, double A, int lstart, int lend, std::vector< double > &r, int nr, int start, int end, TwoIndex< double > &F)
Definition: ecpint.cc:510
double calcG(int l, int m) const
Worker functions for calculating terms in the USP to spherical transformation coefficients.
Definition: ecpint.cc:111
Templated skeleton three index array for convenience.
Definition: multiarr.h:72
void buildU(const GaussianShell &U, int l, int N, GCQuadrature &grid, double *Utab)
Definition: ecpint.cc:398
double B2
Definition: ecpint.h:77