Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 613086 : 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 .GE. MAXVAL(I_diff)
157 : ELSE
158 : acc_check = .TRUE.
159 : END IF
160 :
161 3 : IF (.NOT. acc_check) &
162 0 : CPABORT("Actual error greater than upper bound estimate.")
163 :
164 : END IF
165 : END IF
166 :
167 8 : 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 2 : 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 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: I_test
205 :
206 2 : IF (PRESENT(GG_count)) GG_count = 0
207 2 : IF (PRESENT(GR_count)) GR_count = 0
208 2 : IF (PRESENT(RR_count)) RR_count = 0
209 :
210 2 : ns = 1
211 2 : IF (PRESENT(nsample)) ns = nsample
212 :
213 2 : nzet = SIZE(zet)
214 2 : nR = SIZE(rabc, 2)
215 :
216 2 : IF (PRESENT(potential)) THEN
217 2 : CALL eri_mme_set_potential(param, potential, pot_par)
218 : END IF
219 :
220 2 : IF (param%debug) THEN
221 0 : DO izeta = 1, nzet
222 0 : DO izetb = 1, nzet
223 0 : DO ira = 1, nR
224 0 : DO irb = 1, nR
225 0 : DO ixyz = 1, 3
226 : CALL overlap_dist_expansion_test(l_max, l_max, zet(izeta), zet(izetb), &
227 0 : 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 2 : IF (iw > 0) THEN
236 1 : IF (PRESENT(potential)) THEN
237 2 : SELECT CASE (potential)
238 : CASE (eri_mme_coulomb)
239 1 : WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
240 : CASE (eri_mme_yukawa)
241 0 : WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: Yukawa with a=", pot_par
242 : CASE (eri_mme_longrange)
243 1 : WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: long-range Coulomb with a=", pot_par
244 : END SELECT
245 : ELSE
246 0 : WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
247 : END IF
248 1 : WRITE (iw, '(T2, A)') "ERI_MME| Test results for 3c cpu time"
249 1 : WRITE (iw, '(T11, A)') "la, lb, lc, zeta, zetb, zetc, cpu time"
250 : END IF
251 :
252 10 : ALLOCATE (I_test(ncoset(l_max), ncoset(l_max), ncoset(l_max)))
253 :
254 2 : nintg = 0
255 14 : DO la = 0, l_max
256 86 : DO lb = 0, l_max
257 516 : DO lc = 0, l_max
258 4824 : DO izeta = 1, nzet
259 47952 : DO izetb = 1, nzet
260 479520 : DO izetc = 1, nzet
261 432000 : nintg = nintg + 1
262 475200 : IF (MOD(nintg, ns) .EQ. 0) THEN
263 2145708 : I_test(:, :, :) = 0.0_dp
264 12 : CALL CPU_TIME(t0)
265 24 : DO irep = 1, nrep
266 120 : DO ira = 1, nR
267 876 : DO irb = 1, nR
268 7008 : 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 6912 : GG_count, GR_count, RR_count)
272 : END DO
273 : END DO
274 : END DO
275 : END DO
276 12 : CALL CPU_TIME(t1)
277 12 : time = t1 - t0
278 12 : CALL para_env%sum(time)
279 12 : 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 6 : 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 2 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: C1, C2, C_ol, H1, H2, H_ol
313 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: C_prod_ref, C_prod_test, H_prod_ref, &
314 0 : H_prod_test, h_to_c_1, h_to_c_2, &
315 0 : h_to_c_ol
316 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: E_C, E_H
317 :
318 0 : zetp = zeta + zetb
319 0 : Rp = (zeta*R1 + zetb*R2)/zetp
320 0 : ALLOCATE (C1(0:l_max), H1(0:l_max))
321 0 : ALLOCATE (C2(0:m_max), H2(0:m_max))
322 0 : ALLOCATE (C_ol(0:l_max + m_max))
323 0 : ALLOCATE (H_ol(0:l_max + m_max))
324 0 : ALLOCATE (C_prod_ref(0:l_max, 0:m_max))
325 0 : ALLOCATE (C_prod_test(0:l_max, 0:m_max))
326 0 : ALLOCATE (H_prod_ref(0:l_max, 0:m_max))
327 0 : ALLOCATE (H_prod_test(0:l_max, 0:m_max))
328 :
329 0 : ALLOCATE (E_C(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
330 0 : ALLOCATE (E_H(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
331 0 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, R1, R2, 1, E_C)
332 0 : CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, R1, R2, 2, E_H)
333 0 : 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 0 : DO t = 0, l_max + m_max
338 0 : C_ol(t) = (r - Rp)**t*EXP(-zetp*(r - Rp)**2)
339 : END DO
340 :
341 0 : DO l = 0, l_max
342 0 : C1(l) = (r - R1)**l*EXP(-zeta*(r - R1)**2)
343 : END DO
344 0 : DO m = 0, m_max
345 0 : C2(m) = (r - R2)**m*EXP(-zetb*(r - R2)**2)
346 : END DO
347 :
348 0 : H1(:) = MATMUL(TRANSPOSE(h_to_c_1(0:, 0:)), C1)
349 0 : H2(:) = MATMUL(TRANSPOSE(h_to_c_2(0:, 0:)), C2)
350 0 : H_ol(:) = MATMUL(TRANSPOSE(h_to_c_ol(0:, 0:)), C_ol)
351 :
352 0 : DO m = 0, m_max
353 0 : DO l = 0, l_max
354 0 : C_prod_ref(l, m) = C1(l)*C2(m)
355 0 : H_prod_ref(l, m) = H1(l)*H2(m)
356 0 : C_prod_test(l, m) = 0.0_dp
357 0 : H_prod_test(l, m) = 0.0_dp
358 0 : DO t = 0, l + m
359 0 : C_prod_test(l, m) = C_prod_test(l, m) + E_C(t, l, m)*H_ol(t)
360 0 : 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 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))
366 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))
367 :
368 0 : CPASSERT(C_prod_err .LE. tolerance)
369 0 : CPASSERT(H_prod_err .LE. tolerance)
370 : MARK_USED(tolerance)
371 :
372 0 : END SUBROUTINE overlap_dist_expansion_test
373 :
374 : END MODULE eri_mme_test
|
