Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 :
14 : MODULE 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 :
42 : CONTAINS
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 8 : 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 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: time
77 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: I_diff, I_ref, I_test
78 : REAL(KIND=dp), DIMENSION(3, 3) :: ht
79 :
80 8 : IF (PRESENT(G_count)) G_count = 0
81 8 : IF (PRESENT(R_count)) R_count = 0
82 :
83 8 : nzet = SIZE(zet)
84 8 : nR = SIZE(rabc, 2)
85 :
86 8 : IF (PRESENT(potential)) THEN
87 8 : 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 8 : IF (test_accuracy) THEN
92 78 : ht = twopi*TRANSPOSE(param%h_inv)
93 :
94 36 : ALLOCATE (I_ref(ncoset(l_max), ncoset(l_max), nR, nzet))
95 6 : I_ref(:, :, :, :) = 0.0_dp
96 :
97 30 : DO izet = 1, nzet
98 222 : 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 216 : 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 48 : ALLOCATE (I_test(ncoset(l_max), ncoset(l_max), nR, nzet))
109 40 : ALLOCATE (I_diff(ncoset(l_max), ncoset(l_max), nR, nzet))
110 :
111 32 : ALLOCATE (time(0:l_max, nzet))
112 56 : DO l = 0, l_max
113 320 : DO izet = 1, nzet
114 264 : CALL CPU_TIME(t0)
115 528 : DO irep = 1, nrep
116 2640 : 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 2376 : normalize=.TRUE.)
121 : END DO
122 : END DO
123 264 : CALL CPU_TIME(t1)
124 312 : time(l, izet) = t1 - t0
125 : END DO
126 : END DO
127 :
128 8 : CALL para_env%sum(time)
129 :
130 8 : IF (test_accuracy) THEN
131 613086 : I_diff(:, :, :, :) = ABS(I_test - I_ref)
132 : END IF
133 :
134 8 : IF (iw > 0) THEN
135 4 : WRITE (iw, '(T2, A)') "ERI_MME| Test results for 2c cpu time"
136 4 : WRITE (iw, '(T11, A)') "l, zet, cpu time, accuracy"
137 :
138 28 : DO l = 0, l_max
139 160 : DO izet = 1, nzet
140 132 : IF (test_accuracy) THEN
141 83976 : acc = MAXVAL(I_diff(ncoset(l - 1) + 1:ncoset(l), ncoset(l - 1) + 1:ncoset(l), :, izet))
142 : ELSE
143 60 : acc = 0.0_dp
144 : END IF
145 :
146 : WRITE (iw, '(T11, I1, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2)') &
147 156 : l, zet(izet), time(l, izet)/nrep, acc
148 : END DO
149 : END DO
150 :
151 4 : IF (test_accuracy) THEN
152 3 : WRITE (iw, '(/T2, A, 47X, ES9.2)') "ERI_MME| Maximum error:", &
153 306546 : MAXVAL(I_diff)
154 :
155 3 : IF (param%is_ortho) THEN
156 306543 : acc_check = param%err_mm + param%err_c >= MAXVAL(I_diff)
157 : ELSE
158 : acc_check = .TRUE.
159 : END IF
160 :
161 3 : IF (.NOT. acc_check) THEN
162 0 : CPABORT("Actual error greater than upper bound estimate.")
163 : END IF
164 :
165 : END IF
166 : END IF
167 :
168 8 : END SUBROUTINE eri_mme_2c_perf_acc_test
169 :
170 : ! **************************************************************************************************
171 : !> \brief ...
172 : !> \param param ...
173 : !> \param l_max ...
174 : !> \param zet ...
175 : !> \param rabc ...
176 : !> \param nrep ...
177 : !> \param nsample ...
178 : !> \param para_env ...
179 : !> \param iw ...
180 : !> \param potential ...
181 : !> \param pot_par ...
182 : !> \param GG_count ...
183 : !> \param GR_count ...
184 : !> \param RR_count ...
185 : ! **************************************************************************************************
186 2 : SUBROUTINE eri_mme_3c_perf_acc_test(param, l_max, zet, rabc, nrep, nsample, &
187 : para_env, iw, potential, pot_par, GG_count, GR_count, RR_count)
188 : TYPE(eri_mme_param), INTENT(INOUT) :: param
189 : INTEGER, INTENT(IN) :: l_max
190 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
191 : INTENT(IN) :: zet
192 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
193 : INTENT(IN) :: rabc
194 : INTEGER, INTENT(IN) :: nrep
195 : INTEGER, INTENT(IN), OPTIONAL :: nsample
196 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
197 : INTEGER, INTENT(IN) :: iw
198 : INTEGER, INTENT(IN), OPTIONAL :: potential
199 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par
200 : INTEGER, INTENT(OUT), OPTIONAL :: GG_count, GR_count, RR_count
201 :
202 : INTEGER :: ira, irb, irc, irep, ixyz, izeta, izetb, &
203 : izetc, la, lb, lc, nintg, nR, ns, nzet
204 : REAL(KIND=dp) :: t0, t1, time
205 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: I_test
206 :
207 2 : IF (PRESENT(GG_count)) GG_count = 0
208 2 : IF (PRESENT(GR_count)) GR_count = 0
209 2 : IF (PRESENT(RR_count)) RR_count = 0
210 :
211 2 : ns = 1
212 2 : IF (PRESENT(nsample)) ns = nsample
213 :
214 2 : nzet = SIZE(zet)
215 2 : nR = SIZE(rabc, 2)
216 :
217 2 : IF (PRESENT(potential)) THEN
218 2 : CALL eri_mme_set_potential(param, potential, pot_par)
219 : END IF
220 :
221 2 : IF (param%debug) THEN
222 0 : DO izeta = 1, nzet
223 0 : DO izetb = 1, nzet
224 0 : DO ira = 1, nR
225 0 : DO irb = 1, nR
226 0 : DO ixyz = 1, 3
227 : CALL overlap_dist_expansion_test(l_max, l_max, zet(izeta), zet(izetb), &
228 0 : rabc(ixyz, ira), rabc(ixyz, irb), 0.0_dp, param%debug_delta)
229 : END DO
230 : END DO
231 : END DO
232 : END DO
233 : END DO
234 : END IF
235 :
236 2 : IF (iw > 0) THEN
237 1 : IF (PRESENT(potential)) THEN
238 2 : SELECT CASE (potential)
239 : CASE (eri_mme_coulomb)
240 1 : WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
241 : CASE (eri_mme_yukawa)
242 0 : WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: Yukawa with a=", pot_par
243 : CASE (eri_mme_longrange)
244 1 : WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: long-range Coulomb with a=", pot_par
245 : END SELECT
246 : ELSE
247 0 : WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
248 : END IF
249 1 : WRITE (iw, '(T2, A)') "ERI_MME| Test results for 3c cpu time"
250 1 : WRITE (iw, '(T11, A)') "la, lb, lc, zeta, zetb, zetc, cpu time"
251 : END IF
252 :
253 10 : ALLOCATE (I_test(ncoset(l_max), ncoset(l_max), ncoset(l_max)))
254 :
255 2 : nintg = 0
256 14 : DO la = 0, l_max
257 86 : DO lb = 0, l_max
258 516 : DO lc = 0, l_max
259 4824 : DO izeta = 1, nzet
260 47952 : DO izetb = 1, nzet
261 479520 : DO izetc = 1, nzet
262 432000 : nintg = nintg + 1
263 475200 : IF (MOD(nintg, ns) == 0) THEN
264 12 : I_test(:, :, :) = 0.0_dp
265 12 : CALL CPU_TIME(t0)
266 24 : DO irep = 1, nrep
267 120 : DO ira = 1, nR
268 876 : DO irb = 1, nR
269 7008 : DO irc = 1, nR
270 : CALL eri_mme_3c_integrate(param, 0, la, 0, lb, 0, lc, zet(izeta), zet(izetb), zet(izetc), &
271 : rabc(:, ira), rabc(:, irb), rabc(:, irc), I_test, 0, 0, 0, &
272 6912 : GG_count, GR_count, RR_count)
273 : END DO
274 : END DO
275 : END DO
276 : END DO
277 12 : CALL CPU_TIME(t1)
278 12 : time = t1 - t0
279 12 : CALL para_env%sum(time)
280 12 : IF (iw > 0) THEN
281 : WRITE (iw, '(T11, I1, 1X, I1, 1X, I1, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2)') &
282 6 : la, lb, lc, zet(izeta), zet(izetb), zet(izetc), time/nrep
283 : END IF
284 : END IF
285 : END DO
286 : END DO
287 : END DO
288 : END DO
289 : END DO
290 : END DO
291 :
292 2 : END SUBROUTINE eri_mme_3c_perf_acc_test
293 :
294 : ! **************************************************************************************************
295 : !> \brief check that expanding an overlap distribution of cartesian/hermite Gaussians into a
296 : !> lin combi of single cartesian/hermite Gaussians is correct.
297 : !> \param l_max ...
298 : !> \param m_max ...
299 : !> \param zeta ...
300 : !> \param zetb ...
301 : !> \param R1 ...
302 : !> \param R2 ...
303 : !> \param r ...
304 : !> \param tolerance ...
305 : !> \note STATUS: tested
306 : ! **************************************************************************************************
307 0 : SUBROUTINE overlap_dist_expansion_test(l_max, m_max, zeta, zetb, R1, R2, r, tolerance)
308 : INTEGER, INTENT(IN) :: l_max, m_max
309 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb, R1, R2, r, tolerance
310 :
311 : INTEGER :: l, m, t
312 : REAL(KIND=dp) :: C_prod_err, H_prod_err, Rp, zetp
313 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: C1, C2, C_ol, H1, H2, H_ol
314 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: C_prod_ref, C_prod_test, H_prod_ref, &
315 0 : H_prod_test, h_to_c_1, h_to_c_2, &
316 0 : h_to_c_ol
317 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: E_C, E_H
318 :
319 0 : zetp = zeta + zetb
320 0 : Rp = (zeta*R1 + zetb*R2)/zetp
321 0 : ALLOCATE (C1(0:l_max), H1(0:l_max))
322 0 : ALLOCATE (C2(0:m_max), H2(0:m_max))
323 0 : ALLOCATE (C_ol(0:l_max + m_max))
324 0 : ALLOCATE (H_ol(0:l_max + m_max))
325 0 : ALLOCATE (C_prod_ref(0:l_max, 0:m_max))
326 0 : ALLOCATE (C_prod_test(0:l_max, 0:m_max))
327 0 : ALLOCATE (H_prod_ref(0:l_max, 0:m_max))
328 0 : ALLOCATE (H_prod_test(0:l_max, 0:m_max))
329 :
330 0 : ALLOCATE (E_C(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
331 0 : ALLOCATE (E_H(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
332 0 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, R1, R2, 1, E_C)
333 0 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, R1, R2, 2, E_H)
334 0 : CALL create_hermite_to_cartesian(zetp, l_max + m_max, h_to_c_ol)
335 0 : CALL create_hermite_to_cartesian(zeta, l_max, h_to_c_1)
336 0 : CALL create_hermite_to_cartesian(zetb, m_max, h_to_c_2)
337 :
338 0 : DO t = 0, l_max + m_max
339 0 : C_ol(t) = (r - Rp)**t*EXP(-zetp*(r - Rp)**2)
340 : END DO
341 :
342 0 : DO l = 0, l_max
343 0 : C1(l) = (r - R1)**l*EXP(-zeta*(r - R1)**2)
344 : END DO
345 0 : DO m = 0, m_max
346 0 : C2(m) = (r - R2)**m*EXP(-zetb*(r - R2)**2)
347 : END DO
348 :
349 0 : H1(:) = MATMUL(TRANSPOSE(h_to_c_1(0:, 0:)), C1)
350 0 : H2(:) = MATMUL(TRANSPOSE(h_to_c_2(0:, 0:)), C2)
351 0 : H_ol(:) = MATMUL(TRANSPOSE(h_to_c_ol(0:, 0:)), C_ol)
352 :
353 0 : DO m = 0, m_max
354 0 : DO l = 0, l_max
355 0 : C_prod_ref(l, m) = C1(l)*C2(m)
356 0 : H_prod_ref(l, m) = H1(l)*H2(m)
357 0 : C_prod_test(l, m) = 0.0_dp
358 0 : H_prod_test(l, m) = 0.0_dp
359 0 : DO t = 0, l + m
360 0 : C_prod_test(l, m) = C_prod_test(l, m) + E_C(t, l, m)*H_ol(t)
361 0 : H_prod_test(l, m) = H_prod_test(l, m) + E_H(t, l, m)*H_ol(t)
362 : END DO
363 : END DO
364 : END DO
365 :
366 0 : 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))
367 0 : 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))
368 :
369 0 : CPASSERT(C_prod_err <= tolerance)
370 0 : CPASSERT(H_prod_err <= tolerance)
371 : MARK_USED(tolerance)
372 :
373 0 : END SUBROUTINE overlap_dist_expansion_test
374 :
375 : END MODULE eri_mme_test
|
