Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 424 : 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 424 : 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 424 : 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 424 : 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 424 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
109 424 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
110 : TYPE(section_vals_type), POINTER :: root_section, subsys_section
111 :
112 424 : NULLIFY (analyt_forces, numer_forces, subsys, particles)
113 :
114 424 : root_section => force_env%root_section
115 :
116 424 : 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 424 : "SUBSYS")
119 :
120 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_STRESS_TENSOR", &
121 424 : l_val=debug_stress_tensor)
122 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_FORCES", &
123 424 : l_val=debug_forces)
124 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_DIPOLE", &
125 424 : l_val=debug_dipole)
126 : CALL section_vals_val_get(root_section, "DEBUG%DEBUG_POLARIZABILITY", &
127 424 : l_val=debug_polar)
128 : CALL section_vals_val_get(root_section, "DEBUG%DX", &
129 424 : r_val=dx)
130 : CALL section_vals_val_get(root_section, "DEBUG%DE", &
131 424 : r_val=de)
132 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_DIPOLE_DIRS", &
133 424 : c_val=cval1)
134 424 : dx = ABS(dx)
135 : CALL section_vals_val_get(root_section, "DEBUG%MAX_RELATIVE_ERROR", &
136 424 : r_val=maxerr)
137 : CALL section_vals_val_get(root_section, "DEBUG%EPS_NO_ERROR_CHECK", &
138 424 : r_val=eps_no_error_check)
139 424 : 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 424 : l_val=stop_on_mismatch)
142 :
143 : ! set active DOF
144 424 : CALL uppercase(cval1)
145 424 : do_dof_dipole(1) = (INDEX(cval1, "X") /= 0)
146 424 : do_dof_dipole(2) = (INDEX(cval1, "Y") /= 0)
147 424 : do_dof_dipole(3) = (INDEX(cval1, "Z") /= 0)
148 424 : NULLIFY (cval2)
149 424 : IF (debug_forces) THEN
150 262 : np = subsys%particles%n_els
151 786 : ALLOCATE (do_dof_atom_coor(3, np))
152 262 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", n_rep_val=nrep)
153 262 : IF (nrep > 0) THEN
154 674 : do_dof_atom_coor = .FALSE.
155 102 : DO irep = 1, nrep
156 : CALL section_vals_val_get(root_section, "DEBUG%CHECK_ATOM_FORCE", i_rep_val=irep, &
157 52 : c_vals=cval2)
158 52 : READ (cval2(1), FMT="(I10)") k
159 52 : CALL uppercase(cval2(2))
160 52 : do_dof_atom_coor(1, k) = (INDEX(cval2(2), "X") /= 0)
161 52 : do_dof_atom_coor(2, k) = (INDEX(cval2(2), "Y") /= 0)
162 102 : do_dof_atom_coor(3, k) = (INDEX(cval2(2), "Z") /= 0)
163 : END DO
164 : ELSE
165 3996 : do_dof_atom_coor = .TRUE.
166 : END IF
167 : END IF
168 :
169 424 : logger => cp_get_default_logger()
170 : iw = cp_print_key_unit_nr(logger, root_section, "DEBUG%PROGRAM_RUN_INFO", &
171 424 : extension=".log")
172 424 : IF (debug_stress_tensor) THEN
173 : ! To debug stress tensor the stress tensor calculation must be
174 : ! first enabled..
175 : CALL section_vals_val_get(force_env%force_env_section, "STRESS_TENSOR", &
176 166 : i_val=stress_tensor)
177 166 : skip = .FALSE.
178 0 : SELECT CASE (stress_tensor)
179 : CASE (do_stress_analytical, do_stress_diagonal_anal)
180 : ! OK
181 : CASE (do_stress_numerical, do_stress_diagonal_numer)
182 : ! Nothing to check
183 : CALL cp_warn(__LOCATION__, "Numerical stress tensor was requested in "// &
184 0 : "the FORCE_EVAL section. Nothing to debug!")
185 114 : skip = .TRUE.
186 : CASE DEFAULT
187 : CALL cp_warn(__LOCATION__, "Stress tensor calculation was not enabled in "// &
188 114 : "the FORCE_EVAL section. Nothing to debug!")
189 166 : skip = .TRUE.
190 : END SELECT
191 :
192 : IF (.NOT. skip) THEN
193 :
194 : BLOCK
195 : TYPE(virial_type) :: virial_analytical, virial_numerical
196 : TYPE(virial_type), POINTER :: virial
197 :
198 : ! Compute the analytical stress tensor
199 52 : CALL cp_subsys_get(subsys, virial=virial)
200 52 : CALL virial_set(virial, pv_numer=.FALSE.)
201 : CALL force_env_calc_energy_force(force_env, &
202 : calc_force=.TRUE., &
203 52 : calc_stress_tensor=.TRUE.)
204 :
205 : ! Retrieve the analytical virial
206 52 : virial_analytical = virial
207 :
208 : ! Debug stress tensor (numerical vs analytical)
209 52 : CALL virial_set(virial, pv_numer=.TRUE.)
210 52 : CALL force_env_calc_num_pressure(force_env, dx=dx)
211 :
212 : ! Retrieve the numerical virial
213 52 : CALL cp_subsys_get(subsys, virial=virial)
214 52 : virial_numerical = virial
215 :
216 : ! Print results
217 52 : IF (iw > 0) THEN
218 : WRITE (UNIT=iw, FMT="((T2,A))") &
219 26 : "DEBUG| Numerical pv_virial [a.u.]"
220 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
221 104 : ("DEBUG|", virial_numerical%pv_virial(i, 1:3), i=1, 3)
222 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
223 26 : "DEBUG| Analytical pv_virial [a.u.]"
224 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
225 104 : ("DEBUG|", virial_analytical%pv_virial(i, 1:3), i=1, 3)
226 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
227 26 : "DEBUG| Difference pv_virial [a.u.]"
228 : WRITE (UNIT=iw, FMT="((T2,A,T21,3(1X,F19.12)))") &
229 338 : ("DEBUG|", virial_numerical%pv_virial(i, 1:3) - virial_analytical%pv_virial(i, 1:3), i=1, 3)
230 : WRITE (UNIT=iw, FMT="(T2,A,T61,F20.12)") &
231 26 : "DEBUG| Sum of differences", &
232 364 : SUM(ABS(virial_numerical%pv_virial(:, :) - virial_analytical%pv_virial(:, :)))
233 : END IF
234 :
235 : ! Check and abort on failure
236 52 : check_failed = .FALSE.
237 52 : IF (iw > 0) THEN
238 : WRITE (UNIT=iw, FMT="(/T2,A)") &
239 26 : "DEBUG| Relative error pv_virial"
240 : WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.8)") &
241 26 : "DEBUG| Threshold value for error check [a.u.]", eps_no_error_check
242 : END IF
243 208 : DO i = 1, 3
244 156 : err(:) = 0.0_dp
245 624 : DO k = 1, 3
246 624 : IF (ABS(virial_analytical%pv_virial(i, k)) >= eps_no_error_check) THEN
247 : err(k) = 100.0_dp*(virial_numerical%pv_virial(i, k) - virial_analytical%pv_virial(i, k))/ &
248 424 : virial_analytical%pv_virial(i, k)
249 424 : WRITE (UNIT=line(20*(k - 1) + 1:20*k), FMT="(1X,F17.2,A2)") err(k), " %"
250 : ELSE
251 44 : WRITE (UNIT=line(20*(k - 1) + 1:20*k), FMT="(17X,A3)") "- %"
252 : END IF
253 : END DO
254 156 : IF (iw > 0) THEN
255 : WRITE (UNIT=iw, FMT="(T2,A,T21,A60)") &
256 78 : "DEBUG|", line
257 : END IF
258 676 : IF (ANY(ABS(err(1:3)) > maxerr)) check_failed = .TRUE.
259 : END DO
260 52 : IF (iw > 0) THEN
261 : WRITE (UNIT=iw, FMT="(T2,A,T61,F18.2,A2)") &
262 26 : "DEBUG| Maximum accepted error", maxerr, " %"
263 : END IF
264 23816 : IF (check_failed) THEN
265 : message = "A mismatch between the analytical and the numerical "// &
266 : "stress tensor has been detected. Check the implementation "// &
267 0 : "of the stress tensor"
268 0 : IF (stop_on_mismatch) THEN
269 0 : CPABORT(TRIM(message))
270 : ELSE
271 0 : CPWARN(TRIM(message))
272 : END IF
273 : END IF
274 : END BLOCK
275 : END IF
276 : END IF
277 :
278 424 : IF (debug_forces) THEN
279 : ! Debug forces (numerical vs analytical)
280 262 : particles => subsys%particles%els
281 428 : SELECT CASE (force_env%in_use)
282 : CASE (use_qs_force)
283 166 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
284 166 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
285 : CASE DEFAULT
286 262 : CALL write_fist_particle_coordinates(particles, subsys_section)
287 : END SELECT
288 : ! First evaluate energy and forces
289 : CALL force_env_calc_energy_force(force_env, &
290 : calc_force=.TRUE., &
291 262 : calc_stress_tensor=.FALSE.)
292 : ! Copy forces in array and start the numerical calculation
293 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
294 262 : np = subsys%particles%n_els
295 786 : ALLOCATE (analyt_forces(np, 3))
296 1364 : DO ip = 1, np
297 4670 : analyt_forces(ip, 1:3) = particles(ip)%f(1:3)
298 : END DO
299 : ! Loop on atoms and coordinates
300 : IF (ASSOCIATED(numer_forces)) DEALLOCATE (numer_forces)
301 786 : ALLOCATE (numer_forces(subsys%particles%n_els, 3))
302 1364 : Atom: DO ip = 1, np
303 4408 : Coord: DO k = 1, 3
304 4408 : IF (do_dof_atom_coor(k, ip)) THEN
305 2890 : numer_energy = 0.0_dp
306 2890 : std_value = particles(ip)%r(k)
307 8670 : DO j = 1, 2
308 5780 : particles(ip)%r(k) = std_value - (-1.0_dp)**j*dx
309 8464 : SELECT CASE (force_env%in_use)
310 : CASE (use_qs_force)
311 2684 : CALL get_qs_env(force_env%qs_env, qs_kind_set=qs_kind_set)
312 2684 : CALL write_qs_particle_coordinates(particles, qs_kind_set, subsys_section, "DEBUG")
313 : CASE DEFAULT
314 5780 : CALL write_fist_particle_coordinates(particles, subsys_section)
315 : END SELECT
316 : ! Compute energy
317 : CALL force_env_calc_energy_force(force_env, &
318 : calc_force=.FALSE., &
319 : calc_stress_tensor=.FALSE., &
320 5780 : consistent_energies=.TRUE.)
321 8670 : CALL force_env_get(force_env, potential_energy=numer_energy(j))
322 : END DO
323 2890 : particles(ip)%r(k) = std_value
324 2890 : numer_forces(ip, k) = -0.5_dp*(numer_energy(1) - numer_energy(2))/dx
325 2890 : IF (iw > 0) THEN
326 : WRITE (UNIT=iw, FMT="(/,T2,A,T17,A,F7.4,A,T34,A,F7.4,A,T52,A,T68,A)") &
327 1445 : "DEBUG| Atom", "E("//ACHAR(119 + k)//" +", dx, ")", &
328 1445 : "E("//ACHAR(119 + k)//" -", dx, ")", &
329 2890 : "f(numerical)", "f(analytical)"
330 : WRITE (UNIT=iw, FMT="(T2,A,I5,4(1X,F16.8))") &
331 1445 : "DEBUG|", ip, numer_energy(1:2), numer_forces(ip, k), analyt_forces(ip, k)
332 : END IF
333 : ELSE
334 416 : numer_forces(ip, k) = 0.0_dp
335 : END IF
336 : END DO Coord
337 : ! Check analytical forces using the numerical forces as reference
338 4408 : my_maxerr = maxerr
339 1102 : err(1:3) = 0.0_dp
340 4408 : DO k = 1, 3
341 4408 : IF (do_dof_atom_coor(k, ip)) THEN
342 : ! Calculate percentage but ignore very small force values
343 2890 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
344 2218 : err(k) = 100.0_dp*(numer_forces(ip, k) - analyt_forces(ip, k))/analyt_forces(ip, k)
345 : END IF
346 : ! Increase error tolerance for small force values
347 2890 : IF (ABS(analyt_forces(ip, k)) <= 0.0001_dp) my_maxerr(k) = 5.0_dp*my_maxerr(k)
348 : ELSE
349 416 : err(k) = 0.0_dp
350 : END IF
351 : END DO
352 1102 : IF (iw > 0) THEN
353 755 : IF (ANY(do_dof_atom_coor(1:3, ip))) THEN
354 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
355 499 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
356 : END IF
357 2204 : DO k = 1, 3
358 2204 : IF (do_dof_atom_coor(k, ip)) THEN
359 1445 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
360 1445 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
361 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T71,F10.2)") &
362 1109 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, err(k)
363 : ELSE
364 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T80,A1)") &
365 336 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, "-"
366 : END IF
367 : END IF
368 : END DO
369 : END IF
370 4670 : IF (ANY(ABS(err(1:3)) > my_maxerr(1:3))) THEN
371 : message = "A mismatch between analytical and numerical forces "// &
372 : "has been detected. Check the implementation of the "// &
373 0 : "analytical force calculation"
374 0 : IF (stop_on_mismatch) THEN
375 0 : CPABORT(message)
376 : ELSE
377 0 : CPWARN(message)
378 : END IF
379 : END IF
380 : END DO Atom
381 : ! Print summary
382 262 : IF (iw > 0) THEN
383 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
384 131 : "DEBUG|======================== BEGIN OF SUMMARY ===============================", &
385 262 : "DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
386 131 : sum_of_differences = 0.0_dp
387 131 : errmax = 0.0_dp
388 682 : DO ip = 1, np
389 551 : err(1:3) = 0.0_dp
390 2335 : DO k = 1, 3
391 2204 : IF (do_dof_atom_coor(k, ip)) THEN
392 1445 : difference = analyt_forces(ip, k) - numer_forces(ip, k)
393 1445 : IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
394 1109 : err(k) = 100_dp*(numer_forces(ip, k) - analyt_forces(ip, k))/analyt_forces(ip, k)
395 1109 : errmax = MAX(errmax, ABS(err(k)))
396 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T71,F10.2)") &
397 1109 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, err(k)
398 : ELSE
399 : WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T80,A1)") &
400 336 : "DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, "-"
401 : END IF
402 1445 : sum_of_differences = sum_of_differences + ABS(difference)
403 : END IF
404 : END DO
405 : END DO
406 : WRITE (UNIT=iw, FMT="(T2,A,T57,F12.8,T71,F10.2)") &
407 131 : "DEBUG| Sum of differences:", sum_of_differences, errmax
408 : WRITE (UNIT=iw, FMT="(T2,A)") &
409 131 : "DEBUG|======================== END OF SUMMARY ================================="
410 : END IF
411 : ! Release work storage
412 262 : IF (ASSOCIATED(analyt_forces)) DEALLOCATE (analyt_forces)
413 262 : IF (ASSOCIATED(numer_forces)) DEALLOCATE (numer_forces)
414 262 : DEALLOCATE (do_dof_atom_coor)
415 : END IF
416 :
417 424 : IF (debug_dipole) THEN
418 78 : IF (force_env%in_use == use_qs_force) THEN
419 78 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
420 78 : poldir = (/0.0_dp, 0.0_dp, 1.0_dp/)
421 78 : amplitude = 0.0_dp
422 78 : CALL set_efield(dft_control, amplitude, poldir)
423 78 : CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
424 78 : CALL get_qs_env(force_env%qs_env, results=results)
425 78 : description = "[DIPOLE]"
426 78 : IF (test_for_result(results, description=description)) THEN
427 78 : CALL get_results(results, description=description, values=dipole_moment)
428 : ELSE
429 0 : CALL cp_warn(__LOCATION__, "Debug of dipole moments needs DFT/PRINT/MOMENTS section enabled")
430 0 : CPABORT("DEBUG")
431 : END IF
432 78 : amplitude = de
433 312 : DO k = 1, 3
434 312 : IF (do_dof_dipole(k)) THEN
435 166 : poldir = 0.0_dp
436 166 : poldir(k) = 1.0_dp
437 498 : DO j = 1, 2
438 1328 : poldir = -1.0_dp*poldir
439 332 : CALL set_efield(dft_control, amplitude, poldir)
440 332 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
441 498 : CALL force_env_get(force_env, potential_energy=numer_energy(j))
442 : END DO
443 166 : dipole_numer(k) = 0.5_dp*(numer_energy(1) - numer_energy(2))/de
444 : ELSE
445 68 : dipole_numer(k) = 0.0_dp
446 : END IF
447 : END DO
448 78 : IF (iw > 0) THEN
449 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
450 39 : "DEBUG|========================= DIPOLE MOMENTS ================================", &
451 78 : "DEBUG| Coordinate D(numerical) D(analytical) Difference Error [%]"
452 39 : err(1:3) = 0.0_dp
453 156 : DO k = 1, 3
454 156 : IF (do_dof_dipole(k)) THEN
455 83 : dd = dipole_moment(k) - dipole_numer(k)
456 83 : IF (ABS(dipole_moment(k)) > eps_no_error_check) THEN
457 40 : derr = 100._dp*dd/dipole_moment(k)
458 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
459 40 : "DEBUG|", ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd, derr
460 : ELSE
461 43 : derr = 0.0_dp
462 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
463 43 : "DEBUG|", ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd
464 : END IF
465 83 : err(k) = derr
466 : ELSE
467 : WRITE (UNIT=iw, FMT="(T2,A,T13,A1,T21,A16,T38,F16.8)") &
468 34 : "DEBUG|", ACHAR(119 + k), " skipped", dipole_moment(k)
469 : END IF
470 : END DO
471 : WRITE (UNIT=iw, FMT="((T2,A))") &
472 39 : "DEBUG|========================================================================="
473 156 : WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") 'DIPOLE : CheckSum =', SUM(dipole_moment)
474 153 : IF (ANY(ABS(err(1:3)) > maxerr)) THEN
475 : message = "A mismatch between analytical and numerical dipoles "// &
476 : "has been detected. Check the implementation of the "// &
477 1 : "analytical dipole calculation"
478 1 : IF (stop_on_mismatch) THEN
479 0 : CPABORT(message)
480 : ELSE
481 1 : CPWARN(message)
482 : END IF
483 : END IF
484 : END IF
485 :
486 : ELSE
487 0 : CALL cp_warn(__LOCATION__, "Debug of dipole moments only for Quickstep code available")
488 : END IF
489 : END IF
490 :
491 424 : IF (debug_polar) THEN
492 52 : IF (force_env%in_use == use_qs_force) THEN
493 52 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
494 52 : poldir = (/0.0_dp, 0.0_dp, 1.0_dp/)
495 52 : amplitude = 0.0_dp
496 52 : CALL set_efield(dft_control, amplitude, poldir)
497 52 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE.)
498 52 : CALL get_qs_env(force_env%qs_env, results=results)
499 52 : description = "[POLAR]"
500 52 : IF (test_for_result(results, description=description)) THEN
501 52 : CALL get_results(results, description=description, values=pvals)
502 52 : polar_analytic(1:3, 1:3) = RESHAPE(pvals(1:9), (/3, 3/))
503 : ELSE
504 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities needs PROPERTIES/LINRES/POLAR section enabled")
505 0 : CPABORT("DEBUG")
506 : END IF
507 52 : description = "[DIPOLE]"
508 52 : IF (.NOT. test_for_result(results, description=description)) THEN
509 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities need DFT/PRINT/MOMENTS section enabled")
510 0 : CPABORT("DEBUG")
511 : END IF
512 52 : amplitude = de
513 208 : DO k = 1, 3
514 156 : poldir = 0.0_dp
515 156 : poldir(k) = 1.0_dp
516 468 : DO j = 1, 2
517 1248 : poldir = -1.0_dp*poldir
518 312 : CALL set_efield(dft_control, amplitude, poldir)
519 312 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., linres=.TRUE.)
520 468 : CALL get_results(results, description=description, values=dipn(1:3, j))
521 : END DO
522 676 : polar_numeric(1:3, k) = 0.5_dp*(dipn(1:3, 2) - dipn(1:3, 1))/de
523 : END DO
524 52 : IF (iw > 0) THEN
525 : WRITE (UNIT=iw, FMT="(/,(T2,A))") &
526 26 : "DEBUG|========================= POLARIZABILITY ================================", &
527 52 : "DEBUG| Coordinates P(numerical) P(analytical) Difference Error [%]"
528 104 : DO k = 1, 3
529 338 : DO j = 1, 3
530 234 : dd = polar_analytic(k, j) - polar_numeric(k, j)
531 312 : IF (ABS(polar_analytic(k, j)) > eps_no_error_check) THEN
532 170 : derr = 100._dp*dd/polar_analytic(k, j)
533 : WRITE (UNIT=iw, FMT="(T2,A,T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
534 170 : "DEBUG|", ACHAR(119 + k), ACHAR(119 + j), polar_numeric(k, j), polar_analytic(k, j), dd, derr
535 : ELSE
536 : WRITE (UNIT=iw, FMT="(T2,A,T12,A1,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
537 64 : "DEBUG|", ACHAR(119 + k), ACHAR(119 + j), polar_numeric(k, j), polar_analytic(k, j), dd
538 : END IF
539 : END DO
540 : END DO
541 : WRITE (UNIT=iw, FMT="((T2,A))") &
542 26 : "DEBUG|========================================================================="
543 338 : WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") ' POLAR : CheckSum =', SUM(polar_analytic)
544 : END IF
545 : ELSE
546 0 : CALL cp_warn(__LOCATION__, "Debug of polarizabilities only for Quickstep code available")
547 : END IF
548 : END IF
549 :
550 424 : CALL cp_print_key_finished_output(iw, logger, root_section, "DEBUG%PROGRAM_RUN_INFO")
551 :
552 848 : END SUBROUTINE cp2k_debug_energy_and_forces
553 :
554 : ! **************************************************************************************************
555 : !> \brief ...
556 : !> \param dft_control ...
557 : !> \param amplitude ...
558 : !> \param poldir ...
559 : ! **************************************************************************************************
560 774 : SUBROUTINE set_efield(dft_control, amplitude, poldir)
561 : TYPE(dft_control_type), POINTER :: dft_control
562 : REAL(KIND=dp), INTENT(IN) :: amplitude
563 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: poldir
564 :
565 774 : IF (dft_control%apply_efield) THEN
566 746 : dft_control%efield_fields(1)%efield%strength = amplitude
567 2984 : dft_control%efield_fields(1)%efield%polarisation(1:3) = poldir(1:3)
568 28 : ELSEIF (dft_control%apply_period_efield) THEN
569 28 : dft_control%period_efield%strength = amplitude
570 112 : dft_control%period_efield%polarisation(1:3) = poldir(1:3)
571 : ELSE
572 0 : CPABORT("No EFIELD definition available")
573 : END IF
574 :
575 774 : END SUBROUTINE set_efield
576 :
577 : END MODULE cp2k_debug
|
