(git:b195825)
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 ! **************************************************************************************************
12 MODULE free_energy_methods
13  USE colvar_methods, ONLY: colvar_eval_glob_f
14  USE cp_log_handling, ONLY: cp_get_default_logger,&
15  cp_logger_get_default_io_unit,&
16  cp_logger_type,&
17  cp_to_string
18  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
19  cp_print_key_unit_nr
20  USE cp_subsys_types, ONLY: cp_subsys_type
21  USE force_env_types, ONLY: force_env_get,&
22  force_env_type
23  USE fparser, ONLY: evalf,&
24  evalfd,&
25  finalizef,&
26  initf,&
27  parsef
28  USE free_energy_types, ONLY: free_energy_type,&
29  ui_var_type
30  USE input_constants, ONLY: do_fe_ac,&
31  do_fe_ui
32  USE input_section_types, ONLY: section_vals_get_subs_vals,&
33  section_vals_type,&
34  section_vals_val_get
35  USE kinds, ONLY: default_path_length,&
36  default_string_length,&
37  dp
38  USE mathlib, ONLY: diamat_all
39  USE md_environment_types, ONLY: get_md_env,&
40  md_environment_type
41  USE memory_utilities, ONLY: reallocate
42  USE simpar_types, ONLY: simpar_type
43  USE statistical_methods, ONLY: k_test,&
44  min_sample_size,&
45  sw_test,&
46  vn_test
47  USE string_utilities, ONLY: compress
48 #include "../base/base_uses.f90"
49 
50  IMPLICIT NONE
51 
52  PRIVATE
53  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'free_energy_methods'
54  PUBLIC :: free_energy_evaluate
55 
56 CONTAINS
57 
58 ! **************************************************************************************************
59 !> \brief Main driver for free energy calculations
60 !> In this routine we handle specifically biased MD.
61 !> \param md_env ...
62 !> \param converged ...
63 !> \param fe_section ...
64 !> \par History
65 !> Teodoro Laino (01.2007) [tlaino]
66 ! **************************************************************************************************
67  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 
901 END 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
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::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
fparser::parsef
subroutine, public parsef(i, FuncStr, Var)
Parse ith function string FuncStr and compile it into bytecode.
Definition: fparser.F:148
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:376
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