Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2022 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 cp_para_types, ONLY: cp_para_env_type
17 : USE eri_mme_gaussian, ONLY: create_gaussian_overlap_dist_to_hermite,&
18 : create_hermite_to_cartesian
19 : USE eri_mme_integrate, ONLY: eri_mme_2c_integrate,&
20 : eri_mme_3c_integrate
21 : USE eri_mme_types, ONLY: eri_mme_coulomb,&
22 : eri_mme_longrange,&
23 : eri_mme_param,&
24 : eri_mme_set_potential,&
25 : eri_mme_yukawa
26 : USE kinds, ONLY: dp
27 : USE mathconstants, ONLY: twopi
28 : USE message_passing, ONLY: mp_sum
29 : USE orbital_pointers, ONLY: ncoset
30 : #include "../base/base_uses.f90"
31 :
32 : IMPLICIT NONE
33 :
34 : PRIVATE
35 :
36 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
37 :
38 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_test'
39 :
40 : PUBLIC :: eri_mme_2c_perf_acc_test, &
41 : eri_mme_3c_perf_acc_test
42 :
43 : CONTAINS
44 : ! **************************************************************************************************
45 : !> \brief Unit test for performance and accuracy
46 : !> \param param ...
47 : !> \param l_max ...
48 : !> \param zet ...
49 : !> \param rabc ...
50 : !> \param nrep ...
51 : !> \param test_accuracy ...
52 : !> \param para_env ...
53 : !> \param iw ...
54 : !> \param potential ...
55 : !> \param pot_par ...
56 : !> \param G_count ...
57 : !> \param R_count ...
58 : ! **************************************************************************************************
59 8 : SUBROUTINE eri_mme_2c_perf_acc_test(param, l_max, zet, rabc, nrep, test_accuracy, &
60 : para_env, iw, potential, pot_par, G_count, R_count)
61 : TYPE(eri_mme_param), INTENT(INOUT) :: param
62 : INTEGER, INTENT(IN) :: l_max
63 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
64 : INTENT(IN) :: zet
65 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rabc
66 : INTEGER, INTENT(IN) :: nrep
67 : LOGICAL, INTENT(INOUT) :: test_accuracy
68 : TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env
69 : INTEGER, INTENT(IN) :: iw
70 : INTEGER, INTENT(IN), OPTIONAL :: potential
71 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par
72 : INTEGER, INTENT(OUT), OPTIONAL :: G_count, R_count
73 :
74 : INTEGER :: iab, irep, izet, l, nR, nzet
75 : LOGICAL :: acc_check
76 : REAL(KIND=dp) :: acc, t0, t1
77 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: time
78 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: I_diff, I_ref, I_test
79 : REAL(KIND=dp), DIMENSION(3, 3) :: ht
80 :
81 8 : IF (PRESENT(G_count)) G_count = 0
82 8 : IF (PRESENT(R_count)) R_count = 0
83 :
84 8 : nzet = SIZE(zet)
85 8 : nR = SIZE(rabc, 2)
86 :
87 8 : IF (PRESENT(potential)) THEN
88 8 : CALL eri_mme_set_potential(param, potential, pot_par)
89 : END IF
90 :
91 : ! Calculate reference values (Exact expression in G space converged to high precision)
92 8 : IF (test_accuracy) THEN
93 78 : ht = twopi*TRANSPOSE(param%h_inv)
94 :
95 36 : ALLOCATE (I_ref(ncoset(l_max), ncoset(l_max), nR, nzet))
96 613086 : I_ref(:, :, :, :) = 0.0_dp
97 :
98 30 : DO izet = 1, nzet
99 222 : DO iab = 1, nR
100 : CALL eri_mme_2c_integrate(param, 0, l_max, 0, l_max, zet(izet), zet(izet), rabc(:, iab), &
101 : I_ref(:, :, iab, izet), 0, 0, &
102 216 : normalize=.TRUE., potential=potential, pot_par=pot_par)
103 :
104 : END DO
105 : END DO
106 : END IF
107 :
108 : ! test performance and accuracy of MME method
109 48 : ALLOCATE (I_test(ncoset(l_max), ncoset(l_max), nR, nzet))
110 48 : ALLOCATE (I_diff(ncoset(l_max), ncoset(l_max), nR, nzet))
111 :
112 32 : ALLOCATE (time(0:l_max, nzet))
113 56 : DO l = 0, l_max
114 320 : DO izet = 1, nzet
115 264 : CALL CPU_TIME(t0)
116 528 : DO irep = 1, nrep
117 2640 : DO iab = 1, nR
118 : CALL eri_mme_2c_integrate(param, 0, l, 0, l, zet(izet), zet(izet), rabc(:, iab), &
119 : I_test(:, :, iab, izet), 0, 0, &
120 : G_count=G_count, R_count=R_count, &
121 2376 : normalize=.TRUE.)
122 : END DO
123 : END DO
124 264 : CALL CPU_TIME(t1)
125 312 : time(l, izet) = t1 - t0
126 : END DO
127 : END DO
128 :
129 8 : CALL mp_sum(time, para_env%group)
130 :
131 8 : IF (test_accuracy) THEN
132 613086 : I_diff(:, :, :, :) = ABS(I_test - I_ref)
133 : END IF
134 :
135 8 : IF (iw > 0) THEN
136 4 : WRITE (iw, '(T2, A)') "ERI_MME| Test results for 2c cpu time"
137 4 : WRITE (iw, '(T11, A)') "l, zet, cpu time, accuracy"
138 :
139 28 : DO l = 0, l_max
140 160 : DO izet = 1, nzet
141 132 : IF (test_accuracy) THEN
142 83976 : acc = MAXVAL(I_diff(ncoset(l - 1) + 1:ncoset(l), ncoset(l - 1) + 1:ncoset(l), :, izet))
143 : ELSE
144 60 : acc = 0.0_dp
145 : END IF
146 :
147 : WRITE (iw, '(T11, I1, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2)') &
148 156 : l, zet(izet), time(l, izet)/nrep, acc
149 : END DO
150 : END DO
151 :
152 4 : IF (test_accuracy) THEN
153 3 : WRITE (iw, '(/T2, A, 47X, ES9.2)') "ERI_MME| Maximum error:", &
154 306546 : MAXVAL(I_diff)
155 :
156 3 : IF (param%is_ortho) THEN
157 306543 : acc_check = param%err_mm + param%err_c .GE. MAXVAL(I_diff)
158 : ELSE
159 : acc_check = .TRUE.
160 : END IF
161 :
162 3 : IF (.NOT. acc_check) &
163 0 : CPABORT("Actual error greater than upper bound estimate.")
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(cp_para_env_type), INTENT(IN), POINTER :: 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 2 : 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) .EQ. 0) THEN
264 2145708 : 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 mp_sum(time, para_env%group)
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 : CALL create_hermite_to_cartesian(zeta, l_max, h_to_c_1)
336 : 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 .LE. tolerance)
370 0 : CPASSERT(H_prod_err .LE. tolerance)
371 : MARK_USED(tolerance)
372 :
373 0 : END SUBROUTINE overlap_dist_expansion_test
374 :
375 : END MODULE eri_mme_test
|