(git:374b731)
Loading...
Searching...
No Matches
eri_mme_test.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief Methods for testing / debugging.
9!> \par History
10!> 2015 09 created
11!> \author Patrick Seewald
12! **************************************************************************************************
13
14MODULE eri_mme_test
15
16 USE eri_mme_gaussian, ONLY: create_gaussian_overlap_dist_to_hermite,&
17 create_hermite_to_cartesian
18 USE eri_mme_integrate, ONLY: eri_mme_2c_integrate,&
19 eri_mme_3c_integrate
20 USE eri_mme_types, ONLY: eri_mme_coulomb,&
21 eri_mme_longrange,&
22 eri_mme_param,&
23 eri_mme_set_potential,&
24 eri_mme_yukawa
25 USE kinds, ONLY: dp
26 USE mathconstants, ONLY: twopi
27 USE message_passing, ONLY: mp_para_env_type
28 USE orbital_pointers, ONLY: ncoset
29#include "../base/base_uses.f90"
30
31 IMPLICIT NONE
32
33 PRIVATE
34
35 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
36
37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_test'
38
39 PUBLIC :: eri_mme_2c_perf_acc_test, &
40 eri_mme_3c_perf_acc_test
41
42CONTAINS
43! **************************************************************************************************
44!> \brief Unit test for performance and accuracy
45!> \param param ...
46!> \param l_max ...
47!> \param zet ...
48!> \param rabc ...
49!> \param nrep ...
50!> \param test_accuracy ...
51!> \param para_env ...
52!> \param iw ...
53!> \param potential ...
54!> \param pot_par ...
55!> \param G_count ...
56!> \param R_count ...
57! **************************************************************************************************
58 SUBROUTINE eri_mme_2c_perf_acc_test(param, l_max, zet, rabc, nrep, test_accuracy, &
59 para_env, iw, potential, pot_par, G_count, R_count)
60 TYPE(eri_mme_param), INTENT(INOUT) :: param
61 INTEGER, INTENT(IN) :: l_max
62 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
63 INTENT(IN) :: zet
64 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: rabc
65 INTEGER, INTENT(IN) :: nrep
66 LOGICAL, INTENT(INOUT) :: test_accuracy
67 TYPE(mp_para_env_type), INTENT(IN) :: para_env
68 INTEGER, INTENT(IN) :: iw
69 INTEGER, INTENT(IN), OPTIONAL :: potential
70 REAL(kind=dp), INTENT(IN), OPTIONAL :: pot_par
71 INTEGER, INTENT(OUT), OPTIONAL :: g_count, r_count
72
73 INTEGER :: iab, irep, izet, l, nr, nzet
74 LOGICAL :: acc_check
75 REAL(kind=dp) :: acc, t0, t1
76 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: time
77 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: i_diff, i_ref, i_test
78 REAL(kind=dp), DIMENSION(3, 3) :: ht
79
80 IF (PRESENT(g_count)) g_count = 0
81 IF (PRESENT(r_count)) r_count = 0
82
83 nzet = SIZE(zet)
84 nr = SIZE(rabc, 2)
85
86 IF (PRESENT(potential)) THEN
87 CALL eri_mme_set_potential(param, potential, pot_par)
88 END IF
89
90 ! Calculate reference values (Exact expression in G space converged to high precision)
91 IF (test_accuracy) THEN
92 ht = twopi*transpose(param%h_inv)
93
94 ALLOCATE (i_ref(ncoset(l_max), ncoset(l_max), nr, nzet))
95 i_ref(:, :, :, :) = 0.0_dp
96
97 DO izet = 1, nzet
98 DO iab = 1, nr
99 CALL eri_mme_2c_integrate(param, 0, l_max, 0, l_max, zet(izet), zet(izet), rabc(:, iab), &
100 i_ref(:, :, iab, izet), 0, 0, &
101 normalize=.true., potential=potential, pot_par=pot_par)
102
103 END DO
104 END DO
105 END IF
106
107 ! test performance and accuracy of MME method
108 ALLOCATE (i_test(ncoset(l_max), ncoset(l_max), nr, nzet))
109 ALLOCATE (i_diff(ncoset(l_max), ncoset(l_max), nr, nzet))
110
111 ALLOCATE (time(0:l_max, nzet))
112 DO l = 0, l_max
113 DO izet = 1, nzet
114 CALL cpu_time(t0)
115 DO irep = 1, nrep
116 DO iab = 1, nr
117 CALL eri_mme_2c_integrate(param, 0, l, 0, l, zet(izet), zet(izet), rabc(:, iab), &
118 i_test(:, :, iab, izet), 0, 0, &
119 g_count=g_count, r_count=r_count, &
120 normalize=.true.)
121 END DO
122 END DO
123 CALL cpu_time(t1)
124 time(l, izet) = t1 - t0
125 END DO
126 END DO
127
128 CALL para_env%sum(time)
129
130 IF (test_accuracy) THEN
131 i_diff(:, :, :, :) = abs(i_test - i_ref)
132 END IF
133
134 IF (iw > 0) THEN
135 WRITE (iw, '(T2, A)') "ERI_MME| Test results for 2c cpu time"
136 WRITE (iw, '(T11, A)') "l, zet, cpu time, accuracy"
137
138 DO l = 0, l_max
139 DO izet = 1, nzet
140 IF (test_accuracy) THEN
141 acc = maxval(i_diff(ncoset(l - 1) + 1:ncoset(l), ncoset(l - 1) + 1:ncoset(l), :, izet))
142 ELSE
143 acc = 0.0_dp
144 END IF
145
146 WRITE (iw, '(T11, I1, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2)') &
147 l, zet(izet), time(l, izet)/nrep, acc
148 END DO
149 END DO
150
151 IF (test_accuracy) THEN
152 WRITE (iw, '(/T2, A, 47X, ES9.2)') "ERI_MME| Maximum error:", &
153 maxval(i_diff)
154
155 IF (param%is_ortho) THEN
156 acc_check = param%err_mm + param%err_c .GE. maxval(i_diff)
157 ELSE
158 acc_check = .true.
159 END IF
160
161 IF (.NOT. acc_check) &
162 cpabort("Actual error greater than upper bound estimate.")
163
164 END IF
165 END IF
166
167 END SUBROUTINE eri_mme_2c_perf_acc_test
168
169! **************************************************************************************************
170!> \brief ...
171!> \param param ...
172!> \param l_max ...
173!> \param zet ...
174!> \param rabc ...
175!> \param nrep ...
176!> \param nsample ...
177!> \param para_env ...
178!> \param iw ...
179!> \param potential ...
180!> \param pot_par ...
181!> \param GG_count ...
182!> \param GR_count ...
183!> \param RR_count ...
184! **************************************************************************************************
185 SUBROUTINE eri_mme_3c_perf_acc_test(param, l_max, zet, rabc, nrep, nsample, &
186 para_env, iw, potential, pot_par, GG_count, GR_count, RR_count)
187 TYPE(eri_mme_param), INTENT(INOUT) :: param
188 INTEGER, INTENT(IN) :: l_max
189 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
190 INTENT(IN) :: zet
191 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
192 INTENT(IN) :: rabc
193 INTEGER, INTENT(IN) :: nrep
194 INTEGER, INTENT(IN), OPTIONAL :: nsample
195 TYPE(mp_para_env_type), INTENT(IN) :: para_env
196 INTEGER, INTENT(IN) :: iw
197 INTEGER, INTENT(IN), OPTIONAL :: potential
198 REAL(kind=dp), INTENT(IN), OPTIONAL :: pot_par
199 INTEGER, INTENT(OUT), OPTIONAL :: gg_count, gr_count, rr_count
200
201 INTEGER :: ira, irb, irc, irep, ixyz, izeta, izetb, &
202 izetc, la, lb, lc, nintg, nr, ns, nzet
203 REAL(kind=dp) :: t0, t1, time
204 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: i_test
205
206 IF (PRESENT(gg_count)) gg_count = 0
207 IF (PRESENT(gr_count)) gr_count = 0
208 IF (PRESENT(rr_count)) rr_count = 0
209
210 ns = 1
211 IF (PRESENT(nsample)) ns = nsample
212
213 nzet = SIZE(zet)
214 nr = SIZE(rabc, 2)
215
216 IF (PRESENT(potential)) THEN
217 CALL eri_mme_set_potential(param, potential, pot_par)
218 END IF
219
220 IF (param%debug) THEN
221 DO izeta = 1, nzet
222 DO izetb = 1, nzet
223 DO ira = 1, nr
224 DO irb = 1, nr
225 DO ixyz = 1, 3
226 CALL overlap_dist_expansion_test(l_max, l_max, zet(izeta), zet(izetb), &
227 rabc(ixyz, ira), rabc(ixyz, irb), 0.0_dp, param%debug_delta)
228 END DO
229 END DO
230 END DO
231 END DO
232 END DO
233 END IF
234
235 IF (iw > 0) THEN
236 IF (PRESENT(potential)) THEN
237 SELECT CASE (potential)
238 CASE (eri_mme_coulomb)
239 WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
240 CASE (eri_mme_yukawa)
241 WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: Yukawa with a=", pot_par
242 CASE (eri_mme_longrange)
243 WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: long-range Coulomb with a=", pot_par
244 END SELECT
245 ELSE
246 WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
247 END IF
248 WRITE (iw, '(T2, A)') "ERI_MME| Test results for 3c cpu time"
249 WRITE (iw, '(T11, A)') "la, lb, lc, zeta, zetb, zetc, cpu time"
250 END IF
251
252 ALLOCATE (i_test(ncoset(l_max), ncoset(l_max), ncoset(l_max)))
253
254 nintg = 0
255 DO la = 0, l_max
256 DO lb = 0, l_max
257 DO lc = 0, l_max
258 DO izeta = 1, nzet
259 DO izetb = 1, nzet
260 DO izetc = 1, nzet
261 nintg = nintg + 1
262 IF (mod(nintg, ns) .EQ. 0) THEN
263 i_test(:, :, :) = 0.0_dp
264 CALL cpu_time(t0)
265 DO irep = 1, nrep
266 DO ira = 1, nr
267 DO irb = 1, nr
268 DO irc = 1, nr
269 CALL eri_mme_3c_integrate(param, 0, la, 0, lb, 0, lc, zet(izeta), zet(izetb), zet(izetc), &
270 rabc(:, ira), rabc(:, irb), rabc(:, irc), i_test, 0, 0, 0, &
271 gg_count, gr_count, rr_count)
272 END DO
273 END DO
274 END DO
275 END DO
276 CALL cpu_time(t1)
277 time = t1 - t0
278 CALL para_env%sum(time)
279 IF (iw > 0) THEN
280 WRITE (iw, '(T11, I1, 1X, I1, 1X, I1, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2)') &
281 la, lb, lc, zet(izeta), zet(izetb), zet(izetc), time/nrep
282 END IF
283 END IF
284 END DO
285 END DO
286 END DO
287 END DO
288 END DO
289 END DO
290
291 END SUBROUTINE eri_mme_3c_perf_acc_test
292
293! **************************************************************************************************
294!> \brief check that expanding an overlap distribution of cartesian/hermite Gaussians into a
295!> lin combi of single cartesian/hermite Gaussians is correct.
296!> \param l_max ...
297!> \param m_max ...
298!> \param zeta ...
299!> \param zetb ...
300!> \param R1 ...
301!> \param R2 ...
302!> \param r ...
303!> \param tolerance ...
304!> \note STATUS: tested
305! **************************************************************************************************
306 SUBROUTINE overlap_dist_expansion_test(l_max, m_max, zeta, zetb, R1, R2, r, tolerance)
307 INTEGER, INTENT(IN) :: l_max, m_max
308 REAL(kind=dp), INTENT(IN) :: zeta, zetb, r1, r2, r, tolerance
309
310 INTEGER :: l, m, t
311 REAL(kind=dp) :: c_prod_err, h_prod_err, rp, zetp
312 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: c1, c2, c_ol, h1, h2, h_ol
313 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: c_prod_ref, c_prod_test, h_prod_ref, &
314 h_prod_test, h_to_c_1, h_to_c_2, &
315 h_to_c_ol
316 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: e_c, e_h
317
318 zetp = zeta + zetb
319 rp = (zeta*r1 + zetb*r2)/zetp
320 ALLOCATE (c1(0:l_max), h1(0:l_max))
321 ALLOCATE (c2(0:m_max), h2(0:m_max))
322 ALLOCATE (c_ol(0:l_max + m_max))
323 ALLOCATE (h_ol(0:l_max + m_max))
324 ALLOCATE (c_prod_ref(0:l_max, 0:m_max))
325 ALLOCATE (c_prod_test(0:l_max, 0:m_max))
326 ALLOCATE (h_prod_ref(0:l_max, 0:m_max))
327 ALLOCATE (h_prod_test(0:l_max, 0:m_max))
328
329 ALLOCATE (e_c(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
330 ALLOCATE (e_h(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
331 CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, r1, r2, 1, e_c)
332 CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, r1, r2, 2, e_h)
333 CALL create_hermite_to_cartesian(zetp, l_max + m_max, h_to_c_ol)
334 CALL create_hermite_to_cartesian(zeta, l_max, h_to_c_1)
335 CALL create_hermite_to_cartesian(zetb, m_max, h_to_c_2)
336
337 DO t = 0, l_max + m_max
338 c_ol(t) = (r - rp)**t*exp(-zetp*(r - rp)**2)
339 END DO
340
341 DO l = 0, l_max
342 c1(l) = (r - r1)**l*exp(-zeta*(r - r1)**2)
343 END DO
344 DO m = 0, m_max
345 c2(m) = (r - r2)**m*exp(-zetb*(r - r2)**2)
346 END DO
347
348 h1(:) = matmul(transpose(h_to_c_1(0:, 0:)), c1)
349 h2(:) = matmul(transpose(h_to_c_2(0:, 0:)), c2)
350 h_ol(:) = matmul(transpose(h_to_c_ol(0:, 0:)), c_ol)
351
352 DO m = 0, m_max
353 DO l = 0, l_max
354 c_prod_ref(l, m) = c1(l)*c2(m)
355 h_prod_ref(l, m) = h1(l)*h2(m)
356 c_prod_test(l, m) = 0.0_dp
357 h_prod_test(l, m) = 0.0_dp
358 DO t = 0, l + m
359 c_prod_test(l, m) = c_prod_test(l, m) + e_c(t, l, m)*h_ol(t)
360 h_prod_test(l, m) = h_prod_test(l, m) + e_h(t, l, m)*h_ol(t)
361 END DO
362 END DO
363 END DO
364
365 c_prod_err = maxval(abs(c_prod_test - c_prod_ref)/(0.5_dp*(abs(c_prod_test) + abs(c_prod_ref)) + 1.0_dp))
366 h_prod_err = maxval(abs(h_prod_test - h_prod_ref)/(0.5_dp*(abs(c_prod_test) + abs(c_prod_ref)) + 1.0_dp))
367
368 cpassert(c_prod_err .LE. tolerance)
369 cpassert(h_prod_err .LE. tolerance)
370 mark_used(tolerance)
371
372 END SUBROUTINE overlap_dist_expansion_test
373
374END MODULE eri_mme_test
eri_mme_gaussian
Methods related to properties of Hermite and Cartesian Gaussian functions.
Definition eri_mme_gaussian.F:15
eri_mme_gaussian::create_hermite_to_cartesian
pure subroutine, public create_hermite_to_cartesian(zet, l_max, h_to_c)
Create matrix to transform between cartesian and hermite gaussian basis functions.
Definition eri_mme_gaussian.F:51
eri_mme_gaussian::eri_mme_longrange
integer, parameter, public eri_mme_longrange
Definition eri_mme_gaussian.F:29
eri_mme_gaussian::eri_mme_coulomb
integer, parameter, public eri_mme_coulomb
Definition eri_mme_gaussian.F:29
eri_mme_gaussian::create_gaussian_overlap_dist_to_hermite
pure subroutine, public create_gaussian_overlap_dist_to_hermite(l, m, a, b, r1, r2, h_or_c_product, e)
Expand 1d product of cartesian (or hermite) gaussians into single hermite gaussians: Find E_t^{lm} s....
Definition eri_mme_gaussian.F:194
eri_mme_gaussian::eri_mme_yukawa
integer, parameter, public eri_mme_yukawa
Definition eri_mme_gaussian.F:29
eri_mme_integrate
Minimax-Ewald (MME) method for calculating 2-center and 3-center electron repulsion integrals (ERI) o...
Definition eri_mme_integrate.F:23
eri_mme_integrate::eri_mme_3c_integrate
subroutine, public eri_mme_3c_integrate(param, la_min, la_max, lb_min, lb_max, lc_min, lc_max, zeta, zetb, zetc, ra, rb, rc, habc, o1, o2, o3, gg_count, gr_count, rr_count)
Low-level integration routine for 3-center ERIs.
Definition eri_mme_integrate.F:274
eri_mme_integrate::eri_mme_2c_integrate
subroutine, public eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, hab, o1, o2, g_count, r_count, normalize, potential, pot_par)
Low-level integration routine for 2-center ERIs.
Definition eri_mme_integrate.F:72
eri_mme_test
Methods for testing / debugging.
Definition eri_mme_test.F:14
eri_mme_test::eri_mme_3c_perf_acc_test
subroutine, public eri_mme_3c_perf_acc_test(param, l_max, zet, rabc, nrep, nsample, para_env, iw, potential, pot_par, gg_count, gr_count, rr_count)
...
Definition eri_mme_test.F:187
eri_mme_test::eri_mme_2c_perf_acc_test
subroutine, public eri_mme_2c_perf_acc_test(param, l_max, zet, rabc, nrep, test_accuracy, para_env, iw, potential, pot_par, g_count, r_count)
Unit test for performance and accuracy.
Definition eri_mme_test.F:60
eri_mme_types
Types and initialization / release routines for Minimax-Ewald method for electron repulsion integrals...
Definition eri_mme_types.F:16
eri_mme_types::eri_mme_set_potential
subroutine, public eri_mme_set_potential(param, potential, pot_par)
...
Definition eri_mme_types.F:262
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
eri_mme_types::eri_mme_param
Definition eri_mme_types.F:63
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763