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