LCOV - code coverage report
Current view: top level - src/motion - free_energy_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 218 416 52.4 %
Date: 2024-04-20 06:29:22 Functions: 9 13 69.2 %

          Line data    Source code
       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       86270 :    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       43135 :          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       43135 :       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       43135 :       NULLIFY (force_env, istep, subsys, cv, simpar)
      90       86270 :       logger => cp_get_default_logger()
      91       43135 :       CALL timeset(routineN, handle)
      92       43135 :       converged = .FALSE.
      93             :       CALL get_md_env(md_env, force_env=force_env, fe_env=fe_env, simpar=simpar, &
      94       43135 :                       itimes=istep)
      95             :       ! Metadynamics is also a free energy calculation but is handled in a different
      96             :       ! module.
      97       43135 :       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       43135 :       CALL timestop(handle)
     159             : 
     160       43135 :    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          12 :       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 1.15