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