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