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 854 : 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 854 : 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 854 : 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
102 : REAL(KIND=dp), DIMENSION(9) :: pvals
103 854 : 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 854 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
111 854 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
112 : TYPE(section_vals_type), POINTER :: root_section, subsys_section
113 :
114 854 : NULLIFY (analyt_forces, numer_forces, subsys, particles)
115 :
116 854 : root_section => force_env%root_section
117 :
118 854 : 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 854 : "SUBSYS")
121 :
122 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_STRESS_TENSOR", &
123 854 : l_val=debug_stress_tensor)
124 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_FORCES", &
125 854 : l_val=debug_forces)
126 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_DIPOLE", &
127 854 : l_val=debug_dipole)
128 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_POLARIZABILITY", &
129 854 : l_val=debug_polar)
130 : CALL section_vals_val_get(root_section, "DEBUG%DX", &
131 854 : r_val=dx)
132 : CALL section_vals_val_get(root_section, "DEBUG%DE", &
133 854 : r_val=de)
134 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_DIPOLE_DIRS", &
135 854 : c_val=cval1)
136 854 : dx = ABS(dx)
137 : CALL section_vals_val_get(root_section, "DEBUG%MAX_RELATIVE_ERROR", &
138 854 : r_val=maxerr)
139 : CALL section_vals_val_get(root_section, "DEBUG%EPS_NO_ERROR_CHECK", &
140 854 : r_val=eps_no_error_check)
141 854 : 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 854 : l_val=stop_on_mismatch)
144 :
145 : ! set active DOF
146 854 : CALL uppercase(cval1)
147 854 : do_dof_dipole(1) = (INDEX(cval1, "X") /= 0)
148 854 : do_dof_dipole(2) = (INDEX(cval1, "Y") /= 0)
149 854 : do_dof_dipole(3) = (INDEX(cval1, "Z") /= 0)
150 854 : NULLIFY (cval2)
151 854 : IF (debug_forces) THEN
152 566 : np = subsys%particles%n_els
153 1698 : ALLOCATE (do_dof_atom_coor(3, np))
154 566 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", n_rep_val=nrep)
155 566 : IF (nrep > 0) THEN
156 210 : do_dof_atom_coor = .FALSE.
157 426 : DO irep = 1, nrep
158 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", i_rep_val=irep, &
159 216 : c_vals=cval2)
160 216 : READ (cval2(1), FMT="(I10)") k
161 216 : CALL uppercase(cval2(2))
162 216 : do_dof_atom_coor(1, k) = (INDEX(cval2(2), "X") /= 0)
163 216 : do_dof_atom_coor(2, k) = (INDEX(cval2(2), "Y") /= 0)
164 426 : do_dof_atom_coor(3, k) = (INDEX(cval2(2), "Z") /= 0)
165 : END DO
166 : ELSE
167 6252 : do_dof_atom_coor = .TRUE.
168 : END IF
169 : END IF
170 :
171 854 : logger => cp_get_default_logger()
172 : iw = cp_print_key_unit_nr(logger, root_section, "DEBUG%PROGRAM_RUN_INFO", &
173 854 : extension=".log")
174 854 : IF (debug_stress_tensor) THEN
175 1208 : 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 56 : "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 302 : i_val=stress_tensor)
185 302 : 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 120 : skip = .TRUE.
194 : CASE DEFAULT
195 : CALL cp_warn(__LOCATION__, "Stress tensor calculation was not enabled in "// &
196 120 : "the FORCE_EVAL section. Nothing to debug!")
197 302 : 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 182 : CALL cp_subsys_get(subsys, virial=virial)
208 182 : CALL virial_set(virial, pv_numer=.FALSE.)
209 : CALL force_env_calc_energy_force(force_env, &
210 : calc_force=.TRUE., &
211 182 : calc_stress_tensor=.TRUE.)
212 :
213 : ! Retrieve the analytical virial
214 182 : virial_analytical = virial
215 :
216 : ! Debug stress tensor (numerical vs analytical)
217 182 : CALL virial_set(virial, pv_numer=.TRUE.)
218 182 : CALL force_env_calc_num_pressure(force_env, dx=dx)
219 :
220 : ! Retrieve the numerical virial
221 182 : CALL cp_subsys_get(subsys, virial=virial)
222 182 : virial_numerical = virial
223 :
224 : ! Numerical diagonal stress checks only perturb diagonal cell elements.
225 182 : IF (virial_analytical%pv_diagonal .OR. virial_numerical%pv_diagonal) THEN
226 152 : DO i = 1, 3
227 638 : DO k = 1, 3
228 456 : IF (i /= k) THEN
229 228 : virial_analytical%pv_virial(i, k) = 0.0_dp
230 228 : virial_numerical%pv_virial(i, k) = 0.0_dp
231 : END IF
232 : END DO
233 : END DO
234 : END IF
235 :
236 182 : CALL get_cell(cell=cell, periodic=periodic)
237 728 : n_periodic = COUNT(periodic /= 0)
238 2366 : check_stress_element = .TRUE.
239 182 : IF (n_periodic > 0 .AND. n_periodic < 3) THEN
240 48 : check_stress_element = .FALSE.
241 192 : DO i = 1, 3
242 624 : DO k = 1, 3
243 870 : 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 182 : IF (iw > 0) THEN
250 : WRITE (UNIT=iw, FMT="((T2,A))") &
251 91 : "DEBUG| Numerical pv_virial [a.u.]"
252 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
253 364 : ("DEBUG|", virial_numerical%pv_virial(i, 1:3), i=1, 3)
254 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
255 91 : "DEBUG| Analytical pv_virial [a.u.]"
256 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
257 364 : ("DEBUG|", virial_analytical%pv_virial(i, 1:3), i=1, 3)
258 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
259 91 : "DEBUG| Difference pv_virial [a.u.]"
260 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
261 1183 : ("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 91 : "DEBUG| Sum of differences", &
264 1274 : SUM(ABS(virial_numerical%pv_virial(:, :) - virial_analytical%pv_virial(:, :)))
265 91 : IF (n_periodic > 0 .AND. n_periodic < 3) THEN
266 24 : periodic_stress_sum = 0.0_dp
267 96 : DO i = 1, 3
268 312 : DO k = 1, 3
269 288 : 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 69 : 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 24 : "DEBUG| Periodic-subspace sum of differences", periodic_stress_sum
278 : END IF
279 : END IF
280 :
281 : ! Check and abort on failure
282 182 : check_failed = .FALSE.
283 182 : IF (iw > 0) THEN
284 : WRITE (UNIT=iw, FMT="(/T2,A)") &
285 91 : "DEBUG| Relative error pv_virial"
286 : WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.8)") &
287 91 : "DEBUG| Threshold value for error check [a.u.]", eps_no_error_check
288 : END IF
289 728 : DO i = 1, 3
290 546 : err(:) = 0.0_dp
291 2184 : DO k = 1, 3
292 1638 : IF (check_stress_element(i, k) .AND. &
293 546 : 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 976 : virial_analytical%pv_virial(i, k)
296 976 : WRITE (UNIT=line(20*(k - 1) + 1:20*k), FMT="(1X,F17.2,A2)") err(k), " %"
297 : ELSE
298 662 : WRITE (UNIT=line(20*(k - 1) + 1:20*k), FMT="(17X,A3)") "- %"
299 : END IF
300 : END DO
301 546 : IF (iw > 0) THEN
302 : WRITE (UNIT=iw, FMT="(T2,A,T21,A60)") &
303 273 : "DEBUG|", line
304 : END IF
305 2360 : IF (ANY(ABS(err(1:3)) > maxerr)) check_failed = .TRUE.
306 : END DO
307 182 : IF (iw > 0) THEN
308 : WRITE (UNIT=iw, FMT="(T2,A,T61,F18.2,A2)") &
309 91 : "DEBUG| Maximum accepted error", maxerr, " %"
310 : END IF
311 83356 : IF (check_failed) THEN
312 : message = "A mismatch between the analytical and the numerical "// &
313 : "stress tensor has been detected. Check the implementation "// &
314 2 : "of the stress tensor"
315 2 : IF (stop_on_mismatch) THEN
316 0 : CPABORT(TRIM(message))
317 : ELSE
318 2 : CPWARN(TRIM(message))
319 : END IF
320 : END IF
321 : END BLOCK
322 : END IF
323 : END IF
324 :
325 854 : IF (debug_forces) THEN
326 : ! Debug forces (numerical vs analytical)
327 566 : particles => subsys%particles%els
328 1038 : SELECT CASE (force_env%in_use)
329 : CASE (use_qs_force)
330 472 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
331 472 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
332 : CASE DEFAULT
333 566 : 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 566 : calc_stress_tensor=.FALSE.)
339 : ! Copy forces in array and start the numerical calculation
340 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
341 566 : np = subsys%particles%n_els
342 1698 : ALLOCATE (analyt_forces(np, 3))
343 2700 : DO ip = 1, np
344 9102 : 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 1132 : ALLOCATE (numer_forces(subsys%particles%n_els, 3))
349 2700 : Atom: DO ip = 1, np
350 8536 : Coord: DO k = 1, 3
351 8536 : IF (do_dof_atom_coor(k, ip)) THEN
352 4666 : numer_energy = 0.0_dp
353 4666 : std_value = particles(ip)%r(k)
354 13998 : DO j = 1, 2
355 9332 : particles(ip)%r(k) = std_value - (-1.0_dp)**j*dx
356 15640 : SELECT CASE (force_env%in_use)
357 : CASE (use_qs_force)
358 6308 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
359 6308 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
360 : CASE DEFAULT
361 9332 : 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 9332 : consistent_energies=.TRUE.)
368 13998 : CALL force_env_get(force_env, potential_energy=numer_energy(j))
369 : END DO
370 4666 : particles(ip)%r(k) = std_value
371 4666 : numer_forces(ip, k) = -0.5_dp*(numer_energy(1) - numer_energy(2))/dx
372 4666 : 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 2333 : "DEBUG| Atom", "E("//ACHAR(119 + k)//" +", dx, ")", &
375 2333 : "E("//ACHAR(119 + k)//" -", dx, ")", &
376 4666 : "f(numerical)", "f(analytical)"
377 : WRITE (UNIT=iw, FMT="(T2,A,I5,4(1X,F16.8))") &
378 2333 : "DEBUG|", ip, numer_energy(1:2), numer_forces(ip, k), analyt_forces(ip, k)
379 : END IF
380 : ELSE
381 1736 : 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 8536 : my_maxerr = maxerr
386 2134 : err(1:3) = 0.0_dp
387 8536 : DO k = 1, 3
388 8536 : IF (do_dof_atom_coor(k, ip)) THEN
389 : ! Calculate percentage but ignore very small force values
390 4666 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
391 3520 : 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 4666 : IF (ABS(analyt_forces(ip, k)) <= 0.0001_dp) my_maxerr(k) = 5.0_dp*my_maxerr(k)
395 : ELSE
396 1736 : err(k) = 0.0_dp
397 : END IF
398 : END DO
399 2134 : IF (iw > 0) THEN
400 1920 : IF (ANY(do_dof_atom_coor(1:3, ip))) THEN
401 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
402 845 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
403 : END IF
404 4268 : DO k = 1, 3
405 4268 : IF (do_dof_atom_coor(k, ip)) THEN
406 2333 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
407 2333 : 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 1760 : "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 573 : "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 9102 : 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 0 : "analytical force calculation"
421 0 : IF (stop_on_mismatch) THEN
422 0 : CPABORT(message)
423 : ELSE
424 0 : CPWARN(message)
425 : END IF
426 : END IF
427 : END DO Atom
428 : ! Print summary
429 566 : IF (iw > 0) THEN
430 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
431 283 : "DEBUG|======================== BEGIN OF SUMMARY ===============================", &
432 566 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
433 283 : sum_of_differences = 0.0_dp
434 283 : errmax = 0.0_dp
435 1350 : DO ip = 1, np
436 1067 : err(1:3) = 0.0_dp
437 4551 : DO k = 1, 3
438 4268 : IF (do_dof_atom_coor(k, ip)) THEN
439 2333 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
440 2333 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
441 1760 : err(k) = 100_dp*(numer_forces(ip, k) - analyt_forces(ip, k))/analyt_forces(ip, k)
442 1760 : 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 1760 : "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 573 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, "-"
448 : END IF
449 2333 : 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 283 : "DEBUG| Sum of differences:", sum_of_differences, errmax
455 : WRITE (UNIT=iw, FMT="(T2,A)") &
456 283 : "DEBUG|======================== END OF SUMMARY ================================="
457 : END IF
458 : ! Release work storage
459 566 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
460 566 : IF (ASSOCIATED(numer_forces)) DEALLOCATE (numer_forces)
461 566 : DEALLOCATE (do_dof_atom_coor)
462 : END IF
463 :
464 854 : IF (debug_dipole) THEN
465 122 : IF (force_env%in_use == use_qs_force) THEN
466 122 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
467 122 : poldir = [0.0_dp, 0.0_dp, 1.0_dp]
468 122 : amplitude = 0.0_dp
469 122 : CALL set_efield(dft_control, amplitude, poldir)
470 122 : CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
471 122 : CALL get_qs_env(force_env%qs_env, results=results)
472 122 : description = "[DIPOLE]"
473 122 : IF (test_for_result(results, description=description)) THEN
474 122 : 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 122 : amplitude = de
480 488 : DO k = 1, 3
481 488 : IF (do_dof_dipole(k)) THEN
482 242 : poldir = 0.0_dp
483 242 : poldir(k) = 1.0_dp
484 726 : DO j = 1, 2
485 1936 : poldir = -1.0_dp*poldir
486 484 : CALL set_efield(dft_control, amplitude, poldir)
487 484 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
488 726 : CALL force_env_get(force_env, potential_energy=numer_energy(j))
489 : END DO
490 242 : 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 122 : IF (iw > 0) THEN
496 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
497 61 : "DEBUG|========================= DIPOLE MOMENTS ================================", &
498 122 : "DEBUG| Coordinate D(numerical) D(analytical) Difference Error [%]"
499 61 : err(1:3) = 0.0_dp
500 244 : DO k = 1, 3
501 244 : IF (do_dof_dipole(k)) THEN
502 121 : dd = dipole_moment(k) - dipole_numer(k)
503 121 : IF (ABS(dipole_moment(k)) > eps_no_error_check) THEN
504 62 : 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 62 : "DEBUG|", ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd, derr
507 : ELSE
508 59 : derr = 0.0_dp
509 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
510 59 : "DEBUG|", ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd
511 : END IF
512 121 : 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 61 : "DEBUG|========================================================================="
520 244 : WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") 'DIPOLE : CheckSum =', SUM(dipole_moment)
521 241 : 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 854 : IF (debug_polar) THEN
539 58 : IF (force_env%in_use == use_qs_force) THEN
540 58 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
541 58 : poldir = [0.0_dp, 0.0_dp, 1.0_dp]
542 58 : amplitude = 0.0_dp
543 58 : CALL set_efield(dft_control, amplitude, poldir)
544 58 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
545 58 : CALL get_qs_env(force_env%qs_env, results=results)
546 58 : description = "[POLAR]"
547 58 : IF (test_for_result(results, description=description)) THEN
548 58 : CALL get_results(results, description=description, values=pvals)
549 58 : 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 58 : description = "[DIPOLE]"
555 58 : 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 58 : amplitude = de
560 232 : DO k = 1, 3
561 174 : poldir = 0.0_dp
562 174 : poldir(k) = 1.0_dp
563 522 : DO j = 1, 2
564 1392 : poldir = -1.0_dp*poldir
565 348 : CALL set_efield(dft_control, amplitude, poldir)
566 348 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., linres=.TRUE.)
567 522 : CALL get_results(results, description=description, values=dipn(1:3, j))
568 : END DO
569 754 : polar_numeric(1:3, k) = 0.5_dp*(dipn(1:3, 2) - dipn(1:3, 1))/de
570 : END DO
571 58 : IF (iw > 0) THEN
572 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
573 29 : "DEBUG|========================= POLARIZABILITY ================================", &
574 58 : "DEBUG| Coordinates P(numerical) P(analytical) Difference Error [%]"
575 116 : DO k = 1, 3
576 377 : DO j = 1, 3
577 261 : dd = polar_analytic(k, j) - polar_numeric(k, j)
578 348 : IF (ABS(polar_analytic(k, j)) > eps_no_error_check) THEN
579 179 : derr = 100._dp*dd/polar_analytic(k, j)
580 : WRITE (UNIT=iw, FMT="(T2,A,T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
581 179 : "DEBUG|", ACHAR(119 + k), ACHAR(119 + j), polar_numeric(k, j), polar_analytic(k, j), dd, derr
582 : ELSE
583 : WRITE (UNIT=iw, FMT="(T2,A,T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
584 82 : "DEBUG|", ACHAR(119 + k), ACHAR(119 + j), polar_numeric(k, j), polar_analytic(k, j), dd
585 : END IF
586 : END DO
587 : END DO
588 : WRITE (UNIT=iw, FMT="((T2,A))") &
589 29 : "DEBUG|========================================================================="
590 377 : WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") ' POLAR : CheckSum =', SUM(polar_analytic)
591 : END IF
592 : ELSE
593 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities only for Quickstep code available")
594 : END IF
595 : END IF
596 :
597 854 : CALL cp_print_key_finished_output(iw, logger, root_section, "DEBUG%PROGRAM_RUN_INFO")
598 :
599 1708 : END SUBROUTINE cp2k_debug_energy_and_forces
600 :
601 : ! **************************************************************************************************
602 : !> \brief ...
603 : !> \param dft_control ...
604 : !> \param amplitude ...
605 : !> \param poldir ...
606 : ! **************************************************************************************************
607 1012 : SUBROUTINE set_efield(dft_control, amplitude, poldir)
608 : TYPE(dft_control_type), POINTER :: dft_control
609 : REAL(KIND=dp), INTENT(IN) :: amplitude
610 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: poldir
611 :
612 1012 : IF (dft_control%apply_efield) THEN
613 928 : dft_control%efield_fields(1)%efield%strength = amplitude
614 3712 : dft_control%efield_fields(1)%efield%polarisation(1:3) = poldir(1:3)
615 84 : ELSEIF (dft_control%apply_period_efield) THEN
616 84 : dft_control%period_efield%strength = amplitude
617 336 : dft_control%period_efield%polarisation(1:3) = poldir(1:3)
618 : ELSE
619 0 : CPABORT("No EFIELD definition available")
620 : END IF
621 :
622 1012 : END SUBROUTINE set_efield
623 :
624 : END MODULE cp2k_debug
|
