Skip to content

Commit e03492c

Browse files
ifilotifilot
committed
Fixing issue #78 
1 parent 118bb20 commit e03492c

File tree

7 files changed

+275
-54
lines changed

7 files changed

+275
-54
lines changed

‎pyproject.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,7 @@ module-path = "src"
4343
# ---------------------------------------------------------------------------
4444
[project]
4545
name = "pyqint"
46-
version = "1.4.2"
46+
version = "1.4.3"
4747
description = "Python package for evaluating integrals of Gaussian type orbitals"
4848
readme = "README.md"
4949
license = { text = "GPL-3.0-only" }

‎src/pyqint/cgf.cpp‎

Lines changed: 16 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
**************************************************************************/
2121

2222
#include "cgf.h"
23+
#include "gaussian_integrals.h"
2324
#include <cassert>
2425

2526
GTO::GTO(double _c,
@@ -380,29 +381,29 @@ double CGF::get_contraction_norm() const {
380381
}
381382
}
382383

383-
// Standard contraction: all primitives have same (l,m,n)
384-
const int L = l + m + n;
385-
386-
static const double pi_three_half_pow = M_PI * std::sqrt(M_PI);
387-
const double prefactor = pi_three_half_pow *
388-
(l < 1 ? 1.0 : double_factorial(2*l - 1)) *
389-
(m < 1 ? 1.0 : double_factorial(2*m - 1)) *
390-
(n < 1 ? 1.0 : double_factorial(2*n - 1)) /
391-
static_cast<double>(1 << L);
392-
393384
double sum = 0.0;
394385
for (size_t i = 0; i < this->size(); ++i) {
395386
for (size_t j = 0; j < this->size(); ++j) {
396387
// const double a_i = this->get_coefficient_gto(i);
397388
// const double a_j = this->get_coefficient_gto(j);
398389
const double a_i = this->get_norm_gto(i) * this->get_coefficient_gto(i);
399390
const double a_j = this->get_norm_gto(j) * this->get_coefficient_gto(j);
400-
const double alpha_i = this->get_gto(i).get_alpha();
401-
const double alpha_j = this->get_gto(j).get_alpha();
402-
403-
sum += a_i * a_j / std::pow(alpha_i + alpha_j, L + 1.5);
391+
const auto& gto_i = this->get_gto(i);
392+
const auto& gto_j = this->get_gto(j);
393+
394+
sum += a_i * a_j *
395+
integrals::gaussian::overlap_gto(gto_i.get_alpha(),
396+
gto_i.get_l(),
397+
gto_i.get_m(),
398+
gto_i.get_n(),
399+
gto_i.get_position(),
400+
gto_j.get_alpha(),
401+
gto_j.get_l(),
402+
gto_j.get_m(),
403+
gto_j.get_n(),
404+
gto_j.get_position());
404405
}
405406
}
406407

407-
return 1.0 / std::sqrt(prefactor * sum);
408+
return 1.0 / std::sqrt(sum);
408409
}

‎src/pyqint/gaussian_integrals.cpp‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
/**************************************************************************
2+
* This file is part of PyQInt. *
3+
* *
4+
* Author: Ivo Filot <ivo@ivofilot.nl> *
5+
* *
6+
* PyQInt is free software: *
7+
* you can redistribute it and/or modify it under the terms of the *
8+
* GNU General Public License as published by the Free Software *
9+
* Foundation, either version 3 of the License, or (at your option) *
10+
* any later version. *
11+
* *
12+
* PyQInt is distributed in the hope that it will be useful, *
13+
* but WITHOUT ANY WARRANTY; without even the implied warranty *
14+
* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
15+
* See the GNU General Public License for more details. *
16+
* *
17+
* You should have received a copy of the GNU General Public License *
18+
* along with this program. If not, see http://www.gnu.org/licenses/. *
19+
* *
20+
**************************************************************************/
21+
22+
#define _USE_MATH_DEFINES
23+
#include <cmath>
24+
25+
#include "factorials.h"
26+
#include "gaussian_integrals.h"
27+
28+
namespace integrals::gaussian {
29+
30+
/**
31+
* @brief Computes the Gaussian product center of two primitive Gaussians
32+
*
33+
* @param double alpha1 Exponent of first Gaussian
34+
* @param const Vec3& a Center of first Gaussian
35+
* @param double alpha2 Exponent of second Gaussian
36+
* @param const Vec3& b Center of second Gaussian
37+
*
38+
* Calculates the position of the Gaussian product center P.
39+
*
40+
* @return Vec3 Position vector of the Gaussian product center
41+
*/
42+
Vec3 gaussian_product_center(double alpha1, const Vec3& a,
43+
double alpha2, const Vec3& b) {
44+
return (alpha1 * a + alpha2 * b) / (alpha1 + alpha2);
45+
}
46+
47+
/**
48+
* @brief Computes the binomial coefficient
49+
*
50+
* @param int a Upper index
51+
* @param int b Lower index
52+
*
53+
* Returns 1 if indices are outside valid range.
54+
*
55+
* @return double Value of the binomial coefficient
56+
*/
57+
double binomial(int a, int b) {
58+
if ((a < 0) || (b < 0) || (a - b < 0)) {
59+
return 1.0;
60+
}
61+
return factorial(a) / (factorial(b) * factorial(a - b));
62+
}
63+
64+
/**
65+
* @brief Computes binomial prefactor used in Gaussian integral recursion
66+
*
67+
* @param int s Summation index
68+
* @param int ia Angular momentum component of first Gaussian
69+
* @param int ib Angular momentum component of second Gaussian
70+
* @param double xpa Distance P_x - A_x
71+
* @param double xpb Distance P_x - B_x
72+
*
73+
* @return double Value of the binomial prefactor
74+
*/
75+
double binomial_prefactor(int s, int ia, int ib,
76+
double xpa, double xpb) {
77+
double sum = 0.0;
78+
79+
for (int t = 0; t < s + 1; t++) {
80+
if ((s - ia <= t) && (t <= ib)) {
81+
sum += binomial(ia, s - t) *
82+
binomial(ib, t) *
83+
std::pow(xpa, ia - s + t) *
84+
std::pow(xpb, ib - t);
85+
}
86+
}
87+
88+
return sum;
89+
}
90+
91+
/**
92+
* @brief Computes one-dimensional overlap integral component
93+
*
94+
* @param int l1 Angular momentum of first Gaussian along axis
95+
* @param int l2 Angular momentum of second Gaussian along axis
96+
* @param double x1 Distance P_x - A_x
97+
* @param double x2 Distance P_x - B_x
98+
* @param double gamma Sum of Gaussian exponents
99+
*
100+
* @return double Value of the 1D overlap contribution
101+
*/
102+
double overlap_1d(int l1, int l2, double x1, double x2, double gamma) {
103+
double sum = 0.0;
104+
105+
for (int i = 0; i < (1 + std::floor(0.5 * (l1 + l2))); i++) {
106+
sum += binomial_prefactor(2 * i, l1, l2, x1, x2) *
107+
(i == 0 ? 1 : double_factorial(2 * i - 1)) /
108+
std::pow(2 * gamma, i);
109+
}
110+
111+
return sum;
112+
}
113+
114+
/**
115+
* @brief Computes overlap integral between two primitive Gaussian orbitals
116+
*
117+
* @param double alpha1 Exponent of first Gaussian
118+
* @param unsigned int l1,m1,n1 Angular momentum components of first Gaussian
119+
* @param const Vec3& a Center of first Gaussian
120+
* @param double alpha2 Exponent of second Gaussian
121+
* @param unsigned int l2,m2,n2 Angular momentum components of second Gaussian
122+
* @param const Vec3& b Center of second Gaussian
123+
*
124+
* @return double Value of the 3D overlap integral
125+
*/
126+
double overlap_gto(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const Vec3& a,
127+
double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const Vec3& b) {
128+
129+
double rab2 = (a - b).norm2();
130+
double gamma = alpha1 + alpha2;
131+
Vec3 p = gaussian_product_center(alpha1, a, alpha2, b);
132+
133+
double pre = std::pow(M_PI / gamma, 1.5) * std::exp(-alpha1 * alpha2 * rab2 / gamma);
134+
double wx = overlap_1d(l1, l2, p[0] - a[0], p[0] - b[0], gamma);
135+
double wy = overlap_1d(m1, m2, p[1] - a[1], p[1] - b[1], gamma);
136+
double wz = overlap_1d(n1, n2, p[2] - a[2], p[2] - b[2], gamma);
137+
138+
return pre * wx * wy * wz;
139+
}
140+
141+
} // namespace normalization

‎src/pyqint/gaussian_integrals.h‎

Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
/**************************************************************************
2+
* This file is part of PyQInt. *
3+
* *
4+
* Author: Ivo Filot <ivo@ivofilot.nl> *
5+
* *
6+
* PyQInt is free software: *
7+
* you can redistribute it and/or modify it under the terms of the *
8+
* GNU General Public License as published by the Free Software *
9+
* Foundation, either version 3 of the License, or (at your option) *
10+
* any later version. *
11+
* *
12+
* PyQInt is distributed in the hope that it will be useful, *
13+
* but WITHOUT ANY WARRANTY; without even the implied warranty *
14+
* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
15+
* See the GNU General Public License for more details. *
16+
* *
17+
* You should have received a copy of the GNU General Public License *
18+
* along with this program. If not, see http://www.gnu.org/licenses/. *
19+
* *
20+
**************************************************************************/
21+
22+
#pragma once
23+
24+
#include "vec3.h"
25+
26+
namespace integrals::gaussian {
27+
28+
/**
29+
* @brief Computes the Gaussian product center of two primitive Gaussians
30+
*
31+
* @param double alpha1 Exponent of first Gaussian
32+
* @param const Vec3& a Center of first Gaussian
33+
* @param double alpha2 Exponent of second Gaussian
34+
* @param const Vec3& b Center of second Gaussian
35+
*
36+
* Calculates the position of the Gaussian product center P given by:
37+
* P = (alpha1*A + alpha2*B) / (alpha1 + alpha2)
38+
*
39+
* @return Vec3 Position vector of the Gaussian product center
40+
*/
41+
Vec3 gaussian_product_center(double alpha1, const Vec3& a,
42+
double alpha2, const Vec3& b);
43+
44+
/**
45+
* @brief Computes the binomial coefficient
46+
*
47+
* @param int a Upper index
48+
* @param int b Lower index
49+
*
50+
* Calculates the binomial coefficient:
51+
* (a choose b) = a! / (b! (a-b)!)
52+
*
53+
* @return double Value of the binomial coefficient
54+
*/
55+
double binomial(int a, int b);
56+
57+
/**
58+
* @brief Computes binomial prefactor used in Gaussian integral recursion
59+
*
60+
* @param int s Summation index
61+
* @param int ia Angular momentum component of first Gaussian
62+
* @param int ib Angular momentum component of second Gaussian
63+
* @param double xpa Distance P_x - A_x
64+
* @param double xpb Distance P_x - B_x
65+
*
66+
* Evaluates the binomial prefactor appearing in Hermite Gaussian
67+
* expansion formulas for overlap integrals.
68+
*
69+
* @return double Value of the binomial prefactor
70+
*/
71+
double binomial_prefactor(int s, int ia, int ib, double xpa, double xpb);
72+
73+
/**
74+
* @brief Computes one-dimensional overlap integral component
75+
*
76+
* @param int l1 Angular momentum of first Gaussian along axis
77+
* @param int l2 Angular momentum of second Gaussian along axis
78+
* @param double x1 Distance P_x - A_x
79+
* @param double x2 Distance P_x - B_x
80+
* @param double gamma Sum of Gaussian exponents (alpha1 + alpha2)
81+
*
82+
* Computes the 1D overlap term used in separable Gaussian integrals.
83+
*
84+
* @return double Value of the 1D overlap contribution
85+
*/
86+
double overlap_1d(int l1, int l2, double x1, double x2, double gamma);
87+
88+
/**
89+
* @brief Computes overlap integral between two primitive Gaussian orbitals
90+
*
91+
* @param double alpha1 Exponent of first Gaussian
92+
* @param unsigned int l1 Angular momentum in x for first Gaussian
93+
* @param unsigned int m1 Angular momentum in y for first Gaussian
94+
* @param unsigned int n1 Angular momentum in z for first Gaussian
95+
* @param const Vec3& a Center of first Gaussian
96+
* @param double alpha2 Exponent of second Gaussian
97+
* @param unsigned int l2 Angular momentum in x for second Gaussian
98+
* @param unsigned int m2 Angular momentum in y for second Gaussian
99+
* @param unsigned int n2 Angular momentum in z for second Gaussian
100+
* @param const Vec3& b Center of second Gaussian
101+
*
102+
* Calculates the value of < gto1 | gto2 > for primitive Gaussian type orbitals.
103+
*
104+
* @return double Value of the 3D overlap integral
105+
*/
106+
double overlap_gto(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const Vec3& a,
107+
double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const Vec3& b);
108+
109+
} // namespace normalization

‎src/pyqint/integrals.cpp‎

Lines changed: 6 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
**************************************************************************/
2121

2222
#include "integrals.h"
23+
#include "gaussian_integrals.h"
2324

2425
/**
2526
* @brief Construct an Integrator instance and capture build-time information.
@@ -1113,17 +1114,7 @@ double Integrator::repulsion_deriv(const GTO& gto1, const GTO& gto2, const GTO &
11131114
*/
11141115
double Integrator::overlap(double alpha1, unsigned int l1, unsigned int m1, unsigned int n1, const Vec3 &a,
11151116
double alpha2, unsigned int l2, unsigned int m2, unsigned int n2, const Vec3 &b) const {
1116-
1117-
double rab2 = (a-b).norm2();
1118-
double gamma = alpha1 + alpha2;
1119-
Vec3 p = this->gaussian_product_center(alpha1, a, alpha2, b);
1120-
1121-
double pre = std::pow(M_PI / gamma, 1.5) * std::exp(-alpha1 * alpha2 * rab2 / gamma);
1122-
double wx = this->overlap_1D(l1, l2, p[0]-a[0], p[0]-b[0], gamma);
1123-
double wy = this->overlap_1D(m1, m2, p[1]-a[1], p[1]-b[1], gamma);
1124-
double wz = this->overlap_1D(n1, n2, p[2]-a[2], p[2]-b[2], gamma);
1125-
1126-
return pre * wx * wy * wz;
1117+
return integrals::gaussian::overlap_gto(alpha1, l1, m1, n1, a, alpha2, l2, m2, n2, b);
11271118
}
11281119

11291120
/**
@@ -1188,15 +1179,7 @@ double Integrator::dipole(double alpha1, unsigned int l1, unsigned int m1, unsig
11881179
* @return double value of the one dimensional overlap integral
11891180
*/
11901181
double Integrator::overlap_1D(int l1, int l2, double x1, double x2, double gamma) const {
1191-
double sum = 0.0;
1192-
1193-
for(int i=0; i < (1 + std::floor(0.5 * (l1 + l2))); i++) {
1194-
sum += this->binomial_prefactor(2*i, l1, l2, x1, x2) *
1195-
(i == 0 ? 1 : double_factorial(2 * i - 1) ) /
1196-
std::pow(2 * gamma, i);
1197-
}
1198-
1199-
return sum;
1182+
return integrals::gaussian::overlap_1d(l1, l2, x1, x2, gamma);
12001183
}
12011184

12021185
/**
@@ -1212,30 +1195,16 @@ double Integrator::overlap_1D(int l1, int l2, double x1, double x2, double gamma
12121195
*/
12131196
Vec3 Integrator::gaussian_product_center(double alpha1, const Vec3& a,
12141197
double alpha2, const Vec3& b) const {
1215-
return (alpha1 * a + alpha2 * b) / (alpha1 + alpha2);
1198+
return integrals::gaussian::gaussian_product_center(alpha1, a, alpha2, b);
12161199
}
12171200

12181201
double Integrator::binomial_prefactor(int s, int ia, int ib,
12191202
double xpa, double xpb) const {
1220-
double sum = 0.0;
1221-
1222-
for (int t=0; t < s+1; t++) {
1223-
if ((s-ia <= t) && (t <= ib)) {
1224-
sum += this->binomial(ia, s-t) *
1225-
this->binomial(ib, t) *
1226-
std::pow(xpa, ia - s + t) *
1227-
std::pow(xpb, ib - t);
1228-
}
1229-
}
1230-
1231-
return sum;
1203+
return integrals::gaussian::binomial_prefactor(s, ia, ib, xpa, xpb);
12321204
}
12331205

12341206
double Integrator::binomial(int a, int b) const {
1235-
if( (a < 0) || (b < 0) || (a-b < 0) ) {
1236-
return 1.0;
1237-
}
1238-
return factorial(a) / (factorial(b) * factorial(a-b));
1207+
return integrals::gaussian::binomial(a, b);
12391208
}
12401209

12411210
double Integrator::nuclear(const Vec3& a, int l1, int m1, int n1, double alpha1,

‎src/pyqint/meson.build‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ py.extension_module(
2525
'fgamma.cpp',
2626
'hellsing_cache.cpp',
2727
'hellsing_eri.cpp',
28+
'gaussian_integrals.cpp',
2829
],
2930
include_directories: inc,
3031
override_options: ['cython_language=cpp'],

0 commit comments

Comments
 (0)