1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /** @file integrals_1el_single.cc
31
32 @brief Functionality for computing a single 1-electron integral,
33 for a given primitive Gaussian distribution and a given point
34 charge.
35
36 @author: Elias Rudberg <em>responsible</em>
37 */
38
39 #include <stdlib.h>
40 #include <math.h>
41 #include <stdio.h>
42 #include <errno.h>
43 #include <memory.h>
44 #include <time.h>
45 #include <stdarg.h>
46 #include "integrals_1el_single.h"
47 #include "pi.h"
48 #include "boysfunction.h"
49
50 #include "integrals_hermite.h"
51
52
53 static ergo_real
do_1e_repulsion_integral_using_symb_info_h(const DistributionSpecStruct & psi,ergo_real pointCharge,const ergo_real * pointChargeCoords,const IntegralInfo & integralInfo)54 do_1e_repulsion_integral_using_symb_info_h(const DistributionSpecStruct & psi,
55 ergo_real pointCharge,
56 const ergo_real* pointChargeCoords,
57 const IntegralInfo & integralInfo)
58 {
59 // Let the distr be 1 and the pointcharge 2
60 ergo_real alpha1 = psi.exponent;
61 // alpha2 is ~ infinity
62 ergo_real alpha0 = alpha1;
63 // alpham is ~ infinity
64 // g is 1 (not needed)
65
66 int n1 = 0;
67 int n2 = 0;
68 for(int i = 0; i < 3; i++)
69 n1 += psi.monomialInts[i];
70 int n1x = psi.monomialInts[0];
71 int n1y = psi.monomialInts[1];
72 int n1z = psi.monomialInts[2];
73
74 int noOfMonomials_1 = integralInfo.monomial_info.no_of_monomials_list[n1];
75 int noOfMonomials_2 = integralInfo.monomial_info.no_of_monomials_list[n2];
76
77 ergo_real dx0 = pointChargeCoords[0] - psi.centerCoords[0];
78 ergo_real dx1 = pointChargeCoords[1] - psi.centerCoords[1];
79 ergo_real dx2 = pointChargeCoords[2] - psi.centerCoords[2];
80
81 ergo_real resultPreFactor = 2 * pi / alpha1;
82
83 ergo_real primitiveIntegralList_h[noOfMonomials_1*noOfMonomials_2];
84 ergo_real primitiveIntegralList_2[noOfMonomials_1*noOfMonomials_2];
85
86 const JK::ExchWeights CAM_params_not_used;
87
88 get_related_integrals_hermite(integralInfo,
89 CAM_params_not_used,
90 n1, noOfMonomials_1,
91 n2, noOfMonomials_2,
92 dx0,
93 dx1,
94 dx2,
95 alpha0,
96 resultPreFactor,
97 primitiveIntegralList_h);
98
99 integralInfo.multiply_by_hermite_conversion_matrix_from_right(n1,
100 n2,
101 1.0/alpha1,
102 primitiveIntegralList_h,
103 primitiveIntegralList_2);
104 int monomialIndex = integralInfo.monomial_info.monomial_index_list[n1x][n1y][n1z];
105
106 ergo_real result = psi.coeff * pointCharge * primitiveIntegralList_2[monomialIndex];
107
108 return result;
109 }
110
111
112 /* This routine is supposed to compute derivatives of integrals w.r.t. changes in the pointCharge coordinates. */
do_1e_repulsion_integral_derivatives_using_symb_info(const DistributionSpecStruct * psi,ergo_real pointCharge,const ergo_real * pointChargeCoords,const IntegralInfo & integralInfo)113 std::vector<ergo_real> do_1e_repulsion_integral_derivatives_using_symb_info(const DistributionSpecStruct* psi,
114 ergo_real pointCharge,
115 const ergo_real* pointChargeCoords,
116 const IntegralInfo & integralInfo) {
117 // Let the distr be 1 and the pointcharge 2
118 ergo_real alpha1 = psi->exponent;
119 // alpha2 is ~ infinity
120 ergo_real alpha0 = alpha1;
121 // alpham is ~ infinity
122 // g is 1 (not needed)
123
124 int n1 = 0;
125 int n2 = 1;
126 for(int i = 0; i < 3; i++)
127 n1 += psi->monomialInts[i];
128 int n1x = psi->monomialInts[0];
129 int n1y = psi->monomialInts[1];
130 int n1z = psi->monomialInts[2];
131
132 int noOfMonomials_1 = integralInfo.monomial_info.no_of_monomials_list[n1];
133 int noOfMonomials_2 = integralInfo.monomial_info.no_of_monomials_list[n2];
134
135 ergo_real dx0 = pointChargeCoords[0] - psi->centerCoords[0];
136 ergo_real dx1 = pointChargeCoords[1] - psi->centerCoords[1];
137 ergo_real dx2 = pointChargeCoords[2] - psi->centerCoords[2];
138
139 ergo_real resultPreFactor = 2 * pi / alpha1;
140
141 ergo_real primitiveIntegralList_h[noOfMonomials_1*noOfMonomials_2];
142
143 const JK::ExchWeights CAM_params_not_used;
144
145 get_related_integrals_hermite(integralInfo,
146 CAM_params_not_used,
147 n1, noOfMonomials_1,
148 n2, noOfMonomials_2,
149 dx0,
150 dx1,
151 dx2,
152 alpha0,
153 resultPreFactor,
154 primitiveIntegralList_h);
155
156 int n1b = n1;
157 int n2b = 0;
158 ergo_real primitiveIntegralList_h_components[3][noOfMonomials_1];
159 int monomialIndex_x = integralInfo.monomial_info.monomial_index_list[1][0][0];
160 int monomialIndex_y = integralInfo.monomial_info.monomial_index_list[0][1][0];
161 int monomialIndex_z = integralInfo.monomial_info.monomial_index_list[0][0][1];
162 for(int i = 0; i < noOfMonomials_1; i++) {
163 primitiveIntegralList_h_components[0][i] = primitiveIntegralList_h[i*noOfMonomials_2+monomialIndex_x];
164 primitiveIntegralList_h_components[1][i] = primitiveIntegralList_h[i*noOfMonomials_2+monomialIndex_y];
165 primitiveIntegralList_h_components[2][i] = primitiveIntegralList_h[i*noOfMonomials_2+monomialIndex_z];
166 }
167 ergo_real primitiveIntegralList_2_components[3][noOfMonomials_1];
168 for(int i = 0; i < 3; i++)
169 integralInfo.multiply_by_hermite_conversion_matrix_from_right(n1b, n2b, 1.0/alpha1, primitiveIntegralList_h_components[i], primitiveIntegralList_2_components[i]);
170
171 int monomialIndex = integralInfo.monomial_info.monomial_index_list[n1x][n1y][n1z];
172
173 ergo_real result_x = psi->coeff * pointCharge * primitiveIntegralList_2_components[0][monomialIndex];
174 ergo_real result_y = psi->coeff * pointCharge * primitiveIntegralList_2_components[1][monomialIndex];
175 ergo_real result_z = psi->coeff * pointCharge * primitiveIntegralList_2_components[2][monomialIndex];
176
177 std::vector<ergo_real> v(3);
178 v[0] = result_x;
179 v[1] = result_y;
180 v[2] = result_z;
181 return v;
182 }
183
184
185
186
187
188 ergo_real
do_1e_repulsion_integral_using_symb_info(const DistributionSpecStruct & psi,ergo_real pointCharge,const ergo_real * pointChargeCoords,const IntegralInfo & integralInfo)189 do_1e_repulsion_integral_using_symb_info(const DistributionSpecStruct & psi,
190 ergo_real pointCharge,
191 const ergo_real* pointChargeCoords,
192 const IntegralInfo & integralInfo)
193 {
194 return do_1e_repulsion_integral_using_symb_info_h(psi,
195 pointCharge,
196 pointChargeCoords,
197 integralInfo);
198 }
199
200
201
202