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 906 : 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 906 : 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 906 : 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 906 : 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 906 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
111 906 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
112 : TYPE(section_vals_type), POINTER :: root_section, subsys_section
113 :
114 906 : NULLIFY (analyt_forces, numer_forces, subsys, particles)
115 :
116 906 : root_section => force_env%root_section
117 :
118 906 : 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 906 : "SUBSYS")
121 :
122 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_STRESS_TENSOR", &
123 906 : l_val=debug_stress_tensor)
124 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_FORCES", &
125 906 : l_val=debug_forces)
126 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_DIPOLE", &
127 906 : l_val=debug_dipole)
128 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_POLARIZABILITY", &
129 906 : l_val=debug_polar)
130 : CALL section_vals_val_get(root_section, "DEBUG%DX", &
131 906 : r_val=dx)
132 : CALL section_vals_val_get(root_section, "DEBUG%DE", &
133 906 : r_val=de)
134 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_DIPOLE_DIRS", &
135 906 : c_val=cval1)
136 906 : dx = ABS(dx)
137 : CALL section_vals_val_get(root_section, "DEBUG%MAX_RELATIVE_ERROR", &
138 906 : r_val=maxerr)
139 : CALL section_vals_val_get(root_section, "DEBUG%EPS_NO_ERROR_CHECK", &
140 906 : r_val=eps_no_error_check)
141 906 : 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 906 : l_val=stop_on_mismatch)
144 :
145 : ! set active DOF
146 906 : CALL uppercase(cval1)
147 906 : do_dof_dipole(1) = (INDEX(cval1, "X") /= 0)
148 906 : do_dof_dipole(2) = (INDEX(cval1, "Y") /= 0)
149 906 : do_dof_dipole(3) = (INDEX(cval1, "Z") /= 0)
150 906 : NULLIFY (cval2)
151 906 : IF (debug_forces) THEN
152 600 : np = subsys%particles%n_els
153 1800 : ALLOCATE (do_dof_atom_coor(3, np))
154 600 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", n_rep_val=nrep)
155 600 : IF (nrep > 0) THEN
156 238 : do_dof_atom_coor = .FALSE.
157 482 : DO irep = 1, nrep
158 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", i_rep_val=irep, &
159 244 : c_vals=cval2)
160 244 : READ (cval2(1), FMT="(I10)") k
161 244 : CALL uppercase(cval2(2))
162 244 : do_dof_atom_coor(1, k) = (INDEX(cval2(2), "X") /= 0)
163 244 : do_dof_atom_coor(2, k) = (INDEX(cval2(2), "Y") /= 0)
164 482 : do_dof_atom_coor(3, k) = (INDEX(cval2(2), "Z") /= 0)
165 : END DO
166 : ELSE
167 6314 : do_dof_atom_coor = .TRUE.
168 : END IF
169 : END IF
170 :
171 906 : logger => cp_get_default_logger()
172 : iw = cp_print_key_unit_nr(logger, root_section, "DEBUG%PROGRAM_RUN_INFO", &
173 906 : extension=".log")
174 906 : IF (debug_stress_tensor) THEN
175 1280 : 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 320 : i_val=stress_tensor)
185 320 : 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 320 : 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 200 : CALL cp_subsys_get(subsys, virial=virial)
208 200 : CALL virial_set(virial, pv_numer=.FALSE.)
209 : CALL force_env_calc_energy_force(force_env, &
210 : calc_force=.TRUE., &
211 200 : calc_stress_tensor=.TRUE.)
212 :
213 : ! Retrieve the analytical virial
214 200 : virial_analytical = virial
215 :
216 : ! Debug stress tensor (numerical vs analytical)
217 200 : CALL virial_set(virial, pv_numer=.TRUE.)
218 200 : CALL force_env_calc_num_pressure(force_env, dx=dx)
219 :
220 : ! Retrieve the numerical virial
221 200 : CALL cp_subsys_get(subsys, virial=virial)
222 200 : virial_numerical = virial
223 :
224 : ! Numerical diagonal stress checks only perturb diagonal cell elements.
225 200 : IF (virial_analytical%pv_diagonal .OR. virial_numerical%pv_diagonal) THEN
226 168 : DO i = 1, 3
227 704 : 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 200 : CALL get_cell(cell=cell, periodic=periodic)
237 800 : n_periodic = COUNT(periodic /= 0)
238 2600 : check_stress_element = .TRUE.
239 200 : 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 200 : IF (iw > 0) THEN
250 : WRITE (UNIT=iw, FMT="((T2,A))") &
251 100 : "DEBUG| Numerical pv_virial [a.u.]"
252 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
253 400 : ("DEBUG|", virial_numerical%pv_virial(i, 1:3), i=1, 3)
254 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
255 100 : "DEBUG| Analytical pv_virial [a.u.]"
256 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
257 400 : ("DEBUG|", virial_analytical%pv_virial(i, 1:3), i=1, 3)
258 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
259 100 : "DEBUG| Difference pv_virial [a.u.]"
260 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
261 1300 : ("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 100 : "DEBUG| Sum of differences", &
264 1400 : SUM(ABS(virial_numerical%pv_virial(:, :) - virial_analytical%pv_virial(:, :)))
265 100 : 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 200 : check_failed = .FALSE.
283 200 : IF (iw > 0) THEN
284 : WRITE (UNIT=iw, FMT="(/T2,A)") &
285 100 : "DEBUG| Relative error pv_virial"
286 : WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.8)") &
287 100 : "DEBUG| Threshold value for error check [a.u.]", eps_no_error_check
288 : END IF
289 800 : DO i = 1, 3
290 600 : err(:) = 0.0_dp
291 2400 : DO k = 1, 3
292 1800 : IF (check_stress_element(i, k) .AND. &
293 600 : 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 1056 : virial_analytical%pv_virial(i, k)
296 1056 : 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 600 : IF (iw > 0) THEN
302 : WRITE (UNIT=iw, FMT="(T2,A,T21,A60)") &
303 300 : "DEBUG|", line
304 : END IF
305 2594 : IF (ANY(ABS(err(1:3)) > maxerr)) check_failed = .TRUE.
306 : END DO
307 200 : IF (iw > 0) THEN
308 : WRITE (UNIT=iw, FMT="(T2,A,T61,F18.2,A2)") &
309 100 : "DEBUG| Maximum accepted error", maxerr, " %"
310 : END IF
311 91600 : 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 906 : IF (debug_forces) THEN
326 : ! Debug forces (numerical vs analytical)
327 600 : particles => subsys%particles%els
328 1106 : SELECT CASE (force_env%in_use)
329 : CASE (use_qs_force)
330 506 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
331 506 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
332 : CASE DEFAULT
333 600 : 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 600 : calc_stress_tensor=.FALSE.)
339 : ! Copy forces in array and start the numerical calculation
340 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
341 600 : np = subsys%particles%n_els
342 1800 : ALLOCATE (analyt_forces(np, 3))
343 2828 : DO ip = 1, np
344 9512 : 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 1200 : ALLOCATE (numer_forces(subsys%particles%n_els, 3))
349 2828 : Atom: DO ip = 1, np
350 8912 : Coord: DO k = 1, 3
351 8912 : IF (do_dof_atom_coor(k, ip)) THEN
352 4738 : numer_energy = 0.0_dp
353 4738 : std_value = particles(ip)%r(k)
354 14214 : DO j = 1, 2
355 9476 : particles(ip)%r(k) = std_value - (-1.0_dp)**j*dx
356 15928 : SELECT CASE (force_env%in_use)
357 : CASE (use_qs_force)
358 6452 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
359 6452 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
360 : CASE DEFAULT
361 9476 : 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 9476 : consistent_energies=.TRUE.)
368 14214 : CALL force_env_get(force_env, potential_energy=numer_energy(j))
369 : END DO
370 4738 : particles(ip)%r(k) = std_value
371 4738 : numer_forces(ip, k) = -0.5_dp*(numer_energy(1) - numer_energy(2))/dx
372 4738 : 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 2369 : "DEBUG| Atom", "E("//ACHAR(119 + k)//" +", dx, ")", &
375 2369 : "E("//ACHAR(119 + k)//" -", dx, ")", &
376 4738 : "f(numerical)", "f(analytical)"
377 : WRITE (UNIT=iw, FMT="(T2,A,I5,4(1X,F16.8))") &
378 2369 : "DEBUG|", ip, numer_energy(1:2), numer_forces(ip, k), analyt_forces(ip, k)
379 : END IF
380 : ELSE
381 1946 : 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 8912 : my_maxerr = maxerr
386 2228 : err(1:3) = 0.0_dp
387 8912 : DO k = 1, 3
388 8912 : IF (do_dof_atom_coor(k, ip)) THEN
389 : ! Calculate percentage but ignore very small force values
390 4738 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
391 3578 : 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 4738 : IF (ABS(analyt_forces(ip, k)) <= 0.0001_dp) my_maxerr(k) = 5.0_dp*my_maxerr(k)
395 : ELSE
396 1946 : err(k) = 0.0_dp
397 : END IF
398 : END DO
399 2228 : IF (iw > 0) THEN
400 2059 : IF (ANY(do_dof_atom_coor(1:3, ip))) THEN
401 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
402 866 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
403 : END IF
404 4456 : DO k = 1, 3
405 4456 : IF (do_dof_atom_coor(k, ip)) THEN
406 2369 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
407 2369 : 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 1789 : "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 580 : "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 9510 : 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 600 : IF (iw > 0) THEN
430 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
431 300 : "DEBUG|======================== BEGIN OF SUMMARY ===============================", &
432 600 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
433 300 : sum_of_differences = 0.0_dp
434 300 : errmax = 0.0_dp
435 1414 : DO ip = 1, np
436 1114 : err(1:3) = 0.0_dp
437 4756 : DO k = 1, 3
438 4456 : IF (do_dof_atom_coor(k, ip)) THEN
439 2369 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
440 2369 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
441 1789 : err(k) = 100_dp*(numer_forces(ip, k) - analyt_forces(ip, k))/analyt_forces(ip, k)
442 1789 : 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 1789 : "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 580 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, "-"
448 : END IF
449 2369 : 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 300 : "DEBUG| Sum of differences:", sum_of_differences, errmax
455 : WRITE (UNIT=iw, FMT="(T2,A)") &
456 300 : "DEBUG|======================== END OF SUMMARY ================================="
457 : END IF
458 : ! Release work storage
459 600 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
460 600 : IF (ASSOCIATED(numer_forces)) DEALLOCATE (numer_forces)
461 600 : DEALLOCATE (do_dof_atom_coor)
462 : END IF
463 :
464 906 : 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 906 : IF (debug_polar) THEN
539 66 : IF (force_env%in_use == use_qs_force) THEN
540 66 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
541 66 : poldir = [0.0_dp, 0.0_dp, 1.0_dp]
542 66 : amplitude = 0.0_dp
543 66 : CALL set_efield(dft_control, amplitude, poldir)
544 66 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
545 66 : CALL get_qs_env(force_env%qs_env, results=results)
546 66 : description = "[POLAR]"
547 66 : IF (test_for_result(results, description=description)) THEN
548 66 : CALL get_results(results, description=description, values=pvals)
549 66 : 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 66 : description = "[DIPOLE]"
555 66 : 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 66 : amplitude = de
560 264 : DO k = 1, 3
561 198 : poldir = 0.0_dp
562 198 : poldir(k) = 1.0_dp
563 594 : DO j = 1, 2
564 1584 : poldir = -1.0_dp*poldir
565 396 : CALL set_efield(dft_control, amplitude, poldir)
566 396 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., linres=.TRUE.)
567 594 : CALL get_results(results, description=description, values=dipn(1:3, j))
568 : END DO
569 858 : polar_numeric(1:3, k) = 0.5_dp*(dipn(1:3, 2) - dipn(1:3, 1))/de
570 : END DO
571 66 : IF (iw > 0) THEN
572 33 : polerr = 0.0_dp
573 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
574 33 : "DEBUG|========================= POLARIZABILITY ================================", &
575 66 : "DEBUG| Coordinates P(numerical) P(analytical) Difference Error [%]"
576 132 : DO k = 1, 3
577 429 : DO j = 1, 3
578 297 : dd = polar_analytic(k, j) - polar_numeric(k, j)
579 396 : IF (ABS(polar_analytic(k, j)) > eps_no_error_check) THEN
580 142 : 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 142 : "DEBUG|", ACHAR(119 + k), ACHAR(119 + j), polar_numeric(k, j), polar_analytic(k, j), dd, derr
583 142 : 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 155 : "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 33 : "DEBUG|========================================================================="
592 429 : WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") ' POLAR : CheckSum =', SUM(polar_analytic)
593 429 : 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 906 : CALL cp_print_key_finished_output(iw, logger, root_section, "DEBUG%PROGRAM_RUN_INFO")
610 :
611 1812 : END SUBROUTINE cp2k_debug_energy_and_forces
612 :
613 : ! **************************************************************************************************
614 : !> \brief ...
615 : !> \param dft_control ...
616 : !> \param amplitude ...
617 : !> \param poldir ...
618 : ! **************************************************************************************************
619 1068 : 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 1068 : IF (dft_control%apply_efield) THEN
625 984 : dft_control%efield_fields(1)%efield%strength = amplitude
626 3936 : dft_control%efield_fields(1)%efield%polarisation(1:3) = poldir(1:3)
627 84 : ELSE IF (dft_control%apply_period_efield) THEN
628 84 : dft_control%period_efield%strength = amplitude
629 336 : dft_control%period_efield%polarisation(1:3) = poldir(1:3)
630 : ELSE
631 0 : CPABORT("No EFIELD definition available")
632 : END IF
633 :
634 1068 : END SUBROUTINE set_efield
635 :
636 : END MODULE cp2k_debug
|
