(git:374b731)
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-2024 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)) &
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 my_par => force_env%mixed_env%par
144 my_val => force_env%mixed_env%val
145 dx = force_env%mixed_env%dx
146 lerr = force_env%mixed_env%lerr
147 coupling_function = force_env%mixed_env%coupling_function
148 beta = 1/simpar%temp_ext
149 CALL parsef(1, trim(coupling_function), my_par)
150 nforce_eval = SIZE(force_env%sub_force_env)
151 CALL dump_ac_info(my_val, my_par, dx, lerr, fe_section, nforce_eval, &
152 fe_env%covmx, istep, beta)
153 CALL finalizef()
154 CASE DEFAULT
155 ! Do Nothing
156 END SELECT
157 END IF
158 CALL timestop(handle)
159
160 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 SUBROUTINE print_fe_prolog(output_unit)
169 INTEGER, INTENT(IN) :: output_unit
170
171 IF (output_unit > 0) THEN
172 WRITE (output_unit, '(T2,79("*"))')
173 WRITE (output_unit, '(T30,"FREE ENERGY CALCULATION",/)')
174 END IF
175 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 SUBROUTINE print_fe_epilog(output_unit)
184 INTEGER, INTENT(IN) :: output_unit
185
186 IF (output_unit > 0) THEN
187 WRITE (output_unit, '(T2,79("*"),/)')
188 END IF
189 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: wrk
212
213 CALL timeset(routinen, handle)
214 trend_free = .false.
215 test_avg = .true.
216 test_std = .true.
217 ncolvar = fe_env%ncolvar
218 ! Number of coarse grained points
219 IF (output_unit > 0) THEN
220 WRITE (output_unit, *) nr_points, fe_env%conv_par%cg_width
221 END IF
222 ng_points = nr_points/fe_env%conv_par%cg_width
223 my_reject = 0
224 ! Allocate storage
225 CALL create_tmp_data(fe_env, wrk, ng_points, ncolvar)
226 ! Compute the Coarse Grained data set using a reverse cumulative strategy
227 CALL create_csg_data(fe_env, ng_points, output_unit)
228 ! Test on coarse grained average
229 DO j = 1, ncolvar
230 ii = 1
231 DO i = ng_points, 1, -1
232 wrk(ii) = fe_env%cg_data(i)%avg(j)
233 ii = ii + 1
234 END DO
235 DO i = my_reject + 1, ng_points
236 IF ((ng_points - my_reject) .LT. min_sample_size) THEN
237 my_reject = max(0, my_reject - 1)
238 test_avg = .false.
239 EXIT
240 END IF
241 CALL k_test(wrk, my_reject + 1, ng_points, tau, z, prob)
242 print *, prob, fe_env%conv_par%k_conf_lm
243 IF (prob < fe_env%conv_par%k_conf_lm) EXIT
244 my_reject = my_reject + 1
245 END DO
246 my_reject = min(ng_points, my_reject)
247 END DO
248 rejected_points = my_reject*fe_env%conv_par%cg_width
249 ! Print some info
250 IF (output_unit > 0) THEN
251 WRITE (output_unit, *) "Kendall trend test (Average)", test_avg, &
252 "number of points rejected:", rejected_points + fe_env%nr_rejected
253 WRITE (output_unit, *) "Reject Nr.", my_reject, " coarse grained points testing average"
254 END IF
255 ! Test on coarse grained covariance matrix
256 DO j = 1, ncolvar
257 DO k = j, ncolvar
258 ii = 1
259 DO i = ng_points, 1, -1
260 wrk(ii) = fe_env%cg_data(i)%var(j, k)
261 ii = ii + 1
262 END DO
263 DO i = my_reject + 1, ng_points
264 IF ((ng_points - my_reject) .LT. min_sample_size) THEN
265 my_reject = max(0, my_reject - 1)
266 test_std = .false.
267 EXIT
268 END IF
269 CALL k_test(wrk, my_reject + 1, ng_points, tau, z, prob)
270 print *, prob, fe_env%conv_par%k_conf_lm
271 IF (prob < fe_env%conv_par%k_conf_lm) EXIT
272 my_reject = my_reject + 1
273 END DO
274 my_reject = min(ng_points, my_reject)
275 END DO
276 END DO
277 rejected_points = my_reject*fe_env%conv_par%cg_width
278 fe_env%nr_rejected = fe_env%nr_rejected + rejected_points
279 trend_free = test_avg .AND. test_std
280 ! Print some info
281 IF (output_unit > 0) THEN
282 WRITE (output_unit, *) "Kendall trend test (Std. Dev.)", test_std, &
283 "number of points rejected:", fe_env%nr_rejected
284 WRITE (output_unit, *) "Reject Nr.", my_reject, " coarse grained points testing standard dev."
285 WRITE (output_unit, *) "Kendall test passed:", trend_free
286 END IF
287 ! Release storage
288 CALL destroy_tmp_data(fe_env, wrk, ng_points)
289 CALL timestop(handle)
290 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 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 ALLOCATE (fe_env%cg_data(ng_points))
309 DO i = 1, ng_points
310 ALLOCATE (fe_env%cg_data(i)%avg(ncolvar))
311 ALLOCATE (fe_env%cg_data(i)%var(ncolvar, ncolvar))
312 END DO
313 IF (PRESENT(wrk)) THEN
314 ALLOCATE (wrk(ng_points))
315 END IF
316 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 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 DO i = 1, ng_points
334 DEALLOCATE (fe_env%cg_data(i)%avg)
335 DEALLOCATE (fe_env%cg_data(i)%var)
336 END DO
337 DEALLOCATE (fe_env%cg_data)
338 IF (PRESENT(wrk)) THEN
339 DEALLOCATE (wrk)
340 END IF
341 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 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 DO i = 1, ng_points
358 istart = fe_env%nr_points - (i)*fe_env%conv_par%cg_width + 1
359 iend = fe_env%nr_points - (i - 1)*fe_env%conv_par%cg_width
360 IF (output_unit > 0) THEN
361 WRITE (output_unit, *) istart, iend
362 END IF
363 CALL eval_cov_matrix(fe_env, cg_index=i, istart=istart, iend=iend, output_unit=output_unit)
364 END DO
365
366 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 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 CALL timeset(routinen, handle)
388 test_passed = .false.
389 DO WHILE (fe_env%conv_par%cg_width < fe_env%conv_par%max_cg_width)
390 ng_points = nr_points/fe_env%conv_par%cg_width
391 print *, ng_points
392 IF (ng_points < min_sample_size) EXIT
393 CALL ui_check_norm_sc_low(fe_env, nr_points, output_unit)
394 test_passed = fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw
395 IF (test_passed) EXIT
396 fe_env%conv_par%cg_width = fe_env%conv_par%cg_width + 1
397 IF (output_unit > 0) THEN
398 WRITE (output_unit, *) "New coarse grained width:", fe_env%conv_par%cg_width
399 END IF
400 END DO
401 IF (fe_env%conv_par%cg_width == fe_env%conv_par%max_cg_width .AND. (.NOT. (test_passed))) THEN
402 CALL ui_check_norm_sc_low(fe_env, nr_points, output_unit)
403 test_passed = fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw
404 END IF
405 CALL timestop(handle)
406 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: wrk
427
428 CALL timeset(routinen, handle)
429 ncolvar = fe_env%ncolvar
430 ! Compute the Coarse Grained data set using a reverse cumulative strategy
431 fe_env%conv_par%test_sw = .false.
432 fe_env%conv_par%test_vn = .false.
433 ! Number of coarse grained points
434 avg_test_passed = .true.
435 sdv_test_passed = .true.
436 ng_points = nr_points/fe_env%conv_par%cg_width
437 CALL create_tmp_data(fe_env, wrk, ng_points, ncolvar)
438 CALL create_csg_data(fe_env, ng_points, output_unit)
439 ! Testing Averages
440 DO j = 1, ncolvar
441 DO i = 1, ng_points
442 wrk(i) = fe_env%cg_data(i)%avg(j)
443 END DO
444 ! Test of Shapiro - Wilks for normality
445 ! - Average
446 CALL sw_test(wrk, ng_points, w, pw)
447 print *, 1.0_dp - pw, fe_env%conv_par%sw_conf_lm
448 avg_test_passed = (1.0_dp - pw) <= fe_env%conv_par%sw_conf_lm
449 fe_env%conv_par%test_sw = avg_test_passed
450 IF (output_unit > 0) THEN
451 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 CALL vn_test(wrk, ng_points, r, u, prob)
456 print *, prob, fe_env%conv_par%vn_conf_lm
457 avg_test_passed = prob <= fe_env%conv_par%vn_conf_lm
458 fe_env%conv_par%test_vn = avg_test_passed
459 IF (output_unit > 0) THEN
460 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 IF (fe_env%conv_par%test_vn .AND. fe_env%conv_par%test_sw) THEN
465 ! Testing Standard Deviations
466 DO j = 1, ncolvar
467 DO k = j, ncolvar
468 DO i = 1, ng_points
469 wrk(i) = fe_env%cg_data(i)%var(j, k)
470 END DO
471 ! Test of Shapiro - Wilks for normality
472 ! - Standard Deviation
473 CALL sw_test(wrk, ng_points, w, pw)
474 print *, 1.0_dp - pw, fe_env%conv_par%sw_conf_lm
475 sdv_test_passed = (1.0_dp - pw) <= fe_env%conv_par%sw_conf_lm
476 fe_env%conv_par%test_sw = fe_env%conv_par%test_sw .AND. sdv_test_passed
477 IF (output_unit > 0) THEN
478 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 CALL vn_test(wrk, ng_points, r, u, prob)
483 print *, prob, fe_env%conv_par%vn_conf_lm
484 sdv_test_passed = prob <= fe_env%conv_par%vn_conf_lm
485 fe_env%conv_par%test_vn = fe_env%conv_par%test_vn .AND. sdv_test_passed
486 IF (output_unit > 0) THEN
487 WRITE (output_unit, *) "von Neumann serial correlation test (Std. Dev.)", sdv_test_passed
488 END IF
489 END DO
490 END DO
491 CALL destroy_tmp_data(fe_env, wrk, ng_points)
492 ELSE
493 CALL destroy_tmp_data(fe_env, wrk, ng_points)
494 END IF
495 CALL timestop(handle)
496 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: avg_std, avgmx
519 REAL(kind=dp), DIMENSION(:, :), POINTER :: cov_std, covmx
520
521 CALL timeset(routinen, handle)
522 converged = .false.
523 ncolvar = fe_env%ncolvar
524 NULLIFY (avgmx, avg_std, covmx, cov_std)
525 CALL ui_check_norm_sc(fe_env, test_passed, nr_points, output_unit)
526 IF (test_passed) THEN
527 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 CALL create_tmp_data(fe_env, ng_points=ng_points, ncolvar=ncolvar)
531 CALL create_csg_data(fe_env, ng_points, output_unit)
532 ALLOCATE (covmx(ncolvar, ncolvar))
533 ALLOCATE (avgmx(ncolvar))
534 ALLOCATE (cov_std(ncolvar*(ncolvar + 1)/2, ncolvar*(ncolvar + 1)/2))
535 ALLOCATE (avg_std(ncolvar))
536 covmx = 0.0_dp
537 avgmx = 0.0_dp
538 DO i = 1, ng_points
539 covmx = covmx + fe_env%cg_data(i)%var
540 avgmx = avgmx + fe_env%cg_data(i)%avg
541 END DO
542 covmx = covmx/real(ng_points, kind=dp)
543 avgmx = avgmx/real(ng_points, kind=dp)
544
545 ! Compute errors on average and standard deviation
546 CALL compute_avg_std_errors(fe_env, ncolvar, avgmx, covmx, avg_std, cov_std)
547 IF (output_unit > 0) THEN
548 WRITE (output_unit, *) "pippo", avgmx, covmx
549 WRITE (output_unit, *) "pippo", avg_std, cov_std
550 END IF
551 ! Convergence of the averages
552 max_error_avg = sqrt(maxval(abs(avg_std))/real(ng_points, kind=dp))/minval(avgmx)
553 max_error_std = sqrt(maxval(abs(cov_std))/real(ng_points, kind=dp))/minval(covmx)
554 IF (max_error_avg <= fe_env%conv_par%eps_conv .AND. &
555 max_error_std <= fe_env%conv_par%eps_conv) converged = .true.
556
557 IF (output_unit > 0) THEN
558 WRITE (output_unit, '(/,T2,"CG SAMPLING LENGTH = ",I7,20X,"REQUESTED ACCURACY = ",E12.6)') ng_points, &
559 fe_env%conv_par%eps_conv
560 WRITE (output_unit, '(T50,"PRESENT ACCURACY AVG= ",E12.6)') max_error_avg
561 WRITE (output_unit, '(T50,"PRESENT ACCURACY STD= ",E12.6)') max_error_std
562 WRITE (output_unit, '(T50,"CONVERGED FE-DER = ",L12)') converged
563
564 WRITE (output_unit, '(/,T33, "COVARIANCE MATRIX")')
565 WRITE (output_unit, '(T8,'//cp_to_string(ncolvar)//'(3X,I7,6X))') (ic, ic=1, ncolvar)
566 DO ic = 1, ncolvar
567 WRITE (output_unit, '(T2,I6,'//cp_to_string(ncolvar)//'(3X,E12.6,1X))') ic, covmx(ic, :)
568 END DO
569 WRITE (output_unit, '(T33, "ERROR OF COVARIANCE MATRIX")')
570 WRITE (output_unit, '(T8,'//cp_to_string(ncolvar)//'(3X,I7,6X))') (ic, ic=1, ncolvar)
571 DO ic = 1, ncolvar
572 WRITE (output_unit, '(T2,I6,'//cp_to_string(ncolvar)//'(3X,E12.6,1X))') ic, cov_std(ic, :)
573 END DO
574
575 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 (ic, avgmx(ic), sqrt(abs(avg_std(ic))), ic=1, ncolvar)
578 END IF
579 CALL destroy_tmp_data(fe_env, ng_points=ng_points)
580 DEALLOCATE (covmx)
581 DEALLOCATE (avgmx)
582 DEALLOCATE (cov_std)
583 DEALLOCATE (avg_std)
584 END IF
585 CALL timestop(handle)
586 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: awrk, eig, tmp
611 REAL(kind=dp), DIMENSION(:, :), POINTER :: wrk
612
613! Averages
614
615 nvar = ncolvar
616 ALLOCATE (wrk(nvar, nvar))
617 ALLOCATE (eig(nvar))
618 fac = real(SIZE(fe_env%cg_data), kind=dp)
619 wrk = 0.0_dp
620 eig = 0.0_dp
621 DO k = 1, SIZE(fe_env%cg_data)
622 DO j = 1, nvar
623 DO i = j, nvar
624 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 DO j = 1, nvar
629 DO i = j, nvar
630 wrk(i, j) = wrk(i, j) - avgmx(i)*avgmx(j)*fac
631 wrk(j, i) = wrk(i, j)
632 END DO
633 END DO
634 wrk = wrk/(fac - 1.0_dp)
635 ! Diagonalize the covariance matrix and check for the maximum error
636 CALL diamat_all(wrk, eig)
637 DO i = 1, nvar
638 avg_std(i) = eig(i)
639 END DO
640 DEALLOCATE (wrk)
641 DEALLOCATE (eig)
642 ! Standard Deviations
643 nvar = ncolvar*(ncolvar + 1)/2
644 ALLOCATE (wrk(nvar, nvar))
645 ALLOCATE (eig(nvar))
646 ALLOCATE (awrk(nvar))
647 ALLOCATE (tmp(nvar))
648 wrk = 0.0_dp
649 eig = 0.0_dp
650 ind = 0
651 DO i = 1, ncolvar
652 DO j = i, ncolvar
653 ind = ind + 1
654 awrk(ind) = covmx(i, j)
655 END DO
656 END DO
657 DO k = 1, SIZE(fe_env%cg_data)
658 ind = 0
659 DO i = 1, ncolvar
660 DO j = i, ncolvar
661 ind = ind + 1
662 tmp(ind) = fe_env%cg_data(k)%var(i, j)
663 END DO
664 END DO
665 DO i = 1, nvar
666 DO j = i, nvar
667 wrk(i, j) = wrk(i, j) + tmp(i)*tmp(j) - awrk(i)*awrk(j)
668 END DO
669 END DO
670 END DO
671 DO i = 1, nvar
672 DO j = i, nvar
673 wrk(i, j) = wrk(i, j) - fac*awrk(i)*awrk(j)
674 wrk(j, i) = wrk(i, j)
675 END DO
676 END DO
677 wrk = wrk/(fac - 1.0_dp)
678 ! Diagonalize the covariance matrix and check for the maximum error
679 CALL diamat_all(wrk, eig)
680 ind = 0
681 DO i = 1, ncolvar
682 DO j = i, ncolvar
683 ind = ind + 1
684 cov_std(i, j) = eig(ind)
685 cov_std(j, i) = cov_std(i, j)
686 END DO
687 END DO
688 DEALLOCATE (wrk)
689 DEALLOCATE (eig)
690 DEALLOCATE (awrk)
691 DEALLOCATE (tmp)
692
693 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 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 CALL timeset(routinen, handle)
720 ncolvar = fe_env%ncolvar
721 nlength = iend - istart + 1
722 fe_env%cg_data(cg_index)%avg = 0.0_dp
723 fe_env%cg_data(cg_index)%var = 0.0_dp
724 IF (nlength > 1) THEN
725 ! Update the info on averages and variances
726 DO jstep = istart, iend
727 DO ic = 1, ncolvar
728 cv => fe_env%uivar(ic)
729 tmp_ic = cv%ss(jstep)
730 fe_env%cg_data(cg_index)%avg(ic) = fe_env%cg_data(cg_index)%avg(ic) + tmp_ic
731 END DO
732 DO ic = 1, ncolvar
733 cv => fe_env%uivar(ic)
734 tmp_ic = cv%ss(jstep)
735 DO jc = 1, ic
736 cv => fe_env%uivar(jc)
737 tmp_jc = cv%ss(jstep)
738 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 fe_env%cg_data(cg_index)%var = fe_env%cg_data(cg_index)%var/real(nlength - 1, kind=dp)
745 fe_env%cg_data(cg_index)%avg = fe_env%cg_data(cg_index)%avg/real(nlength, kind=dp)
746 ! Compute the covariance matrix
747 DO ic = 1, ncolvar
748 tmp_ic = fe_env%cg_data(cg_index)%avg(ic)
749 DO jc = 1, ic
750 tmp_jc = fe_env%cg_data(cg_index)%avg(jc)*real(nlength, kind=dp)/real(nlength - 1, kind=dp)
751 fe_env%cg_data(cg_index)%var(jc, ic) = fe_env%cg_data(cg_index)%var(jc, ic) - tmp_ic*tmp_jc
752 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 IF (output_unit > 0) THEN
756 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 IF (PRESENT(covmx)) covmx = fe_env%cg_data(cg_index)%var
759 IF (PRESENT(avgs)) avgs = fe_env%cg_data(cg_index)%avg
760 END IF
761 CALL timestop(handle)
762 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 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 logger => cp_get_default_logger()
800 alch_section => section_vals_get_subs_vals(fe_section, "ALCHEMICAL_CHANGE")
801 CALL section_vals_val_get(alch_section, "PARAMETER", c_val=par)
802 DO i = 1, SIZE(my_par)
803 IF (my_par(i) == par) EXIT
804 END DO
805 cpassert(i <= SIZE(my_par))
806 ipar = i
807 dedf = evalfd(1, ipar, my_val, dx, err)
808 IF (abs(err) > lerr) THEN
809 WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
810 WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
811 CALL compress(this_error, .true.)
812 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 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 CALL section_vals_val_get(alch_section, "WEIGHTING_FUNCTION", c_val=coupling_function)
823 CALL compress(coupling_function, full=.true.)
824 CALL parsef(2, trim(coupling_function), my_par)
825 ene_w = evalf(2, my_val)
826 d_ene_w = evalfd(2, ipar, my_val, dx, err)
827 IF (abs(err) > lerr) THEN
828 WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
829 WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
830 CALL compress(this_error, .true.)
831 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 trim(def_error)//' .')
836 END IF
837 CALL section_vals_val_get(alch_section, "NEQUIL_STEPS", i_val=nequilstep)
838 ! Store results
839 IF (istep > nequilstep) THEN
840 isize = SIZE(cum_res, 2) + 1
841 CALL reallocate(cum_res, 1, 3, 1, isize)
842 cum_res(1, isize) = dedf
843 cum_res(2, isize) = dedf - d_ene_w
844 cum_res(3, isize) = ene_w
845 ! Compute derivative of biased and total energy
846 ! Total Free Energy
847 avg_det = sum(cum_res(1, 1:isize))/real(isize, kind=dp)
848 std_det = sum(cum_res(1, 1:isize)**2)/real(isize, kind=dp)
849 ! Unbiased Free Energy
850 avg_bp = sum(cum_res(3, 1:isize))/real(isize, kind=dp)
851 wfac = 0.0_dp
852 DO j = 1, isize
853 wfac = wfac + exp(beta*(cum_res(3, j) - avg_bp))
854 END DO
855 avg_due = 0.0_dp
856 std_due = 0.0_dp
857 DO j = 1, isize
858 tmp = cum_res(2, j)
859 tmp2 = exp(beta*(cum_res(3, j) - avg_bp))/wfac
860 avg_due = avg_due + tmp*tmp2
861 std_due = std_due + tmp**2*tmp2
862 END DO
863 IF (isize > 1) THEN
864 err_due = sqrt(std_due - avg_due**2)/sqrt(real(isize - 1, kind=dp))
865 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 extension=".free_energy")
870 IF (iw > 0) THEN
871 WRITE (iw, '(T2,79("-"),T37," oOo ")')
872 DO iforce_eval = 1, nforce_eval
873 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| FORCE_EVAL Nr.",I5,T48,"ENERGY (Hartree)= ",F15.9)') &
874 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 trim(par), dedf
878 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE OF BIASED ENERGY [ PARAMETER (",A,") ]",T66,F15.9)') &
879 trim(par), dedf - d_ene_w
880 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| BIASING UMBRELLA POTENTIAL ",T66,F15.9)') &
881 ene_w
882
883 IF (isize > 1) THEN
884 WRITE (iw, '(/,T2,"ALCHEMICAL CHANGE| DERIVATIVE TOTAL FREE ENERGY ",T50,F15.9,1X,"+/-",1X,F11.9)') &
885 avg_det, err_det
886 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE UNBIASED FREE ENERGY ",T50,F15.9,1X,"+/-",1X,F11.9)') &
887 avg_due, err_due
888 ELSE
889 WRITE (iw, '(/,T2,"ALCHEMICAL CHANGE| DERIVATIVE TOTAL FREE ENERGY ",T50,F15.9,1X,"+/-",1X,T76,A)') &
890 avg_det, "UNDEF"
891 WRITE (iw, '(T2,"ALCHEMICAL CHANGE| DERIVATIVE UNBIASED FREE ENERGY ",T50,F15.9,1X,"+/-",1X,T76,A)') &
892 avg_due, "UNDEF"
893 END IF
894 WRITE (iw, '(T2,79("-"))')
895 END IF
896 END IF
897 CALL cp_print_key_finished_output(iw, logger, fe_section, "FREE_ENERGY_INFO")
898
899 END SUBROUTINE dump_ac_info
900
901END MODULE free_energy_methods
902
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:48
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:1726
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:515
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:803
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:1013
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)
returns various attributes about the force environment
Definition force_env_types.F:329
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:148
fparser::evalf
real(rn) function, public evalf(i, val)
...
Definition fparser.F:180
fparser::evalfd
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
Definition fparser.F:976
fparser::finalizef
subroutine, public finalizef()
...
Definition fparser.F:101
fparser::initf
subroutine, public initf(n)
...
Definition fparser.F:130
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:17
input_constants::do_fe_ac
integer, parameter, public do_fe_ac
Definition input_constants.F:362
input_constants::do_fe_ui
integer, parameter, public do_fe_ui
Definition input_constants.F:362
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:724
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:1040
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:372
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:3190
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:139
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:126
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