LCOV - code coverage report
Current view: top level - src/motion - free_energy_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 52.4 % 416 218
Test Date: 2025-07-25 12:55:17 Functions: 69.2 % 13 9

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

Generated by: LCOV version 2.0-1