1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5! **************************************************************************************************
6!> \brief Methods for testing / debugging.
7!> \par History
8!>       2015 09 created
9!> \author Patrick Seewald
10! **************************************************************************************************
11
12MODULE eri_mme_test
13
14   USE cp_para_types,                   ONLY: cp_para_env_type
15   USE eri_mme_gaussian,                ONLY: create_gaussian_overlap_dist_to_hermite,&
16                                              create_hermite_to_cartesian
17   USE eri_mme_integrate,               ONLY: eri_mme_2c_integrate,&
18                                              eri_mme_3c_integrate
19   USE eri_mme_types,                   ONLY: eri_mme_coulomb,&
20                                              eri_mme_longrange,&
21                                              eri_mme_param,&
22                                              eri_mme_set_potential,&
23                                              eri_mme_yukawa
24   USE kinds,                           ONLY: dp
25   USE mathconstants,                   ONLY: twopi
26   USE message_passing,                 ONLY: mp_sum
27   USE orbital_pointers,                ONLY: ncoset
28#include "../base/base_uses.f90"
29
30   IMPLICIT NONE
31
32   PRIVATE
33
34   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
35
36   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_test'
37
38   PUBLIC :: eri_mme_2c_perf_acc_test, &
39             eri_mme_3c_perf_acc_test
40
41CONTAINS
42! **************************************************************************************************
43!> \brief Unit test for performance and accuracy
44!> \param param ...
45!> \param l_max ...
46!> \param zet ...
47!> \param rabc ...
48!> \param nrep ...
49!> \param test_accuracy ...
50!> \param para_env ...
51!> \param iw ...
52!> \param potential ...
53!> \param pot_par ...
54!> \param G_count ...
55!> \param R_count ...
56! **************************************************************************************************
57   SUBROUTINE eri_mme_2c_perf_acc_test(param, l_max, zet, rabc, nrep, test_accuracy, &
58                                       para_env, iw, potential, pot_par, G_count, R_count)
59      TYPE(eri_mme_param), INTENT(INOUT)                 :: param
60      INTEGER, INTENT(IN)                                :: l_max
61      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
62         INTENT(IN)                                      :: zet
63      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: rabc
64      INTEGER, INTENT(IN)                                :: nrep
65      LOGICAL, INTENT(INOUT)                             :: test_accuracy
66      TYPE(cp_para_env_type), INTENT(IN), POINTER        :: para_env
67      INTEGER, INTENT(IN)                                :: iw
68      INTEGER, INTENT(IN), OPTIONAL                      :: potential
69      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pot_par
70      INTEGER, INTENT(OUT), OPTIONAL                     :: G_count, R_count
71
72      CHARACTER(len=*), PARAMETER :: routineN = 'eri_mme_2c_perf_acc_test', &
73         routineP = moduleN//':'//routineN
74
75      INTEGER                                            :: iab, irep, izet, l, nR, nzet
76      LOGICAL                                            :: acc_check
77      REAL(KIND=dp)                                      :: acc, t0, t1
78      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: time
79      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: I_diff, I_ref, I_test
80      REAL(KIND=dp), DIMENSION(3, 3)                     :: ht
81
82      IF (PRESENT(G_count)) G_count = 0
83      IF (PRESENT(R_count)) R_count = 0
84
85      nzet = SIZE(zet)
86      nR = SIZE(rabc, 2)
87
88      IF (PRESENT(potential)) THEN
89         CALL eri_mme_set_potential(param, potential, pot_par)
90      ENDIF
91
92      ! Calculate reference values (Exact expression in G space converged to high precision)
93      IF (test_accuracy) THEN
94         ht = twopi*TRANSPOSE(param%h_inv)
95
96         ALLOCATE (I_ref(ncoset(l_max), ncoset(l_max), nR, nzet))
97         I_ref(:, :, :, :) = 0.0_dp
98
99         DO izet = 1, nzet
100            DO iab = 1, nR
101               CALL eri_mme_2c_integrate(param, 0, l_max, 0, l_max, zet(izet), zet(izet), rabc(:, iab), &
102                                         I_ref(:, :, iab, izet), 0, 0, &
103                                         normalize=.TRUE., potential=potential, pot_par=pot_par)
104
105            ENDDO
106         ENDDO
107      ENDIF
108
109      ! test performance and accuracy of MME method
110      ALLOCATE (I_test(ncoset(l_max), ncoset(l_max), nR, nzet))
111      ALLOCATE (I_diff(ncoset(l_max), ncoset(l_max), nR, nzet))
112
113      ALLOCATE (time(0:l_max, nzet))
114      DO l = 0, l_max
115         DO izet = 1, nzet
116            CALL CPU_TIME(t0)
117            DO irep = 1, nrep
118               DO iab = 1, nR
119                  CALL eri_mme_2c_integrate(param, 0, l, 0, l, zet(izet), zet(izet), rabc(:, iab), &
120                                            I_test(:, :, iab, izet), 0, 0, &
121                                            G_count=G_count, R_count=R_count, &
122                                            normalize=.TRUE.)
123               ENDDO
124            ENDDO
125            CALL CPU_TIME(t1)
126            time(l, izet) = t1 - t0
127         ENDDO
128      ENDDO
129
130      CALL mp_sum(time, para_env%group)
131
132      IF (test_accuracy) THEN
133         I_diff(:, :, :, :) = ABS(I_test - I_ref)
134      ENDIF
135
136      IF (iw > 0) THEN
137         WRITE (iw, '(T2, A)') "ERI_MME| Test results for 2c cpu time"
138         WRITE (iw, '(T11, A)') "l, zet, cpu time, accuracy"
139
140         DO l = 0, l_max
141            DO izet = 1, nzet
142               IF (test_accuracy) THEN
143                  acc = MAXVAL(I_diff(ncoset(l - 1) + 1:ncoset(l), ncoset(l - 1) + 1:ncoset(l), :, izet))
144               ELSE
145                  acc = 0.0_dp
146               ENDIF
147
148               WRITE (iw, '(T11, I1, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2)') &
149                  l, zet(izet), time(l, izet)/nrep, acc
150            ENDDO
151         ENDDO
152
153         IF (test_accuracy) THEN
154            WRITE (iw, '(/T2, A, 47X, ES9.2)') "ERI_MME| Maximum error:", &
155               MAXVAL(I_diff)
156
157            IF (param%is_ortho) THEN
158               acc_check = param%err_mm + param%err_c .GE. MAXVAL(I_diff)
159            ELSE
160               acc_check = .TRUE.
161            ENDIF
162
163            IF (.NOT. acc_check) &
164               CPABORT("Actual error greater than upper bound estimate.")
165
166         ENDIF
167      ENDIF
168
169   END SUBROUTINE eri_mme_2c_perf_acc_test
170
171! **************************************************************************************************
172!> \brief ...
173!> \param param ...
174!> \param l_max ...
175!> \param zet ...
176!> \param rabc ...
177!> \param nrep ...
178!> \param nsample ...
179!> \param para_env ...
180!> \param iw ...
181!> \param potential ...
182!> \param pot_par ...
183!> \param GG_count ...
184!> \param GR_count ...
185!> \param RR_count ...
186! **************************************************************************************************
187   SUBROUTINE eri_mme_3c_perf_acc_test(param, l_max, zet, rabc, nrep, nsample, &
188                                       para_env, iw, potential, pot_par, GG_count, GR_count, RR_count)
189      TYPE(eri_mme_param), INTENT(INOUT)                 :: param
190      INTEGER, INTENT(IN)                                :: l_max
191      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
192         INTENT(IN)                                      :: zet
193      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
194         INTENT(IN)                                      :: rabc
195      INTEGER, INTENT(IN)                                :: nrep
196      INTEGER, INTENT(IN), OPTIONAL                      :: nsample
197      TYPE(cp_para_env_type), INTENT(IN), POINTER        :: para_env
198      INTEGER, INTENT(IN)                                :: iw
199      INTEGER, INTENT(IN), OPTIONAL                      :: potential
200      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pot_par
201      INTEGER, INTENT(OUT), OPTIONAL                     :: GG_count, GR_count, RR_count
202
203      CHARACTER(len=*), PARAMETER :: routineN = 'eri_mme_3c_perf_acc_test', &
204         routineP = moduleN//':'//routineN
205
206      INTEGER                                            :: ira, irb, irc, irep, ixyz, izeta, izetb, &
207                                                            izetc, la, lb, lc, nintg, nR, ns, nzet
208      REAL(KIND=dp)                                      :: t0, t1, time
209      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: I_test
210
211      IF (PRESENT(GG_count)) GG_count = 0
212      IF (PRESENT(GR_count)) GR_count = 0
213      IF (PRESENT(RR_count)) RR_count = 0
214
215      ns = 1
216      IF (PRESENT(nsample)) ns = nsample
217
218      nzet = SIZE(zet)
219      nR = SIZE(rabc, 2)
220
221      IF (PRESENT(potential)) THEN
222         CALL eri_mme_set_potential(param, potential, pot_par)
223      ENDIF
224
225      IF (param%debug) THEN
226         DO izeta = 1, nzet
227         DO izetb = 1, nzet
228            DO ira = 1, nR
229            DO irb = 1, nR
230               DO ixyz = 1, 3
231                  CALL overlap_dist_expansion_test(l_max, l_max, zet(izeta), zet(izetb), &
232                                                   rabc(ixyz, ira), rabc(ixyz, irb), 0.0_dp, param%debug_delta)
233               ENDDO
234            ENDDO
235            ENDDO
236         ENDDO
237         ENDDO
238      ENDIF
239
240      IF (iw > 0) THEN
241         IF (PRESENT(potential)) THEN
242            SELECT CASE (potential)
243            CASE (eri_mme_coulomb)
244               WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
245            CASE (eri_mme_yukawa)
246               WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: Yukawa with a=", pot_par
247            CASE (eri_mme_longrange)
248               WRITE (iw, '(/T2, A, ES9.2)') "ERI_MME| Potential: long-range Coulomb with a=", pot_par
249            END SELECT
250         ELSE
251            WRITE (iw, '(/T2, A)') "ERI_MME| Potential: Coulomb"
252         ENDIF
253         WRITE (iw, '(T2, A)') "ERI_MME| Test results for 3c cpu time"
254         WRITE (iw, '(T11, A)') "la, lb, lc, zeta, zetb, zetc, cpu time"
255      ENDIF
256
257      ALLOCATE (I_test(ncoset(l_max), ncoset(l_max), ncoset(l_max)))
258
259      nintg = 0
260      DO la = 0, l_max
261      DO lb = 0, l_max
262      DO lc = 0, l_max
263         DO izeta = 1, nzet
264         DO izetb = 1, nzet
265         DO izetc = 1, nzet
266            nintg = nintg + 1
267            IF (MOD(nintg, ns) .EQ. 0) THEN
268               I_test(:, :, :) = 0.0_dp
269               CALL CPU_TIME(t0)
270               DO irep = 1, nrep
271                  DO ira = 1, nR
272                  DO irb = 1, nR
273                  DO irc = 1, nR
274                     CALL eri_mme_3c_integrate(param, 0, la, 0, lb, 0, lc, zet(izeta), zet(izetb), zet(izetc), &
275                                               rabc(:, ira), rabc(:, irb), rabc(:, irc), I_test, 0, 0, 0, &
276                                               GG_count, GR_count, RR_count)
277                  ENDDO
278                  ENDDO
279                  ENDDO
280               ENDDO
281               CALL CPU_TIME(t1)
282               time = t1 - t0
283               CALL mp_sum(time, para_env%group)
284               IF (iw > 0) THEN
285                  WRITE (iw, '(T11, I1, 1X, I1, 1X, I1, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2, 1X, ES9.2)') &
286                     la, lb, lc, zet(izeta), zet(izetb), zet(izetc), time/nrep
287               ENDIF
288            ENDIF
289         ENDDO
290         ENDDO
291         ENDDO
292      ENDDO
293      ENDDO
294      ENDDO
295
296   END SUBROUTINE eri_mme_3c_perf_acc_test
297
298! **************************************************************************************************
299!> \brief check that expanding an overlap distribution of cartesian/hermite Gaussians into a
300!>        lin combi of single cartesian/hermite Gaussians is correct.
301!> \param l_max ...
302!> \param m_max ...
303!> \param zeta ...
304!> \param zetb ...
305!> \param R1 ...
306!> \param R2 ...
307!> \param r ...
308!> \param tolerance ...
309!> \note STATUS: tested
310! **************************************************************************************************
311   SUBROUTINE overlap_dist_expansion_test(l_max, m_max, zeta, zetb, R1, R2, r, tolerance)
312      INTEGER, INTENT(IN)                                :: l_max, m_max
313      REAL(KIND=dp), INTENT(IN)                          :: zeta, zetb, R1, R2, r, tolerance
314
315      INTEGER                                            :: l, m, t
316      REAL(KIND=dp)                                      :: C_prod_err, H_prod_err, Rp, zetp
317      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: C1, C2, C_ol, H1, H2, H_ol
318      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: C_prod_ref, C_prod_test, H_prod_ref, &
319                                                            H_prod_test, h_to_c_1, h_to_c_2, &
320                                                            h_to_c_ol
321      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: E_C, E_H
322
323      zetp = zeta + zetb
324      Rp = (zeta*R1 + zetb*R2)/zetp
325      ALLOCATE (C1(0:l_max), H1(0:l_max))
326      ALLOCATE (C2(0:m_max), H2(0:m_max))
327      ALLOCATE (C_ol(0:l_max + m_max))
328      ALLOCATE (H_ol(0:l_max + m_max))
329      ALLOCATE (C_prod_ref(0:l_max, 0:m_max))
330      ALLOCATE (C_prod_test(0:l_max, 0:m_max))
331      ALLOCATE (H_prod_ref(0:l_max, 0:m_max))
332      ALLOCATE (H_prod_test(0:l_max, 0:m_max))
333
334      ALLOCATE (E_C(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
335      ALLOCATE (E_H(-1:l_max + m_max + 1, -1:l_max, -1:m_max))
336      CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, R1, R2, 1, E_C)
337      CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, R1, R2, 2, E_H)
338      CALL create_hermite_to_cartesian(zetp, l_max + m_max, h_to_c_ol)
339      CALL create_hermite_to_cartesian(zeta, l_max, h_to_c_1)
340      CALL create_hermite_to_cartesian(zetb, m_max, h_to_c_2)
341
342      DO t = 0, l_max + m_max
343         C_ol(t) = (r - Rp)**t*EXP(-zetp*(r - Rp)**2)
344      ENDDO
345
346      DO l = 0, l_max
347         C1(l) = (r - R1)**l*EXP(-zeta*(r - R1)**2)
348      ENDDO
349      DO m = 0, m_max
350         C2(m) = (r - R2)**m*EXP(-zetb*(r - R2)**2)
351      ENDDO
352
353      H1(:) = MATMUL(TRANSPOSE(h_to_c_1(0:, 0:)), C1)
354      H2(:) = MATMUL(TRANSPOSE(h_to_c_2(0:, 0:)), C2)
355      H_ol(:) = MATMUL(TRANSPOSE(h_to_c_ol(0:, 0:)), C_ol)
356
357      DO m = 0, m_max
358         DO l = 0, l_max
359            C_prod_ref(l, m) = C1(l)*C2(m)
360            H_prod_ref(l, m) = H1(l)*H2(m)
361            C_prod_test(l, m) = 0.0_dp
362            H_prod_test(l, m) = 0.0_dp
363            DO t = 0, l + m
364               C_prod_test(l, m) = C_prod_test(l, m) + E_C(t, l, m)*H_ol(t)
365               H_prod_test(l, m) = H_prod_test(l, m) + E_H(t, l, m)*H_ol(t)
366            ENDDO
367         ENDDO
368      ENDDO
369
370      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))
371      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))
372
373      CPASSERT(C_prod_err .LE. tolerance)
374      CPASSERT(H_prod_err .LE. tolerance)
375
376   END SUBROUTINE overlap_dist_expansion_test
377
378END MODULE eri_mme_test
379