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 : ! **************************************************************************************************
9 : !> \brief Debug the derivatives of the the rotational matrices
10 : !>
11 : !> \param sepi ...
12 : !> \param sepj ...
13 : !> \param rjiv ...
14 : !> \param ij_matrix ...
15 : !> \param do_invert ...
16 : !> \date 04.2008 [tlaino]
17 : !> \author Teodoro Laino [tlaino] - University of Zurich
18 : ! **************************************************************************************************
19 0 : SUBROUTINE check_rotmat_der(sepi, sepj, rjiv, ij_matrix, do_invert)
20 :
21 : USE kinds, ONLY: dp
22 : USE semi_empirical_int_utils, ONLY: rotmat
23 : USE semi_empirical_types, ONLY: rotmat_create, &
24 : rotmat_release, &
25 : rotmat_type, &
26 : semi_empirical_type, &
27 : se_int_control_type
28 : #include "./base/base_uses.f90"
29 : IMPLICIT NONE
30 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
31 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rjiv
32 : TYPE(rotmat_type), POINTER :: ij_matrix
33 : LOGICAL, INTENT(IN) :: do_invert
34 :
35 : CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
36 : routineN = 'check_rotmat_der', routineP = moduleN//':'//routineN
37 :
38 : REAL(KIND=dp) :: dx, r, r0(3), x(3)
39 : TYPE(rotmat_type), POINTER :: matrix, matrix_m, matrix_n, &
40 : matrix_p
41 : INTEGER :: imap(3), i, j, k, l
42 :
43 : INTERFACE
44 : FUNCTION check_value(num, ana, minval, thrs) RESULT(passed)
45 : USE kinds, ONLY: dp
46 : IMPLICIT NONE
47 : REAL(KIND=dp) :: num, ana, thrs, minval
48 : LOGICAL :: passed
49 : END FUNCTION check_value
50 : END INTERFACE
51 :
52 0 : NULLIFY (matrix_m, matrix_p, matrix_n, matrix)
53 0 : CALL rotmat_create(matrix_p)
54 0 : CALL rotmat_create(matrix_m)
55 0 : CALL rotmat_create(matrix_n)
56 0 : dx = 1.0E-6_dp
57 0 : imap(1) = 1
58 0 : imap(2) = 2
59 0 : imap(3) = 3
60 0 : IF (do_invert) THEN
61 0 : imap(1) = 3
62 : imap(2) = 2
63 0 : imap(3) = 1
64 : END IF
65 : ! Check derivatives: analytical VS numerical
66 0 : WRITE (*, *) "DEBUG::"//routineP
67 0 : DO j = 1, 3
68 0 : x = 0.0_dp
69 0 : x(imap(j)) = dx
70 0 : DO i = 1, 2
71 0 : IF (i == 1) matrix => matrix_p
72 0 : IF (i == 2) matrix => matrix_m
73 0 : r0 = rjiv + (-1.0_dp)**(i - 1)*x
74 0 : r = SQRT(DOT_PRODUCT(r0, r0))
75 0 : CALL rotmat(sepi, sepj, r0, r, matrix, do_derivatives=.FALSE., debug_invert=do_invert)
76 : END DO
77 : ! SP
78 0 : matrix_n%sp_d(j, :, :) = (matrix_p%sp - matrix_m%sp)/(2.0_dp*dx)
79 0 : DO i = 1, 3
80 0 : DO k = 1, 3
81 0 : IF (.NOT. check_value(matrix_n%sp_d(j, k, i), ij_matrix%sp_d(j, k, i), dx, 0.1_dp)) THEN
82 0 : WRITE (*, *) "ERROR for SP rotation matrix derivative SP(j,k,i), j,k,i::", j, k, i
83 0 : CPABORT("")
84 : END IF
85 : END DO
86 : END DO
87 : ! PP
88 0 : matrix_n%pp_d(j, :, :, :) = (matrix_p%pp - matrix_m%pp)/(2.0_dp*dx)
89 0 : DO i = 1, 3
90 0 : DO k = 1, 3
91 0 : DO l = 1, 6
92 0 : IF (.NOT. check_value(matrix_n%pp_d(j, l, k, i), ij_matrix%pp_d(j, l, k, i), dx, 0.1_dp)) THEN
93 0 : WRITE (*, *) "ERROR for PP rotation matrix derivative PP(j,l,k,i), j,l,k,i::", j, l, k, i
94 0 : CPABORT("")
95 : END IF
96 : END DO
97 : END DO
98 : END DO
99 : ! d-orbitals debug
100 0 : IF (sepi%dorb .OR. sepj%dorb) THEN
101 : ! SD
102 0 : matrix_n%sd_d(j, :, :) = (matrix_p%sd - matrix_m%sd)/(2.0_dp*dx)
103 0 : DO i = 1, 5
104 0 : DO k = 1, 5
105 0 : IF (.NOT. check_value(matrix_n%sd_d(j, k, i), ij_matrix%sd_d(j, k, i), dx, 0.1_dp)) THEN
106 0 : WRITE (*, *) "ERROR for SD rotation matrix derivative SD(j,k,i), j,k,i::", j, k, i
107 0 : CPABORT("")
108 : END IF
109 : END DO
110 : END DO
111 : ! DP
112 0 : matrix_n%pd_d(j, :, :, :) = (matrix_p%pd - matrix_m%pd)/(2.0_dp*dx)
113 0 : DO i = 1, 3
114 0 : DO k = 1, 5
115 0 : DO l = 1, 15
116 0 : IF (.NOT. check_value(matrix_n%pd_d(j, l, k, i), ij_matrix%pd_d(j, l, k, i), dx, 0.1_dp)) THEN
117 0 : WRITE (*, *) "ERROR for DP rotation matrix derivative DP(j,l,k,i), j,l,k,i::", j, l, k, i
118 0 : CPABORT("")
119 : END IF
120 : END DO
121 : END DO
122 : END DO
123 : ! DD
124 0 : matrix_n%dd_d(j, :, :, :) = (matrix_p%dd - matrix_m%dd)/(2.0_dp*dx)
125 0 : DO i = 1, 5
126 0 : DO k = 1, 5
127 0 : DO l = 1, 15
128 0 : IF (.NOT. check_value(matrix_n%dd_d(j, l, k, i), ij_matrix%dd_d(j, l, k, i), dx, 0.1_dp)) THEN
129 0 : WRITE (*, *) "ERROR for DD rotation matrix derivative DD(j,l,k,i), j,l,k,i::", j, l, k, i
130 0 : CPABORT("")
131 : END IF
132 : END DO
133 : END DO
134 : END DO
135 : END IF
136 : END DO
137 0 : CALL rotmat_release(matrix_p)
138 0 : CALL rotmat_release(matrix_m)
139 0 : CALL rotmat_release(matrix_n)
140 0 : END SUBROUTINE check_rotmat_der
141 :
142 : ! **************************************************************************************************
143 : !> \brief Check Numerical Vs Analytical
144 : !> \param sepi ...
145 : !> \param sepj ...
146 : !> \param rijv ...
147 : !> \param se_int_control ...
148 : !> \param se_taper ...
149 : !> \param invert ...
150 : !> \param ii ...
151 : !> \param kk ...
152 : !> \param v_d ...
153 : !> \par History
154 : !> 04.2008 created [tlaino]
155 : !> \author Teodoro Laino - Zurich University
156 : !> \note
157 : !> Debug routine
158 : ! **************************************************************************************************
159 0 : SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, v_d)
160 :
161 : USE kinds, ONLY: dp
162 : USE semi_empirical_int_utils, ONLY: rotmat
163 : USE semi_empirical_int_arrays, ONLY: indexb
164 : USE semi_empirical_int_num, ONLY: terep_num
165 : USE semi_empirical_types, ONLY: semi_empirical_type, &
166 : rotmat_type, &
167 : rotmat_create, &
168 : rotmat_release, &
169 : se_int_control_type, &
170 : se_taper_type
171 : USE semi_empirical_int_utils, ONLY: rot_2el_2c_first
172 : #include "./base/base_uses.f90"
173 : IMPLICIT NONE
174 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
175 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rijv
176 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
177 : LOGICAL, INTENT(IN) :: invert
178 : TYPE(se_taper_type), POINTER :: se_taper
179 : INTEGER, INTENT(IN) :: ii, kk
180 : REAL(KIND=dp), DIMENSION(3, 45, 45), &
181 : INTENT(IN) :: v_d
182 :
183 : CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
184 : routineN = 'rot_2el_2c_first', routineP = moduleN//':'//routineN
185 :
186 : REAL(KIND=dp), DIMENSION(491) :: rep
187 : LOGICAL, DIMENSION(45, 45) :: logv
188 : REAL(KIND=dp), DIMENSION(45, 45) :: v_n, v_p, v_m
189 :
190 : REAL(KIND=dp) :: dx, r, r0(3), x(3)
191 : TYPE(rotmat_type), POINTER :: matrix
192 : INTEGER :: imap(3), i, j, k, limkl
193 :
194 : INTERFACE
195 : FUNCTION check_value(num, ana, minval, thrs) RESULT(passed)
196 : USE kinds, ONLY: dp
197 : IMPLICIT NONE
198 : REAL(KIND=dp) :: num, ana, thrs, minval
199 : LOGICAL :: passed
200 : END FUNCTION check_value
201 : END INTERFACE
202 :
203 0 : NULLIFY (matrix)
204 0 : dx = 1.0E-6_dp
205 0 : imap(1) = 1
206 0 : imap(2) = 2
207 0 : imap(3) = 3
208 0 : IF (invert) THEN
209 0 : imap(1) = 3
210 : imap(2) = 2
211 0 : imap(3) = 1
212 : END IF
213 0 : limkl = indexb(kk, kk)
214 : ! Check derivatives: analytical VS numerical
215 0 : WRITE (*, *) "DEBUG::"//routineP
216 0 : DO j = 1, 3
217 0 : x = 0.0_dp
218 0 : x(imap(j)) = dx
219 0 : DO i = 1, 2
220 0 : r0 = rijv + (-1.0_dp)**(i - 1)*x
221 0 : r = SQRT(DOT_PRODUCT(r0, r0))
222 :
223 0 : CALL rotmat_create(matrix)
224 0 : CALL rotmat(sepi, sepj, r0, r, matrix, do_derivatives=.FALSE., debug_invert=invert)
225 :
226 : ! Compute integrals in diatomic frame
227 0 : CALL terep_num(sepi, sepj, r, rep, se_taper=se_taper, se_int_control=se_int_control)
228 :
229 0 : IF (i == 1) THEN
230 : CALL rot_2el_2c_first(sepi, sepj, r0, se_int_control, se_taper, invert, ii, kk, rep, logv, matrix, &
231 0 : v_p, lgrad=.FALSE.)
232 : END IF
233 0 : IF (i == 2) THEN
234 : CALL rot_2el_2c_first(sepi, sepj, r0, se_int_control, se_taper, invert, ii, kk, rep, logv, matrix, &
235 0 : v_m, lgrad=.FALSE.)
236 : END IF
237 0 : CALL rotmat_release(matrix)
238 : END DO
239 : ! Check numerical VS analytical
240 0 : DO i = 1, 45
241 0 : DO k = 1, limkl
242 : ! Compute the numerical derivative
243 0 : v_n(i, k) = (v_p(i, k) - v_m(i, k))/(2.0_dp*dx)
244 0 : IF (.NOT. check_value(v_d(j, i, k), v_n(i, k), dx, 0.1_dp)) THEN
245 0 : WRITE (*, *) "ERROR for rot_2el_2c_first derivative V_D(j,i,k), j,i,k::", j, i, k
246 0 : CPABORT("")
247 : END IF
248 : END DO
249 : END DO
250 : END DO
251 0 : END SUBROUTINE rot_2el_2c_first_debug
252 :
253 : ! **************************************************************************************************
254 : !> \brief Check Numerical Vs Analytical ssss
255 : !> \param sepi ...
256 : !> \param sepj ...
257 : !> \param r ...
258 : !> \param dssss ...
259 : !> \param itype ...
260 : !> \param se_int_control ...
261 : !> \param se_taper ...
262 : !> \par History
263 : !> 04.2008 created [tlaino]
264 : !> \author Teodoro Laino - Zurich University
265 : !> \note
266 : !> Debug routine
267 : ! **************************************************************************************************
268 0 : SUBROUTINE check_dssss_nucint_ana(sepi, sepj, r, dssss, itype, se_int_control, se_taper)
269 :
270 : USE kinds, ONLY: dp
271 : USE semi_empirical_int_num, ONLY: ssss_nucint_num
272 : USE semi_empirical_types, ONLY: semi_empirical_type, &
273 : se_int_control_type, &
274 : se_taper_type
275 : #include "./base/base_uses.f90"
276 : IMPLICIT NONE
277 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
278 : REAL(dp), INTENT(IN) :: r
279 : REAL(dp), INTENT(IN) :: dssss
280 : INTEGER, INTENT(IN) :: itype
281 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
282 : TYPE(se_taper_type), POINTER :: se_taper
283 :
284 : CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
285 : routineN = 'check_dssss_nucint_ana', routineP = moduleN//':'//routineN
286 :
287 : REAL(dp) :: delta, nssss, od, rn, ssssm, ssssp
288 :
289 : INTERFACE
290 : FUNCTION check_value(num, ana, minval, thrs) RESULT(passed)
291 : USE kinds, ONLY: dp
292 : IMPLICIT NONE
293 : REAL(KIND=dp) :: num, ana, thrs, minval
294 : LOGICAL :: passed
295 : END FUNCTION check_value
296 : END INTERFACE
297 :
298 0 : delta = 1.0E-8_dp
299 0 : od = 0.5_dp/delta
300 0 : rn = r + delta
301 0 : CALL ssss_nucint_num(sepi, sepj, rn, ssssp, itype, se_taper, se_int_control)
302 0 : rn = r - delta
303 0 : CALL ssss_nucint_num(sepi, sepj, rn, ssssm, itype, se_taper, se_int_control)
304 0 : nssss = od*(ssssp - ssssm)
305 : ! check
306 0 : WRITE (*, *) "DEBUG::"//routineP
307 0 : IF (.NOT. check_value(nssss, dssss, delta, 0.1_dp)) THEN
308 0 : WRITE (*, *) "ERROR for SSSS derivative SSSS"
309 0 : CPABORT("")
310 : END IF
311 0 : END SUBROUTINE check_dssss_nucint_ana
312 :
313 : ! **************************************************************************************************
314 : !> \brief Check Numerical Vs Analytical core
315 : !> \param sepi ...
316 : !> \param sepj ...
317 : !> \param r ...
318 : !> \param dcore ...
319 : !> \param itype ...
320 : !> \param se_int_control ...
321 : !> \param se_taper ...
322 : !> \par History
323 : !> 04.2008 created [tlaino]
324 : !> \author Teodoro Laino - Zurich University
325 : !> \note
326 : !> Debug routine
327 : ! **************************************************************************************************
328 0 : SUBROUTINE check_dcore_nucint_ana(sepi, sepj, r, dcore, itype, se_int_control, se_taper)
329 :
330 : USE kinds, ONLY: dp
331 : USE semi_empirical_int_num, ONLY: core_nucint_num
332 : USE semi_empirical_types, ONLY: semi_empirical_type, &
333 : se_int_control_type, &
334 : se_taper_type
335 : #include "./base/base_uses.f90"
336 : IMPLICIT NONE
337 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
338 : REAL(dp), INTENT(IN) :: r
339 : REAL(dp), DIMENSION(10, 2), INTENT(IN) :: dcore
340 : INTEGER, INTENT(IN) :: itype
341 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
342 : TYPE(se_taper_type), POINTER :: se_taper
343 :
344 : CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
345 : routineN = 'check_dcore_nucint_ana', routineP = moduleN//':'//routineN
346 :
347 : INTEGER :: i, j
348 : REAL(dp) :: delta, od, rn
349 : REAL(dp), DIMENSION(10, 2) :: corem, corep, ncore
350 :
351 : INTERFACE
352 : FUNCTION check_value(num, ana, minval, thrs) RESULT(passed)
353 : USE kinds, ONLY: dp
354 : IMPLICIT NONE
355 : REAL(KIND=dp) :: num, ana, thrs, minval
356 : LOGICAL :: passed
357 : END FUNCTION check_value
358 : END INTERFACE
359 :
360 0 : delta = 1.0E-8_dp
361 0 : od = 0.5_dp/delta
362 0 : rn = r + delta
363 0 : CALL core_nucint_num(sepi, sepj, rn, corep, itype, se_taper, se_int_control)
364 0 : rn = r - delta
365 0 : CALL core_nucint_num(sepi, sepj, rn, corem, itype, se_taper, se_int_control)
366 0 : ncore = od*(corep - corem)
367 : ! check
368 0 : WRITE (*, *) "DEBUG::"//routineP
369 0 : DO i = 1, 2
370 0 : DO j = 1, 10
371 0 : IF (.NOT. check_value(ncore(j, i), dcore(j, i), delta, 0.1_dp)) THEN
372 0 : WRITE (*, *) "ERROR for CORE derivative CORE(j,i), j,i::", j, i
373 0 : CPABORT("")
374 : END IF
375 : END DO
376 : END DO
377 0 : END SUBROUTINE check_dcore_nucint_ana
378 :
379 : ! **************************************************************************************************
380 : !> \brief Low level comparison function between numberical and analytical values
381 : !> \param num ...
382 : !> \param ana ...
383 : !> \param minval ...
384 : !> \param thrs ...
385 : !> \return ...
386 : !> \par History
387 : !> 04.2008 created [tlaino]
388 : !> \author Teodoro Laino - Zurich University
389 : !> \note
390 : !> Debug routine
391 : ! **************************************************************************************************
392 0 : FUNCTION check_value(num, ana, minval, thrs) RESULT(passed)
393 : USE kinds, ONLY: dp
394 : IMPLICIT NONE
395 : REAL(KIND=dp) :: num, ana, thrs, minval
396 : LOGICAL :: passed
397 :
398 0 : passed = .TRUE.
399 :
400 0 : IF ((ABS(num) < minval) .AND. (ABS(ana) > minval)) THEN
401 0 : WRITE (*, *) "WARNING ---> ", num, ana, thrs
402 0 : ELSE IF ((ABS(num) > minval) .AND. (ABS(ana) < minval)) THEN
403 0 : WRITE (*, *) "WARNING ---> ", num, ana, thrs
404 0 : ELSE IF ((ABS(num) < minval) .AND. (ABS(ana) < minval)) THEN
405 : ! skip..
406 : RETURN
407 : END IF
408 0 : IF (ABS((num - ana)/num*100._dp) > thrs) THEN
409 0 : WRITE (*, *) ABS(num - ana)/num*100._dp, thrs
410 : passed = .FALSE.
411 : END IF
412 : IF (.NOT. passed) THEN
413 0 : WRITE (*, *) "ANALYTICAL ::", ana
414 0 : WRITE (*, *) "NUMERICAL ::", num
415 : END IF
416 : END FUNCTION check_value
417 :
418 : ! **************************************************************************************************
419 : !> \brief Check Numerical Vs Analytical
420 : !> \param sepi ...
421 : !> \param sepj ...
422 : !> \param rijv ...
423 : !> \param itype ...
424 : !> \param se_int_control ...
425 : !> \param se_taper ...
426 : !> \param e1b ...
427 : !> \param e2a ...
428 : !> \param de1b ...
429 : !> \param de2a ...
430 : !> \par History
431 : !> 04.2008 created [tlaino]
432 : !> \author Teodoro Laino - Zurich University
433 : !> \note
434 : !> Debug routine
435 : ! **************************************************************************************************
436 0 : SUBROUTINE check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, e1b, e2a, de1b, de2a)
437 :
438 : USE kinds, ONLY: dp
439 : USE semi_empirical_int_num, ONLY: rotnuc_num, &
440 : drotnuc_num
441 : USE semi_empirical_types, ONLY: semi_empirical_type, &
442 : se_int_control_type, &
443 : se_taper_type
444 : #include "./base/base_uses.f90"
445 : IMPLICIT NONE
446 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
447 : REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
448 : INTEGER, INTENT(IN) :: itype
449 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
450 : TYPE(se_taper_type), POINTER :: se_taper
451 : REAL(dp), DIMENSION(45), INTENT(IN), &
452 : OPTIONAL :: e1b, e2a
453 : REAL(dp), DIMENSION(3, 45), &
454 : INTENT(IN), OPTIONAL :: de1b, de2a
455 :
456 : CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
457 : routineN = 'check_drotnuc_ana', routineP = moduleN//':'//routineN
458 :
459 : INTEGER :: i, j
460 : LOGICAL :: l_de1b, l_de2a, l_e1b, l_e2a, &
461 : lgrad
462 : REAL(dp) :: delta
463 : REAL(KIND=dp), DIMENSION(45) :: e1b2, e2a2
464 : REAL(KIND=dp), DIMENSION(3, 45) :: de1b2, de2a2
465 :
466 : INTERFACE
467 : FUNCTION check_value(num, ana, minval, thrs) RESULT(passed)
468 : USE kinds, ONLY: dp
469 : IMPLICIT NONE
470 : REAL(KIND=dp) :: num, ana, thrs, minval
471 : LOGICAL :: passed
472 : END FUNCTION check_value
473 : END INTERFACE
474 :
475 0 : l_e1b = PRESENT(e1b)
476 0 : l_e2a = PRESENT(e2a)
477 0 : l_de1b = PRESENT(de1b)
478 0 : l_de2a = PRESENT(de2a)
479 0 : lgrad = l_de1b .OR. l_de2a
480 0 : delta = 1.0E-5_dp
481 : ! Check value of integrals
482 0 : WRITE (*, *) "DEBUG::"//routineP
483 0 : CALL rotnuc_num(sepi, sepj, rijv, e1b2, e2a2, itype, se_int_control, se_taper=se_taper)
484 0 : IF (l_e1b) THEN
485 0 : DO j = 1, 45
486 0 : IF (.NOT. check_value(e1b2(j), e1b(j), delta, 0.1_dp)) THEN
487 0 : WRITE (*, *) "ERROR for E1B value E1B(j), j::", j
488 0 : CPABORT("")
489 : END IF
490 : END DO
491 : END IF
492 0 : IF (l_e2a) THEN
493 0 : DO J = 1, 45
494 0 : IF (.NOT. check_value(e2a2(j), e2a(j), delta, 0.1_dp)) THEN
495 0 : WRITE (*, *) "ERROR for E2A value E2A(j), j::", j
496 0 : CPABORT("")
497 : END IF
498 : END DO
499 : END IF
500 :
501 : ! Check derivatives
502 0 : IF (lgrad) THEN
503 : CALL drotnuc_num(sepi, sepj, rijv, de1b2, de2a2, itype, delta=delta, &
504 0 : se_int_control=se_int_control, se_taper=se_taper)
505 0 : IF (l_de1b) THEN
506 0 : DO i = 1, 3
507 0 : DO j = 1, 45
508 : ! Additional check on the value of the integral before checking derivatives
509 0 : IF (ABS(e1b2(j)) > delta) THEN
510 0 : IF (.NOT. check_value(de1b2(i, j), de1b(i, j), delta, 0.1_dp)) THEN
511 0 : WRITE (*, *) "ERROR for derivative of E1B. DE1B(i,j), i,j::", i, j
512 0 : CPABORT("")
513 : END IF
514 : END IF
515 : END DO
516 : END DO
517 : END IF
518 0 : IF (l_de2a) THEN
519 0 : DO i = 1, 3
520 0 : DO j = 1, 45
521 : ! Additional check on the value of the integral before checking derivatives
522 0 : IF (ABS(e2a2(j)) > delta) THEN
523 0 : IF (.NOT. check_value(de2a2(i, j), de2a(i, j), delta, 0.1_dp)) THEN
524 0 : WRITE (*, *) "ERROR for derivative of E2A. DE2A(i,j), i,j::", i, j
525 0 : CPABORT("")
526 : END IF
527 : END IF
528 : END DO
529 : END DO
530 : END IF
531 : END IF
532 0 : END SUBROUTINE check_drotnuc_ana
533 :
534 : ! **************************************************************************************************
535 : !> \brief Check Numerical Vs Analytical CORE-CORE
536 : !> \param sepi ...
537 : !> \param sepj ...
538 : !> \param rijv ...
539 : !> \param itype ...
540 : !> \param se_int_control ...
541 : !> \param se_taper ...
542 : !> \param enuc ...
543 : !> \param denuc ...
544 : !> \par History
545 : !> 04.2007 created [tlaino]
546 : !> \author Teodoro Laino - Zurich University
547 : !> \note
548 : !> Debug routine
549 : ! **************************************************************************************************
550 0 : SUBROUTINE check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, enuc, denuc)
551 :
552 : USE kinds, ONLY: dp
553 : USE semi_empirical_int_num, ONLY: corecore_num, &
554 : dcorecore_num
555 : USE semi_empirical_types, ONLY: semi_empirical_type, &
556 : se_int_control_type, &
557 : se_taper_type
558 : #include "./base/base_uses.f90"
559 : IMPLICIT NONE
560 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
561 : REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
562 : INTEGER, INTENT(IN) :: itype
563 : REAL(dp), INTENT(IN), OPTIONAL :: enuc
564 : REAL(dp), DIMENSION(3), INTENT(IN), &
565 : OPTIONAL :: denuc
566 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
567 : TYPE(se_taper_type), POINTER :: se_taper
568 :
569 : CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
570 : routineN = 'check_dcorecore_ana', routineP = moduleN//':'//routineN
571 :
572 : INTEGER :: j
573 : REAL(dp) :: enuc_num, delta
574 : REAL(dp), DIMENSION(3) :: denuc_num
575 :
576 : INTERFACE
577 : FUNCTION check_value(num, ana, minval, thrs) RESULT(passed)
578 : USE kinds, ONLY: dp
579 : IMPLICIT NONE
580 : REAL(KIND=dp) :: num, ana, thrs, minval
581 : LOGICAL :: passed
582 : END FUNCTION check_value
583 : END INTERFACE
584 :
585 0 : WRITE (*, *) "DEBUG::"//routineP
586 0 : delta = 1.0E-7_dp
587 : ! check
588 0 : IF (PRESENT(enuc)) THEN
589 0 : CALL corecore_num(sepi, sepj, rijv, enuc_num, itype, se_int_control, se_taper)
590 0 : IF (.NOT. check_value(enuc, enuc_num, delta, 0.001_dp)) THEN
591 0 : WRITE (*, *) "ERROR for CORE-CORE energy value (numerical different from analytical)!"
592 0 : CPABORT("")
593 : END IF
594 : END IF
595 0 : IF (PRESENT(denuc)) THEN
596 0 : CALL dcorecore_num(sepi, sepj, rijv, denuc_num, itype, delta, se_int_control, se_taper)
597 0 : DO j = 1, 3
598 0 : IF (.NOT. check_value(denuc(j), denuc_num(j), delta, 0.001_dp)) THEN
599 0 : WRITE (*, *) "ERROR for CORE-CORE energy derivative value (numerical different from analytical). DENUC(j), j::", j
600 0 : CPABORT("")
601 : END IF
602 : END DO
603 : END IF
604 0 : END SUBROUTINE check_dcorecore_ana
605 :
606 : ! **************************************************************************************************
607 : !> \brief Check Numerical Vs Analytical DTEREP_ANA
608 : !> \param sepi ...
609 : !> \param sepj ...
610 : !> \param r ...
611 : !> \param ri ...
612 : !> \param dri ...
613 : !> \param se_int_control ...
614 : !> \param se_taper ...
615 : !> \param lgrad ...
616 : !> \par History
617 : !> 04.2007 created [tlaino]
618 : !> \author Teodoro Laino - Zurich University
619 : !> \note
620 : !> Debug routine
621 : ! **************************************************************************************************
622 0 : SUBROUTINE check_dterep_ana(sepi, sepj, r, ri, dri, se_int_control, se_taper, lgrad)
623 :
624 : USE kinds, ONLY: dp
625 : USE semi_empirical_int_num, ONLY: terep_num
626 : USE semi_empirical_types, ONLY: semi_empirical_type, &
627 : se_int_control_type, &
628 : se_taper_type
629 : #include "./base/base_uses.f90"
630 : IMPLICIT NONE
631 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
632 : REAL(dp), INTENT(IN) :: r
633 : REAL(dp), DIMENSION(491), INTENT(IN) :: ri, dri
634 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
635 : LOGICAL, INTENT(IN) :: lgrad
636 : TYPE(se_taper_type), POINTER :: se_taper
637 :
638 : CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
639 : routineN = 'check_dterep_ana', routineP = moduleN//':'//routineN
640 :
641 : INTEGER :: j, i
642 : REAL(dp) :: delta, od, rn
643 : REAL(dp), DIMENSION(491) :: nri, ri0, rim, rip
644 :
645 : INTERFACE
646 : FUNCTION check_value(num, ana, minval, thrs) RESULT(passed)
647 : USE kinds, ONLY: dp
648 : IMPLICIT NONE
649 : REAL(KIND=dp) :: num, ana, thrs, minval
650 : LOGICAL :: passed
651 : END FUNCTION check_value
652 : END INTERFACE
653 :
654 0 : delta = 1.0E-8_dp
655 0 : od = 0.5_dp/delta
656 0 : rn = r
657 0 : CALL terep_num(sepi, sepj, rn, ri0, se_taper, se_int_control)
658 0 : IF (lgrad) THEN
659 0 : rn = r + delta
660 0 : CALL terep_num(sepi, sepj, rn, rip, se_taper, se_int_control)
661 0 : rn = r - delta
662 0 : CALL terep_num(sepi, sepj, rn, rim, se_taper, se_int_control)
663 0 : nri = od*(rip - rim)
664 : END IF
665 : ! check
666 0 : WRITE (*, *) "DEBUG::"//routineP
667 0 : DO j = 1, 491
668 0 : IF (ABS(ri(j) - ri0(j)) > EPSILON(0.0_dp)) THEN
669 0 : WRITE (*, *) "Error in value of the integral RI", j, ri(j), ri0(j)
670 0 : CPABORT("")
671 : END IF
672 0 : IF (lgrad) THEN
673 0 : IF (.NOT. check_value(dri(j), nri(j), delta*100.0_dp, 0.1_dp)) THEN
674 0 : WRITE (*, *) "ERROR for derivative of RI integral, RI(j), j::", j
675 0 : WRITE (*, *) "FULL SET OF INTEGRALS: INDX ANAL NUM DIFF"
676 0 : DO i = 1, 491
677 0 : WRITE (*, '(I5,3F15.9)') i, dri(i), nri(i), nri(i) - dri(i)
678 : END DO
679 0 : CPABORT("")
680 : END IF
681 : END IF
682 : END DO
683 0 : END SUBROUTINE check_dterep_ana
684 :
685 : ! **************************************************************************************************
686 : !> \brief Check Numerical Vs Analytical ROTINT_ANA
687 : !> \param sepi ...
688 : !> \param sepj ...
689 : !> \param rijv ...
690 : !> \param w ...
691 : !> \param dw ...
692 : !> \param se_int_control ...
693 : !> \param se_taper ...
694 : !> \par History
695 : !> 04.2008 created [tlaino]
696 : !> \author Teodoro Laino - Zurich University
697 : !> \note
698 : !> Debug routine
699 : ! **************************************************************************************************
700 0 : SUBROUTINE check_rotint_ana(sepi, sepj, rijv, w, dw, se_int_control, se_taper)
701 :
702 : USE kinds, ONLY: dp
703 : USE semi_empirical_int_num, ONLY: rotint_num, &
704 : drotint_num
705 : USE semi_empirical_types, ONLY: semi_empirical_type, &
706 : se_int_control_type, &
707 : se_taper_type
708 : #include "./base/base_uses.f90"
709 : IMPLICIT NONE
710 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
711 : REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
712 : REAL(dp), DIMENSION(2025), INTENT(IN), &
713 : OPTIONAL :: w
714 : REAL(dp), DIMENSION(3, 2025), &
715 : INTENT(IN), OPTIONAL :: dw
716 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
717 : TYPE(se_taper_type), POINTER :: se_taper
718 :
719 : CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
720 : routineN = 'rotint_ana', routineP = moduleN//':'//routineN
721 :
722 : REAL(dp), DIMENSION(2025) :: w2
723 : REAL(dp), DIMENSION(3, 2025) :: dw2
724 : INTEGER :: j, i
725 : REAL(KIND=dp) :: delta
726 :
727 : INTERFACE
728 : FUNCTION check_value(num, ana, minval, thrs) RESULT(passed)
729 : USE kinds, ONLY: dp
730 : IMPLICIT NONE
731 : REAL(KIND=dp) :: num, ana, thrs, minval
732 : LOGICAL :: passed
733 : END FUNCTION check_value
734 : END INTERFACE
735 :
736 0 : delta = 1.0E-6_dp
737 0 : WRITE (*, *) "DEBUG::"//routineP
738 0 : IF (PRESENT(w)) THEN
739 0 : w2 = 0.0_dp
740 0 : CALL rotint_num(sepi, sepj, rijv, w2, se_int_control, se_taper=se_taper)
741 0 : DO j = 1, 2025
742 0 : IF (.NOT. check_value(w(j), w2(j), delta, 0.1_dp)) THEN
743 0 : WRITE (*, *) "ERROR for integral value W(j), j::", j
744 0 : CPABORT("")
745 : END IF
746 : END DO
747 : END IF
748 0 : IF (PRESENT(dw)) THEN
749 : ! Numerical derivatives are obviosly a big problem..
750 : ! First of all let's decide if the value we get for delta is compatible
751 : ! with a reasonable value of the integral.. (compatible if the value of the
752 : ! integral is greater than 1.0E-6)
753 0 : CALL drotint_num(sepi, sepj, rijv, dw2, delta=delta, se_int_control=se_int_control, se_taper=se_taper)
754 0 : CALL rotint_num(sepi, sepj, rijv, w2, se_int_control=se_int_control, se_taper=se_taper)
755 0 : DO i = 1, 3
756 0 : DO j = 1, 2025
757 0 : IF ((ABS(w2(j)) > delta) .AND. (ABS(dw2(i, j)) > delta*10)) THEN
758 0 : IF (.NOT. check_value(dw(i, j), dw2(i, j), delta, 0.1_dp)) THEN
759 0 : WRITE (*, *) "ERROR for derivative of the integral value W(j). DW(i,j) i,j::", i, j
760 0 : CPABORT("")
761 : END IF
762 : END IF
763 : END DO
764 : END DO
765 : END IF
766 0 : END SUBROUTINE check_rotint_ana
|
