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 : USE cp_control_types, ONLY: dft_control_type
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
32 : cp_print_key_unit_nr
33 : USE cp_result_methods, ONLY: get_results,&
34 : test_for_result
35 : USE cp_result_types, ONLY: cp_result_type
36 : USE cp_subsys_types, ONLY: cp_subsys_get,&
37 : cp_subsys_type
38 : USE force_env_methods, ONLY: force_env_calc_energy_force,&
39 : force_env_calc_num_pressure
40 : USE force_env_types, ONLY: force_env_get,&
41 : force_env_type,&
42 : use_qs_force
43 : USE input_constants, ONLY: do_stress_analytical,&
44 : do_stress_diagonal_anal,&
45 : do_stress_diagonal_numer,&
46 : do_stress_numerical
47 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
48 : section_vals_type,&
49 : section_vals_val_get
50 : USE kinds, ONLY: default_string_length,&
51 : dp
52 : USE message_passing, ONLY: mp_para_env_type
53 : USE particle_methods, ONLY: write_fist_particle_coordinates,&
54 : write_qs_particle_coordinates
55 : USE particle_types, ONLY: particle_type
56 : USE qs_environment_types, ONLY: get_qs_env
57 : USE qs_kind_types, ONLY: qs_kind_type
58 : USE string_utilities, ONLY: uppercase
59 : USE virial_types, ONLY: virial_set,&
60 : virial_type
61 : #include "./base/base_uses.f90"
62 :
63 : IMPLICIT NONE
64 : PRIVATE
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp2k_debug'
66 :
67 : PUBLIC :: cp2k_debug_energy_and_forces
68 :
69 : CONTAINS
70 :
71 : ! **************************************************************************************************
72 : !> \brief ...
73 : !> \param force_env ...
74 : ! **************************************************************************************************
75 736 : SUBROUTINE cp2k_debug_energy_and_forces(force_env)
76 :
77 : TYPE(force_env_type), POINTER :: force_env
78 :
79 : CHARACTER(LEN=3) :: cval1
80 : CHARACTER(LEN=3*default_string_length) :: message
81 : CHARACTER(LEN=60) :: line
82 736 : CHARACTER(LEN=80), DIMENSION(:), POINTER :: cval2
83 : CHARACTER(LEN=default_string_length) :: description
84 : INTEGER :: i, ip, irep, iw, j, k, np, nrep, &
85 : stress_tensor
86 : LOGICAL :: check_failed, debug_dipole, &
87 : debug_forces, debug_polar, &
88 : debug_stress_tensor, skip, &
89 : stop_on_mismatch
90 736 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: do_dof_atom_coor
91 : LOGICAL, DIMENSION(3) :: do_dof_dipole
92 : REAL(KIND=dp) :: amplitude, dd, de, derr, difference, dx, &
93 : eps_no_error_check, errmax, maxerr, &
94 : std_value, sum_of_differences
95 : REAL(KIND=dp), DIMENSION(2) :: numer_energy
96 : REAL(KIND=dp), DIMENSION(3) :: dipole_moment, dipole_numer, err, &
97 : my_maxerr, poldir
98 : REAL(KIND=dp), DIMENSION(3, 2) :: dipn
99 : REAL(KIND=dp), DIMENSION(3, 3) :: polar_analytic, polar_numeric
100 : REAL(KIND=dp), DIMENSION(9) :: pvals
101 736 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: analyt_forces, numer_forces
102 : TYPE(cell_type), POINTER :: cell
103 : TYPE(cp_logger_type), POINTER :: logger
104 : TYPE(cp_result_type), POINTER :: results
105 : TYPE(cp_subsys_type), POINTER :: subsys
106 : TYPE(dft_control_type), POINTER :: dft_control
107 : TYPE(mp_para_env_type), POINTER :: para_env
108 736 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
109 736 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
110 : TYPE(section_vals_type), POINTER :: root_section, subsys_section
111 :
112 736 : NULLIFY (analyt_forces, numer_forces, subsys, particles)
113 :
114 736 : root_section => force_env%root_section
115 :
116 736 : CALL force_env_get(force_env, para_env=para_env, subsys=subsys, cell=cell)
117 : subsys_section => section_vals_get_subs_vals(force_env%force_env_section, &
118 736 : "SUBSYS")
119 :
120 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_STRESS_TENSOR", &
121 736 : l_val=debug_stress_tensor)
122 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_FORCES", &
123 736 : l_val=debug_forces)
124 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_DIPOLE", &
125 736 : l_val=debug_dipole)
126 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_POLARIZABILITY", &
127 736 : l_val=debug_polar)
128 : CALL section_vals_val_get(root_section, "DEBUG%DX", &
129 736 : r_val=dx)
130 : CALL section_vals_val_get(root_section, "DEBUG%DE", &
131 736 : r_val=de)
132 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_DIPOLE_DIRS", &
133 736 : c_val=cval1)
134 736 : dx = ABS(dx)
135 : CALL section_vals_val_get(root_section, "DEBUG%MAX_RELATIVE_ERROR", &
136 736 : r_val=maxerr)
137 : CALL section_vals_val_get(root_section, "DEBUG%EPS_NO_ERROR_CHECK", &
138 736 : r_val=eps_no_error_check)
139 736 : eps_no_error_check = MAX(eps_no_error_check, EPSILON(0.0_dp))
140 : CALL section_vals_val_get(root_section, "DEBUG%STOP_ON_MISMATCH", &
141 736 : l_val=stop_on_mismatch)
142 :
143 : ! set active DOF
144 736 : CALL uppercase(cval1)
145 736 : do_dof_dipole(1) = (INDEX(cval1, "X") /= 0)
146 736 : do_dof_dipole(2) = (INDEX(cval1, "Y") /= 0)
147 736 : do_dof_dipole(3) = (INDEX(cval1, "Z") /= 0)
148 736 : NULLIFY (cval2)
149 736 : IF (debug_forces) THEN
150 478 : np = subsys%particles%n_els
151 1434 : ALLOCATE (do_dof_atom_coor(3, np))
152 478 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", n_rep_val=nrep)
153 478 : IF (nrep > 0) THEN
154 2674 : do_dof_atom_coor = .FALSE.
155 394 : DO irep = 1, nrep
156 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", i_rep_val=irep, &
157 200 : c_vals=cval2)
158 200 : READ (cval2(1), FMT="(I10)") k
159 200 : CALL uppercase(cval2(2))
160 200 : do_dof_atom_coor(1, k) = (INDEX(cval2(2), "X") /= 0)
161 200 : do_dof_atom_coor(2, k) = (INDEX(cval2(2), "Y") /= 0)
162 394 : do_dof_atom_coor(3, k) = (INDEX(cval2(2), "Z") /= 0)
163 : END DO
164 : ELSE
165 5276 : do_dof_atom_coor = .TRUE.
166 : END IF
167 : END IF
168 :
169 736 : logger => cp_get_default_logger()
170 : iw = cp_print_key_unit_nr(logger, root_section, "DEBUG%PROGRAM_RUN_INFO", &
171 736 : extension=".log")
172 736 : IF (debug_stress_tensor) THEN
173 896 : IF (SUM(cell%perd) == 0) THEN
174 : CALL cp_warn(__LOCATION__, &
175 : "DEBUG_STRESS_TENSOR requested for PERIODIC NONE. "// &
176 : "The isolated-system virial is useful for finite-difference diagnostics, "// &
177 40 : "but it is not a physically meaningful bulk stress.")
178 : END IF
179 : ! To debug stress tensor the stress tensor calculation must be
180 : ! first enabled..
181 : CALL section_vals_val_get(force_env%force_env_section, "STRESS_TENSOR", &
182 224 : i_val=stress_tensor)
183 224 : skip = .FALSE.
184 0 : SELECT CASE (stress_tensor)
185 : CASE (do_stress_analytical, do_stress_diagonal_anal)
186 : ! OK
187 : CASE (do_stress_numerical, do_stress_diagonal_numer)
188 : ! Nothing to check
189 : CALL cp_warn(__LOCATION__, "Numerical stress tensor was requested in "// &
190 0 : "the FORCE_EVAL section. Nothing to debug!")
191 120 : skip = .TRUE.
192 : CASE DEFAULT
193 : CALL cp_warn(__LOCATION__, "Stress tensor calculation was not enabled in "// &
194 120 : "the FORCE_EVAL section. Nothing to debug!")
195 224 : skip = .TRUE.
196 : END SELECT
197 :
198 : IF (.NOT. skip) THEN
199 :
200 : BLOCK
201 : TYPE(virial_type) :: virial_analytical, virial_numerical
202 : TYPE(virial_type), POINTER :: virial
203 :
204 : ! Compute the analytical stress tensor
205 104 : CALL cp_subsys_get(subsys, virial=virial)
206 104 : CALL virial_set(virial, pv_numer=.FALSE.)
207 : CALL force_env_calc_energy_force(force_env, &
208 : calc_force=.TRUE., &
209 104 : calc_stress_tensor=.TRUE.)
210 :
211 : ! Retrieve the analytical virial
212 104 : virial_analytical = virial
213 :
214 : ! Debug stress tensor (numerical vs analytical)
215 104 : CALL virial_set(virial, pv_numer=.TRUE.)
216 104 : CALL force_env_calc_num_pressure(force_env, dx=dx)
217 :
218 : ! Retrieve the numerical virial
219 104 : CALL cp_subsys_get(subsys, virial=virial)
220 104 : virial_numerical = virial
221 :
222 : ! Print results
223 104 : IF (iw > 0) THEN
224 : WRITE (UNIT=iw, FMT="((T2,A))") &
225 52 : "DEBUG| Numerical pv_virial [a.u.]"
226 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
227 208 : ("DEBUG|", virial_numerical%pv_virial(i, 1:3), i=1, 3)
228 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
229 52 : "DEBUG| Analytical pv_virial [a.u.]"
230 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
231 208 : ("DEBUG|", virial_analytical%pv_virial(i, 1:3), i=1, 3)
232 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
233 52 : "DEBUG| Difference pv_virial [a.u.]"
234 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
235 676 : ("DEBUG|", virial_numerical%pv_virial(i, 1:3) - virial_analytical%pv_virial(i, 1:3), i=1, 3)
236 : WRITE (UNIT=iw, FMT="(T2,A,T61,F20.12)") &
237 52 : "DEBUG| Sum of differences", &
238 728 : SUM(ABS(virial_numerical%pv_virial(:, :) - virial_analytical%pv_virial(:, :)))
239 : END IF
240 :
241 : ! Check and abort on failure
242 104 : check_failed = .FALSE.
243 104 : IF (iw > 0) THEN
244 : WRITE (UNIT=iw, FMT="(/T2,A)") &
245 52 : "DEBUG| Relative error pv_virial"
246 : WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.8)") &
247 52 : "DEBUG| Threshold value for error check [a.u.]", eps_no_error_check
248 : END IF
249 416 : DO i = 1, 3
250 312 : err(:) = 0.0_dp
251 1248 : DO k = 1, 3
252 1248 : IF (ABS(virial_analytical%pv_virial(i, k)) >= eps_no_error_check) THEN
253 : err(k) = 100.0_dp*(virial_numerical%pv_virial(i, k) - virial_analytical%pv_virial(i, k))/ &
254 700 : virial_analytical%pv_virial(i, k)
255 700 : WRITE (UNIT=line(20*(k - 1) + 1:20*k), FMT="(1X,F17.2,A2)") err(k), " %"
256 : ELSE
257 236 : WRITE (UNIT=line(20*(k - 1) + 1:20*k), FMT="(17X,A3)") "- %"
258 : END IF
259 : END DO
260 312 : IF (iw > 0) THEN
261 : WRITE (UNIT=iw, FMT="(T2,A,T21,A60)") &
262 156 : "DEBUG|", line
263 : END IF
264 1352 : IF (ANY(ABS(err(1:3)) > maxerr)) check_failed = .TRUE.
265 : END DO
266 104 : IF (iw > 0) THEN
267 : WRITE (UNIT=iw, FMT="(T2,A,T61,F18.2,A2)") &
268 52 : "DEBUG| Maximum accepted error", maxerr, " %"
269 : END IF
270 47632 : IF (check_failed) THEN
271 : message = "A mismatch between the analytical and the numerical "// &
272 : "stress tensor has been detected. Check the implementation "// &
273 0 : "of the stress tensor"
274 0 : IF (stop_on_mismatch) THEN
275 0 : CPABORT(TRIM(message))
276 : ELSE
277 0 : CPWARN(TRIM(message))
278 : END IF
279 : END IF
280 : END BLOCK
281 : END IF
282 : END IF
283 :
284 736 : IF (debug_forces) THEN
285 : ! Debug forces (numerical vs analytical)
286 478 : particles => subsys%particles%els
287 862 : SELECT CASE (force_env%in_use)
288 : CASE (use_qs_force)
289 384 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
290 384 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
291 : CASE DEFAULT
292 478 : CALL write_fist_particle_coordinates(particles, subsys_section)
293 : END SELECT
294 : ! First evaluate energy and forces
295 : CALL force_env_calc_energy_force(force_env, &
296 : calc_force=.TRUE., &
297 478 : calc_stress_tensor=.FALSE.)
298 : ! Copy forces in array and start the numerical calculation
299 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
300 478 : np = subsys%particles%n_els
301 1434 : ALLOCATE (analyt_forces(np, 3))
302 2346 : DO ip = 1, np
303 7950 : analyt_forces(ip, 1:3) = particles(ip)%f(1:3)
304 : END DO
305 : ! Loop on atoms and coordinates
306 : IF (ASSOCIATED(numer_forces)) DEALLOCATE (numer_forces)
307 956 : ALLOCATE (numer_forces(subsys%particles%n_els, 3))
308 2346 : Atom: DO ip = 1, np
309 7472 : Coord: DO k = 1, 3
310 7472 : IF (do_dof_atom_coor(k, ip)) THEN
311 3972 : numer_energy = 0.0_dp
312 3972 : std_value = particles(ip)%r(k)
313 11916 : DO j = 1, 2
314 7944 : particles(ip)%r(k) = std_value - (-1.0_dp)**j*dx
315 12864 : SELECT CASE (force_env%in_use)
316 : CASE (use_qs_force)
317 4920 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
318 4920 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
319 : CASE DEFAULT
320 7944 : CALL write_fist_particle_coordinates(particles, subsys_section)
321 : END SELECT
322 : ! Compute energy
323 : CALL force_env_calc_energy_force(force_env, &
324 : calc_force=.FALSE., &
325 : calc_stress_tensor=.FALSE., &
326 7944 : consistent_energies=.TRUE.)
327 11916 : CALL force_env_get(force_env, potential_energy=numer_energy(j))
328 : END DO
329 3972 : particles(ip)%r(k) = std_value
330 3972 : numer_forces(ip, k) = -0.5_dp*(numer_energy(1) - numer_energy(2))/dx
331 3972 : IF (iw > 0) THEN
332 : WRITE (UNIT=iw, FMT="(/,T2,A,T17,A,F7.4,A,T34,A,F7.4,A,T52,A,T68,A)") &
333 1986 : "DEBUG| Atom", "E("//ACHAR(119 + k)//" +", dx, ")", &
334 1986 : "E("//ACHAR(119 + k)//" -", dx, ")", &
335 3972 : "f(numerical)", "f(analytical)"
336 : WRITE (UNIT=iw, FMT="(T2,A,I5,4(1X,F16.8))") &
337 1986 : "DEBUG|", ip, numer_energy(1:2), numer_forces(ip, k), analyt_forces(ip, k)
338 : END IF
339 : ELSE
340 1632 : numer_forces(ip, k) = 0.0_dp
341 : END IF
342 : END DO Coord
343 : ! Check analytical forces using the numerical forces as reference
344 7472 : my_maxerr = maxerr
345 1868 : err(1:3) = 0.0_dp
346 7472 : DO k = 1, 3
347 7472 : IF (do_dof_atom_coor(k, ip)) THEN
348 : ! Calculate percentage but ignore very small force values
349 3972 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
350 2998 : err(k) = 100.0_dp*(numer_forces(ip, k) - analyt_forces(ip, k))/analyt_forces(ip, k)
351 : END IF
352 : ! Increase error tolerance for small force values
353 3972 : IF (ABS(analyt_forces(ip, k)) <= 0.0001_dp) my_maxerr(k) = 5.0_dp*my_maxerr(k)
354 : ELSE
355 1632 : err(k) = 0.0_dp
356 : END IF
357 : END DO
358 1868 : IF (iw > 0) THEN
359 1737 : IF (ANY(do_dof_atom_coor(1:3, ip))) THEN
360 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
361 724 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
362 : END IF
363 3736 : DO k = 1, 3
364 3736 : IF (do_dof_atom_coor(k, ip)) THEN
365 1986 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
366 1986 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
367 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T71,F10.2)") &
368 1499 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, err(k)
369 : ELSE
370 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T80,A1)") &
371 487 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, "-"
372 : END IF
373 : END IF
374 : END DO
375 : END IF
376 7950 : IF (ANY(ABS(err(1:3)) > my_maxerr(1:3))) THEN
377 : message = "A mismatch between analytical and numerical forces "// &
378 : "has been detected. Check the implementation of the "// &
379 0 : "analytical force calculation"
380 0 : IF (stop_on_mismatch) THEN
381 0 : CPABORT(message)
382 : ELSE
383 0 : CPWARN(message)
384 : END IF
385 : END IF
386 : END DO Atom
387 : ! Print summary
388 478 : IF (iw > 0) THEN
389 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
390 239 : "DEBUG|======================== BEGIN OF SUMMARY ===============================", &
391 478 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
392 239 : sum_of_differences = 0.0_dp
393 239 : errmax = 0.0_dp
394 1173 : DO ip = 1, np
395 934 : err(1:3) = 0.0_dp
396 3975 : DO k = 1, 3
397 3736 : IF (do_dof_atom_coor(k, ip)) THEN
398 1986 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
399 1986 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
400 1499 : err(k) = 100_dp*(numer_forces(ip, k) - analyt_forces(ip, k))/analyt_forces(ip, k)
401 1499 : errmax = MAX(errmax, ABS(err(k)))
402 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T71,F10.2)") &
403 1499 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, err(k)
404 : ELSE
405 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T80,A1)") &
406 487 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, "-"
407 : END IF
408 1986 : sum_of_differences = sum_of_differences + ABS(difference)
409 : END IF
410 : END DO
411 : END DO
412 : WRITE (UNIT=iw, FMT="(T2,A,T57,F12.8,T71,F10.2)") &
413 239 : "DEBUG| Sum of differences:", sum_of_differences, errmax
414 : WRITE (UNIT=iw, FMT="(T2,A)") &
415 239 : "DEBUG|======================== END OF SUMMARY ================================="
416 : END IF
417 : ! Release work storage
418 478 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
419 478 : IF (ASSOCIATED(numer_forces)) DEALLOCATE (numer_forces)
420 478 : DEALLOCATE (do_dof_atom_coor)
421 : END IF
422 :
423 736 : IF (debug_dipole) THEN
424 122 : IF (force_env%in_use == use_qs_force) THEN
425 122 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
426 122 : poldir = [0.0_dp, 0.0_dp, 1.0_dp]
427 122 : amplitude = 0.0_dp
428 122 : CALL set_efield(dft_control, amplitude, poldir)
429 122 : CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
430 122 : CALL get_qs_env(force_env%qs_env, results=results)
431 122 : description = "[DIPOLE]"
432 122 : IF (test_for_result(results, description=description)) THEN
433 122 : CALL get_results(results, description=description, values=dipole_moment)
434 : ELSE
435 0 : CALL cp_warn(__LOCATION__, "Debug of dipole moments needs DFT/PRINT/MOMENTS section enabled")
436 0 : CPABORT("DEBUG")
437 : END IF
438 122 : amplitude = de
439 488 : DO k = 1, 3
440 488 : IF (do_dof_dipole(k)) THEN
441 242 : poldir = 0.0_dp
442 242 : poldir(k) = 1.0_dp
443 726 : DO j = 1, 2
444 1936 : poldir = -1.0_dp*poldir
445 484 : CALL set_efield(dft_control, amplitude, poldir)
446 484 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
447 726 : CALL force_env_get(force_env, potential_energy=numer_energy(j))
448 : END DO
449 242 : dipole_numer(k) = 0.5_dp*(numer_energy(1) - numer_energy(2))/de
450 : ELSE
451 124 : dipole_numer(k) = 0.0_dp
452 : END IF
453 : END DO
454 122 : IF (iw > 0) THEN
455 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
456 61 : "DEBUG|========================= DIPOLE MOMENTS ================================", &
457 122 : "DEBUG| Coordinate D(numerical) D(analytical) Difference Error [%]"
458 61 : err(1:3) = 0.0_dp
459 244 : DO k = 1, 3
460 244 : IF (do_dof_dipole(k)) THEN
461 121 : dd = dipole_moment(k) - dipole_numer(k)
462 121 : IF (ABS(dipole_moment(k)) > eps_no_error_check) THEN
463 62 : derr = 100._dp*dd/dipole_moment(k)
464 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
465 62 : "DEBUG|", ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd, derr
466 : ELSE
467 59 : derr = 0.0_dp
468 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
469 59 : "DEBUG|", ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd
470 : END IF
471 121 : err(k) = derr
472 : ELSE
473 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,A16,T38,F16.8)") &
474 62 : "DEBUG|", ACHAR(119 + k), " skipped", dipole_moment(k)
475 : END IF
476 : END DO
477 : WRITE (UNIT=iw, FMT="((T2,A))") &
478 61 : "DEBUG|========================================================================="
479 244 : WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") 'DIPOLE : CheckSum =', SUM(dipole_moment)
480 241 : IF (ANY(ABS(err(1:3)) > maxerr)) THEN
481 : message = "A mismatch between analytical and numerical dipoles "// &
482 : "has been detected. Check the implementation of the "// &
483 1 : "analytical dipole calculation"
484 1 : IF (stop_on_mismatch) THEN
485 0 : CPABORT(message)
486 : ELSE
487 1 : CPWARN(message)
488 : END IF
489 : END IF
490 : END IF
491 :
492 : ELSE
493 0 : CALL cp_warn(__LOCATION__, "Debug of dipole moments only for Quickstep code available")
494 : END IF
495 : END IF
496 :
497 736 : IF (debug_polar) THEN
498 58 : IF (force_env%in_use == use_qs_force) THEN
499 58 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
500 58 : poldir = [0.0_dp, 0.0_dp, 1.0_dp]
501 58 : amplitude = 0.0_dp
502 58 : CALL set_efield(dft_control, amplitude, poldir)
503 58 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
504 58 : CALL get_qs_env(force_env%qs_env, results=results)
505 58 : description = "[POLAR]"
506 58 : IF (test_for_result(results, description=description)) THEN
507 58 : CALL get_results(results, description=description, values=pvals)
508 58 : polar_analytic(1:3, 1:3) = RESHAPE(pvals(1:9), [3, 3])
509 : ELSE
510 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities needs PROPERTIES/LINRES/POLAR section enabled")
511 0 : CPABORT("DEBUG")
512 : END IF
513 58 : description = "[DIPOLE]"
514 58 : IF (.NOT. test_for_result(results, description=description)) THEN
515 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities need DFT/PRINT/MOMENTS section enabled")
516 0 : CPABORT("DEBUG")
517 : END IF
518 58 : amplitude = de
519 232 : DO k = 1, 3
520 174 : poldir = 0.0_dp
521 174 : poldir(k) = 1.0_dp
522 522 : DO j = 1, 2
523 1392 : poldir = -1.0_dp*poldir
524 348 : CALL set_efield(dft_control, amplitude, poldir)
525 348 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., linres=.TRUE.)
526 522 : CALL get_results(results, description=description, values=dipn(1:3, j))
527 : END DO
528 754 : polar_numeric(1:3, k) = 0.5_dp*(dipn(1:3, 2) - dipn(1:3, 1))/de
529 : END DO
530 58 : IF (iw > 0) THEN
531 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
532 29 : "DEBUG|========================= POLARIZABILITY ================================", &
533 58 : "DEBUG| Coordinates P(numerical) P(analytical) Difference Error [%]"
534 116 : DO k = 1, 3
535 377 : DO j = 1, 3
536 261 : dd = polar_analytic(k, j) - polar_numeric(k, j)
537 348 : IF (ABS(polar_analytic(k, j)) > eps_no_error_check) THEN
538 179 : derr = 100._dp*dd/polar_analytic(k, j)
539 : WRITE (UNIT=iw, FMT="(T2,A,T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
540 179 : "DEBUG|", ACHAR(119 + k), ACHAR(119 + j), polar_numeric(k, j), polar_analytic(k, j), dd, derr
541 : ELSE
542 : WRITE (UNIT=iw, FMT="(T2,A,T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
543 82 : "DEBUG|", ACHAR(119 + k), ACHAR(119 + j), polar_numeric(k, j), polar_analytic(k, j), dd
544 : END IF
545 : END DO
546 : END DO
547 : WRITE (UNIT=iw, FMT="((T2,A))") &
548 29 : "DEBUG|========================================================================="
549 377 : WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") ' POLAR : CheckSum =', SUM(polar_analytic)
550 : END IF
551 : ELSE
552 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities only for Quickstep code available")
553 : END IF
554 : END IF
555 :
556 736 : CALL cp_print_key_finished_output(iw, logger, root_section, "DEBUG%PROGRAM_RUN_INFO")
557 :
558 1472 : END SUBROUTINE cp2k_debug_energy_and_forces
559 :
560 : ! **************************************************************************************************
561 : !> \brief ...
562 : !> \param dft_control ...
563 : !> \param amplitude ...
564 : !> \param poldir ...
565 : ! **************************************************************************************************
566 1012 : SUBROUTINE set_efield(dft_control, amplitude, poldir)
567 : TYPE(dft_control_type), POINTER :: dft_control
568 : REAL(KIND=dp), INTENT(IN) :: amplitude
569 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: poldir
570 :
571 1012 : IF (dft_control%apply_efield) THEN
572 928 : dft_control%efield_fields(1)%efield%strength = amplitude
573 3712 : dft_control%efield_fields(1)%efield%polarisation(1:3) = poldir(1:3)
574 84 : ELSEIF (dft_control%apply_period_efield) THEN
575 84 : dft_control%period_efield%strength = amplitude
576 336 : dft_control%period_efield%polarisation(1:3) = poldir(1:3)
577 : ELSE
578 0 : CPABORT("No EFIELD definition available")
579 : END IF
580 :
581 1012 : END SUBROUTINE set_efield
582 :
583 : END MODULE cp2k_debug
|
