(git:e7e05ae)
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 
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  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 
374 END 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_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::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_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_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_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_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