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 
46 #include "psi4/pragma.h"
47 #include "psi4/libmints/multiarr.h"
49 #include "psi4/libmints/typedefs.h"
50 #include "psi4/libmints/onebody.h"
52 #include "psi4/libmints/bessel.h"
53 
54 namespace psi {
55 
56 class BasisSet;
57 class GaussianShell;
58 class IntegralFactory;
59 class SphericalTransform;
60 
68 static TwoIndex<double> realSphericalHarmonics(int lmax, double x, double phi);
69 
75 struct ShellPairData {
76  int LA, LB, maxLBasis;
77  int ncartA, ncartB;
78  double A[3], B[3];
79  double A2, Am, B2, Bm, RAB2, RABm;
80 };
81 
91  private:
93  int LB, LE;
95  int wDim, maxL;
96 
101 
103  double calcG(int l, int m) const;
104  double calcH1(int i, int j, int l, int m) const;
105  double calcH2(int i, int j, int k, int m) const;
106 
113  ThreeIndex<double> uklm(int lam, int mu) const;
119  ThreeIndex<double> Pijk(int maxI) const;
120 
130  void makeW(FiveIndex<double> &U);
136 
137  public:
139  AngularIntegral();
141  AngularIntegral(int LB, int LE);
147  void init(int LB, int LE);
151  void compute();
152 
154  void clear();
155 
165  double getIntegral(int k, int l, int m, int lam, int mu) const;
177  double getIntegral(int k, int l, int m, int lam, int mu, int rho, int sigma) const;
178 
180  bool isZero(int k, int l, int m, int lam, int mu, double tolerance) const;
182  bool isZero(int k, int l, int m, int lam, int mu, int rho, int sigma, double tolerance) const;
183 };
184 
194  private:
201 
204 
206  double tolerance;
207 
209  static double integrand(double r, double *p, int ix);
210 
219  void buildBessel(std::vector<double> &r, int nr, int maxL, TwoIndex<double> &values, double weight = 1.0);
220 
221  double calcKij(double Na, double Nb, double zeta_a, double zeta_b, double R2) const;
222 
231  void buildU(const GaussianShell &U, int l, int N, GCQuadrature &grid, double *Utab);
232 
244  void buildF(const GaussianShell &shell, double A, int lstart, int lend, std::vector<double> &r, int nr, int start,
245  int end, TwoIndex<double> &F);
246 
257  int integrate(int maxL, int gridSize, TwoIndex<double> &intValues, GCQuadrature &grid, std::vector<double> &values,
258  int offset = 0, int skip = 1);
259 
260  public:
262  RadialIntegral();
263 
273  void init(int maxL, double tol = 1e-15, int small = 256, int large = 1024);
274 
282  void buildParameters(const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data);
283 
296  void type1(int maxL, int N, int offset, const GaussianShell &U, const GaussianShell &shellA,
297  const GaussianShell &shellB, ShellPairData &data, TwoIndex<double> &values);
298 
314  void type2(int lam, int l1start, int l1end, int l2start, int l2end, int N, const GaussianShell &U,
315  const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data, TwoIndex<double> &values);
316 };
317 
326 class ECPInt : public OneBodyAOInt {
327  private:
332 
334  double calcC(int a, int m, double A) const;
335  void makeC(FiveIndex<double> &C, int L, double *A);
336 
338  void type1(const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB, ShellPairData &data,
341  void type2(int l, const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB,
343 
345  void compute_pair(const GaussianShell &shellA, const GaussianShell &shellB) override;
346 
348  void compute_shell_pair(const GaussianShell &U, const GaussianShell &shellA, const GaussianShell &shellB,
349  TwoIndex<double> &values, int shiftA = 0, int shiftB = 0);
350 
351  public:
357  ECPInt(std::vector<SphericalTransform> &, std::shared_ptr<BasisSet>, std::shared_ptr<BasisSet>, int deriv = 0);
358  ~ECPInt() override;
359 };
360 
361 class ECPSOInt : public OneBodySOInt {
362  int natom_;
363 
364  public:
365  ECPSOInt(const std::shared_ptr<OneBodyAOInt> &, const std::shared_ptr<IntegralFactory> &);
366  ECPSOInt(const std::shared_ptr<OneBodyAOInt> &, const IntegralFactory *);
367 };
368 
369 } // namespace psi
370 #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:93
RadialIntegral radInts
The interface to the radial integral calculation.
Definition: ecpint.h:329
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:79
double A2
Definition: ecpint.h:79
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:64
Stores the (shifted) angular momenta, number of cartesians in a shell pair, and shifted centers...
Definition: ecpint.h:75
Definition: pointgrp.h:104
FiveIndex< double > W
Stores the type 1 angular integrals.
Definition: ecpint.h:98
double calcKij(double Na, double Nb, double zeta_a, double zeta_b, double R2) const
Definition: ecpint.cc:363
int ncartB
Definition: ecpint.h:77
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:59
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:203
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:198
int maxL
Definition: ecpint.h:95
int ncartA
Definition: ecpint.h:77
RadialIntegral()
Default constructor creates an empty object.
Definition: ecpint.cc:343
double A[3]
Definition: ecpint.h:78
int maxLBasis
Definition: ecpint.h:76
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:100
Definition: sointegral_onebody.h:40
~ECPInt() override
Definition: ecpint.cc:656
Definition: ecpint.h:326
double RABm
Definition: ecpint.h:79
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:79
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:206
void makeOmega(FiveIndex< double > &U)
Definition: ecpint.cc:265
Definition: ecpint.h:361
int LB
Definition: ecpint.h:76
int wDim
Limits for the w-integral indices, and angular momentum indices.
Definition: ecpint.h:95
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:193
double Am
Definition: ecpint.h:79
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:90
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:93
int natom_
Definition: ecpint.h:362
Definition: integral.h:374
TwoIndex< double > P2
Definition: ecpint.h:203
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:331
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:203
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:203
int LA
Definition: ecpint.h:76
GCQuadrature bigGrid
The larger integration grid for type 1 integrals, and for when the smaller grid fails for type 2 inte...
Definition: ecpint.h:196
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:200
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:78
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
void compute_pair(const GaussianShell &shellA, const GaussianShell &shellB) override
Overridden shell-pair integral calculation over all ECP centers.
Definition: ecpint.cc:915
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:79