(git:fc9bb57)
Loading...
Searching...
No Matches
free_energy_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
56CONTAINS
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 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 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 NULLIFY (force_env, istep, subsys, cv, simpar)
90 logger => cp_get_default_logger()
91 CALL timeset(routinen, handle)
92 converged = .false.
93 CALL get_md_env(md_env, force_env=force_env, fe_env=fe_env, simpar=simpar, &
94 itimes=istep)
95 ! Metadynamics is also a free energy calculation but is handled in a different
96 ! module.
97 IF (.NOT. ASSOCIATED(force_env%meta_env) .AND. ASSOCIATED(fe_env)) THEN
98 SELECT CASE (fe_env%type)
99 CASE (do_fe_ui)
100 ! Umbrella Integration..
101 CALL force_env_get(force_env, subsys=subsys)
102 fe_env%nr_points = fe_env%nr_points + 1
103 output_unit = cp_logger_get_default_io_unit(logger)
104 DO ic = 1, fe_env%ncolvar
105 cv => fe_env%uivar(ic)
106 icolvar = cv%icolvar
107 CALL colvar_eval_glob_f(icolvar, force_env)
108 CALL reallocate(cv%ss, 1, fe_env%nr_points)
109 cv%ss(fe_env%nr_points) = subsys%colvar_p(icolvar)%colvar%ss
110 IF (output_unit > 0) THEN
111 WRITE (output_unit, *) "COLVAR::", cv%ss(fe_env%nr_points)
112 END IF
113 END DO
114 stat_sign_points = fe_env%nr_points - fe_env%nr_rejected
115 IF (output_unit > 0) THEN
116 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 IF ((fe_env%conv_par%cg_width*fe_env%conv_par%cg_points <= stat_sign_points) .AND. &
120 (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 extension=".FreeEnergyLog", log_filename=.false.)
123 CALL print_fe_prolog(output_unit)
124 ! Trend test.. recomputes the number of statistically significant points..
125 CALL ui_check_trend(fe_env, fe_env%conv_par%test_k, stat_sign_points, output_unit)
126 stat_sign_points = fe_env%nr_points - fe_env%nr_rejected
127 ! Normality and serial correlation tests..
128 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 CALL ui_check_convergence(fe_env, converged, stat_sign_points, output_unit)
132 END IF
133 CALL print_fe_epilog(output_unit)
134 CALL cp_print_key_finished_output(output_unit, logger, fe_section, "FREE_ENERGY_INFO")
135 END IF
136 CASE (do_fe_ac)
137 CALL initf(2)
138 ! Alchemical Changes
139 IF (.NOT. ASSOCIATED(force_env%mixed_env)) THEN
140 CALL cp_abort(__location__, &
141 'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
142 ' Free Energy calculations require the definition of a mixed env!')
143 END IF
144 my_par => force_env%mixed_env%par
145 my_val => force_env%mixed_env%val
146 dx = force_env%mixed_env%dx
147 lerr = force_env%mixed_env%lerr
148 coupling_function = force_env%mixed_env%coupling_function
149 beta = 1/simpar%temp_ext
150 CALL parsef(1, trim(coupling_function), my_par)
151 nforce_eval = SIZE(force_env%sub_force_env)
152 CALL dump_ac_info(my_val, my_par, dx, lerr, fe_section, nforce_eval, &
153 fe_env%covmx, istep, beta)
154 CALL finalizef()
155 CASE DEFAULT
156 ! Do Nothing
157 END SELECT
158 END IF
159 CALL timestop(handle)
160
161 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 SUBROUTINE print_fe_prolog(output_unit)
170 INTEGER, INTENT(IN) :: output_unit
171
172 IF (output_unit > 0) THEN
173 WRITE (output_unit, '(T2,79("*"))')
174 WRITE (output_unit, '(T30,"FREE ENERGY CALCULATION",/)')
175 END IF
176 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 SUBROUTINE print_fe_epilog(output_unit)
185 INTEGER, INTENT(IN) :: output_unit
186
187 IF (output_unit > 0) THEN
188 WRITE (output_unit, '(T2,79("*"),/)')
189 END IF
190 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: wrk
213
214 CALL timeset(routinen, handle)
215 trend_free = .false.
216 test_avg = .true.
217 test_std = .true.
218 ncolvar = fe_env%ncolvar
219 ! Number of coarse grained points
220 IF (output_unit > 0) THEN
221 WRITE (output_unit, *) nr_points, fe_env%conv_par%cg_width
222 END IF
223 ng_points = nr_points/fe_env%conv_par%cg_width
224 my_reject = 0
225 ! Allocate storage
226 CALL create_tmp_data(fe_env, wrk, ng_points, ncolvar)
227 ! Compute the Coarse Grained data set using a reverse cumulative strategy
228 CALL create_csg_data(fe_env, ng_points, output_unit)
229 ! Test on coarse grained average
230 DO j = 1, ncolvar
231 ii = 1
232 DO i = ng_points, 1, -1
233 wrk(ii) = fe_env%cg_data(i)%avg(j)
234 ii = ii + 1
235 END DO
236 DO i = my_reject + 1, ng_points
237 IF ((ng_points - my_reject) < min_sample_size) THEN
238 my_reject = max(0, my_reject - 1)
239 test_avg = .false.
240 EXIT
241 END IF
242 CALL k_test(wrk, my_reject + 1, ng_points, tau, z, prob)
243 print *, prob, fe_env%conv_par%k_conf_lm
244 IF (prob < fe_env%conv_par%k_conf_lm) EXIT
245 my_reject = my_reject + 1
246 END DO
247 my_reject = min(ng_points, my_reject)
248 END DO
249 rejected_points = my_reject*fe_env%conv_par%cg_width
250 ! Print some info
251 IF (output_unit > 0) THEN
252 WRITE (output_unit, *) "Kendall trend test (Average)", test_avg, &
253 "number of points rejected:", rejected_points + fe_env%nr_rejected
254 WRITE (output_unit, *) "Reject Nr.", my_reject, " coarse grained points testing average"
255 END IF
256 ! Test on coarse grained covariance matrix
257 DO j = 1, ncolvar
258 DO k = j, ncolvar
259 ii = 1
260 DO i = ng_points, 1, -1
261 wrk(ii) = fe_env%cg_data(i)%var(j, k)
262 ii = ii + 1
263 END DO
264 DO i = my_reject + 1, ng_points
265 IF ((ng_points - my_reject) < min_sample_size) THEN
266 my_reject = max(0, my_reject - 1)
267 test_std = .false.
268 EXIT
269 END IF
270 CALL k_test(wrk, my_reject + 1, ng_points, tau, z, prob)
271 print *, prob, fe_env%conv_par%k_conf_lm
272 IF (prob < fe_env%conv_par%k_conf_lm) EXIT
273 my_reject = my_reject + 1
274 END DO
275 my_reject = min(ng_points, my_reject)
276 END DO
277 END DO
278 rejected_points = my_reject*fe_env%conv_par%cg_width
279 fe_env%nr_rejected = fe_env%nr_rejected + rejected_points
280 trend_free = test_avg .AND. test_std
281 ! Print some info
282 IF (output_unit > 0) THEN
283 WRITE (output_unit, *) "Kendall trend test (Std. Dev.)", test_std, &
284 "number of points rejected:", fe_env%nr_rejected
285 WRITE (output_unit, *) "Reject Nr.", my_reject, " coarse grained points testing standard dev."
286 WRITE (output_unit, *) "Kendall test passed:", trend_free
287 END IF
288 ! Release storage
289 CALL destroy_tmp_data(fe_env, wrk, ng_points)
290 CALL timestop(handle)
291 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 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 ALLOCATE (fe_env%cg_data(ng_points))
310 DO i = 1, ng_points
311 ALLOCATE (fe_env%cg_data(i)%avg(ncolvar))
312 ALLOCATE (fe_env%cg_data(i)%var(ncolvar, ncolvar))
313 END DO
314 IF (PRESENT(wrk)) THEN
315 ALLOCATE (wrk(ng_points))
316 END IF
317 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 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 DO i = 1, ng_points
335 DEALLOCATE (fe_env%cg_data(i)%avg)
336 DEALLOCATE (fe_env%cg_data(i)%var)
337 END DO
338 DEALLOCATE (fe_env%cg_data)
339 IF (PRESENT(wrk)) THEN
340 DEALLOCATE (wrk)
341 END IF
342 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 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 DO i = 1, ng_points
359 istart = fe_env%nr_points - (i)*fe_env%conv_par%cg_width + 1
360 iend = fe_env%nr_points - (i - 1)*fe_env%conv_par%cg_width
361 IF (output_unit > 0) THEN
362 WRITE (output_unit, *) istart, iend
363 END IF
364 CALL eval_cov_matrix(fe_env, cg_index=i, istart=istart, iend=iend, output_unit=output_unit)
365 END DO
366
367 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 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 CALL timeset(routinen, handle)
389 test_passed = .false.
390 DO WHILE (fe_env%conv_par%cg_width < fe_env%conv_par%max_cg_width)
391 ng_points = nr_points/fe_env%conv_par%cg_width
392 print *, ng_points
393 IF (ng_points < min_sample_size) EXIT
394 CALL ui_check_norm_sc_low(fe_env, nr_points, output_unit)
395 test_passed = fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw
396 IF (test_passed) EXIT
397 fe_env%conv_par%cg_width = fe_env%conv_par%cg_width + 1
398 IF (output_unit > 0) THEN
399 WRITE (output_unit, *) "New coarse grained width:", fe_env%conv_par%cg_width
400 END IF
401 END DO
402 IF (fe_env%conv_par%cg_width == fe_env%conv_par%max_cg_width .AND. (.NOT. (test_passed))) THEN
403 CALL ui_check_norm_sc_low(fe_env, nr_points, output_unit)
404 test_passed = fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw
405 END IF
406 CALL timestop(handle)
407 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: wrk
428
429 CALL timeset(routinen, handle)
430 ncolvar = fe_env%ncolvar
431 ! Compute the Coarse Grained data set using a reverse cumulative strategy
432 fe_env%conv_par%test_sw = .false.
433 fe_env%conv_par%test_vn = .false.
434 ! Number of coarse grained points
435 avg_test_passed = .true.
436 sdv_test_passed = .true.
437 ng_points = nr_points/fe_env%conv_par%cg_width
438 CALL create_tmp_data(fe_env, wrk, ng_points, ncolvar)
439 CALL create_csg_data(fe_env, ng_points, output_unit)
440 ! Testing Averages
441 DO j = 1, ncolvar
442 DO i = 1, ng_points
443 wrk(i) = fe_env%cg_data(i)%avg(j)
444 END DO
445 ! Test of Shapiro - Wilks for normality
446 ! - Average
447 CALL sw_test(wrk, ng_points, w, pw)
448 print *, 1.0_dp - pw, fe_env%conv_par%sw_conf_lm
449 avg_test_passed = (1.0_dp - pw) <= fe_env%conv_par%sw_conf_lm
450 fe_env%conv_par%test_sw = avg_test_passed
451 IF (output_unit > 0) THEN
452 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 CALL vn_test(wrk, ng_points, r, u, prob)
457 print *, prob, fe_env%conv_par%vn_conf_lm
458 avg_test_passed = prob <= fe_env%conv_par%vn_conf_lm
459 fe_env%conv_par%test_vn = avg_test_passed
460 IF (output_unit > 0) THEN
461 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 IF (fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw) THEN
466 ! Testing Standard Deviations
467 DO j = 1, ncolvar
468 DO k = j, ncolvar
469 DO i = 1, ng_points
470 wrk(i) = fe_env%cg_data(i)%var(j, k)
471 END DO
472 ! Test of Shapiro - Wilks for normality
473 ! - Standard Deviation
474 CALL sw_test(wrk, ng_points, w, pw)
475 print *, 1.0_dp - pw, fe_env%conv_par%sw_conf_lm
476 sdv_test_passed = (1.0_dp - pw) <= fe_env%conv_par%sw_conf_lm
477 fe_env%conv_par%test_sw = fe_env%conv_par%test_sw .AND. sdv_test_passed
478 IF (output_unit > 0) THEN
479 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 CALL vn_test(wrk, ng_points, r, u, prob)
484 print *, prob, fe_env%conv_par%vn_conf_lm
485 sdv_test_passed = prob <= fe_env%conv_par%vn_conf_lm
486 fe_env%conv_par%test_vn = fe_env%conv_par%test_vn .AND. sdv_test_passed
487 IF (output_unit > 0) THEN
488 WRITE (output_unit, *) "von Neumann serial correlation test (Std. Dev.)", sdv_test_passed
489 END IF
490 END DO
491 END DO
492 CALL destroy_tmp_data(fe_env, wrk, ng_points)
493 ELSE
494 CALL destroy_tmp_data(fe_env, wrk, ng_points)
495 END IF
496 CALL timestop(handle)
497 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 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 CALL timeset(routinen, handle)
523 converged = .false.
524 ncolvar = fe_env%ncolvar
525 NULLIFY (avgmx, avg_std, covmx, cov_std)
526 CALL ui_check_norm_sc(fe_env, test_passed, nr_points, output_unit)
527 IF (test_passed) THEN
528 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 CALL create_tmp_data(fe_env, ng_points=ng_points, ncolvar=ncolvar)
532 CALL create_csg_data(fe_env, ng_points, output_unit)
533 ALLOCATE (covmx(ncolvar, ncolvar))
534 ALLOCATE (avgmx(ncolvar))
535 ALLOCATE (cov_std(ncolvar*(ncolvar + 1)/2, ncolvar*(ncolvar + 1)/2))
536 ALLOCATE (avg_std(ncolvar))
537 covmx = 0.0_dp
538 avgmx = 0.0_dp
539 DO i = 1, ng_points
540 covmx = covmx + fe_env%cg_data(i)%var
541 avgmx = avgmx + fe_env%cg_data(i)%avg
542 END DO
543 covmx = covmx/real(ng_points, kind=dp)
544 avgmx = avgmx/real(ng_points, kind=dp)
545
546 ! Compute errors on average and standard deviation
547 CALL compute_avg_std_errors(fe_env, ncolvar, avgmx, covmx, avg_std, cov_std)
548 IF (output_unit > 0) THEN
549 WRITE (output_unit, *) "pippo", avgmx, covmx
550 WRITE (output_unit, *) "pippo", avg_std, cov_std
551 END IF
552 ! Convergence of the averages
553 max_error_avg = sqrt(maxval(abs(avg_std))/real(ng_points, kind=dp))/minval(avgmx)
554 max_error_std = sqrt(maxval(abs(cov_std))/real(ng_points, kind=dp))/minval(covmx)
555 IF (max_error_avg <= fe_env%conv_par%eps_conv .AND. &
556 max_error_std <= fe_env%conv_par%eps_conv) converged = .true.
557
558 IF (output_unit > 0) THEN
559 WRITE (output_unit, '(/,T2,"CG SAMPLING LENGTH = ",I7,20X,"REQUESTED ACCURACY = ",E12.6)') ng_points, &
560 fe_env%conv_par%eps_conv
561 WRITE (output_unit, '(T50,"PRESENT ACCURACY AVG= ",E12.6)') max_error_avg
562 WRITE (output_unit, '(T50,"PRESENT ACCURACY STD= ",E12.6)') max_error_std
563 WRITE (output_unit, '(T50,"CONVERGED FE-DER = ",L12)') converged
564
565 WRITE (output_unit, '(/,T33, "COVARIANCE MATRIX")')
566 WRITE (output_unit, '(T8,'//cp_to_string(ncolvar)//'(3X,I7,6X))') (ic, ic=1, ncolvar)
567 DO ic = 1, ncolvar
568 WRITE (output_unit, '(T2,I6,'//cp_to_string(ncolvar)//'(3X,E12.6,1X))') ic, covmx(ic, :)
569 END DO
570 WRITE (output_unit, '(T33, "ERROR OF COVARIANCE MATRIX")')
571 WRITE (output_unit, '(T8,'//cp_to_string(ncolvar)//'(3X,I7,6X))') (ic, ic=1, ncolvar)
572 DO ic = 1, ncolvar
573 WRITE (output_unit, '(T2,I6,'//cp_to_string(ncolvar)//'(3X,E12.6,1X))') ic, cov_std(ic, :)
574 END DO
575
576 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 (ic, avgmx(ic), sqrt(abs(avg_std(ic))), ic=1, ncolvar)
579 END IF
580 CALL destroy_tmp_data(fe_env, ng_points=ng_points)
581 DEALLOCATE (covmx)
582 DEALLOCATE (avgmx)
583 DEALLOCATE (cov_std)
584 DEALLOCATE (avg_std)
585 END IF
586 CALL timestop(handle)
587 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: awrk, eig, tmp
612 REAL(kind=dp), DIMENSION(:, :), POINTER :: wrk
613
614! Averages
615
616 nvar = ncolvar
617 ALLOCATE (wrk(nvar, nvar))
618 ALLOCATE (eig(nvar))
619 fac = real(SIZE(fe_env%cg_data), kind=dp)
620 wrk = 0.0_dp
621 eig = 0.0_dp
622 DO k = 1, SIZE(fe_env%cg_data)
623 DO j = 1, nvar
624 DO i = j, nvar
625 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 DO j = 1, nvar
630 DO i = j, nvar
631 wrk(i, j) = wrk(i, j) - avgmx(i)*avgmx(j)*fac
632 wrk(j, i) = wrk(i, j)
633 END DO
634 END DO
635 wrk = wrk/(fac - 1.0_dp)
636 ! Diagonalize the covariance matrix and check for the maximum error
637 CALL diamat_all(wrk, eig)
638 DO i = 1, nvar
639 avg_std(i) = eig(i)
640 END DO
641 DEALLOCATE (wrk)
642 DEALLOCATE (eig)
643 ! Standard Deviations
644 nvar = ncolvar*(ncolvar + 1)/2
645 ALLOCATE (wrk(nvar, nvar))
646 ALLOCATE (eig(nvar))
647 ALLOCATE (awrk(nvar))
648 ALLOCATE (tmp(nvar))
649 wrk = 0.0_dp
650 eig = 0.0_dp
651 ind = 0
652 DO i = 1, ncolvar
653 DO j = i, ncolvar
654 ind = ind + 1
655 awrk(ind) = covmx(i, j)
656 END DO
657 END DO
658 DO k = 1, SIZE(fe_env%cg_data)
659 ind = 0
660 DO i = 1, ncolvar
661 DO j = i, ncolvar
662 ind = ind + 1
663 tmp(ind) = fe_env%cg_data(k)%var(i, j)
664 END DO
665 END DO
666 DO i = 1, nvar
667 DO j = i, nvar
668 wrk(i, j) = wrk(i, j) + tmp(i)*tmp(j) - awrk(i)*awrk(j)
669 END DO
670 END DO
671 END DO
672 DO i = 1, nvar
673 DO j = i, nvar
674 wrk(i, j) = wrk(i, j) - fac*awrk(i)*awrk(j)
675 wrk(j, i) = wrk(i, j)
676 END DO
677 END DO
678 wrk = wrk/(fac - 1.0_dp)
679 ! Diagonalize the covariance matrix and check for the maximum error
680 CALL diamat_all(wrk, eig)
681 ind = 0
682 DO i = 1, ncolvar
683 DO j = i, ncolvar
684 ind = ind + 1
685 cov_std(i, j) = eig(ind)
686 cov_std(j, i) = cov_std(i, j)
687 END DO
688 END DO
689 DEALLOCATE (wrk)
690 DEALLOCATE (eig)
691 DEALLOCATE (awrk)
692 DEALLOCATE (tmp)
693
694 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 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 CALL timeset(routinen, handle)
721 ncolvar = fe_env%ncolvar
722 nlength = iend - istart + 1
723 fe_env%cg_data(cg_index)%avg = 0.0_dp
724 fe_env%cg_data(cg_index)%var = 0.0_dp
725 IF (nlength > 1) THEN
726 ! Update the info on averages and variances
727 DO jstep = istart, iend
728 DO ic = 1, ncolvar
729 cv => fe_env%uivar(ic)
730 tmp_ic = cv%ss(jstep)
731 fe_env%cg_data(cg_index)%avg(ic) = fe_env%cg_data(cg_index)%avg(ic) + tmp_ic
732 END DO
733 DO ic = 1, ncolvar
734 cv => fe_env%uivar(ic)
735 tmp_ic = cv%ss(jstep)
736 DO jc = 1, ic
737 cv => fe_env%uivar(jc)
738 tmp_jc = cv%ss(jstep)
739 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 fe_env%cg_data(cg_index)%var = fe_env%cg_data(cg_index)%var/real(nlength - 1, kind=dp)
746 fe_env%cg_data(cg_index)%avg = fe_env%cg_data(cg_index)%avg/real(nlength, kind=dp)
747 ! Compute the covariance matrix
748 DO ic = 1, ncolvar
749 tmp_ic = fe_env%cg_data(cg_index)%avg(ic)
750 DO jc = 1, ic
751 tmp_jc = fe_env%cg_data(cg_index)%avg(jc)*real(nlength, kind=dp)/real(nlength - 1, kind=dp)
752 fe_env%cg_data(cg_index)%var(jc, ic) = fe_env%cg_data(cg_index)%var(jc, ic) - tmp_ic*tmp_jc
753 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 IF (output_unit > 0) THEN
757 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 IF (PRESENT(covmx)) covmx = fe_env%cg_data(cg_index)%var
760 IF (PRESENT(avgs)) avgs = fe_env%cg_data(cg_index)%avg
761 END IF
762 CALL timestop(handle)
763 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 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 logger => cp_get_default_logger()
801 alch_section => section_vals_get_subs_vals(fe_section, "ALCHEMICAL_CHANGE")
802 CALL section_vals_val_get(alch_section, "PARAMETER", c_val=par)
803 DO i = 1, SIZE(my_par)
804 IF (my_par(i) == par) EXIT
805 END DO
806 cpassert(i <= SIZE(my_par))
807 ipar = i
808 dedf = evalfd(1, ipar, my_val, dx, err)
809 IF (abs(err) > lerr) THEN
810 WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
811 WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
812 CALL compress(this_error, .true.)
813 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 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 CALL section_vals_val_get(alch_section, "WEIGHTING_FUNCTION", c_val=coupling_function)
824 CALL compress(coupling_function, full=.true.)
825 CALL parsef(2, trim(coupling_function), my_par)
826 ene_w = evalf(2, my_val)
827 d_ene_w = evalfd(2, ipar, my_val, dx, err)
828 IF (abs(err) > lerr) THEN
829 WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
830 WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
831 CALL compress(this_error, .true.)
832 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 trim(def_error)//' .')
837 END IF
838 CALL section_vals_val_get(alch_section, "NEQUIL_STEPS", i_val=nequilstep)
839 ! Store results
840 IF (istep > nequilstep) THEN
841 isize = SIZE(cum_res, 2) + 1
842 CALL reallocate(cum_res, 1, 3, 1, isize)
843 cum_res(1, isize) = dedf
844 cum_res(2, isize) = dedf - d_ene_w
845 cum_res(3, isize) = ene_w
846 ! Compute derivative of biased and total energy
847 ! Total Free Energy
848 avg_det = sum(cum_res(1, 1:isize))/real(isize, kind=dp)
849 std_det = sum(cum_res(1, 1:isize)**2)/real(isize, kind=dp)
850 ! Unbiased Free Energy
851 avg_bp = sum(cum_res(3, 1:isize))/real(isize, kind=dp)
852 wfac = 0.0_dp
853 DO j = 1, isize
854 wfac = wfac + exp(beta*(cum_res(3, j) - avg_bp))
855 END DO
856 avg_due = 0.0_dp
857 std_due = 0.0_dp
858 DO j = 1, isize
859 tmp = cum_res(2, j)
860 tmp2 = exp(beta*(cum_res(3, j) - avg_bp))/wfac
861 avg_due = avg_due + tmp*tmp2
862 std_due = std_due + tmp**2*tmp2
863 END DO
864 IF (isize > 1) THEN
865 err_due = sqrt(std_due - avg_due**2)/sqrt(real(isize - 1, kind=dp))
866 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 extension=".free_energy")
871 IF (iw > 0) THEN
872 WRITE (iw, '(T2,79("-"),T37," oOo ")')
873 DO iforce_eval = 1, nforce_eval
874 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| FORCE_EVAL Nr.",I5,T48,"ENERGY (Hartree)= ",F15.9)') &
875 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 trim(par), dedf
879 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE OF BIASED ENERGY [ PARAMETER (",A,") ]",T66,F15.9)') &
880 trim(par), dedf - d_ene_w
881 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| BIASING UMBRELLA POTENTIAL ",T66,F15.9)') &
882 ene_w
883
884 IF (isize > 1) THEN
885 WRITE (iw, '(/,T2,"ALCHEMICAL CHANGE| DERIVATIVE TOTAL FREE ENERGY ",T50,F15.9,1X,"+/-",1X,F11.9)') &
886 avg_det, err_det
887 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE UNBIASED FREE ENERGY ",T50,F15.9,1X,"+/-",1X,F11.9)') &
888 avg_due, err_due
889 ELSE
890 WRITE (iw, '(/,T2,"ALCHEMICAL CHANGE| DERIVATIVE TOTAL FREE ENERGY ",T50,F15.9,1X,"+/-",1X,T76,A)') &
891 avg_det, "UNDEF"
892 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE UNBIASED FREE ENERGY ",T50,F15.9,1X,"+/-",1X,T76,A)') &
893 avg_due, "UNDEF"
894 END IF
895 WRITE (iw, '(T2,79("-"))')
896 END IF
897 END IF
898 CALL cp_print_key_finished_output(iw, logger, fe_section, "FREE_ENERGY_INFO")
899
900 END SUBROUTINE dump_ac_info
901
902END MODULE free_energy_methods
903
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:56
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
memory_utilities::reallocate
Definition memory_utilities.F:27
colvar_methods
defines collective variables s({R}) and the derivative of this variable wrt R these can then be used ...
Definition colvar_methods.F:16
colvar_methods::colvar_eval_glob_f
subroutine, public colvar_eval_glob_f(icolvar, force_env)
evaluates the derivatives (dsdr) given and due to the given colvar
Definition colvar_methods.F:1757
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:344
fparser
This public domain function parser module is intended for applications where a set of mathematical ex...
Definition fparser.F:17
fparser::parsef
subroutine, public parsef(i, funcstr, var)
Parse ith function string FuncStr and compile it into bytecode.
Definition fparser.F:174
fparser::evalf
real(rn) function, public evalf(i, val)
...
Definition fparser.F:206
fparser::evalfd
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
Definition fparser.F:1097
fparser::finalizef
subroutine, public finalizef()
...
Definition fparser.F:127
fparser::initf
subroutine, public initf(n)
...
Definition fparser.F:156
free_energy_methods
Methods to perform free energy and free energy derivatives calculations.
Definition free_energy_methods.F:12
free_energy_methods::free_energy_evaluate
subroutine, public free_energy_evaluate(md_env, converged, fe_section)
Main driver for free energy calculations In this routine we handle specifically biased MD.
Definition free_energy_methods.F:68
free_energy_types
defines types for metadynamics calculation
Definition free_energy_types.F:13
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_fe_ac
integer, parameter, public do_fe_ac
Definition input_constants.F:377
input_constants::do_fe_ui
integer, parameter, public do_fe_ui
Definition input_constants.F:377
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:381
md_environment_types
Definition md_environment_types.F:15
md_environment_types::get_md_env
subroutine, public get_md_env(md_env, itimes, constant, used_time, cell, simpar, npt, force_env, para_env, reftraj, t, init, first_time, fe_env, thermostats, barostat, thermostat_coeff, thermostat_part, thermostat_shell, thermostat_baro, thermostat_fast, thermostat_slow, md_ener, averages, thermal_regions, ehrenfest_md)
get components of MD environment type
Definition md_environment_types.F:195
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
simpar_types
Type for storing MD parameters.
Definition simpar_types.F:14
statistical_methods
Methods to perform on the fly statistical analysis of data -) Schiferl and Wallace,...
Definition statistical_methods.F:16
statistical_methods::vn_test
subroutine, public vn_test(xdata, n, r, u, prob)
Von Neumann test for serial correlation.
Definition statistical_methods.F:398
statistical_methods::min_sample_size
integer, parameter, public min_sample_size
Definition statistical_methods.F:28
statistical_methods::k_test
subroutine, public k_test(xdata, istart, n, tau, z, prob)
Kandall's test for correlation.
Definition statistical_methods.F:356
statistical_methods::sw_test
subroutine, public sw_test(ix, n, w, pw)
Shapiro - Wilk's test or W-statistic to test normality of a distribution R94 APPL....
Definition statistical_methods.F:47
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::compress
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Definition string_utilities.F:3191
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
free_energy_types::free_energy_type
Definition free_energy_types.F:64
free_energy_types::ui_var_type
Definition free_energy_types.F:37
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
md_environment_types::md_environment_type
Definition md_environment_types.F:53
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30