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_2el_explicit.cc
31
32 @brief Code for explicit computation of 4-index 2-electron
33 integrals.
34
35 @author: Elias Rudberg <em>responsible</em>
36 */
37
38 /* Written by Elias Rudberg, KTH, Stockholm */
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_2el_explicit.h"
47 #include "memorymanag.h"
48 #include "pi.h"
49 #include "output.h"
50 #include "utilities.h"
51 #include "boysfunction.h"
52 #include "integral_info.h"
53 #include "integrals_general.h"
54 #include "integrals_2el_single.h"
55
56
57
58 typedef struct
59 {
60 int a, b, c, d;
61 int poly_ab_index;
62 int poly_cd_index;
63 } abcd_struct;
64
65 #define set_abcd_list_item_macro(i,A,B,C,D) \
66 list[i].a = A; list[i].b = B; list[i].c = C; list[i].d = D;
67
68
69
70
71 static int globalCount = 0;
72
73
74 ergo_real
do_2e_integral(int mu,int nu,int la,int si,const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo)75 do_2e_integral(int mu,
76 int nu,
77 int la,
78 int si,
79 const BasisInfoStruct & basisInfo,
80 const IntegralInfo & integralInfo) {
81 return do_2e_integral_general(mu,
82 nu,
83 la,
84 si,
85 basisInfo,
86 basisInfo,
87 basisInfo,
88 basisInfo,
89 integralInfo);
90 }
91
92 ergo_real
do_2e_integral_general(int mu,int nu,int la,int si,const BasisInfoStruct & basisInfo_mu,const BasisInfoStruct & basisInfo_nu,const BasisInfoStruct & basisInfo_la,const BasisInfoStruct & basisInfo_si,const IntegralInfo & integralInfo)93 do_2e_integral_general(int mu,
94 int nu,
95 int la,
96 int si,
97 const BasisInfoStruct & basisInfo_mu,
98 const BasisInfoStruct & basisInfo_nu,
99 const BasisInfoStruct & basisInfo_la,
100 const BasisInfoStruct & basisInfo_si,
101 const IntegralInfo & integralInfo)
102 {
103 int n_psi2, i, j;
104 ergo_real sum, currIntegral;
105 const int maxCount = 1000;
106 DistributionSpecStruct list_psi1[maxCount];
107 DistributionSpecStruct list_psi2[maxCount];
108
109 /* form product of basisfuncs mu and nu, store product in psi1 */
110 int n_psi1 = get_product_simple_primitives(basisInfo_mu, mu,
111 basisInfo_nu, nu,
112 list_psi1,
113 maxCount,
114 0);
115 if(n_psi1 <= 0)
116 {
117 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_product_simple_primitives\n");
118 exit(EXIT_FAILURE);
119 return 0;
120 }
121
122 /* form product of basisfuncs la and si, store product in psi2 */
123 n_psi2 = get_product_simple_primitives(basisInfo_la, la,
124 basisInfo_si, si,
125 list_psi2,
126 maxCount,
127 0);
128 if(n_psi2 <= 0)
129 {
130 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in get_product_simple_primitives\n");
131 exit(EXIT_FAILURE);
132 return 0;
133 }
134
135 const JK::ExchWeights CAM_params_not_used;
136
137 sum = 0;
138 for(i = 0; i < n_psi1; i++)
139 {
140 DistributionSpecStruct* prim_psi1 = &list_psi1[i];
141 for(j = 0; j < n_psi2; j++)
142 {
143 DistributionSpecStruct* prim_psi2 = &list_psi2[j];
144 globalCount++;
145 currIntegral = do_2e_integral_using_symb_info(CAM_params_not_used, prim_psi1, prim_psi2, integralInfo);
146 sum += currIntegral;
147 } /* END FOR j */
148 } /* END FOR i */
149
150 return sum;
151 }
152
153
154 /** compute_2e_matrix_simple computes the 2el matrix in the simplest
155 possible way. It assumes that the matrix is computed for closed
156 shell. The weight of the HF exchange is controlled by @param hf_weight
157 which is equal 1 for ordinary Hartree-Fock calculation. No
158 assumption are made regarding symmetry of the density matrix
159 @param dens . The computed two-electron part of the Fock matrix is
160 returned in @param result .
161 @param basisInfo info about the used basis set.
162 @param integralInfo info needed for evaluation of integrals of Gaussian functions.
163 */
164 int
compute_2e_matrix_simple(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,ergo_real hf_weight,ergo_real * result,const ergo_real * dens)165 compute_2e_matrix_simple(const BasisInfoStruct & basisInfo,
166 const IntegralInfo & integralInfo,
167 ergo_real hf_weight,
168 ergo_real* result,
169 const ergo_real* dens)
170 {
171 int mu, nu, sigma, lambda;
172 int nbast = basisInfo.noOfBasisFuncs;
173 ergo_real munusila, mulasinu, sum;
174 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "entering compute_2e_matrix HF_WEIGHT=%f", (double)hf_weight);
175
176 for(mu = 0; mu < nbast; mu++)
177 {
178 for(nu = 0; nu < nbast; nu++)
179 {
180 sum = 0;
181 for(lambda = 0; lambda < nbast; lambda++)
182 {
183 for(sigma = 0; sigma < nbast; sigma++)
184 {
185 munusila = do_2e_integral(mu, nu, sigma, lambda, basisInfo, integralInfo);
186 mulasinu = do_2e_integral(mu, lambda, sigma, nu, basisInfo, integralInfo);
187 sum +=
188 dens[lambda*nbast+sigma] *
189 (munusila - 0.5 * hf_weight * mulasinu);
190 } /* END FOR sigma */
191 } /* END FOR lambda */
192 result[mu*nbast+nu] = sum;
193 } /* END FOR nu */
194 } /* END FOR mu */
195
196 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_2e_matrix ending OK\n");
197 return 0;
198 }
199
200
201
202
203 static int
compute_J_and_K_integraldriven(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,ergo_real * J,ergo_real * K,ergo_real * dens)204 compute_J_and_K_integraldriven(const BasisInfoStruct & basisInfo,
205 const IntegralInfo & integralInfo,
206 ergo_real* J,
207 ergo_real* K,
208 ergo_real* dens)
209 {
210 int n, nBytes, i, j, count, a, b, c, d;
211
212 Util::TimeMeter timeMeter;
213
214 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "entering compute_J_and_K_integraldriven");
215 n = basisInfo.noOfBasisFuncs;
216 nBytes = n * n * sizeof(ergo_real);
217 memset(J, 0, nBytes);
218 memset(K, 0, nBytes);
219
220 count = 0;
221 a = 0;
222 b = 0;
223 c = 0;
224 d = 0;
225 while(d < n)
226 {
227 abcd_struct list[8];
228
229 /* compute integral */
230 ergo_real integralValue = do_2e_integral(a, b, c, d, basisInfo, integralInfo);
231
232 count++;
233
234 /* determine unique configurations */
235 set_abcd_list_item_macro(0, a, b, c, d);
236 set_abcd_list_item_macro(1, a, b, d, c);
237 set_abcd_list_item_macro(2, b, a, c, d);
238 set_abcd_list_item_macro(3, b, a, d, c);
239 set_abcd_list_item_macro(4, c, d, a, b);
240 set_abcd_list_item_macro(5, d, c, a, b);
241 set_abcd_list_item_macro(6, c, d, b, a);
242 set_abcd_list_item_macro(7, d, c, b, a);
243
244 for(i = 0; i < 8; i++)
245 {
246 abcd_struct* abcd = &list[i];
247 int aa, bb, cc, dd;
248
249 /* check if this is a new unique configuration */
250 int unique = 1;
251 for(j = 0; j < i; j++)
252 {
253 if(abcd->a == list[j].a &&
254 abcd->b == list[j].b &&
255 abcd->c == list[j].c &&
256 abcd->d == list[j].d)
257 unique = 0;
258 }
259 if(unique == 0)
260 continue;
261 /* now we know that this configuration is unique. */
262 aa = abcd->a;
263 bb = abcd->b;
264 cc = abcd->c;
265 dd = abcd->d;
266
267 #if 1
268 /* add contribution to coulomb matrix */
269 J[aa*n+bb] += dens[cc*n+dd] * integralValue;
270
271 /* add contribution to exchange matrix */
272 K[aa*n+dd] += -0.5 * dens[bb*n+cc] * integralValue;
273 #endif
274 } /* END FOR i go through 8 configurations */
275
276 /* now get numbers for next unique integral */
277 d++;
278 if(d < n)
279 continue;
280
281 /* d has hit the roof */
282 c++;
283 if(c < n)
284 {
285 d = c;
286 continue;
287 }
288
289 /* c has hit roof */
290 b++;
291 if(b < n)
292 {
293 c = a;
294 d = b;
295 continue;
296 }
297
298 /* b has hit roof */
299 a++;
300 if(a < n)
301 {
302 b = a;
303 c = a;
304 d = a;
305 continue;
306 }
307
308 /* a has hit roof. This means that we are done. */
309 break;
310
311 } /* END WHILE more unique integrals */
312
313 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_J_and_K_integraldriven ending OK");
314 timeMeter.print(LOG_AREA_INTEGRALS, "compute_J_and_K_integraldriven");
315 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "number of unique integrals computed: %i", count);
316
317 return 0;
318 }
319
320
321
322 int
compute_2e_matrix_list_explicit(const BasisInfoStruct & basisInfo,const IntegralInfo & integralInfo,ergo_real ** resultList,ergo_real ** densList,int noOfMatrices,ergo_real threshold)323 compute_2e_matrix_list_explicit(const BasisInfoStruct & basisInfo,
324 const IntegralInfo & integralInfo,
325 ergo_real** resultList,
326 ergo_real** densList,
327 int noOfMatrices,
328 ergo_real threshold)
329 {
330 ergo_real* J;
331 int n, i, j;
332
333 if(noOfMatrices != 1)
334 do_output(LOG_CAT_INFO, LOG_AREA_INTEGRALS, "compute_2e_matrix_list_explicit: (noOfMatrices != 1), will take some time");
335
336 n = basisInfo.noOfBasisFuncs;
337 J = (ergo_real*)ergo_malloc(n*n*sizeof(ergo_real));
338 for(j = 0; j < noOfMatrices; j++)
339 {
340 if(compute_J_and_K_integraldriven(basisInfo,
341 integralInfo,
342 J,
343 resultList[j],
344 densList[j]) != 0)
345 {
346 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_J_and_K_integraldriven");
347 return -1;
348 }
349 for(i = 0; i < n*n; i++)
350 {
351 resultList[j][i] += J[i];
352 } // END FOR i
353 } // END FOR j
354 ergo_free(J);
355 return 0;
356 }
357
358
359