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 Methods to perform free energy and free energy derivatives calculations
10 : !> \author Teodoro Laino (01.2007) [tlaino]
11 : ! **************************************************************************************************
12 : MODULE free_energy_methods
13 : USE colvar_methods, ONLY: colvar_eval_glob_f
14 : USE cp_log_handling, ONLY: cp_get_default_logger,&
15 : cp_logger_get_default_io_unit,&
16 : cp_logger_type,&
17 : cp_to_string
18 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
19 : cp_print_key_unit_nr
20 : USE cp_subsys_types, ONLY: cp_subsys_type
21 : USE force_env_types, ONLY: force_env_get,&
22 : force_env_type
23 : USE fparser, ONLY: evalf,&
24 : evalfd,&
25 : finalizef,&
26 : initf,&
27 : parsef
28 : USE free_energy_types, ONLY: free_energy_type,&
29 : ui_var_type
30 : USE input_constants, ONLY: do_fe_ac,&
31 : do_fe_ui
32 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
33 : section_vals_type,&
34 : section_vals_val_get
35 : USE kinds, ONLY: default_path_length,&
36 : default_string_length,&
37 : dp
38 : USE mathlib, ONLY: diamat_all
39 : USE md_environment_types, ONLY: get_md_env,&
40 : md_environment_type
41 : USE memory_utilities, ONLY: reallocate
42 : USE simpar_types, ONLY: simpar_type
43 : USE statistical_methods, ONLY: k_test,&
44 : min_sample_size,&
45 : sw_test,&
46 : vn_test
47 : USE string_utilities, ONLY: compress
48 : #include "../base/base_uses.f90"
49 :
50 : IMPLICIT NONE
51 :
52 : PRIVATE
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'free_energy_methods'
54 : PUBLIC :: free_energy_evaluate
55 :
56 : CONTAINS
57 :
58 : ! **************************************************************************************************
59 : !> \brief Main driver for free energy calculations
60 : !> In this routine we handle specifically biased MD.
61 : !> \param md_env ...
62 : !> \param converged ...
63 : !> \param fe_section ...
64 : !> \par History
65 : !> Teodoro Laino (01.2007) [tlaino]
66 : ! **************************************************************************************************
67 79038 : SUBROUTINE free_energy_evaluate(md_env, converged, fe_section)
68 : TYPE(md_environment_type), POINTER :: md_env
69 : LOGICAL, INTENT(OUT) :: converged
70 : TYPE(section_vals_type), POINTER :: fe_section
71 :
72 : CHARACTER(LEN=*), PARAMETER :: routineN = 'free_energy_evaluate'
73 :
74 : INTEGER :: handle, ic, icolvar, nforce_eval, &
75 : output_unit, stat_sign_points
76 : INTEGER, POINTER :: istep
77 : REAL(KIND=dp) :: beta
78 : TYPE(cp_logger_type), POINTER :: logger
79 : TYPE(cp_subsys_type), POINTER :: subsys
80 : TYPE(force_env_type), POINTER :: force_env
81 : TYPE(free_energy_type), POINTER :: fe_env
82 : TYPE(simpar_type), POINTER :: simpar
83 : TYPE(ui_var_type), POINTER :: cv
84 :
85 39519 : NULLIFY (force_env, istep, subsys, cv, simpar)
86 39519 : logger => cp_get_default_logger()
87 39519 : CALL timeset(routineN, handle)
88 39519 : converged = .FALSE.
89 : CALL get_md_env(md_env, force_env=force_env, fe_env=fe_env, simpar=simpar, &
90 39519 : itimes=istep)
91 : ! Metadynamics is also a free energy calculation but is handled in a different
92 : ! module.
93 39519 : IF (.NOT. ASSOCIATED(force_env%meta_env) .AND. ASSOCIATED(fe_env)) THEN
94 210 : SELECT CASE (fe_env%type)
95 : CASE (do_fe_ui)
96 : ! Umbrella Integration..
97 20 : CALL force_env_get(force_env, subsys=subsys)
98 20 : fe_env%nr_points = fe_env%nr_points + 1
99 20 : output_unit = cp_logger_get_default_io_unit(logger)
100 40 : DO ic = 1, fe_env%ncolvar
101 20 : cv => fe_env%uivar(ic)
102 20 : icolvar = cv%icolvar
103 20 : CALL colvar_eval_glob_f(icolvar, force_env)
104 20 : CALL reallocate(cv%ss, 1, fe_env%nr_points)
105 20 : cv%ss(fe_env%nr_points) = subsys%colvar_p(icolvar)%colvar%ss
106 40 : IF (output_unit > 0) THEN
107 10 : WRITE (output_unit, *) "COLVAR::", cv%ss(fe_env%nr_points)
108 : END IF
109 : END DO
110 20 : stat_sign_points = fe_env%nr_points - fe_env%nr_rejected
111 20 : IF (output_unit > 0) THEN
112 10 : WRITE (output_unit, *) fe_env%nr_points, stat_sign_points
113 : END IF
114 : ! Start statistical analysis when enough CG data points have been collected
115 20 : IF ((fe_env%conv_par%cg_width*fe_env%conv_par%cg_points <= stat_sign_points) .AND. &
116 170 : (MOD(stat_sign_points, fe_env%conv_par%cg_width) == 0)) THEN
117 : output_unit = cp_print_key_unit_nr(logger, fe_section, "FREE_ENERGY_INFO", &
118 4 : extension=".FreeEnergyLog", log_filename=.FALSE.)
119 4 : CALL print_fe_prolog(output_unit)
120 : ! Trend test.. recomputes the number of statistically significant points..
121 4 : CALL ui_check_trend(fe_env, fe_env%conv_par%test_k, stat_sign_points, output_unit)
122 4 : stat_sign_points = fe_env%nr_points - fe_env%nr_rejected
123 : ! Normality and serial correlation tests..
124 4 : IF (fe_env%conv_par%cg_width*fe_env%conv_par%cg_points <= stat_sign_points .AND. &
125 : fe_env%conv_par%test_k) THEN
126 : ! Statistical tests
127 0 : CALL ui_check_convergence(fe_env, converged, stat_sign_points, output_unit)
128 : END IF
129 4 : CALL print_fe_epilog(output_unit)
130 4 : CALL cp_print_key_finished_output(output_unit, logger, fe_section, "FREE_ENERGY_INFO")
131 : END IF
132 : CASE (do_fe_ac)
133 : ! Alchemical Changes
134 170 : IF (.NOT. ASSOCIATED(force_env%mixed_env)) THEN
135 : CALL cp_abort(__LOCATION__, &
136 : "Free Energy calculations with the alchemical method require the "// &
137 : "definition of a mixed force env with the FORCE_EVAL/MIXED/GENERIC "// &
138 0 : "section, but none is found!")
139 : END IF
140 170 : beta = 1.0_dp/simpar%temp_ext
141 170 : nforce_eval = SIZE(force_env%sub_force_env)
142 : CALL dump_ac_info(force_env%mixed_env%val, force_env%mixed_env%par, &
143 : force_env%mixed_env%coupling_function, &
144 : force_env%mixed_env%dx, force_env%mixed_env%lerr, &
145 360 : fe_section, nforce_eval, fe_env%covmx, istep, beta)
146 : CASE DEFAULT
147 : ! Do Nothing
148 : END SELECT
149 : END IF
150 39519 : CALL timestop(handle)
151 :
152 39519 : END SUBROUTINE free_energy_evaluate
153 :
154 : ! **************************************************************************************************
155 : !> \brief Print prolog of free energy output section
156 : !> \param output_unit which unit to print to
157 : !> \par History
158 : !> Teodoro Laino (02.2007) [tlaino]
159 : ! **************************************************************************************************
160 4 : SUBROUTINE print_fe_prolog(output_unit)
161 : INTEGER, INTENT(IN) :: output_unit
162 :
163 4 : IF (output_unit > 0) THEN
164 2 : WRITE (output_unit, '(T2,79("*"))')
165 2 : WRITE (output_unit, '(T30,"FREE ENERGY CALCULATION",/)')
166 : END IF
167 4 : END SUBROUTINE print_fe_prolog
168 :
169 : ! **************************************************************************************************
170 : !> \brief Print epilog of free energy output section
171 : !> \param output_unit which unit to print to
172 : !> \par History
173 : !> Teodoro Laino (02.2007) [tlaino]
174 : ! **************************************************************************************************
175 4 : SUBROUTINE print_fe_epilog(output_unit)
176 : INTEGER, INTENT(IN) :: output_unit
177 :
178 4 : IF (output_unit > 0) THEN
179 2 : WRITE (output_unit, '(T2,79("*"),/)')
180 : END IF
181 4 : END SUBROUTINE print_fe_epilog
182 :
183 : ! **************************************************************************************************
184 : !> \brief Test for trend in coarse grained data set
185 : !> \param fe_env ...
186 : !> \param trend_free ...
187 : !> \param nr_points ...
188 : !> \param output_unit which unit to print to
189 : !> \par History
190 : !> Teodoro Laino (01.2007) [tlaino]
191 : ! **************************************************************************************************
192 4 : SUBROUTINE ui_check_trend(fe_env, trend_free, nr_points, output_unit)
193 : TYPE(free_energy_type), POINTER :: fe_env
194 : LOGICAL, INTENT(OUT) :: trend_free
195 : INTEGER, INTENT(IN) :: nr_points, output_unit
196 :
197 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ui_check_trend'
198 :
199 : INTEGER :: handle, i, ii, j, k, my_reject, ncolvar, &
200 : ng_points, rejected_points
201 : LOGICAL :: test_avg, test_std
202 : REAL(KIND=dp) :: prob, tau, z
203 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: wrk
204 :
205 4 : CALL timeset(routineN, handle)
206 4 : trend_free = .FALSE.
207 4 : test_avg = .TRUE.
208 4 : test_std = .TRUE.
209 4 : ncolvar = fe_env%ncolvar
210 : ! Number of coarse grained points
211 4 : IF (output_unit > 0) THEN
212 2 : WRITE (output_unit, *) nr_points, fe_env%conv_par%cg_width
213 : END IF
214 4 : ng_points = nr_points/fe_env%conv_par%cg_width
215 4 : my_reject = 0
216 : ! Allocate storage
217 4 : CALL create_tmp_data(fe_env, wrk, ng_points, ncolvar)
218 : ! Compute the Coarse Grained data set using a reverse cumulative strategy
219 4 : CALL create_csg_data(fe_env, ng_points, output_unit)
220 : ! Test on coarse grained average
221 8 : DO j = 1, ncolvar
222 : ii = 1
223 14 : DO i = ng_points, 1, -1
224 10 : wrk(ii) = fe_env%cg_data(i)%avg(j)
225 14 : ii = ii + 1
226 : END DO
227 4 : DO i = my_reject + 1, ng_points
228 4 : IF ((ng_points - my_reject) < min_sample_size) THEN
229 4 : my_reject = MAX(0, my_reject - 1)
230 4 : test_avg = .FALSE.
231 4 : EXIT
232 : END IF
233 0 : CALL k_test(wrk, my_reject + 1, ng_points, tau, z, prob)
234 0 : PRINT *, prob, fe_env%conv_par%k_conf_lm
235 0 : IF (prob < fe_env%conv_par%k_conf_lm) EXIT
236 0 : my_reject = my_reject + 1
237 : END DO
238 8 : my_reject = MIN(ng_points, my_reject)
239 : END DO
240 4 : rejected_points = my_reject*fe_env%conv_par%cg_width
241 : ! Print some info
242 4 : IF (output_unit > 0) THEN
243 2 : WRITE (output_unit, *) "Kendall trend test (Average)", test_avg, &
244 4 : "number of points rejected:", rejected_points + fe_env%nr_rejected
245 2 : WRITE (output_unit, *) "Reject Nr.", my_reject, " coarse grained points testing average"
246 : END IF
247 : ! Test on coarse grained covariance matrix
248 8 : DO j = 1, ncolvar
249 12 : DO k = j, ncolvar
250 : ii = 1
251 14 : DO i = ng_points, 1, -1
252 10 : wrk(ii) = fe_env%cg_data(i)%var(j, k)
253 14 : ii = ii + 1
254 : END DO
255 4 : DO i = my_reject + 1, ng_points
256 4 : IF ((ng_points - my_reject) < min_sample_size) THEN
257 4 : my_reject = MAX(0, my_reject - 1)
258 4 : test_std = .FALSE.
259 4 : EXIT
260 : END IF
261 0 : CALL k_test(wrk, my_reject + 1, ng_points, tau, z, prob)
262 0 : PRINT *, prob, fe_env%conv_par%k_conf_lm
263 0 : IF (prob < fe_env%conv_par%k_conf_lm) EXIT
264 0 : my_reject = my_reject + 1
265 : END DO
266 8 : my_reject = MIN(ng_points, my_reject)
267 : END DO
268 : END DO
269 4 : rejected_points = my_reject*fe_env%conv_par%cg_width
270 4 : fe_env%nr_rejected = fe_env%nr_rejected + rejected_points
271 4 : trend_free = test_avg .AND. test_std
272 : ! Print some info
273 4 : IF (output_unit > 0) THEN
274 2 : WRITE (output_unit, *) "Kendall trend test (Std. Dev.)", test_std, &
275 4 : "number of points rejected:", fe_env%nr_rejected
276 2 : WRITE (output_unit, *) "Reject Nr.", my_reject, " coarse grained points testing standard dev."
277 2 : WRITE (output_unit, *) "Kendall test passed:", trend_free
278 : END IF
279 : ! Release storage
280 4 : CALL destroy_tmp_data(fe_env, wrk, ng_points)
281 4 : CALL timestop(handle)
282 4 : END SUBROUTINE ui_check_trend
283 :
284 : ! **************************************************************************************************
285 : !> \brief Creates temporary data structures
286 : !> \param fe_env ...
287 : !> \param wrk ...
288 : !> \param ng_points ...
289 : !> \param ncolvar ...
290 : !> \par History
291 : !> Teodoro Laino (02.2007) [tlaino]
292 : ! **************************************************************************************************
293 4 : SUBROUTINE create_tmp_data(fe_env, wrk, ng_points, ncolvar)
294 : TYPE(free_energy_type), POINTER :: fe_env
295 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: wrk
296 : INTEGER, INTENT(IN) :: ng_points, ncolvar
297 :
298 : INTEGER :: i
299 :
300 22 : ALLOCATE (fe_env%cg_data(ng_points))
301 14 : DO i = 1, ng_points
302 30 : ALLOCATE (fe_env%cg_data(i)%avg(ncolvar))
303 44 : ALLOCATE (fe_env%cg_data(i)%var(ncolvar, ncolvar))
304 : END DO
305 4 : IF (PRESENT(wrk)) THEN
306 12 : ALLOCATE (wrk(ng_points))
307 : END IF
308 4 : END SUBROUTINE create_tmp_data
309 :
310 : ! **************************************************************************************************
311 : !> \brief Destroys temporary data structures
312 : !> \param fe_env ...
313 : !> \param wrk ...
314 : !> \param ng_points ...
315 : !> \par History
316 : !> Teodoro Laino (02.2007) [tlaino]
317 : ! **************************************************************************************************
318 4 : SUBROUTINE destroy_tmp_data(fe_env, wrk, ng_points)
319 : TYPE(free_energy_type), POINTER :: fe_env
320 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: wrk
321 : INTEGER, INTENT(IN) :: ng_points
322 :
323 : INTEGER :: i
324 :
325 14 : DO i = 1, ng_points
326 10 : DEALLOCATE (fe_env%cg_data(i)%avg)
327 14 : DEALLOCATE (fe_env%cg_data(i)%var)
328 : END DO
329 4 : DEALLOCATE (fe_env%cg_data)
330 4 : IF (PRESENT(wrk)) THEN
331 4 : DEALLOCATE (wrk)
332 : END IF
333 4 : END SUBROUTINE destroy_tmp_data
334 :
335 : ! **************************************************************************************************
336 : !> \brief Fills in temporary arrays with coarse grained data
337 : !> \param fe_env ...
338 : !> \param ng_points ...
339 : !> \param output_unit which unit to print to
340 : !> \par History
341 : !> Teodoro Laino (02.2007) [tlaino]
342 : ! **************************************************************************************************
343 4 : SUBROUTINE create_csg_data(fe_env, ng_points, output_unit)
344 : TYPE(free_energy_type), POINTER :: fe_env
345 : INTEGER, INTENT(IN) :: ng_points, output_unit
346 :
347 : INTEGER :: i, iend, istart
348 :
349 14 : DO i = 1, ng_points
350 10 : istart = fe_env%nr_points - (i)*fe_env%conv_par%cg_width + 1
351 10 : iend = fe_env%nr_points - (i - 1)*fe_env%conv_par%cg_width
352 10 : IF (output_unit > 0) THEN
353 5 : WRITE (output_unit, *) istart, iend
354 : END IF
355 14 : CALL eval_cov_matrix(fe_env, cg_index=i, istart=istart, iend=iend, output_unit=output_unit)
356 : END DO
357 :
358 4 : END SUBROUTINE create_csg_data
359 :
360 : ! **************************************************************************************************
361 : !> \brief Checks Normality of the distribution and Serial Correlation of
362 : !> coarse grained data
363 : !> \param fe_env ...
364 : !> \param test_passed ...
365 : !> \param nr_points ...
366 : !> \param output_unit which unit to print to
367 : !> \par History
368 : !> Teodoro Laino (02.2007) [tlaino]
369 : ! **************************************************************************************************
370 0 : SUBROUTINE ui_check_norm_sc(fe_env, test_passed, nr_points, output_unit)
371 : TYPE(free_energy_type), POINTER :: fe_env
372 : LOGICAL, INTENT(OUT) :: test_passed
373 : INTEGER, INTENT(IN) :: nr_points, output_unit
374 :
375 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ui_check_norm_sc'
376 :
377 : INTEGER :: handle, ng_points
378 :
379 0 : CALL timeset(routineN, handle)
380 0 : test_passed = .FALSE.
381 0 : DO WHILE (fe_env%conv_par%cg_width < fe_env%conv_par%max_cg_width)
382 0 : ng_points = nr_points/fe_env%conv_par%cg_width
383 0 : PRINT *, ng_points
384 0 : IF (ng_points < min_sample_size) EXIT
385 0 : CALL ui_check_norm_sc_low(fe_env, nr_points, output_unit)
386 0 : test_passed = fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw
387 0 : IF (test_passed) EXIT
388 0 : fe_env%conv_par%cg_width = fe_env%conv_par%cg_width + 1
389 0 : IF (output_unit > 0) THEN
390 0 : WRITE (output_unit, *) "New coarse grained width:", fe_env%conv_par%cg_width
391 : END IF
392 : END DO
393 0 : IF (fe_env%conv_par%cg_width == fe_env%conv_par%max_cg_width .AND. (.NOT. (test_passed))) THEN
394 0 : CALL ui_check_norm_sc_low(fe_env, nr_points, output_unit)
395 0 : test_passed = fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw
396 : END IF
397 0 : CALL timestop(handle)
398 0 : END SUBROUTINE ui_check_norm_sc
399 :
400 : ! **************************************************************************************************
401 : !> \brief Checks Normality of the distribution and Serial Correlation of
402 : !> coarse grained data - Low Level routine
403 : !> \param fe_env ...
404 : !> \param nr_points ...
405 : !> \param output_unit which unit to print to
406 : !> \par History
407 : !> Teodoro Laino (02.2007) [tlaino]
408 : ! **************************************************************************************************
409 0 : SUBROUTINE ui_check_norm_sc_low(fe_env, nr_points, output_unit)
410 : TYPE(free_energy_type), POINTER :: fe_env
411 : INTEGER, INTENT(IN) :: nr_points, output_unit
412 :
413 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ui_check_norm_sc_low'
414 :
415 : INTEGER :: handle, i, j, k, ncolvar, ng_points
416 : LOGICAL :: avg_test_passed, sdv_test_passed
417 : REAL(KIND=dp) :: prob, pw, r, u, w
418 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: wrk
419 :
420 0 : CALL timeset(routineN, handle)
421 0 : ncolvar = fe_env%ncolvar
422 : ! Compute the Coarse Grained data set using a reverse cumulative strategy
423 0 : fe_env%conv_par%test_sw = .FALSE.
424 0 : fe_env%conv_par%test_vn = .FALSE.
425 : ! Number of coarse grained points
426 : avg_test_passed = .TRUE.
427 : sdv_test_passed = .TRUE.
428 0 : ng_points = nr_points/fe_env%conv_par%cg_width
429 0 : CALL create_tmp_data(fe_env, wrk, ng_points, ncolvar)
430 0 : CALL create_csg_data(fe_env, ng_points, output_unit)
431 : ! Testing Averages
432 0 : DO j = 1, ncolvar
433 0 : DO i = 1, ng_points
434 0 : wrk(i) = fe_env%cg_data(i)%avg(j)
435 : END DO
436 : ! Test of Shapiro - Wilks for normality
437 : ! - Average
438 0 : CALL sw_test(wrk, ng_points, w, pw)
439 0 : PRINT *, 1.0_dp - pw, fe_env%conv_par%sw_conf_lm
440 0 : avg_test_passed = (1.0_dp - pw) <= fe_env%conv_par%sw_conf_lm
441 0 : fe_env%conv_par%test_sw = avg_test_passed
442 0 : IF (output_unit > 0) THEN
443 0 : WRITE (output_unit, *) "Shapiro-Wilks normality test (Avg)", avg_test_passed
444 : END IF
445 : ! Test of von Neumann for serial correlation
446 : ! - Average
447 0 : CALL vn_test(wrk, ng_points, r, u, prob)
448 0 : PRINT *, prob, fe_env%conv_par%vn_conf_lm
449 0 : avg_test_passed = prob <= fe_env%conv_par%vn_conf_lm
450 0 : fe_env%conv_par%test_vn = avg_test_passed
451 0 : IF (output_unit > 0) THEN
452 0 : WRITE (output_unit, *) "von Neumann serial correlation test (Avg)", avg_test_passed
453 : END IF
454 : END DO
455 : ! If tests on average are ok let's proceed with Standard Deviation
456 0 : IF (fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw) THEN
457 : ! Testing Standard Deviations
458 0 : DO j = 1, ncolvar
459 0 : DO k = j, ncolvar
460 0 : DO i = 1, ng_points
461 0 : wrk(i) = fe_env%cg_data(i)%var(j, k)
462 : END DO
463 : ! Test of Shapiro - Wilks for normality
464 : ! - Standard Deviation
465 0 : CALL sw_test(wrk, ng_points, w, pw)
466 0 : PRINT *, 1.0_dp - pw, fe_env%conv_par%sw_conf_lm
467 0 : sdv_test_passed = (1.0_dp - pw) <= fe_env%conv_par%sw_conf_lm
468 0 : fe_env%conv_par%test_sw = fe_env%conv_par%test_sw .AND. sdv_test_passed
469 0 : IF (output_unit > 0) THEN
470 0 : WRITE (output_unit, *) "Shapiro-Wilks normality test (Std. Dev.)", sdv_test_passed
471 : END IF
472 : ! Test of von Neumann for serial correlation
473 : ! - Standard Deviation
474 0 : CALL vn_test(wrk, ng_points, r, u, prob)
475 0 : PRINT *, prob, fe_env%conv_par%vn_conf_lm
476 0 : sdv_test_passed = prob <= fe_env%conv_par%vn_conf_lm
477 0 : fe_env%conv_par%test_vn = fe_env%conv_par%test_vn .AND. sdv_test_passed
478 0 : IF (output_unit > 0) THEN
479 0 : WRITE (output_unit, *) "von Neumann serial correlation test (Std. Dev.)", sdv_test_passed
480 : END IF
481 : END DO
482 : END DO
483 0 : CALL destroy_tmp_data(fe_env, wrk, ng_points)
484 : ELSE
485 0 : CALL destroy_tmp_data(fe_env, wrk, ng_points)
486 : END IF
487 0 : CALL timestop(handle)
488 0 : END SUBROUTINE ui_check_norm_sc_low
489 :
490 : ! **************************************************************************************************
491 : !> \brief Convergence criteria (Error on average and covariance matrix)
492 : !> for free energy method
493 : !> \param fe_env ...
494 : !> \param converged ...
495 : !> \param nr_points ...
496 : !> \param output_unit which unit to print to
497 : !> \par History
498 : !> Teodoro Laino (01.2007) [tlaino]
499 : ! **************************************************************************************************
500 0 : SUBROUTINE ui_check_convergence(fe_env, converged, nr_points, output_unit)
501 : TYPE(free_energy_type), POINTER :: fe_env
502 : LOGICAL, INTENT(OUT) :: converged
503 : INTEGER, INTENT(IN) :: nr_points, output_unit
504 :
505 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ui_check_convergence'
506 :
507 : INTEGER :: handle, i, ic, ncolvar, ng_points
508 : LOGICAL :: test_passed
509 : REAL(KIND=dp) :: max_error_avg, max_error_std
510 : REAL(KIND=dp), DIMENSION(:), POINTER :: avg_std, avgmx
511 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cov_std, covmx
512 :
513 0 : CALL timeset(routineN, handle)
514 0 : converged = .FALSE.
515 0 : ncolvar = fe_env%ncolvar
516 : NULLIFY (avgmx, avg_std, covmx, cov_std)
517 0 : CALL ui_check_norm_sc(fe_env, test_passed, nr_points, output_unit)
518 0 : IF (test_passed) THEN
519 0 : ng_points = nr_points/fe_env%conv_par%cg_width
520 : ! We can finally compute the error on average and covariance matrix
521 : ! and check if we converged..
522 0 : CALL create_tmp_data(fe_env, ng_points=ng_points, ncolvar=ncolvar)
523 0 : CALL create_csg_data(fe_env, ng_points, output_unit)
524 0 : ALLOCATE (covmx(ncolvar, ncolvar))
525 0 : ALLOCATE (avgmx(ncolvar))
526 0 : ALLOCATE (cov_std(ncolvar*(ncolvar + 1)/2, ncolvar*(ncolvar + 1)/2))
527 0 : ALLOCATE (avg_std(ncolvar))
528 0 : covmx = 0.0_dp
529 0 : avgmx = 0.0_dp
530 0 : DO i = 1, ng_points
531 0 : covmx = covmx + fe_env%cg_data(i)%var
532 0 : avgmx = avgmx + fe_env%cg_data(i)%avg
533 : END DO
534 0 : covmx = covmx/REAL(ng_points, KIND=dp)
535 0 : avgmx = avgmx/REAL(ng_points, KIND=dp)
536 :
537 : ! Compute errors on average and standard deviation
538 0 : CALL compute_avg_std_errors(fe_env, ncolvar, avgmx, covmx, avg_std, cov_std)
539 0 : IF (output_unit > 0) THEN
540 0 : WRITE (output_unit, *) "pippo", avgmx, covmx
541 0 : WRITE (output_unit, *) "pippo", avg_std, cov_std
542 : END IF
543 : ! Convergence of the averages
544 0 : max_error_avg = SQRT(MAXVAL(ABS(avg_std))/REAL(ng_points, KIND=dp))/MINVAL(avgmx)
545 0 : max_error_std = SQRT(MAXVAL(ABS(cov_std))/REAL(ng_points, KIND=dp))/MINVAL(covmx)
546 0 : IF (max_error_avg <= fe_env%conv_par%eps_conv .AND. &
547 0 : max_error_std <= fe_env%conv_par%eps_conv) converged = .TRUE.
548 :
549 0 : IF (output_unit > 0) THEN
550 0 : WRITE (output_unit, '(/,T2,"CG SAMPLING LENGTH = ",I7,20X,"REQUESTED ACCURACY = ",E12.6)') ng_points, &
551 0 : fe_env%conv_par%eps_conv
552 0 : WRITE (output_unit, '(T50,"PRESENT ACCURACY AVG= ",E12.6)') max_error_avg
553 0 : WRITE (output_unit, '(T50,"PRESENT ACCURACY STD= ",E12.6)') max_error_std
554 0 : WRITE (output_unit, '(T50,"CONVERGED FE-DER = ",L12)') converged
555 :
556 0 : WRITE (output_unit, '(/,T33, "COVARIANCE MATRIX")')
557 0 : WRITE (output_unit, '(T8,'//cp_to_string(ncolvar)//'(3X,I7,6X))') (ic, ic=1, ncolvar)
558 0 : DO ic = 1, ncolvar
559 0 : WRITE (output_unit, '(T2,I6,'//cp_to_string(ncolvar)//'(3X,E12.6,1X))') ic, covmx(ic, :)
560 : END DO
561 0 : WRITE (output_unit, '(T33, "ERROR OF COVARIANCE MATRIX")')
562 0 : WRITE (output_unit, '(T8,'//cp_to_string(ncolvar)//'(3X,I7,6X))') (ic, ic=1, ncolvar)
563 0 : DO ic = 1, ncolvar
564 0 : WRITE (output_unit, '(T2,I6,'//cp_to_string(ncolvar)//'(3X,E12.6,1X))') ic, cov_std(ic, :)
565 : END DO
566 :
567 0 : WRITE (output_unit, '(/,T2,"COLVAR Nr.",18X,13X,"AVERAGE",13X,"STANDARD DEVIATION")')
568 : WRITE (output_unit, '(T2,"CV",I8,21X,7X,E12.6,14X,E12.6)') &
569 0 : (ic, avgmx(ic), SQRT(ABS(avg_std(ic))), ic=1, ncolvar)
570 : END IF
571 0 : CALL destroy_tmp_data(fe_env, ng_points=ng_points)
572 0 : DEALLOCATE (covmx)
573 0 : DEALLOCATE (avgmx)
574 0 : DEALLOCATE (cov_std)
575 0 : DEALLOCATE (avg_std)
576 : END IF
577 0 : CALL timestop(handle)
578 0 : END SUBROUTINE ui_check_convergence
579 :
580 : ! **************************************************************************************************
581 : !> \brief Computes the errors on averages and standard deviations for a
582 : !> correlation-independent coarse grained data set
583 : !> \param fe_env ...
584 : !> \param ncolvar ...
585 : !> \param avgmx ...
586 : !> \param covmx ...
587 : !> \param avg_std ...
588 : !> \param cov_std ...
589 : !> \par History
590 : !> Teodoro Laino (02.2007) [tlaino]
591 : ! **************************************************************************************************
592 0 : SUBROUTINE compute_avg_std_errors(fe_env, ncolvar, avgmx, covmx, avg_std, cov_std)
593 : TYPE(free_energy_type), POINTER :: fe_env
594 : INTEGER, INTENT(IN) :: ncolvar
595 : REAL(KIND=dp), DIMENSION(:), POINTER :: avgmx
596 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: covmx
597 : REAL(KIND=dp), DIMENSION(:), POINTER :: avg_std
598 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cov_std
599 :
600 : INTEGER :: i, ind, j, k, nvar
601 : REAL(KIND=dp) :: fac
602 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: awrk, eig, tmp
603 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: wrk
604 :
605 : ! Averages
606 :
607 0 : nvar = ncolvar
608 0 : ALLOCATE (wrk(nvar, nvar))
609 0 : ALLOCATE (eig(nvar))
610 0 : fac = REAL(SIZE(fe_env%cg_data), KIND=dp)
611 0 : wrk = 0.0_dp
612 0 : eig = 0.0_dp
613 0 : DO k = 1, SIZE(fe_env%cg_data)
614 0 : DO j = 1, nvar
615 0 : DO i = j, nvar
616 0 : wrk(i, j) = wrk(i, j) + fe_env%cg_data(k)%avg(i)*fe_env%cg_data(k)%avg(j)
617 : END DO
618 : END DO
619 : END DO
620 0 : DO j = 1, nvar
621 0 : DO i = j, nvar
622 0 : wrk(i, j) = wrk(i, j) - avgmx(i)*avgmx(j)*fac
623 0 : wrk(j, i) = wrk(i, j)
624 : END DO
625 : END DO
626 0 : wrk = wrk/(fac - 1.0_dp)
627 : ! Diagonalize the covariance matrix and check for the maximum error
628 0 : CALL diamat_all(wrk, eig)
629 0 : DO i = 1, nvar
630 0 : avg_std(i) = eig(i)
631 : END DO
632 0 : DEALLOCATE (wrk)
633 0 : DEALLOCATE (eig)
634 : ! Standard Deviations
635 0 : nvar = ncolvar*(ncolvar + 1)/2
636 0 : ALLOCATE (wrk(nvar, nvar))
637 0 : ALLOCATE (eig(nvar))
638 0 : ALLOCATE (awrk(nvar))
639 0 : ALLOCATE (tmp(nvar))
640 0 : wrk = 0.0_dp
641 0 : eig = 0.0_dp
642 : ind = 0
643 0 : DO i = 1, ncolvar
644 0 : DO j = i, ncolvar
645 0 : ind = ind + 1
646 0 : awrk(ind) = covmx(i, j)
647 : END DO
648 : END DO
649 0 : DO k = 1, SIZE(fe_env%cg_data)
650 : ind = 0
651 0 : DO i = 1, ncolvar
652 0 : DO j = i, ncolvar
653 0 : ind = ind + 1
654 0 : tmp(ind) = fe_env%cg_data(k)%var(i, j)
655 : END DO
656 : END DO
657 0 : DO i = 1, nvar
658 0 : DO j = i, nvar
659 0 : wrk(i, j) = wrk(i, j) + tmp(i)*tmp(j) - awrk(i)*awrk(j)
660 : END DO
661 : END DO
662 : END DO
663 0 : DO i = 1, nvar
664 0 : DO j = i, nvar
665 0 : wrk(i, j) = wrk(i, j) - fac*awrk(i)*awrk(j)
666 0 : wrk(j, i) = wrk(i, j)
667 : END DO
668 : END DO
669 0 : wrk = wrk/(fac - 1.0_dp)
670 : ! Diagonalize the covariance matrix and check for the maximum error
671 0 : CALL diamat_all(wrk, eig)
672 0 : ind = 0
673 0 : DO i = 1, ncolvar
674 0 : DO j = i, ncolvar
675 0 : ind = ind + 1
676 0 : cov_std(i, j) = eig(ind)
677 0 : cov_std(j, i) = cov_std(i, j)
678 : END DO
679 : END DO
680 0 : DEALLOCATE (wrk)
681 0 : DEALLOCATE (eig)
682 0 : DEALLOCATE (awrk)
683 0 : DEALLOCATE (tmp)
684 :
685 0 : END SUBROUTINE compute_avg_std_errors
686 :
687 : ! **************************************************************************************************
688 : !> \brief Computes the covariance matrix
689 : !> \param fe_env ...
690 : !> \param cg_index ...
691 : !> \param istart ...
692 : !> \param iend ...
693 : !> \param output_unit which unit to print to
694 : !> \param covmx ...
695 : !> \param avgs ...
696 : !> \par History
697 : !> Teodoro Laino (01.2007) [tlaino]
698 : ! **************************************************************************************************
699 10 : SUBROUTINE eval_cov_matrix(fe_env, cg_index, istart, iend, output_unit, covmx, avgs)
700 : TYPE(free_energy_type), POINTER :: fe_env
701 : INTEGER, INTENT(IN) :: cg_index, istart, iend, output_unit
702 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: covmx
703 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: avgs
704 :
705 : CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_cov_matrix'
706 :
707 : INTEGER :: handle, ic, jc, jstep, ncolvar, nlength
708 : REAL(KIND=dp) :: tmp_ic, tmp_jc
709 : TYPE(ui_var_type), POINTER :: cv
710 :
711 10 : CALL timeset(routineN, handle)
712 10 : ncolvar = fe_env%ncolvar
713 10 : nlength = iend - istart + 1
714 20 : fe_env%cg_data(cg_index)%avg = 0.0_dp
715 30 : fe_env%cg_data(cg_index)%var = 0.0_dp
716 10 : IF (nlength > 1) THEN
717 : ! Update the info on averages and variances
718 40 : DO jstep = istart, iend
719 60 : DO ic = 1, ncolvar
720 30 : cv => fe_env%uivar(ic)
721 30 : tmp_ic = cv%ss(jstep)
722 60 : fe_env%cg_data(cg_index)%avg(ic) = fe_env%cg_data(cg_index)%avg(ic) + tmp_ic
723 : END DO
724 70 : DO ic = 1, ncolvar
725 30 : cv => fe_env%uivar(ic)
726 30 : tmp_ic = cv%ss(jstep)
727 90 : DO jc = 1, ic
728 30 : cv => fe_env%uivar(jc)
729 30 : tmp_jc = cv%ss(jstep)
730 60 : fe_env%cg_data(cg_index)%var(jc, ic) = fe_env%cg_data(cg_index)%var(jc, ic) + tmp_ic*tmp_jc
731 : END DO
732 : END DO
733 : END DO
734 : ! Normalized the variances and the averages
735 : ! Unbiased estimator
736 30 : fe_env%cg_data(cg_index)%var = fe_env%cg_data(cg_index)%var/REAL(nlength - 1, KIND=dp)
737 20 : fe_env%cg_data(cg_index)%avg = fe_env%cg_data(cg_index)%avg/REAL(nlength, KIND=dp)
738 : ! Compute the covariance matrix
739 20 : DO ic = 1, ncolvar
740 10 : tmp_ic = fe_env%cg_data(cg_index)%avg(ic)
741 30 : DO jc = 1, ic
742 10 : tmp_jc = fe_env%cg_data(cg_index)%avg(jc)*REAL(nlength, KIND=dp)/REAL(nlength - 1, KIND=dp)
743 10 : fe_env%cg_data(cg_index)%var(jc, ic) = fe_env%cg_data(cg_index)%var(jc, ic) - tmp_ic*tmp_jc
744 20 : fe_env%cg_data(cg_index)%var(ic, jc) = fe_env%cg_data(cg_index)%var(jc, ic)
745 : END DO
746 : END DO
747 10 : IF (output_unit > 0) THEN
748 20 : WRITE (output_unit, *) "eval_cov_matrix", istart, iend, fe_env%cg_data(cg_index)%avg, fe_env%cg_data(cg_index)%var
749 : END IF
750 10 : IF (PRESENT(covmx)) covmx = fe_env%cg_data(cg_index)%var
751 10 : IF (PRESENT(avgs)) avgs = fe_env%cg_data(cg_index)%avg
752 : END IF
753 10 : CALL timestop(handle)
754 10 : END SUBROUTINE eval_cov_matrix
755 :
756 : ! **************************************************************************************************
757 : !> \brief Dumps information when performing an alchemical change run
758 : !> \param my_val ...
759 : !> \param my_par ...
760 : !> \param coupling_function ...
761 : !> \param dx ...
762 : !> \param lerr ...
763 : !> \param fe_section ...
764 : !> \param nforce_eval ...
765 : !> \param cum_res ...
766 : !> \param istep ...
767 : !> \param beta ...
768 : !> \author Teodoro Laino - University of Zurich [tlaino] - 05.2007
769 : ! **************************************************************************************************
770 340 : SUBROUTINE dump_ac_info(my_val, my_par, coupling_function, dx, lerr, fe_section, &
771 : nforce_eval, cum_res, istep, beta)
772 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
773 : INTENT(INOUT) :: my_val
774 : CHARACTER(LEN=default_string_length), &
775 : ALLOCATABLE, DIMENSION(:), INTENT(IN) :: my_par
776 : CHARACTER(LEN=default_path_length), INTENT(IN) :: coupling_function
777 : REAL(KIND=dp), INTENT(IN) :: dx, lerr
778 : TYPE(section_vals_type), POINTER :: fe_section
779 : INTEGER, INTENT(IN) :: nforce_eval
780 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cum_res
781 : INTEGER, POINTER :: istep
782 : REAL(KIND=dp), INTENT(IN) :: beta
783 :
784 : CHARACTER(LEN=default_path_length) :: weighting_function
785 : CHARACTER(LEN=default_string_length) :: def_error, par, this_error
786 : INTEGER :: i, iforce_eval, ipar, isize, iw, j, &
787 : NEquilStep
788 : REAL(KIND=dp) :: avg_BP, avg_DET, avg_DUE, d_ene_w, dedf, &
789 : ene_w, err, Err_DET, Err_DUE, std_DET, &
790 : std_DUE, tmp, tmp2, wfac
791 : TYPE(cp_logger_type), POINTER :: logger
792 : TYPE(section_vals_type), POINTER :: alch_section
793 :
794 170 : logger => cp_get_default_logger()
795 170 : alch_section => section_vals_get_subs_vals(fe_section, "ALCHEMICAL_CHANGE")
796 170 : CALL section_vals_val_get(alch_section, "PARAMETER", c_val=par)
797 570 : DO i = 1, SIZE(my_par)
798 570 : IF (my_par(i) == par) EXIT
799 : END DO
800 170 : CPASSERT(i <= SIZE(my_par))
801 170 : ipar = i
802 :
803 : ! Two functions should be defined in separate input sections:
804 : ! 1 for the FORCE_EVAL/MIXED/GENERIC/MIXING_FUNCTION
805 : ! 2 for the MOTION/FREE_ENERGY/ALCHEMICAL_CHANGE/WEIGHTING_FUNCTION
806 : ! The definitions for parameters, variables, dx and error_limit all
807 : ! follow from 1, with which 2 is assumed consistent.
808 : ! MOTION/FREE_ENERGY/ALCHEMICAL_CHANGE/PARAMETER should be one of the
809 : ! symbols defined in FORCE_EVAL/MIXED/GENERIC/PARAMETERS.
810 170 : CALL initf(2)
811 170 : CALL parsef(1, TRIM(coupling_function), my_par)
812 170 : dedf = evalfd(1, ipar, my_val, dx, err)
813 170 : IF (ABS(err) > lerr) THEN
814 0 : WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
815 0 : WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
816 0 : CALL compress(this_error, .TRUE.)
817 0 : CALL compress(def_error, .TRUE.)
818 : CALL cp_warn(__LOCATION__, &
819 : 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
820 : ' Error '//TRIM(this_error)//' in computing numerical derivatives larger than'// &
821 0 : TRIM(def_error)//' .')
822 : END IF
823 :
824 : ! We must print now the energy of the biased system, the weigthing energy
825 : ! and the derivative w.r.t.the coupling parameter of the biased energy
826 : ! Retrieve the expression of the weighting function:
827 170 : CALL section_vals_val_get(alch_section, "WEIGHTING_FUNCTION", c_val=weighting_function)
828 170 : CALL compress(weighting_function, full=.TRUE.)
829 170 : CALL parsef(2, TRIM(weighting_function), my_par)
830 170 : ene_w = evalf(2, my_val)
831 170 : d_ene_w = evalfd(2, ipar, my_val, dx, err)
832 170 : IF (ABS(err) > lerr) THEN
833 0 : WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
834 0 : WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
835 0 : CALL compress(this_error, .TRUE.)
836 0 : CALL compress(def_error, .TRUE.)
837 : CALL cp_warn(__LOCATION__, &
838 : 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
839 : ' Error '//TRIM(this_error)//' in computing numerical derivatives larger than'// &
840 0 : TRIM(def_error)//' .')
841 : END IF
842 170 : CALL finalizef()
843 :
844 170 : CALL section_vals_val_get(alch_section, "NEQUIL_STEPS", i_val=NEquilStep)
845 : ! Store results
846 170 : IF (istep > NEquilStep) THEN
847 170 : isize = SIZE(cum_res, 2) + 1
848 170 : CALL reallocate(cum_res, 1, 3, 1, isize)
849 170 : cum_res(1, isize) = dedf
850 170 : cum_res(2, isize) = dedf - d_ene_w
851 170 : cum_res(3, isize) = ene_w
852 : ! Compute derivative of biased and total energy
853 : ! Total Free Energy
854 1680 : avg_DET = SUM(cum_res(1, 1:isize))/REAL(isize, KIND=dp)
855 1680 : std_DET = SUM(cum_res(1, 1:isize)**2)/REAL(isize, KIND=dp)
856 : ! Unbiased Free Energy
857 1680 : avg_BP = SUM(cum_res(3, 1:isize))/REAL(isize, KIND=dp)
858 170 : wfac = 0.0_dp
859 1680 : DO j = 1, isize
860 1680 : wfac = wfac + EXP(beta*(cum_res(3, j) - avg_BP))
861 : END DO
862 170 : avg_DUE = 0.0_dp
863 170 : std_DUE = 0.0_dp
864 1680 : DO j = 1, isize
865 1510 : tmp = cum_res(2, j)
866 1510 : tmp2 = EXP(beta*(cum_res(3, j) - avg_BP))/wfac
867 1510 : avg_DUE = avg_DUE + tmp*tmp2
868 1680 : std_DUE = std_DUE + tmp**2*tmp2
869 : END DO
870 170 : IF (isize > 1) THEN
871 158 : Err_DUE = SQRT(std_DUE - avg_DUE**2)/SQRT(REAL(isize - 1, KIND=dp))
872 158 : Err_DET = SQRT(std_DET - avg_DET**2)/SQRT(REAL(isize - 1, KIND=dp))
873 : END IF
874 : ! Print info
875 : iw = cp_print_key_unit_nr(logger, fe_section, "FREE_ENERGY_INFO", &
876 170 : extension=".free_energy")
877 170 : IF (iw > 0) THEN
878 85 : WRITE (iw, '(T2,79("-"),T37," oOo ")')
879 285 : DO iforce_eval = 1, nforce_eval
880 : WRITE (iw, '(T2,"ALCHEMICAL CHANGE| FORCE_EVAL Nr.",I5,T48,"ENERGY (Hartree)= ",F15.9)') &
881 285 : iforce_eval, my_val(iforce_eval)
882 : END DO
883 : WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE OF TOTAL ENERGY [ PARAMETER (",A,") ]",T66,F15.9)') &
884 85 : TRIM(par), dedf
885 : WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE OF BIASED ENERGY [ PARAMETER (",A,") ]",T66,F15.9)') &
886 85 : TRIM(par), dedf - d_ene_w
887 : WRITE (iw, '(T2,"ALCHEMICAL CHANGE| BIASING UMBRELLA POTENTIAL ",T66,F15.9)') &
888 85 : ene_w
889 :
890 85 : IF (isize > 1) THEN
891 : WRITE (iw, '(/,T2,"ALCHEMICAL CHANGE| DERIVATIVE TOTAL FREE ENERGY ",T50,F15.9,1X,"+/-",1X,F11.9)') &
892 79 : avg_DET, Err_DET
893 : WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE UNBIASED FREE ENERGY ",T50,F15.9,1X,"+/-",1X,F11.9)') &
894 79 : avg_DUE, Err_DUE
895 : ELSE
896 : WRITE (iw, '(/,T2,"ALCHEMICAL CHANGE| DERIVATIVE TOTAL FREE ENERGY ",T50,F15.9,1X,"+/-",1X,T76,A)') &
897 6 : avg_DET, "UNDEF"
898 : WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE UNBIASED FREE ENERGY ",T50,F15.9,1X,"+/-",1X,T76,A)') &
899 6 : avg_DUE, "UNDEF"
900 : END IF
901 85 : WRITE (iw, '(T2,79("-"))')
902 : END IF
903 : END IF
904 170 : CALL cp_print_key_finished_output(iw, logger, fe_section, "FREE_ENERGY_INFO")
905 :
906 170 : END SUBROUTINE dump_ac_info
907 :
908 : END MODULE free_energy_methods
909 :
|
