Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Debug energy and derivatives w.r.t. finite differences
10 : !> \note
11 : !> Use INTERPOLATION USE_GUESS, in order to perform force and energy
12 : !> calculations with the same density. This is not compulsory when iterating
13 : !> to selfconsistency, but essential in the non-selfconsistent case [08.2005,TdK].
14 : !> \par History
15 : !> 12.2004 created [tlaino]
16 : !> 08.2005 consistent_energies option added, to allow FD calculations
17 : !> with the correct energies in the non-selfconsistent case, but
18 : !> keep in mind, that the QS energies and forces are then NOT
19 : !> consistent to each other [TdK].
20 : !> 08.2005 In case the Harris functional is used, consistent_energies is
21 : !> et to .FALSE., otherwise the QS energies are spuriously used [TdK].
22 : !> 01.2015 Remove Harris functional option
23 : !> - Revised (20.11.2013,MK)
24 : !> \author Teodoro Laino
25 : ! **************************************************************************************************
26 : MODULE cp2k_debug
27 : USE cell_types, ONLY: cell_type,&
28 : get_cell
29 : USE cp_control_types, ONLY: dft_control_type
30 : USE cp_log_handling, ONLY: cp_get_default_logger,&
31 : cp_logger_type
32 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
33 : cp_print_key_unit_nr
34 : USE cp_result_methods, ONLY: get_results,&
35 : test_for_result
36 : USE cp_result_types, ONLY: cp_result_type
37 : USE cp_subsys_types, ONLY: cp_subsys_get,&
38 : cp_subsys_type
39 : USE force_env_methods, ONLY: force_env_calc_energy_force,&
40 : force_env_calc_num_pressure
41 : USE force_env_types, ONLY: force_env_get,&
42 : force_env_type,&
43 : use_qs_force
44 : USE input_constants, ONLY: do_stress_analytical,&
45 : do_stress_diagonal_anal,&
46 : do_stress_diagonal_numer,&
47 : do_stress_numerical
48 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
49 : section_vals_type,&
50 : section_vals_val_get
51 : USE kinds, ONLY: default_string_length,&
52 : dp
53 : USE message_passing, ONLY: mp_para_env_type
54 : USE particle_methods, ONLY: write_fist_particle_coordinates,&
55 : write_qs_particle_coordinates
56 : USE particle_types, ONLY: particle_type
57 : USE qs_environment_types, ONLY: get_qs_env
58 : USE qs_kind_types, ONLY: qs_kind_type
59 : USE string_utilities, ONLY: uppercase
60 : USE virial_types, ONLY: virial_set,&
61 : virial_type
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 : PRIVATE
66 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp2k_debug'
67 :
68 : PUBLIC :: cp2k_debug_energy_and_forces
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief ...
74 : !> \param force_env ...
75 : ! **************************************************************************************************
76 970 : SUBROUTINE cp2k_debug_energy_and_forces(force_env)
77 :
78 : TYPE(force_env_type), POINTER :: force_env
79 :
80 : CHARACTER(LEN=3) :: cval1
81 : CHARACTER(LEN=3*default_string_length) :: message
82 : CHARACTER(LEN=60) :: line
83 970 : CHARACTER(LEN=80), DIMENSION(:), POINTER :: cval2
84 : CHARACTER(LEN=default_string_length) :: description
85 : INTEGER :: i, ip, irep, iw, j, k, n_periodic, np, &
86 : nrep, stress_tensor
87 : INTEGER, DIMENSION(3) :: periodic
88 : LOGICAL :: check_failed, debug_dipole, &
89 : debug_forces, debug_polar, &
90 : debug_stress_tensor, skip, &
91 : stop_on_mismatch
92 970 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: do_dof_atom_coor
93 : LOGICAL, DIMENSION(3) :: do_dof_dipole
94 : LOGICAL, DIMENSION(3, 3) :: check_stress_element
95 : REAL(KIND=dp) :: amplitude, dd, de, derr, difference, dx, eps_no_error_check, errmax, &
96 : maxerr, periodic_stress_sum, std_value, sum_of_differences
97 : REAL(KIND=dp), DIMENSION(2) :: numer_energy
98 : REAL(KIND=dp), DIMENSION(3) :: dipole_moment, dipole_numer, err, &
99 : my_maxerr, poldir
100 : REAL(KIND=dp), DIMENSION(3, 2) :: dipn
101 : REAL(KIND=dp), DIMENSION(3, 3) :: polar_analytic, polar_numeric, polerr
102 : REAL(KIND=dp), DIMENSION(9) :: pvals
103 970 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: analyt_forces, numer_forces
104 : TYPE(cell_type), POINTER :: cell
105 : TYPE(cp_logger_type), POINTER :: logger
106 : TYPE(cp_result_type), POINTER :: results
107 : TYPE(cp_subsys_type), POINTER :: subsys
108 : TYPE(dft_control_type), POINTER :: dft_control
109 : TYPE(mp_para_env_type), POINTER :: para_env
110 970 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
111 970 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
112 : TYPE(section_vals_type), POINTER :: root_section, subsys_section
113 :
114 970 : NULLIFY (analyt_forces, numer_forces, subsys, particles)
115 :
116 970 : root_section => force_env%root_section
117 :
118 970 : CALL force_env_get(force_env, para_env=para_env, subsys=subsys, cell=cell)
119 : subsys_section => section_vals_get_subs_vals(force_env%force_env_section, &
120 970 : "SUBSYS")
121 :
122 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_STRESS_TENSOR", &
123 970 : l_val=debug_stress_tensor)
124 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_FORCES", &
125 970 : l_val=debug_forces)
126 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_DIPOLE", &
127 970 : l_val=debug_dipole)
128 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_POLARIZABILITY", &
129 970 : l_val=debug_polar)
130 : CALL section_vals_val_get(root_section, "DEBUG%DX", &
131 970 : r_val=dx)
132 : CALL section_vals_val_get(root_section, "DEBUG%DE", &
133 970 : r_val=de)
134 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_DIPOLE_DIRS", &
135 970 : c_val=cval1)
136 970 : dx = ABS(dx)
137 : CALL section_vals_val_get(root_section, "DEBUG%MAX_RELATIVE_ERROR", &
138 970 : r_val=maxerr)
139 : CALL section_vals_val_get(root_section, "DEBUG%EPS_NO_ERROR_CHECK", &
140 970 : r_val=eps_no_error_check)
141 970 : eps_no_error_check = MAX(eps_no_error_check, EPSILON(0.0_dp))
142 : CALL section_vals_val_get(root_section, "DEBUG%STOP_ON_MISMATCH", &
143 970 : l_val=stop_on_mismatch)
144 :
145 : ! set active DOF
146 970 : CALL uppercase(cval1)
147 970 : do_dof_dipole(1) = (INDEX(cval1, "X") /= 0)
148 970 : do_dof_dipole(2) = (INDEX(cval1, "Y") /= 0)
149 970 : do_dof_dipole(3) = (INDEX(cval1, "Z") /= 0)
150 970 : NULLIFY (cval2)
151 970 : IF (debug_forces) THEN
152 652 : np = subsys%particles%n_els
153 1956 : ALLOCATE (do_dof_atom_coor(3, np))
154 652 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", n_rep_val=nrep)
155 652 : IF (nrep > 0) THEN
156 268 : do_dof_atom_coor = .FALSE.
157 546 : DO irep = 1, nrep
158 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", i_rep_val=irep, &
159 278 : c_vals=cval2)
160 278 : READ (cval2(1), FMT="(I10)") k
161 278 : CALL uppercase(cval2(2))
162 278 : do_dof_atom_coor(1, k) = (INDEX(cval2(2), "X") /= 0)
163 278 : do_dof_atom_coor(2, k) = (INDEX(cval2(2), "Y") /= 0)
164 546 : do_dof_atom_coor(3, k) = (INDEX(cval2(2), "Z") /= 0)
165 : END DO
166 : ELSE
167 6680 : do_dof_atom_coor = .TRUE.
168 : END IF
169 : END IF
170 :
171 970 : logger => cp_get_default_logger()
172 : iw = cp_print_key_unit_nr(logger, root_section, "DEBUG%PROGRAM_RUN_INFO", &
173 970 : extension=".log")
174 970 : IF (debug_stress_tensor) THEN
175 1376 : IF (SUM(cell%perd) == 0) THEN
176 : CALL cp_warn(__LOCATION__, &
177 : "DEBUG_STRESS_TENSOR requested for PERIODIC NONE. "// &
178 : "The isolated-system virial is useful for finite-difference diagnostics, "// &
179 58 : "but it is not a physically meaningful bulk stress.")
180 : END IF
181 : ! To debug stress tensor the stress tensor calculation must be
182 : ! first enabled..
183 : CALL section_vals_val_get(force_env%force_env_section, "STRESS_TENSOR", &
184 344 : i_val=stress_tensor)
185 344 : skip = .FALSE.
186 0 : SELECT CASE (stress_tensor)
187 : CASE (do_stress_analytical, do_stress_diagonal_anal)
188 : ! OK
189 : CASE (do_stress_numerical, do_stress_diagonal_numer)
190 : ! Nothing to check
191 : CALL cp_warn(__LOCATION__, "Numerical stress tensor was requested in "// &
192 0 : "the FORCE_EVAL section. Nothing to debug!")
193 146 : skip = .TRUE.
194 : CASE DEFAULT
195 : CALL cp_warn(__LOCATION__, "Stress tensor calculation was not enabled in "// &
196 146 : "the FORCE_EVAL section. Nothing to debug!")
197 344 : skip = .TRUE.
198 : END SELECT
199 :
200 : IF (.NOT. skip) THEN
201 :
202 : BLOCK
203 : TYPE(virial_type) :: virial_analytical, virial_numerical
204 : TYPE(virial_type), POINTER :: virial
205 :
206 : ! Compute the analytical stress tensor
207 198 : CALL cp_subsys_get(subsys, virial=virial)
208 198 : CALL virial_set(virial, pv_numer=.FALSE.)
209 : CALL force_env_calc_energy_force(force_env, &
210 : calc_force=.TRUE., &
211 198 : calc_stress_tensor=.TRUE.)
212 :
213 : ! Retrieve the analytical virial
214 198 : virial_analytical = virial
215 :
216 : ! Debug stress tensor (numerical vs analytical)
217 198 : CALL virial_set(virial, pv_numer=.TRUE.)
218 198 : CALL force_env_calc_num_pressure(force_env, dx=dx)
219 :
220 : ! Retrieve the numerical virial
221 198 : CALL cp_subsys_get(subsys, virial=virial)
222 198 : virial_numerical = virial
223 :
224 : ! Numerical diagonal stress checks only perturb diagonal cell elements.
225 198 : IF (virial_analytical%pv_diagonal .OR. virial_numerical%pv_diagonal) THEN
226 168 : DO i = 1, 3
227 702 : DO k = 1, 3
228 504 : IF (i /= k) THEN
229 252 : virial_analytical%pv_virial(i, k) = 0.0_dp
230 252 : virial_numerical%pv_virial(i, k) = 0.0_dp
231 : END IF
232 : END DO
233 : END DO
234 : END IF
235 :
236 198 : CALL get_cell(cell=cell, periodic=periodic)
237 792 : n_periodic = COUNT(periodic /= 0)
238 2574 : check_stress_element = .TRUE.
239 198 : IF (n_periodic > 0 .AND. n_periodic < 3) THEN
240 50 : check_stress_element = .FALSE.
241 200 : DO i = 1, 3
242 650 : DO k = 1, 3
243 904 : check_stress_element(i, k) = periodic(i) /= 0 .AND. periodic(k) /= 0
244 : END DO
245 : END DO
246 : END IF
247 :
248 : ! Print results
249 198 : IF (iw > 0) THEN
250 : WRITE (UNIT=iw, FMT="((T2,A))") &
251 99 : "DEBUG| Numerical pv_virial [a.u.]"
252 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
253 396 : ("DEBUG|", virial_numerical%pv_virial(i, 1:3), i=1, 3)
254 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
255 99 : "DEBUG| Analytical pv_virial [a.u.]"
256 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
257 396 : ("DEBUG|", virial_analytical%pv_virial(i, 1:3), i=1, 3)
258 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
259 99 : "DEBUG| Difference pv_virial [a.u.]"
260 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
261 1287 : ("DEBUG|", virial_numerical%pv_virial(i, 1:3) - virial_analytical%pv_virial(i, 1:3), i=1, 3)
262 : WRITE (UNIT=iw, FMT="(T2,A,T61,F20.12)") &
263 99 : "DEBUG| Sum of differences", &
264 1386 : SUM(ABS(virial_numerical%pv_virial(:, :) - virial_analytical%pv_virial(:, :)))
265 99 : IF (n_periodic > 0 .AND. n_periodic < 3) THEN
266 25 : periodic_stress_sum = 0.0_dp
267 100 : DO i = 1, 3
268 325 : DO k = 1, 3
269 300 : IF (periodic(i) /= 0 .AND. periodic(k) /= 0) THEN
270 : periodic_stress_sum = periodic_stress_sum + &
271 : ABS(virial_numerical%pv_virial(i, k) - &
272 73 : virial_analytical%pv_virial(i, k))
273 : END IF
274 : END DO
275 : END DO
276 : WRITE (UNIT=iw, FMT="(T2,A,T61,F20.12)") &
277 25 : "DEBUG| Periodic-subspace sum of differences", periodic_stress_sum
278 : END IF
279 : END IF
280 :
281 : ! Check and abort on failure
282 198 : check_failed = .FALSE.
283 198 : IF (iw > 0) THEN
284 : WRITE (UNIT=iw, FMT="(/T2,A)") &
285 99 : "DEBUG| Relative error pv_virial"
286 : WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.8)") &
287 99 : "DEBUG| Threshold value for error check [a.u.]", eps_no_error_check
288 : END IF
289 792 : DO i = 1, 3
290 594 : err(:) = 0.0_dp
291 2376 : DO k = 1, 3
292 1782 : IF (check_stress_element(i, k) .AND. &
293 594 : ABS(virial_analytical%pv_virial(i, k)) >= eps_no_error_check) THEN
294 : err(k) = 100.0_dp*(virial_numerical%pv_virial(i, k) - virial_analytical%pv_virial(i, k))/ &
295 1038 : virial_analytical%pv_virial(i, k)
296 1038 : WRITE (UNIT=line(20*(k - 1) + 1:20*k), FMT="(1X,F17.2,A2)") err(k), " %"
297 : ELSE
298 744 : WRITE (UNIT=line(20*(k - 1) + 1:20*k), FMT="(17X,A3)") "- %"
299 : END IF
300 : END DO
301 594 : IF (iw > 0) THEN
302 : WRITE (UNIT=iw, FMT="(T2,A,T21,A60)") &
303 297 : "DEBUG|", line
304 : END IF
305 2574 : IF (ANY(ABS(err(1:3)) > maxerr)) check_failed = .TRUE.
306 : END DO
307 198 : IF (iw > 0) THEN
308 : WRITE (UNIT=iw, FMT="(T2,A,T61,F18.2,A2)") &
309 99 : "DEBUG| Maximum accepted error", maxerr, " %"
310 : END IF
311 90684 : IF (check_failed) THEN
312 : message = "A mismatch between the analytical and the numerical "// &
313 : "stress tensor has been detected. Check the implementation "// &
314 0 : "of the stress tensor"
315 0 : IF (stop_on_mismatch) THEN
316 0 : CPABORT(TRIM(message))
317 : ELSE
318 0 : CPWARN(TRIM(message))
319 : END IF
320 : END IF
321 : END BLOCK
322 : END IF
323 : END IF
324 :
325 970 : IF (debug_forces) THEN
326 : ! Debug forces (numerical vs analytical)
327 652 : particles => subsys%particles%els
328 1204 : SELECT CASE (force_env%in_use)
329 : CASE (use_qs_force)
330 552 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
331 552 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
332 : CASE DEFAULT
333 652 : CALL write_fist_particle_coordinates(particles, subsys_section)
334 : END SELECT
335 : ! First evaluate energy and forces
336 : CALL force_env_calc_energy_force(force_env, &
337 : calc_force=.TRUE., &
338 652 : calc_stress_tensor=.FALSE.)
339 : ! Copy forces in array and start the numerical calculation
340 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
341 652 : np = subsys%particles%n_els
342 1956 : ALLOCATE (analyt_forces(np, 3))
343 3060 : DO ip = 1, np
344 10284 : analyt_forces(ip, 1:3) = particles(ip)%f(1:3)
345 : END DO
346 : ! Loop on atoms and coordinates
347 : IF (ASSOCIATED(numer_forces)) DEALLOCATE (numer_forces)
348 1304 : ALLOCATE (numer_forces(subsys%particles%n_els, 3))
349 3060 : Atom: DO ip = 1, np
350 9632 : Coord: DO k = 1, 3
351 9632 : IF (do_dof_atom_coor(k, ip)) THEN
352 5042 : numer_energy = 0.0_dp
353 5042 : std_value = particles(ip)%r(k)
354 15126 : DO j = 1, 2
355 10084 : particles(ip)%r(k) = std_value - (-1.0_dp)**j*dx
356 16892 : SELECT CASE (force_env%in_use)
357 : CASE (use_qs_force)
358 6808 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
359 6808 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
360 : CASE DEFAULT
361 10084 : CALL write_fist_particle_coordinates(particles, subsys_section)
362 : END SELECT
363 : ! Compute energy
364 : CALL force_env_calc_energy_force(force_env, &
365 : calc_force=.FALSE., &
366 : calc_stress_tensor=.FALSE., &
367 10084 : consistent_energies=.TRUE.)
368 15126 : CALL force_env_get(force_env, potential_energy=numer_energy(j))
369 : END DO
370 5042 : particles(ip)%r(k) = std_value
371 5042 : numer_forces(ip, k) = -0.5_dp*(numer_energy(1) - numer_energy(2))/dx
372 5042 : IF (iw > 0) THEN
373 : WRITE (UNIT=iw, FMT="(/,T2,A,T17,A,F7.4,A,T34,A,F7.4,A,T52,A,T68,A)") &
374 2521 : "DEBUG| Atom", "E("//ACHAR(119 + k)//" +", dx, ")", &
375 2521 : "E("//ACHAR(119 + k)//" -", dx, ")", &
376 5042 : "f(numerical)", "f(analytical)"
377 : WRITE (UNIT=iw, FMT="(T2,A,I5,4(1X,F16.8))") &
378 2521 : "DEBUG|", ip, numer_energy(1:2), numer_forces(ip, k), analyt_forces(ip, k)
379 : END IF
380 : ELSE
381 2182 : numer_forces(ip, k) = 0.0_dp
382 : END IF
383 : END DO Coord
384 : ! Check analytical forces using the numerical forces as reference
385 9632 : my_maxerr = maxerr
386 2408 : err(1:3) = 0.0_dp
387 9632 : DO k = 1, 3
388 9632 : IF (do_dof_atom_coor(k, ip)) THEN
389 : ! Calculate percentage but ignore very small force values
390 5042 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
391 3786 : err(k) = 100.0_dp*(numer_forces(ip, k) - analyt_forces(ip, k))/analyt_forces(ip, k)
392 : END IF
393 : ! Increase error tolerance for small force values
394 5042 : IF (ABS(analyt_forces(ip, k)) <= 0.0001_dp) my_maxerr(k) = 5.0_dp*my_maxerr(k)
395 : ELSE
396 2182 : err(k) = 0.0_dp
397 : END IF
398 : END DO
399 2408 : IF (iw > 0) THEN
400 2260 : IF (ANY(do_dof_atom_coor(1:3, ip))) THEN
401 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
402 926 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
403 : END IF
404 4816 : DO k = 1, 3
405 4816 : IF (do_dof_atom_coor(k, ip)) THEN
406 2521 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
407 2521 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
408 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T71,F10.2)") &
409 1893 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, err(k)
410 : ELSE
411 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T80,A1)") &
412 628 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, "-"
413 : END IF
414 : END IF
415 : END DO
416 : END IF
417 10282 : IF (ANY(ABS(err(1:3)) > my_maxerr(1:3))) THEN
418 : message = "A mismatch between analytical and numerical forces "// &
419 : "has been detected. Check the implementation of the "// &
420 2 : "analytical force calculation"
421 2 : IF (stop_on_mismatch) THEN
422 0 : CPABORT(message)
423 : ELSE
424 2 : CPWARN(message)
425 : END IF
426 : END IF
427 : END DO Atom
428 : ! Print summary
429 652 : IF (iw > 0) THEN
430 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
431 326 : "DEBUG|======================== BEGIN OF SUMMARY ===============================", &
432 652 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
433 326 : sum_of_differences = 0.0_dp
434 326 : errmax = 0.0_dp
435 1530 : DO ip = 1, np
436 1204 : err(1:3) = 0.0_dp
437 5142 : DO k = 1, 3
438 4816 : IF (do_dof_atom_coor(k, ip)) THEN
439 2521 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
440 2521 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
441 1893 : err(k) = 100_dp*(numer_forces(ip, k) - analyt_forces(ip, k))/analyt_forces(ip, k)
442 1893 : errmax = MAX(errmax, ABS(err(k)))
443 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T71,F10.2)") &
444 1893 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, err(k)
445 : ELSE
446 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T80,A1)") &
447 628 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, "-"
448 : END IF
449 2521 : sum_of_differences = sum_of_differences + ABS(difference)
450 : END IF
451 : END DO
452 : END DO
453 : WRITE (UNIT=iw, FMT="(T2,A,T57,F12.8,T71,F10.2)") &
454 326 : "DEBUG| Sum of differences:", sum_of_differences, errmax
455 : WRITE (UNIT=iw, FMT="(T2,A)") &
456 326 : "DEBUG|======================== END OF SUMMARY ================================="
457 : END IF
458 : ! Release work storage
459 652 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
460 652 : IF (ASSOCIATED(numer_forces)) DEALLOCATE (numer_forces)
461 652 : DEALLOCATE (do_dof_atom_coor)
462 : END IF
463 :
464 970 : IF (debug_dipole) THEN
465 116 : IF (force_env%in_use == use_qs_force) THEN
466 116 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
467 116 : poldir = [0.0_dp, 0.0_dp, 1.0_dp]
468 116 : amplitude = 0.0_dp
469 116 : CALL set_efield(dft_control, amplitude, poldir)
470 116 : CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
471 116 : CALL get_qs_env(force_env%qs_env, results=results)
472 116 : description = "[DIPOLE]"
473 116 : IF (test_for_result(results, description=description)) THEN
474 116 : CALL get_results(results, description=description, values=dipole_moment)
475 : ELSE
476 0 : CALL cp_warn(__LOCATION__, "Debug of dipole moments needs DFT/PRINT/MOMENTS section enabled")
477 0 : CPABORT("DEBUG")
478 : END IF
479 116 : amplitude = de
480 464 : DO k = 1, 3
481 464 : IF (do_dof_dipole(k)) THEN
482 224 : poldir = 0.0_dp
483 224 : poldir(k) = 1.0_dp
484 672 : DO j = 1, 2
485 1792 : poldir = -1.0_dp*poldir
486 448 : CALL set_efield(dft_control, amplitude, poldir)
487 448 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
488 672 : CALL force_env_get(force_env, potential_energy=numer_energy(j))
489 : END DO
490 224 : dipole_numer(k) = 0.5_dp*(numer_energy(1) - numer_energy(2))/de
491 : ELSE
492 124 : dipole_numer(k) = 0.0_dp
493 : END IF
494 : END DO
495 116 : IF (iw > 0) THEN
496 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
497 58 : "DEBUG|========================= DIPOLE MOMENTS ================================", &
498 116 : "DEBUG| Coordinate D(numerical) D(analytical) Difference Error [%]"
499 58 : err(1:3) = 0.0_dp
500 232 : DO k = 1, 3
501 232 : IF (do_dof_dipole(k)) THEN
502 112 : dd = dipole_moment(k) - dipole_numer(k)
503 112 : IF (ABS(dipole_moment(k)) > eps_no_error_check) THEN
504 59 : derr = 100._dp*dd/dipole_moment(k)
505 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
506 59 : "DEBUG|", ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd, derr
507 : ELSE
508 53 : derr = 0.0_dp
509 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
510 53 : "DEBUG|", ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd
511 : END IF
512 112 : err(k) = derr
513 : ELSE
514 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,A16,T38,F16.8)") &
515 62 : "DEBUG|", ACHAR(119 + k), " skipped", dipole_moment(k)
516 : END IF
517 : END DO
518 : WRITE (UNIT=iw, FMT="((T2,A))") &
519 58 : "DEBUG|========================================================================="
520 232 : WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") 'DIPOLE : CheckSum =', SUM(dipole_moment)
521 229 : IF (ANY(ABS(err(1:3)) > maxerr)) THEN
522 : message = "A mismatch between analytical and numerical dipoles "// &
523 : "has been detected. Check the implementation of the "// &
524 1 : "analytical dipole calculation"
525 1 : IF (stop_on_mismatch) THEN
526 0 : CPABORT(message)
527 : ELSE
528 1 : CPWARN(message)
529 : END IF
530 : END IF
531 : END IF
532 :
533 : ELSE
534 0 : CALL cp_warn(__LOCATION__, "Debug of dipole moments only for Quickstep code available")
535 : END IF
536 : END IF
537 :
538 970 : IF (debug_polar) THEN
539 82 : IF (force_env%in_use == use_qs_force) THEN
540 82 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
541 82 : poldir = [0.0_dp, 0.0_dp, 1.0_dp]
542 82 : amplitude = 0.0_dp
543 82 : CALL set_efield(dft_control, amplitude, poldir)
544 82 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
545 82 : CALL get_qs_env(force_env%qs_env, results=results)
546 82 : description = "[POLAR]"
547 82 : IF (test_for_result(results, description=description)) THEN
548 82 : CALL get_results(results, description=description, values=pvals)
549 82 : polar_analytic(1:3, 1:3) = RESHAPE(pvals(1:9), [3, 3])
550 : ELSE
551 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities needs PROPERTIES/LINRES/POLAR section enabled")
552 0 : CPABORT("DEBUG")
553 : END IF
554 82 : description = "[DIPOLE]"
555 82 : IF (.NOT. test_for_result(results, description=description)) THEN
556 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities need DFT/PRINT/MOMENTS section enabled")
557 0 : CPABORT("DEBUG")
558 : END IF
559 82 : amplitude = de
560 328 : DO k = 1, 3
561 246 : poldir = 0.0_dp
562 246 : poldir(k) = 1.0_dp
563 738 : DO j = 1, 2
564 1968 : poldir = -1.0_dp*poldir
565 492 : CALL set_efield(dft_control, amplitude, poldir)
566 492 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., linres=.TRUE.)
567 738 : CALL get_results(results, description=description, values=dipn(1:3, j))
568 : END DO
569 1066 : polar_numeric(1:3, k) = 0.5_dp*(dipn(1:3, 2) - dipn(1:3, 1))/de
570 : END DO
571 82 : IF (iw > 0) THEN
572 41 : polerr = 0.0_dp
573 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
574 41 : "DEBUG|========================= POLARIZABILITY ================================", &
575 82 : "DEBUG| Coordinates P(numerical) P(analytical) Difference Error [%]"
576 164 : DO k = 1, 3
577 533 : DO j = 1, 3
578 369 : dd = polar_analytic(k, j) - polar_numeric(k, j)
579 492 : IF (ABS(polar_analytic(k, j)) > eps_no_error_check) THEN
580 176 : derr = 100._dp*dd/polar_analytic(k, j)
581 : WRITE (UNIT=iw, FMT="(T2,A,T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
582 176 : "DEBUG|", ACHAR(119 + k), ACHAR(119 + j), polar_numeric(k, j), polar_analytic(k, j), dd, derr
583 176 : polerr(k, j) = derr
584 : ELSE
585 : WRITE (UNIT=iw, FMT="(T2,A,T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
586 193 : "DEBUG|", ACHAR(119 + k), ACHAR(119 + j), polar_numeric(k, j), polar_analytic(k, j), dd
587 : END IF
588 : END DO
589 : END DO
590 : WRITE (UNIT=iw, FMT="((T2,A))") &
591 41 : "DEBUG|========================================================================="
592 533 : WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") ' POLAR : CheckSum =', SUM(polar_analytic)
593 533 : IF (ANY(ABS(polerr(1:3, 1:3)) > maxerr)) THEN
594 : message = "A mismatch between analytical and numerical polarizabilities "// &
595 : "has been detected. Check the implementation of the "// &
596 0 : "analytical polarizabilitie calculation"
597 0 : IF (stop_on_mismatch) THEN
598 0 : CPABORT(message)
599 : ELSE
600 0 : CPWARN(message)
601 : END IF
602 : END IF
603 : END IF
604 : ELSE
605 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities only for Quickstep code available")
606 : END IF
607 : END IF
608 :
609 970 : CALL cp_print_key_finished_output(iw, logger, root_section, "DEBUG%PROGRAM_RUN_INFO")
610 :
611 1940 : END SUBROUTINE cp2k_debug_energy_and_forces
612 :
613 : ! **************************************************************************************************
614 : !> \brief ...
615 : !> \param dft_control ...
616 : !> \param amplitude ...
617 : !> \param poldir ...
618 : ! **************************************************************************************************
619 1138 : SUBROUTINE set_efield(dft_control, amplitude, poldir)
620 : TYPE(dft_control_type), POINTER :: dft_control
621 : REAL(KIND=dp), INTENT(IN) :: amplitude
622 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: poldir
623 :
624 1138 : IF (dft_control%apply_efield) THEN
625 1068 : dft_control%efield_fields(1)%efield%strength = amplitude
626 4272 : dft_control%efield_fields(1)%efield%polarisation(1:3) = poldir(1:3)
627 70 : ELSE IF (dft_control%apply_period_efield) THEN
628 70 : dft_control%period_efield%strength = amplitude
629 280 : dft_control%period_efield%polarisation(1:3) = poldir(1:3)
630 : ELSE
631 0 : CPABORT("No EFIELD definition available")
632 : END IF
633 :
634 1138 : END SUBROUTINE set_efield
635 :
636 : END MODULE cp2k_debug
|
