tmc_tree_acceptance.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief tree nodes acceptance
!>        code is separated in 3 parts, first the acceptance criteria,
!>        second the tree node acceptance handling, searching etc. and
!>        than the acceptance probability handling
!> \par History
!>      11.2012 created [Mandes Schoenherr]
!> \author Mandes
! **************************************************************************************************
 
MODULE tmc_tree_acceptance
   USE cp_log_handling,                 ONLY: cp_to_string
   USE kinds,                           ONLY: dp
   USE physcon,                         ONLY: boltzmann,&
                                              joule
   USE tmc_calculations,                ONLY: compute_estimated_prob,&
                                              get_scaled_cell
   USE tmc_dot_tree,                    ONLY: create_dot_color,&
                                              create_global_tree_dot_color
   USE tmc_file_io,                     ONLY: write_result_list_element
   USE tmc_move_handle,                 ONLY: prob_update
   USE tmc_move_types,                  ONLY: mv_type_md,&
                                              mv_type_volume_move
   USE tmc_stati,                       ONLY: task_type_gaussian_adaptation
   USE tmc_tree_build,                  ONLY: remove_unused_g_tree
   USE tmc_tree_search,                 ONLY: get_subtree_elements_to_check,&
                                              search_canceling_elements,&
                                              search_next_gt_element_to_check,&
                                              search_parent_element
   USE tmc_tree_types,                  ONLY: &
        add_to_list, global_tree_type, gt_elem_list_type, status_accepted, status_accepted_result, &
        status_calc_approx_ener, status_calculate_md, status_calculate_nmc_steps, &
        status_calculate_energy, status_calculated, status_cancel_ener, status_cancel_nmc, &
        status_canceled_ener, status_canceled_nmc, status_created, status_deleted, &
        status_deleted_result, status_rejected, status_rejected_result, tree_type
   USE tmc_types,                       ONLY: tmc_env_type,&
                                              tmc_param_type
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_tree_acceptance'
 
   PUBLIC :: acceptance_check
   PUBLIC :: check_acceptance_of_depending_subtree_nodes, &
             check_elements_for_acc_prob_update
   PUBLIC :: tree_update
 
   INTEGER, PARAMETER :: DEBUG = 0
 
CONTAINS
 
   !============================================================================
   ! acceptance criteria calculations
   !============================================================================
! **************************************************************************************************
!> \brief standard Monte Carlo and 2 potential acceptance check
!>        acceptance check of move from old(last accepted) to new configuration
!>        the sum of kinetic and potential energy is used
!>        acc(o->n)=min(1,exp( -beta*(H(n)-H(o)) ))
!> \param tree_element new/actual configuration
!> \param parent_element last accepted configuration
!> \param tmc_params TMC global parameters
!> \param temperature actual temperature configuration should be checked with
!> \param diff_pot_check 2potential check or not?
!> \param accept result (configuration accepted of rejected)
!> \param rnd_nr random number for acceptance check
!> \param approx_ener for NMC the approximated energies schould be used
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE acceptance_check(tree_element, parent_element, tmc_params, &
                               temperature, diff_pot_check, accept, rnd_nr, &
                               approx_ener)
      TYPE(tree_type), POINTER                           :: tree_element, parent_element
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      REAL(kind=dp)                                      :: temperature
      LOGICAL                                            :: diff_pot_check, accept
      REAL(kind=dp)                                      :: rnd_nr
      LOGICAL, OPTIONAL                                  :: approx_ener
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'acceptance_check'
 
      INTEGER                                            :: handle
      REAL(kind=dp)                                      :: ekin_last_acc, elem_ener, kb, parent_ener
 
      cpassert(ASSOCIATED(tree_element))
      cpassert(ASSOCIATED(parent_element))
      cpassert(ASSOCIATED(tmc_params))
      cpassert(temperature > 0.0_dp)
      cpassert(rnd_nr >= 0.0_dp .AND. rnd_nr <= 1.0_dp)
 
      kb = boltzmann/joule
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      IF (tmc_params%task_type == task_type_gaussian_adaptation) THEN
         cpabort("task_type_gaussian_adaptation not yet implemented")
!TODO       CALL acc_check_gauss_adapt(f=tree_element%potential, ct=tree_element%ekin_before_md, acc=accept)
!DO NOT DO       RETURN
      END IF
 
      !-- using 2 different potentials, the energy difference between these potentials
      !   and the two configurations have to be regarded.
      !   The ensamble should be in equilibrium state of the approximate potential
      IF (diff_pot_check .AND. (tmc_params%NMC_inp_file /= "")) THEN
         IF ((tree_element%potential == huge(tree_element%potential)) .OR. &
             tree_element%e_pot_approx == huge(tree_element%e_pot_approx)) THEN
            elem_ener = huge(elem_ener)
         ELSE
            !for different potentials we have to regard the differences in energy
            ! min(1,e^{-\beta*[(E_{exact}(n)-E_{approx}(n))-(E_{exact}(o)-E_{approx}(o))]})
            elem_ener = 1.0_dp/(kb*temperature)*tree_element%potential &
                        - 1.0_dp/(kb*tmc_params%Temp(tree_element%temp_created)) &
                        *tree_element%e_pot_approx
         END IF
         parent_ener = 1.0_dp/(kb*temperature)*parent_element%potential &
                       - 1.0_dp/(kb*tmc_params%Temp(tree_element%temp_created)) &
                       *parent_element%e_pot_approx
 
         !-- always accepted if new energy is smaller than old energy
         IF (elem_ener <= parent_ener) THEN
            accept = .true.
         ELSE
            !-- gaussian distributed acceptance if new energy is greater than old energy
            IF (rnd_nr < &
                exp(-(elem_ener - parent_ener))) THEN
               accept = .true.
            ELSE
               accept = .false.
            END IF
         END IF
      ELSE
         !-- standard MC check, but also using the kinetic energy for Hybrid Monte Carlo
         IF (tree_element%move_type == mv_type_md) THEN
            ekin_last_acc = tree_element%ekin_before_md
            ! using the kinetic energy before velocity change
         ELSE
            ekin_last_acc = parent_element%ekin
         END IF
         ! comparing aproximated energies
         IF (PRESENT(approx_ener)) THEN
            elem_ener = tree_element%e_pot_approx &
                        + tree_element%ekin
            parent_ener = parent_element%e_pot_approx &
                          + ekin_last_acc
         ELSE
            elem_ener = tree_element%potential &
                        + tree_element%ekin
            parent_ener = parent_element%potential &
                          + ekin_last_acc
         END IF
 
         !-- always accepted if new energy is smaller than old energy
         IF (elem_ener <= parent_ener) THEN
            accept = .true.
         ELSE
            !-- gaussian distributed acceptance if new energy is greater than old energy
            IF (rnd_nr < &
                exp(-1.0_dp/(kb*temperature)*(elem_ener - parent_ener))) THEN
               accept = .true.
            ELSE
               accept = .false.
            END IF
         END IF
      END IF
 
      ! update the estimated energy acceptance probability distribution
      IF (diff_pot_check) THEN
         cpassert(ASSOCIATED(tmc_params%prior_NMC_acc))
         tmc_params%prior_NMC_acc%counter = tmc_params%prior_NMC_acc%counter + 1
         tmc_params%prior_NMC_acc%aver = (tmc_params%prior_NMC_acc%aver*(tmc_params%prior_NMC_acc%counter - 1) + &
                                          ((elem_ener - parent_ener)))/real(tmc_params%prior_NMC_acc%counter, kind=dp)
         tmc_params%prior_NMC_acc%aver_2 = (tmc_params%prior_NMC_acc%aver_2*(tmc_params%prior_NMC_acc%counter - 1) + &
                                            (elem_ener - parent_ener)**2)/real(tmc_params%prior_NMC_acc%counter, kind=dp)
      END IF
 
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE acceptance_check
 
! **************************************************************************************************
!> \brief parallel tempering swap acceptance check,
!>        check swapped configurations of different temperatures
!>        using sum of potential and kinetic energy
!>        always accepted if energy of configuration with higher temperature
!>        is smaller than energy of configuration with lower temperature
!>        acc(o->n)=min(1,exp((beta_i-beta_j)*(U_i-U_j)))
!> \param tree_elem actual global tree element
!> \param conf1 sub tree element of lower temperature
!> \param conf2 sub tree element of higher temperature
!> \param tmc_params TMC environment parameters
!> \param accept acceptance of configurational change
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE swap_acceptance_check(tree_elem, conf1, conf2, tmc_params, accept)
      TYPE(global_tree_type), POINTER                    :: tree_elem
      TYPE(tree_type), POINTER                           :: conf1, conf2
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      LOGICAL                                            :: accept
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'swap_acceptance_check'
 
      INTEGER                                            :: handle
      REAL(kind=dp)                                      :: delta, elem1_ener, elem2_ener, kb, vol1, &
                                                            vol2
 
      kb = boltzmann/joule
 
      cpassert(ASSOCIATED(tree_elem))
      cpassert(tree_elem%rnd_nr >= 0.0_dp)
      cpassert(ASSOCIATED(conf1))
      cpassert(ASSOCIATED(conf2))
      cpassert(ASSOCIATED(tmc_params))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      IF (tmc_params%pressure > 0.0_dp) THEN
         ! pt-NVT
         elem1_ener = conf1%potential &
                      + conf1%ekin
         elem2_ener = conf2%potential &
                      + conf2%ekin
         ! the swap is done with prob: exp((\beta_i-\beta_j)(U_i-U_j)),
         ! BUT because they are already swaped we exchange the energies.
         IF (tree_elem%rnd_nr < exp((1.0_dp/(kb*tmc_params%Temp(tree_elem%mv_conf)) - &
                                     1.0_dp/(kb*tmc_params%Temp(tree_elem%mv_conf + 1))) &
                                    *(elem2_ener - elem1_ener))) THEN
            accept = .true.
         ELSE
            accept = .false.
         END IF
      ELSE
         ! pt-NpT (parallel Tempering with constant pressure, temperature and num. particle)
         CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf1%box_scale, &
                              vol=vol1)
         CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf2%box_scale, &
                              vol=vol2)
         ! delta= (beta_m-beta_n)(U_n-U_m) + (beta_m*P_m-beta_n*P_n)(V_n-V_m)
         delta = (1.0_dp/(kb*tmc_params%Temp(tree_elem%mv_conf)) &
                  - 1.0_dp/(kb*tmc_params%Temp(tree_elem%mv_conf + 1)))* &
                 ((conf2%potential + conf2%ekin) - (conf1%potential + conf1%ekin)) + &
                 (1.0_dp/(kb*tmc_params%Temp(tree_elem%mv_conf))*tmc_params%pressure &
                  - 1.0_dp/(kb*tmc_params%Temp(tree_elem%mv_conf + 1))*tmc_params%pressure)* &
                 (vol2 - vol1)
         IF (tree_elem%rnd_nr < exp(delta)) THEN
            accept = .true.
         ELSE
            accept = .false.
         END IF
      END IF
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE swap_acceptance_check
 
! **************************************************************************************************
!> \brief volume move acceptance check
!>        sampling NPT, we need to do volume moves. Their acceptance are
!>        checked regarding the old and new energy, the volume difference
!>        and the constant pressure
!> \param parent_elem old tree element with old box size
!> \param new_elem actual tree element with new box size
!> \param temperature actual temperature
!> \param rnd_nr random number to check with
!> \param tmc_params TMC environment parameters
!> \param accept Monte Carlo move acceptance (result)
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE volume_acceptance_check(parent_elem, new_elem, temperature, &
                                      rnd_nr, tmc_params, accept)
      TYPE(tree_type), POINTER                           :: parent_elem, new_elem
      REAL(kind=dp)                                      :: temperature, rnd_nr
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      LOGICAL                                            :: accept
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'volume_acceptance_check'
 
      INTEGER                                            :: handle
      REAL(kind=dp)                                      :: d_enthalpy, kb, n_vol, p_vol
 
      kb = boltzmann/joule
 
      cpassert(ASSOCIATED(parent_elem))
      cpassert(ASSOCIATED(new_elem))
      cpassert(temperature > 0.0_dp)
      cpassert(rnd_nr > 0.0_dp)
      cpassert(ASSOCIATED(tmc_params))
      cpassert(tmc_params%pressure >= 0.0_dp)
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      CALL get_scaled_cell(cell=tmc_params%cell, box_scale=parent_elem%box_scale, &
                           vol=p_vol)
      CALL get_scaled_cell(cell=tmc_params%cell, box_scale=new_elem%box_scale, &
                           vol=n_vol)
      ! H=enthalpy, U=energy, P=pressure, V=volume, kB=Boltzmann constant, T=temperature, N=amount particle
      ! delta_H  =      delta_U                               + P*delta_V           - kB*T*N*ln(V_n/V_p)
      IF (.false.) THEN
         ! the volume move in volume space (dV)
         d_enthalpy = (new_elem%potential - parent_elem%potential) + &
                      tmc_params%pressure*(n_vol - p_vol) - &
                      kb*temperature*(SIZE(new_elem%pos)/ &
                                      tmc_params%dim_per_elem)* &
                      log(n_vol/p_vol)
      ELSE
         IF (tmc_params%v_isotropic) THEN
            d_enthalpy = (new_elem%potential - parent_elem%potential) + &
                         tmc_params%pressure*(n_vol - p_vol) - &
                         kb*temperature*((SIZE(new_elem%pos)/ &
                                          tmc_params%dim_per_elem) + 2/real(3, kind=dp))* &
                         log(n_vol/p_vol)
         ELSE
            d_enthalpy = (new_elem%potential - parent_elem%potential) + &
                         tmc_params%pressure*(n_vol - p_vol) - &
                         kb*temperature*(SIZE(new_elem%pos)/ &
                                         tmc_params%dim_per_elem)* &
                         log(n_vol/p_vol)
         END IF
      END IF
      ! acc(o->n) = min(1, exp(-beta*delta_H))
      IF (d_enthalpy <= 0.0_dp) THEN
         accept = .true.
      ELSE
         IF (rnd_nr < exp(-1.0_dp/(kb*temperature)*d_enthalpy)) THEN
            accept = .true.
         ELSE
            accept = .false.
         END IF
      END IF
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE volume_acceptance_check
 
   !============================================================================
   ! tree node acceptance
   !============================================================================
! **************************************************************************************************
!> \brief check acceptance of energy calculated element and related childs,
!>        when ready
!> \param tree_elem actual tree element with calculated energy
!> \param tmc_env TMC environment parameters
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE check_acceptance_of_depending_subtree_nodes(tree_elem, tmc_env)
      TYPE(tree_type), POINTER                           :: tree_elem
      TYPE(tmc_env_type), POINTER                        :: tmc_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'check_acceptance_of_depending_subtree_nodes'
 
      INTEGER                                            :: handle
      LOGICAL                                            :: parent_ready
      TYPE(tree_type), POINTER                           :: act_elem, parent_elem
 
      NULLIFY (parent_elem, act_elem)
 
      cpassert(ASSOCIATED(tmc_env))
      cpassert(ASSOCIATED(tree_elem))
      cpassert(tree_elem%stat == status_calculated)
      cpassert(ASSOCIATED(tree_elem%parent))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      act_elem => tree_elem
 
      ! ------------------------------------------------------
      ! check this element
      parent_elem => search_parent_element(act_elem)
      cpassert(.NOT. ASSOCIATED(act_elem, parent_elem))
 
      ! check status of parent element
      SELECT CASE (parent_elem%stat)
      CASE (status_created, status_calculate_energy, status_calculate_md, &
            status_calculate_nmc_steps, status_canceled_ener, &
            status_canceled_nmc, status_cancel_nmc, status_cancel_ener)
         parent_ready = .false.
      CASE (status_accepted_result, status_rejected_result, &
            status_accepted, status_rejected, status_calculated)
         parent_ready = .true.
      CASE DEFAULT
         cpabort("unknown parent element status"//cp_to_string(parent_elem%stat))
      END SELECT
 
      ! ready ?
      IF (parent_ready) THEN
         ! acceptance check
         CALL check_and_change_status_of_subtree_elem(act_elem=act_elem, &
                                                      parent_elem=parent_elem, tmc_env=tmc_env)
      END IF
      !------------------------------------------------------
      ! check acc child
      parent_elem => tree_elem
      IF (ASSOCIATED(tree_elem%acc)) THEN
         act_elem => tree_elem%acc
         IF (act_elem%stat == status_calculated) THEN
            CALL check_and_change_status_of_subtree_elem(act_elem=act_elem, &
                                                         parent_elem=parent_elem, tmc_env=tmc_env)
         END IF
         !------------------------------------------------------
         ! check all nacc childs
         nacc_loop: DO
            IF (.NOT. ASSOCIATED(act_elem%nacc)) EXIT nacc_loop
            act_elem => act_elem%nacc
            IF (act_elem%stat == status_calculated) THEN
               CALL check_and_change_status_of_subtree_elem(act_elem=act_elem, &
                                                            parent_elem=parent_elem, tmc_env=tmc_env)
            END IF
         END DO nacc_loop
      END IF
 
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE check_acceptance_of_depending_subtree_nodes
 
! **************************************************************************************************
!> \brief checking the elements and change the status
!> \param act_elem actual tree element (new configuration)
!> \param parent_elem parent tree element (old configuration)
!> \param tmc_env TMC environment parameters
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE check_and_change_status_of_subtree_elem(act_elem, parent_elem, tmc_env)
      TYPE(tree_type), POINTER                           :: act_elem, parent_elem
      TYPE(tmc_env_type), POINTER                        :: tmc_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'check_and_change_status_of_subtree_elem'
 
      INTEGER                                            :: handle
      LOGICAL                                            :: flag, result_acc
      TYPE(global_tree_type), POINTER                    :: tmp_gt_ptr
      TYPE(gt_elem_list_type), POINTER                   :: tmp_gt_list_ptr
 
      NULLIFY (tmp_gt_list_ptr, tmp_gt_ptr)
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      cpassert(ASSOCIATED(tmc_env))
      cpassert(ASSOCIATED(tmc_env%params))
      cpassert(ASSOCIATED(act_elem))
      cpassert(ASSOCIATED(parent_elem))
      cpassert(act_elem%stat == status_calculated)
      mark_used(parent_elem)
 
      flag = .false.
 
      tmp_gt_list_ptr => act_elem%gt_nodes_references
      ! check for all global tree elements referring to this subtree element
      ! same subtree element can be accepted at a certain temperature and rejected at another one
      DO WHILE (ASSOCIATED(tmp_gt_list_ptr))
         cpassert(tmp_gt_list_ptr%gt_elem%stat /= status_accepted_result)
         cpassert(tmp_gt_list_ptr%gt_elem%stat /= status_rejected_result)
 
         CALL check_elements(gt_act_elem=tmp_gt_list_ptr%gt_elem, tmc_env=tmc_env, &
                             check_done=flag, result_acc=result_acc)
         IF (flag) THEN
            ! probability update
            CALL prob_update(move_types=tmc_env%params%move_types, elem=act_elem, &
                             acc=result_acc, prob_opt=tmc_env%params%esimate_acc_prob)
 
            ! change status
            ! accepted case: delete (and cancel) not accepted branch
            IF (result_acc .EQV. .true.) THEN
               tmp_gt_list_ptr%gt_elem%stat = status_accepted
               IF (ASSOCIATED(tmp_gt_list_ptr%gt_elem%nacc)) THEN
                  tmp_gt_ptr => tmp_gt_list_ptr%gt_elem%nacc
                  CALL remove_unused_g_tree(begin_ptr=tmp_gt_ptr, &
                                            end_ptr=tmp_gt_list_ptr%gt_elem, removed=flag, &
                                            tmc_env=tmc_env)
               END IF
            ELSE
               ! not accepted case: delete (and cancel) accepted branch
               tmp_gt_list_ptr%gt_elem%stat = status_rejected
               IF (ASSOCIATED(tmp_gt_list_ptr%gt_elem%acc)) THEN
                  tmp_gt_ptr => tmp_gt_list_ptr%gt_elem%acc
                  CALL remove_unused_g_tree(begin_ptr=tmp_gt_ptr, &
                                            end_ptr=tmp_gt_list_ptr%gt_elem, removed=flag, &
                                            tmc_env=tmc_env)
               END IF
            END IF
 
            IF (tmc_env%params%DRAW_TREE) THEN
               CALL create_global_tree_dot_color(gt_tree_element=tmp_gt_list_ptr%gt_elem, &
                                                 tmc_params=tmc_env%params)
            END IF
         END IF
         tmp_gt_list_ptr => tmp_gt_list_ptr%next
      END DO
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE check_and_change_status_of_subtree_elem
 
! **************************************************************************************************
!> \brief change status flag of all subtree elements
!>        most important to draw the right subtrees
!> \param gt_ptr global tree pointer, the related/changed sub tree element
!>        should change status
!> \param stat the status of the global tree element
!> \param tmc_params TMC emvironment parameters for drawing
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE subtree_configuration_stat_change(gt_ptr, stat, tmc_params)
      TYPE(global_tree_type), POINTER                    :: gt_ptr
      INTEGER                                            :: stat
      TYPE(tmc_param_type), POINTER                      :: tmc_params
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'subtree_configuration_stat_change'
 
      INTEGER                                            :: handle
      TYPE(tree_type), POINTER                           :: ptr
 
      NULLIFY (ptr)
 
      cpassert(ASSOCIATED(gt_ptr))
      cpassert(ASSOCIATED(tmc_params))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      ! don't respect swaped nodes (subtree nodes are already in the list,
      !   and we dont want them in marked in subtrees)
      IF (.NOT. gt_ptr%swaped) THEN
         ptr => gt_ptr%conf(gt_ptr%mv_conf)%elem
         ! dependent on the status (don't map each status 1:1,
         !   e.g. not the result ones)
         SELECT CASE (stat)
         CASE (status_accepted_result)
            ptr%stat = status_accepted
         CASE (status_rejected_result)
            ptr%stat = status_rejected
         CASE (status_accepted, status_rejected)
            ptr%stat = stat
         CASE DEFAULT
            CALL cp_abort(__location__, &
                          "unknown global tree status"// &
                          cp_to_string(stat))
         END SELECT
 
         IF (tmc_params%DRAW_TREE) THEN
            CALL create_dot_color(tree_element=ptr, tmc_params=tmc_params)
         END IF
      END IF
 
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE subtree_configuration_stat_change
 
! **************************************************************************************************
!> \brief checks if the element is ready for an acceptance check
!> \param elem sub tree element
!> \param ready return value
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE elem_ready_to_check(elem, ready)
      TYPE(tree_type), POINTER                           :: elem
      LOGICAL                                            :: ready
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'elem_ready_to_check'
 
      INTEGER                                            :: handle
 
      cpassert(ASSOCIATED(elem))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      ready = .false.
 
      SELECT CASE (elem%stat)
      CASE (status_created, status_calculate_energy, &
            status_calculate_md, status_calculate_nmc_steps, status_calc_approx_ener, &
            status_cancel_nmc, status_cancel_ener, status_canceled_nmc, status_canceled_ener)
         ready = .false.
      CASE (status_calculated, status_accepted_result, status_accepted, status_rejected, status_rejected_result)
         ready = .true.
      CASE DEFAULT
         CALL cp_abort(__location__, &
                       "status of actual tree node is unknown" &
                       //cp_to_string(elem%stat))
      END SELECT
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE elem_ready_to_check
 
! **************************************************************************************************
!> \brief checking the elements (find element and do acceptance check)
!> \param gt_act_elem actual global tree element
!> \param tmc_env TMC environment
!> \param check_done successful checked? (result)
!> \param result_acc checked configuration accepted? (result)
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE check_elements(gt_act_elem, tmc_env, check_done, result_acc)
      TYPE(global_tree_type), POINTER                    :: gt_act_elem
      TYPE(tmc_env_type), POINTER                        :: tmc_env
      LOGICAL                                            :: check_done, result_acc
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'check_elements'
 
      INTEGER                                            :: handle
      LOGICAL                                            :: act_ready, tmp_ready
      TYPE(tree_type), POINTER                           :: act_element, tmp_element
 
      NULLIFY (act_element, tmp_element)
      check_done = .false.
 
      cpassert(ASSOCIATED(gt_act_elem))
      cpassert(ASSOCIATED(tmc_env))
      cpassert(ASSOCIATED(tmc_env%m_env))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      !====================================================================================
      ! check if global tree element is already checked
      SELECT CASE (gt_act_elem%stat)
      CASE (status_accepted, status_rejected)
         check_done = .true.
         IF (gt_act_elem%stat == status_accepted) THEN
            result_acc = .true.
         ELSE IF (gt_act_elem%stat == status_rejected) THEN
            result_acc = .false.
         ELSE
            CALL cp_abort(__location__, &
                          "undefinite status of already checked elements:" &
                          //cp_to_string(gt_act_elem%stat))
         END IF
      CASE DEFAULT
      END SELECT
 
      IF (.NOT. check_done) THEN
         !====================================================================================
         ! get elements related to global tree element to check
         CALL get_subtree_elements_to_check(gt_act_elem=gt_act_elem, elem1=act_element, &
                                            elem2=tmp_element)
         cpassert(ASSOCIATED(act_element))
         cpassert(ASSOCIATED(tmp_element))
 
         ! check status of both elements (if they are already calculated and hence ready to check)
         CALL elem_ready_to_check(elem=act_element, ready=act_ready)
         CALL elem_ready_to_check(elem=tmp_element, ready=tmp_ready)
 
         ! both ready ? check
         IF (act_ready .AND. tmp_ready) THEN
            ! 2 temperature (swap) check
            IF (gt_act_elem%swaped) THEN
               CALL swap_acceptance_check(tree_elem=gt_act_elem, conf1=act_element, &
                                          conf2=tmp_element, tmc_params=tmc_env%params, &
                                          accept=result_acc)
               ! volume move check
            ELSE IF (act_element%move_type == mv_type_volume_move) THEN
               CALL volume_acceptance_check(parent_elem=tmp_element, new_elem=act_element, &
                                            temperature=tmc_env%params%Temp(gt_act_elem%mv_conf), &
                                            rnd_nr=gt_act_elem%rnd_nr, &
                                            tmc_params=tmc_env%params, accept=result_acc)
            ELSE
               IF (tmc_env%m_env%temp_decrease /= 1.0_dp) THEN
                  CALL acceptance_check(tree_element=act_element, parent_element=tmp_element, &
                                        tmc_params=tmc_env%params, temperature=gt_act_elem%Temp, &
                                        diff_pot_check=.true., accept=result_acc, &
                                        rnd_nr=gt_act_elem%rnd_nr)
               ELSE
                  CALL acceptance_check(tree_element=act_element, parent_element=tmp_element, &
                                        tmc_params=tmc_env%params, &
                                        temperature=tmc_env%params%Temp(gt_act_elem%mv_conf), &
                                        diff_pot_check=.true., accept=result_acc, &
                                        rnd_nr=gt_act_elem%rnd_nr)
               END IF
            END IF
            check_done = .true.
         ELSE
            ! if element is not ready, update status of global tree element
            !TODO maybe update this part (how much status we need in global tree from sub tree??)
            SELECT CASE (gt_act_elem%stat)
               !-- check if at least the MD move is already done
               !-- hence new configurations can be created on basis of this configuration
            CASE (status_calculate_md, status_calculate_nmc_steps, status_calculate_energy, &
                  status_created, status_calc_approx_ener)
               IF (gt_act_elem%conf(gt_act_elem%mv_conf)%elem%stat /= gt_act_elem%stat) THEN
                  gt_act_elem%stat = gt_act_elem%conf(gt_act_elem%mv_conf)%elem%stat
               END IF
            CASE (status_calculated)
            CASE DEFAULT
               CALL cp_abort(__location__, &
                             "status of actual checked node is unknown" &
                             //cp_to_string(gt_act_elem%stat))
            END SELECT
         END IF
      END IF
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE check_elements
 
! **************************************************************************************************
!> \brief searching tree nodes to check for Markov Chain, elements are marked
!>        and stored in lists ... (main entry point)
!> \param tmc_env TMC environment
!> \param result_acc checked configuration accepted? (result)
!> \param something_updated ...
!> \param
!> \param
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE tree_update(tmc_env, result_acc, something_updated)
      TYPE(tmc_env_type), POINTER                        :: tmc_env
      LOGICAL                                            :: result_acc, something_updated
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'tree_update'
 
      INTEGER                                            :: handle, itmp
      LOGICAL                                            :: found
      TYPE(global_tree_type), POINTER                    :: gt_act_elem
      TYPE(tree_type), POINTER                           :: act_element, tmp_element
 
      NULLIFY (gt_act_elem, tmp_element, act_element)
 
      cpassert(ASSOCIATED(tmc_env))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      result_acc = .false.
      something_updated = .false.
 
      gt_act_elem => tmc_env%m_env%gt_act
      search_calculated_element_loop: DO
         NULLIFY (act_element, tmp_element)
         ! search for element to check
         CALL search_next_gt_element_to_check(ptr=gt_act_elem, found=found)
         IF (.NOT. found) EXIT search_calculated_element_loop
 
         ! check the elements status end, if possible, do acceptance check
         CALL check_elements(gt_act_elem=gt_act_elem, tmc_env=tmc_env, &
                             check_done=found, result_acc=result_acc)
         ! if check was not possible, update the status of the global tree element and return
         IF (.NOT. found) EXIT search_calculated_element_loop
 
         ! get elements related to global tree element, which were checked
         CALL get_subtree_elements_to_check(gt_act_elem=gt_act_elem, &
                                            elem1=act_element, elem2=tmp_element)
 
         !========================================================================
         !-- set result counters
         ! counter for certain temperature
         tmc_env%m_env%result_count(gt_act_elem%mv_conf) = &
            tmc_env%m_env%result_count(gt_act_elem%mv_conf) + 1
         ! in case of swapped also count the result for
         !  the other swapped temperature
         IF (gt_act_elem%swaped) THEN
            tmc_env%m_env%result_count(gt_act_elem%mv_conf + 1) = &
               tmc_env%m_env%result_count(gt_act_elem%mv_conf + 1) + 1
         END IF
         ! count also for global tree Markov Chain
         tmc_env%m_env%result_count(0) = tmc_env%m_env%result_count(0) + 1
 
         ! flag for doing tree cleaning with canceling and deallocation
         something_updated = .true.
 
         !IF((.NOT.gt_act_elem%swaped) .AND. act_element%temp_created/=0 .AND. &
         !   act_element%temp_created/=gt_act_elem%mv_conf)&
         !  count_wrong_temp_move(gt_act_elem%mv_conf) = &
         !       count_wrong_temp_move(gt_act_elem%mv_conf)+1
 
         !========================================================================
         !-- sort element in lists
         !-- case NOT ACCEPTED
         IF (.NOT. result_acc) THEN !result NOT accepted
            IF (gt_act_elem%prob_acc > 0.0_dp) THEN
               IF (gt_act_elem%prob_acc >= 0.5_dp) THEN
                  ! wrong estimate (estimated accepted)
                  tmc_env%m_env%estim_corr_wrong(4) = tmc_env%m_env%estim_corr_wrong(4) + 1
                  IF (debug > 0) WRITE (tmc_env%m_env%io_unit, *) &
                     "Wrong guess for NACC (elem/estim acc prob)", gt_act_elem%nr, gt_act_elem%prob_acc
               ELSE
                  tmc_env%m_env%estim_corr_wrong(3) = tmc_env%m_env%estim_corr_wrong(3) + 1
               END IF
            END IF
            gt_act_elem%stat = status_rejected_result
            gt_act_elem%prob_acc = 0.0_dp
            itmp = gt_act_elem%conf(gt_act_elem%mv_conf)%elem%sub_tree_nr
 
            !-- case ACCEPTED
         ELSE
            IF (gt_act_elem%prob_acc > 0.0_dp) THEN
               IF (gt_act_elem%prob_acc <= 0.5_dp) THEN
                  ! wrong estimate (estimated NOT accepted)
                  tmc_env%m_env%estim_corr_wrong(2) = tmc_env%m_env%estim_corr_wrong(2) + 1
                  IF (debug > 0) WRITE (tmc_env%m_env%io_unit, *) &
                     "wrong guess for ACC  (elem/estim acc prob)", gt_act_elem%nr, gt_act_elem%prob_acc
               ELSE
                  tmc_env%m_env%estim_corr_wrong(1) = tmc_env%m_env%estim_corr_wrong(1) + 1
               END IF
            END IF
            gt_act_elem%stat = status_accepted_result
            gt_act_elem%prob_acc = 1.0_dp
 
            ! save result in the result list
            IF (.NOT. gt_act_elem%swaped) THEN
               !-- saving moved/changed configuration of result node
               tmc_env%m_env%result_list(gt_act_elem%mv_conf)%elem => &
                  gt_act_elem%conf(gt_act_elem%mv_conf)%elem
            ELSE
               ! ATTENTION: act_element != gt_act_elem%conf(gt_act_elem%mv_conf),
               !   because we take the last accepted conf
               tmc_env%m_env%result_list(gt_act_elem%mv_conf)%elem => act_element
               tmc_env%m_env%result_list(gt_act_elem%mv_conf + 1)%elem => tmp_element
            END IF
            tmc_env%m_env%gt_act => gt_act_elem
         END IF ! result acceptance check
         !======================================================================
         !-- set status accepted or rejected of certain subtree elements
         IF (.NOT. gt_act_elem%swaped) THEN
            CALL subtree_configuration_stat_change(gt_ptr=gt_act_elem, &
                                                   stat=gt_act_elem%stat, &
                                                   tmc_params=tmc_env%params)
         END IF
 
         IF (tmc_env%params%DRAW_TREE) THEN
            CALL create_global_tree_dot_color(gt_tree_element=gt_act_elem, &
                                              tmc_params=tmc_env%params)
         END IF
 
         ! probability update
         CALL prob_update(move_types=tmc_env%params%move_types, &
                          pt_el=gt_act_elem, &
                          prob_opt=tmc_env%params%esimate_acc_prob)
 
         !writes only different configurations with repetition file if possible
         CALL write_result_list_element(result_list=tmc_env%m_env%result_list, &
                                        result_count=tmc_env%m_env%result_count, &
                                        conf_updated=gt_act_elem%mv_conf, accepted=result_acc, &
                                        tmc_params=tmc_env%params)
         IF (gt_act_elem%swaped) THEN
            CALL write_result_list_element(result_list=tmc_env%m_env%result_list, &
                                           result_count=tmc_env%m_env%result_count, &
                                           conf_updated=gt_act_elem%mv_conf + 1, accepted=result_acc, &
                                           tmc_params=tmc_env%params)
         END IF
 
         ! save for analysis
         IF (tmc_env%tmc_comp_set%para_env_m_ana%num_pe > 1 .AND. result_acc) THEN
            CALL add_to_list(elem=tmc_env%m_env%result_list(gt_act_elem%mv_conf)%elem, &
                             list=tmc_env%m_env%analysis_list, &
                             temp_ind=gt_act_elem%mv_conf, &
                             nr=tmc_env%m_env%result_count(gt_act_elem%mv_conf))
            IF (gt_act_elem%swaped) THEN
               CALL add_to_list(elem=tmc_env%m_env%result_list(gt_act_elem%mv_conf + 1)%elem, &
                                list=tmc_env%m_env%analysis_list, &
                                temp_ind=gt_act_elem%mv_conf + 1, &
                                nr=tmc_env%m_env%result_count(gt_act_elem%mv_conf + 1))
            END IF
         END IF
      END DO search_calculated_element_loop
 
      ! end the timing
      CALL timestop(handle)
 
   END SUBROUTINE tree_update
 
   !============================================================================
   ! tree node acceptance probability
   !============================================================================
! **************************************************************************************************
!> \brief checks if element is ready for accaptance probability update
!>        checks status and the amount of already received intermediate energies
!> \param elem sub tree element to update
!> \return return value
!> \author Mandes 12.2012
! **************************************************************************************************
   FUNCTION ready_for_update_acc_prob(elem) RESULT(ready)
      TYPE(tree_type), POINTER                           :: elem
      LOGICAL                                            :: ready
 
      cpassert(ASSOCIATED(elem))
      ready = .false.
      IF ((elem%scf_energies_count >= 4) &
          .AND. (elem%stat /= status_deleted) .AND. (elem%stat /= status_deleted_result) &
          .AND. (elem%stat /= status_canceled_ener)) THEN
         ready = .true.
      END IF
   END FUNCTION ready_for_update_acc_prob
 
! **************************************************************************************************
!> \brief updates the subtree acceptance probability
!>        the swap probabilities are handled within the certain checks of the
!>        global tree elements (pt references)
!> \param tree_elem sub tree element to update
!> \param tmc_env TMC environment
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE check_elements_for_acc_prob_update(tree_elem, tmc_env)
      TYPE(tree_type), POINTER                           :: tree_elem
      TYPE(tmc_env_type), POINTER                        :: tmc_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'check_elements_for_acc_prob_update'
 
      INTEGER                                            :: handle
      TYPE(tree_type), POINTER                           :: act_elem, parent_elem
 
      NULLIFY (parent_elem, act_elem)
 
      cpassert(ASSOCIATED(tree_elem))
      cpassert(ASSOCIATED(tmc_env))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      act_elem => tree_elem
      !-- nothing to do if option is disabled or element not ready
      IF (tmc_env%params%esimate_acc_prob .AND. &
          ready_for_update_acc_prob(act_elem)) THEN
         ! ------------------------------------------------------
         ! check this element
         !   for usual moves and swapping
         ! get parent subtree elment for case of not swapped configurations
         parent_elem => search_parent_element(act_elem)
         ! update the prob of accaptance
         CALL update_prob_gt_node_list(reference_list=act_elem%gt_nodes_references, &
                                       act_elem=act_elem, parent_elem=parent_elem, act_parent=.true., &
                                       tmc_env=tmc_env)
 
         !------------------------------------------------------
         ! check the childs of this element
         parent_elem => tree_elem
 
         ! check ACC child (parent element is the entered/updated element)
         IF (ASSOCIATED(tree_elem%acc)) THEN
            act_elem => tree_elem%acc
            ! update the prob of accaptance
            CALL update_prob_gt_node_list(reference_list=act_elem%gt_nodes_references, &
                                          act_elem=act_elem, parent_elem=parent_elem, act_parent=.false., &
                                          tmc_env=tmc_env)
         END IF
 
         ! check all NACC childs of next accepted one
         nacc_loop: DO
            IF (.NOT. ASSOCIATED(act_elem%nacc)) EXIT nacc_loop
            act_elem => act_elem%nacc
            CALL update_prob_gt_node_list(reference_list=act_elem%gt_nodes_references, &
                                          act_elem=act_elem, parent_elem=parent_elem, act_parent=.false., &
                                          tmc_env=tmc_env)
         END DO nacc_loop
      END IF
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE check_elements_for_acc_prob_update
 
! **************************************************************************************************
!> \brief search all pt nodes in reference list and update the estimated
!>        acceptance probabilities
!> \param reference_list list of global tree element referring to the updated
!>        subtree element
!> \param act_elem actual sub tree element updated or child of the updated one
!> \param parent_elem parent or the actual updated subtree element
!> \param act_parent flag if updated element is the actual (TRUE) or
!>        the parent (FALSE) element
!> \param tmc_env TMC environment
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE update_prob_gt_node_list(reference_list, act_elem, parent_elem, act_parent, tmc_env)
      TYPE(gt_elem_list_type), POINTER                   :: reference_list
      TYPE(tree_type), POINTER                           :: act_elem, parent_elem
      LOGICAL                                            :: act_parent
      TYPE(tmc_env_type), POINTER                        :: tmc_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'update_prob_gt_node_list'
 
      INTEGER                                            :: handle, integration_steps
      REAL(kind=dp)                                      :: kb, tmp_prob
      TYPE(gt_elem_list_type), POINTER                   :: tmp_pt_ptr
      TYPE(tree_type), POINTER                           :: tmp_elem, tmp_parent_elem
 
      kb = boltzmann/joule
 
      NULLIFY (tmp_pt_ptr, tmp_elem, tmp_parent_elem)
 
      IF (.NOT. ASSOCIATED(reference_list)) RETURN ! could be canceled already
 
      cpassert(ASSOCIATED(reference_list%gt_elem))
      cpassert(ASSOCIATED(act_elem))
      cpassert(ASSOCIATED(parent_elem))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      ! check if element is ready
      IF (ready_for_update_acc_prob(act_elem)) THEN
         ! set the number of integration steps used for 3 point approximation
         ! of the exact energy, using the sub step energies
         ! still fixed value
         integration_steps = 100
 
         tmp_pt_ptr => reference_list
         ! do for all references using the updated subtree node
         reference_loop: DO WHILE (ASSOCIATED(tmp_pt_ptr))
            tmp_prob = -1.0_dp
            IF (tmp_pt_ptr%gt_elem%swaped) THEN
               NULLIFY (tmp_elem, tmp_parent_elem)
               ! in case of swap use the other swaped configuration as related one
               CALL get_subtree_elements_to_check(gt_act_elem=tmp_pt_ptr%gt_elem, &
                                                  elem1=tmp_elem, elem2=tmp_parent_elem)
               ! NOT if parent is the updated one, and check for correct elements ready
               IF (act_parent .AND. ASSOCIATED(act_elem, tmp_elem) .AND. &
                   ready_for_update_acc_prob(elem=tmp_parent_elem)) THEN
                  ! using ln(rnd)/-(beta_i-beta_j) < U_j-U_i)
                  tmp_prob = compute_estimated_prob(elem_old=tmp_parent_elem, &
                                                    elem_new=act_elem, &
                                                    e_classical_diff=0.0_dp, &
                                                    rnd_nr=tmp_pt_ptr%gt_elem%rnd_nr, &
                                                    beta=1.0_dp/(kb*(tmc_env%params%Temp(tmp_pt_ptr%gt_elem%mv_conf) - &
                                                                     tmc_env%params%Temp(tmp_pt_ptr%gt_elem%mv_conf + 1))), &
                                                    tmc_params=tmc_env%params)
               ELSE
                  tmp_pt_ptr => tmp_pt_ptr%next
                  cycle reference_loop
               END IF
            ELSE
               ! if no swap, use the parent configuration as related one
               IF (ready_for_update_acc_prob(parent_elem)) THEN
                  tmp_prob = compute_estimated_prob(elem_old=parent_elem, &
                                                    elem_new=act_elem, &
                                                    e_classical_diff=act_elem%e_pot_approx - parent_elem%e_pot_approx, &
                                                    rnd_nr=tmp_pt_ptr%gt_elem%rnd_nr, &
                                                    beta=1.0_dp/(kb*tmc_env%params%Temp(tmp_pt_ptr%gt_elem%mv_conf)), &
                                                    tmc_params=tmc_env%params)
               END IF
            END IF
            !successful probability update
            IF (tmp_prob >= 0.0_dp) THEN
               tmp_pt_ptr%gt_elem%prob_acc = tmp_prob
               !-- speculative canceling for the related direction
               IF (tmc_env%params%SPECULATIVE_CANCELING) THEN
                  CALL search_canceling_elements(pt_elem_in=tmp_pt_ptr%gt_elem, &
                                                 prob=tmp_pt_ptr%gt_elem%prob_acc, tmc_env=tmc_env)
               END IF
            END IF
 
            ! get next related global tree pointer
            tmp_pt_ptr => tmp_pt_ptr%next
         END DO reference_loop ! global tree references of actual subtree element
      END IF
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE update_prob_gt_node_list
 
END MODULE tmc_tree_acceptance