LCOV - code coverage report
Current view: top level - src/tmc - tmc_tree_acceptance.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:9f1e18e) Lines: 227 303 74.9 %
Date: 2021-09-22 20:58:48 Functions: 9 12 75.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2021 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief tree nodes acceptance
      10             : !>        code is separated in 3 parts, first the acceptance criteria,
      11             : !>        second the tree node acceptance handling, searching etc. and
      12             : !>        than the acceptance probability handling
      13             : !> \par History
      14             : !>      11.2012 created [Mandes Schoenherr]
      15             : !> \author Mandes
      16             : ! **************************************************************************************************
      17             : 
      18             : MODULE tmc_tree_acceptance
      19             :    USE cp_log_handling,                 ONLY: cp_to_string
      20             :    USE kinds,                           ONLY: dp
      21             :    USE physcon,                         ONLY: boltzmann,&
      22             :                                               joule
      23             :    USE tmc_calculations,                ONLY: compute_estimated_prob,&
      24             :                                               get_scaled_cell
      25             :    USE tmc_dot_tree,                    ONLY: create_dot_color,&
      26             :                                               create_global_tree_dot_color
      27             :    USE tmc_file_io,                     ONLY: write_result_list_element
      28             :    USE tmc_move_handle,                 ONLY: prob_update
      29             :    USE tmc_move_types,                  ONLY: mv_type_MD,&
      30             :                                               mv_type_volume_move
      31             :    USE tmc_stati,                       ONLY: task_type_gaussian_adaptation
      32             :    USE tmc_tree_build,                  ONLY: remove_unused_g_tree
      33             :    USE tmc_tree_search,                 ONLY: get_subtree_elements_to_check,&
      34             :                                               search_canceling_elements,&
      35             :                                               search_next_gt_element_to_check,&
      36             :                                               search_parent_element
      37             :    USE tmc_tree_types,                  ONLY: &
      38             :         add_to_list, global_tree_type, gt_elem_list_type, status_accepted, status_accepted_result, &
      39             :         status_calc_approx_ener, status_calculate_MD, status_calculate_NMC_steps, &
      40             :         status_calculate_energy, status_calculated, status_cancel_ener, status_cancel_nmc, &
      41             :         status_canceled_ener, status_canceled_nmc, status_created, status_deleted, &
      42             :         status_deleted_result, status_rejected, status_rejected_result, tree_type
      43             :    USE tmc_types,                       ONLY: tmc_env_type,&
      44             :                                               tmc_param_type
      45             : #include "../base/base_uses.f90"
      46             : 
      47             :    IMPLICIT NONE
      48             : 
      49             :    PRIVATE
      50             : 
      51             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_tree_acceptance'
      52             : 
      53             :    PUBLIC :: acceptance_check
      54             :    PUBLIC :: check_acceptance_of_depending_subtree_nodes, &
      55             :              check_elements_for_acc_prob_update
      56             :    PUBLIC :: tree_update
      57             : 
      58             :    INTEGER, PARAMETER :: DEBUG = 0
      59             : 
      60             : CONTAINS
      61             : 
      62             :    !============================================================================
      63             :    ! acceptance criteria calculations
      64             :    !============================================================================
      65             : ! **************************************************************************************************
      66             : !> \brief standard Monte Carlo and 2 potential acceptance check
      67             : !>        acceptance check of move from old(last accepted) to new configuration
      68             : !>        the sum of kinetic and potential energy is used
      69             : !>        acc(o->n)=min(1,exp( -beta*(H(n)-H(o)) ))
      70             : !> \param tree_element new/actual configuration
      71             : !> \param parent_element last accepted configuration
      72             : !> \param tmc_params TMC global parameters
      73             : !> \param temperature actual temperature configuration should be checked with
      74             : !> \param diff_pot_check 2potential check or not?
      75             : !> \param accept result (configuration accepted of rejected)
      76             : !> \param rnd_nr random number for acceptance check
      77             : !> \param approx_ener for NMC the approximated energies schould be used
      78             : !> \author Mandes 12.2012
      79             : ! **************************************************************************************************
      80        3844 :    SUBROUTINE acceptance_check(tree_element, parent_element, tmc_params, &
      81             :                                temperature, diff_pot_check, accept, rnd_nr, &
      82             :                                approx_ener)
      83             :       TYPE(tree_type), POINTER                           :: tree_element, parent_element
      84             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
      85             :       REAL(KIND=dp)                                      :: temperature
      86             :       LOGICAL                                            :: diff_pot_check, accept
      87             :       REAL(KIND=dp)                                      :: rnd_nr
      88             :       LOGICAL, OPTIONAL                                  :: approx_ener
      89             : 
      90             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'acceptance_check'
      91             : 
      92             :       INTEGER                                            :: handle
      93             :       REAL(KIND=dp)                                      :: ekin_last_acc, elem_ener, kB, parent_ener
      94             : 
      95        1922 :       CPASSERT(ASSOCIATED(tree_element))
      96        1922 :       CPASSERT(ASSOCIATED(parent_element))
      97        1922 :       CPASSERT(ASSOCIATED(tmc_params))
      98        1922 :       CPASSERT(temperature .GT. 0.0_dp)
      99        1922 :       CPASSERT(rnd_nr .GE. 0.0_dp .AND. rnd_nr .LE. 1.0_dp)
     100             : 
     101        1922 :       kB = boltzmann/joule
     102             : 
     103             :       ! start the timing
     104        1922 :       CALL timeset(routineN, handle)
     105             : 
     106        1922 :       IF (tmc_params%task_type .EQ. task_type_gaussian_adaptation) THEN
     107           0 :          CPABORT("")
     108             : !TODO       CALL acc_check_gauss_adapt(f=tree_element%potential, ct=tree_element%ekin_before_md, acc=accept)
     109             : !DO NOT DO       RETURN
     110             :       END IF
     111             : 
     112             :       !-- using 2 different potentials, the energy difference between these potentials
     113             :       !   and the two configurations have to be regarded.
     114             :       !   The ensamble should be in equilibrium state of the approximate potential
     115        1922 :       IF (diff_pot_check .AND. (tmc_params%NMC_inp_file .NE. "")) THEN
     116          66 :          IF ((tree_element%potential .EQ. HUGE(tree_element%potential)) .OR. &
     117             :              tree_element%e_pot_approx .EQ. HUGE(tree_element%e_pot_approx)) THEN
     118             :             elem_ener = HUGE(elem_ener)
     119             :          ELSE
     120             :             !for different potentials we have to regard the differences in energy
     121             :             ! min(1,e^{-\beta*[(E_{exact}(n)-E_{approx}(n))-(E_{exact}(o)-E_{approx}(o))]})
     122             :             elem_ener = 1.0_dp/(kB*temperature)*tree_element%potential &
     123             :                         - 1.0_dp/(kB*tmc_params%Temp(tree_element%temp_created)) &
     124          66 :                         *tree_element%e_pot_approx
     125             :          END IF
     126             :          parent_ener = 1.0_dp/(kB*temperature)*parent_element%potential &
     127             :                        - 1.0_dp/(kB*tmc_params%Temp(tree_element%temp_created)) &
     128          66 :                        *parent_element%e_pot_approx
     129             : 
     130             :          !-- always accepted if new energy is smaller than old energy
     131          66 :          IF (elem_ener .LE. parent_ener) THEN
     132          27 :             accept = .TRUE.
     133             :          ELSE
     134             :             !-- gaussian distributed acceptance if new energy is greater than old energy
     135          39 :             IF (rnd_nr .LT. &
     136             :                 EXP(-(elem_ener - parent_ener))) THEN
     137           8 :                accept = .TRUE.
     138             :             ELSE
     139          31 :                accept = .FALSE.
     140             :             END IF
     141             :          END IF
     142             :       ELSE
     143             :          !-- standard MC check, but also using the kinetic energy for Hybrid Monte Carlo
     144        1856 :          IF (tree_element%move_type .EQ. mv_type_MD) THEN
     145           0 :             ekin_last_acc = tree_element%ekin_before_md
     146             :             ! using the kinetic energy before velocity change
     147             :          ELSE
     148        1856 :             ekin_last_acc = parent_element%ekin
     149             :          END IF
     150             :          ! comparing aproximated energies
     151        1856 :          IF (PRESENT(approx_ener)) THEN
     152             :             elem_ener = tree_element%e_pot_approx &
     153         126 :                         + tree_element%ekin
     154             :             parent_ener = parent_element%e_pot_approx &
     155         126 :                           + ekin_last_acc
     156             :          ELSE
     157             :             elem_ener = tree_element%potential &
     158        1730 :                         + tree_element%ekin
     159             :             parent_ener = parent_element%potential &
     160        1730 :                           + ekin_last_acc
     161             :          END IF
     162             : 
     163             :          !-- always accepted if new energy is smaller than old energy
     164        1856 :          IF (elem_ener .LE. parent_ener) THEN
     165         365 :             accept = .TRUE.
     166             :          ELSE
     167             :             !-- gaussian distributed acceptance if new energy is greater than old energy
     168        1491 :             IF (rnd_nr .LT. &
     169             :                 EXP(-1.0_dp/(kB*temperature)*(elem_ener - parent_ener))) THEN
     170         177 :                accept = .TRUE.
     171             :             ELSE
     172        1314 :                accept = .FALSE.
     173             :             END IF
     174             :          END IF
     175             :       END IF
     176             : 
     177             :       ! update the estimated energy acceptance probability distribution
     178        1922 :       IF (diff_pot_check) THEN
     179        1796 :          CPASSERT(ASSOCIATED(tmc_params%prior_NMC_acc))
     180        1796 :          tmc_params%prior_NMC_acc%counter = tmc_params%prior_NMC_acc%counter + 1
     181             :          tmc_params%prior_NMC_acc%aver = (tmc_params%prior_NMC_acc%aver*(tmc_params%prior_NMC_acc%counter - 1) + &
     182        1796 :                                           ((elem_ener - parent_ener)))/REAL(tmc_params%prior_NMC_acc%counter, KIND=dp)
     183             :          tmc_params%prior_NMC_acc%aver_2 = (tmc_params%prior_NMC_acc%aver_2*(tmc_params%prior_NMC_acc%counter - 1) + &
     184        1796 :                                             (elem_ener - parent_ener)**2)/REAL(tmc_params%prior_NMC_acc%counter, KIND=dp)
     185             :       END IF
     186             : 
     187             :       ! end the timing
     188        1922 :       CALL timestop(handle)
     189        1922 :    END SUBROUTINE acceptance_check
     190             : 
     191             : ! **************************************************************************************************
     192             : !> \brief parallel tempering swap acceptance check,
     193             : !>        check swapped configurations of different temperatures
     194             : !>        using sum of potential and kinetic energy
     195             : !>        always accepted if energy of configuration with higher temperature
     196             : !>        is smaller than energy of configuration with lower temperature
     197             : !>        acc(o->n)=min(1,exp((beta_i-beta_j)*(U_i-U_j)))
     198             : !> \param tree_elem actual global tree element
     199             : !> \param conf1 sub tree element of lower temperature
     200             : !> \param conf2 sub tree element of higher temperature
     201             : !> \param tmc_params TMC environment parameters
     202             : !> \param accept acceptance of configurational change
     203             : !> \author Mandes 12.2012
     204             : ! **************************************************************************************************
     205         336 :    SUBROUTINE swap_acceptance_check(tree_elem, conf1, conf2, tmc_params, accept)
     206             :       TYPE(global_tree_type), POINTER                    :: tree_elem
     207             :       TYPE(tree_type), POINTER                           :: conf1, conf2
     208             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     209             :       LOGICAL                                            :: accept
     210             : 
     211             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'swap_acceptance_check'
     212             : 
     213             :       INTEGER                                            :: handle
     214             :       REAL(KIND=dp)                                      :: delta, elem1_ener, elem2_ener, kB, vol1, &
     215             :                                                             vol2
     216             : 
     217         168 :       kB = boltzmann/joule
     218             : 
     219         168 :       CPASSERT(ASSOCIATED(tree_elem))
     220         168 :       CPASSERT(tree_elem%rnd_nr .GE. 0.0_dp)
     221         168 :       CPASSERT(ASSOCIATED(conf1))
     222         168 :       CPASSERT(ASSOCIATED(conf2))
     223         168 :       CPASSERT(ASSOCIATED(tmc_params))
     224             : 
     225             :       ! start the timing
     226         168 :       CALL timeset(routineN, handle)
     227             : 
     228         168 :       IF (tmc_params%pressure .GT. 0.0_dp) THEN
     229             :          ! pt-NVT
     230             :          elem1_ener = conf1%potential &
     231          18 :                       + conf1%ekin
     232             :          elem2_ener = conf2%potential &
     233          18 :                       + conf2%ekin
     234             :          ! the swap is done with prob: exp((\beta_i-\beta_j)(U_i-U_j)),
     235             :          ! BUT because they are already swaped we exchange the energies.
     236          18 :          IF (tree_elem%rnd_nr .LT. EXP((1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf)) - &
     237             :                                         1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf + 1))) &
     238             :                                        *(elem2_ener - elem1_ener))) THEN
     239          18 :             accept = .TRUE.
     240             :          ELSE
     241           0 :             accept = .FALSE.
     242             :          END IF
     243             :       ELSE
     244             :          ! pt-NpT (parallel Tempering with constant pressure, temperature and num. particle)
     245             :          CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf1%box_scale, &
     246         150 :                               vol=vol1)
     247             :          CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf2%box_scale, &
     248         150 :                               vol=vol2)
     249             :          ! delta= (beta_m-beta_n)(U_n-U_m) + (beta_m*P_m-beta_n*P_n)(V_n-V_m)
     250             :          delta = (1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf)) &
     251             :                   - 1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf + 1)))* &
     252             :                  ((conf2%potential + conf2%ekin) - (conf1%potential + conf1%ekin)) + &
     253             :                  (1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf))*tmc_params%pressure &
     254             :                   - 1.0_dp/(kB*tmc_params%Temp(tree_elem%mv_conf + 1))*tmc_params%pressure)* &
     255         150 :                  (vol2 - vol1)
     256         150 :          IF (tree_elem%rnd_nr .LT. EXP(delta)) THEN
     257         112 :             accept = .TRUE.
     258             :          ELSE
     259          38 :             accept = .FALSE.
     260             :          END IF
     261             :       END IF
     262             :       ! end the timing
     263         168 :       CALL timestop(handle)
     264         168 :    END SUBROUTINE swap_acceptance_check
     265             : 
     266             : ! **************************************************************************************************
     267             : !> \brief volume move acceptance check
     268             : !>        sampling NPT, we need to do volume moves. Their acceptance are
     269             : !>        checked regarding the old and new energy, the volume difference
     270             : !>        and the constant pressure
     271             : !> \param parent_elem old tree element with old box size
     272             : !> \param new_elem actual tree element with new box size
     273             : !> \param temperature actual temperature
     274             : !> \param rnd_nr random number to check with
     275             : !> \param tmc_params TMC environment parameters
     276             : !> \param accept Monte Carlo move acceptance (result)
     277             : !> \author Mandes 12.2012
     278             : ! **************************************************************************************************
     279         340 :    SUBROUTINE volume_acceptance_check(parent_elem, new_elem, temperature, &
     280             :                                       rnd_nr, tmc_params, accept)
     281             :       TYPE(tree_type), POINTER                           :: parent_elem, new_elem
     282             :       REAL(KIND=dp)                                      :: temperature, rnd_nr
     283             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     284             :       LOGICAL                                            :: accept
     285             : 
     286             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'volume_acceptance_check'
     287             : 
     288             :       INTEGER                                            :: handle
     289             :       REAL(KIND=dp)                                      :: d_enthalpy, kB, n_vol, p_vol
     290             : 
     291         170 :       kB = boltzmann/joule
     292             : 
     293         170 :       CPASSERT(ASSOCIATED(parent_elem))
     294         170 :       CPASSERT(ASSOCIATED(new_elem))
     295         170 :       CPASSERT(temperature .GT. 0.0_dp)
     296         170 :       CPASSERT(rnd_nr .GT. 0.0_dp)
     297         170 :       CPASSERT(ASSOCIATED(tmc_params))
     298         170 :       CPASSERT(tmc_params%pressure .GE. 0.0_dp)
     299             : 
     300             :       ! start the timing
     301         170 :       CALL timeset(routineN, handle)
     302             : 
     303             :       CALL get_scaled_cell(cell=tmc_params%cell, box_scale=parent_elem%box_scale, &
     304         170 :                            vol=p_vol)
     305             :       CALL get_scaled_cell(cell=tmc_params%cell, box_scale=new_elem%box_scale, &
     306         170 :                            vol=n_vol)
     307             :       ! H=enthalpy, U=energy, P=pressure, V=volume, kB=Boltzmann constant, T=temperature, N=amount particle
     308             :       ! delta_H  =      delta_U                               + P*delta_V           - kB*T*N*ln(V_n/V_p)
     309             :       IF (.FALSE.) THEN
     310             :          ! the volume move in volume space (dV)
     311             :          d_enthalpy = (new_elem%potential - parent_elem%potential) + &
     312             :                       tmc_params%pressure*(n_vol - p_vol) - &
     313             :                       kB*temperature*(SIZE(new_elem%pos)/ &
     314             :                                       tmc_params%dim_per_elem)* &
     315             :                       LOG(n_vol/p_vol)
     316             :       ELSE
     317         170 :          IF (tmc_params%v_isotropic) THEN
     318             :             d_enthalpy = (new_elem%potential - parent_elem%potential) + &
     319             :                          tmc_params%pressure*(n_vol - p_vol) - &
     320             :                          kB*temperature*((SIZE(new_elem%pos)/ &
     321             :                                           tmc_params%dim_per_elem) + 2/REAL(3, KIND=dp))* &
     322         170 :                          LOG(n_vol/p_vol)
     323             :          ELSE
     324             :             d_enthalpy = (new_elem%potential - parent_elem%potential) + &
     325             :                          tmc_params%pressure*(n_vol - p_vol) - &
     326             :                          kB*temperature*(SIZE(new_elem%pos)/ &
     327             :                                          tmc_params%dim_per_elem)* &
     328           0 :                          LOG(n_vol/p_vol)
     329             :          END IF
     330             :       END IF
     331             :       ! acc(o->n) = min(1, exp(-beta*delta_H))
     332         170 :       IF (d_enthalpy .LE. 0.0_dp) THEN
     333          53 :          accept = .TRUE.
     334             :       ELSE
     335         117 :          IF (rnd_nr .LT. EXP(-1.0_dp/(kB*temperature)*d_enthalpy)) THEN
     336           9 :             accept = .TRUE.
     337             :          ELSE
     338         108 :             accept = .FALSE.
     339             :          END IF
     340             :       END IF
     341             :       ! end the timing
     342         170 :       CALL timestop(handle)
     343         170 :    END SUBROUTINE volume_acceptance_check
     344             : 
     345             :    !============================================================================
     346             :    ! tree node acceptance
     347             :    !============================================================================
     348             : ! **************************************************************************************************
     349             : !> \brief check acceptance of energy calculated element and related childs,
     350             : !>        when ready
     351             : !> \param tree_elem actual tree element with calculated energy
     352             : !> \param tmc_env TMC environment parameters
     353             : !> \author Mandes 12.2012
     354             : ! **************************************************************************************************
     355        3932 :    SUBROUTINE check_acceptance_of_depending_subtree_nodes(tree_elem, tmc_env)
     356             :       TYPE(tree_type), POINTER                           :: tree_elem
     357             :       TYPE(tmc_env_type), POINTER                        :: tmc_env
     358             : 
     359             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'check_acceptance_of_depending_subtree_nodes'
     360             : 
     361             :       INTEGER                                            :: handle
     362             :       LOGICAL                                            :: parent_ready
     363             :       TYPE(tree_type), POINTER                           :: act_elem, parent_elem
     364             : 
     365        1966 :       NULLIFY (parent_elem, act_elem)
     366             : 
     367        1966 :       CPASSERT(ASSOCIATED(tmc_env))
     368        1966 :       CPASSERT(ASSOCIATED(tree_elem))
     369        1966 :       CPASSERT(tree_elem%stat .EQ. status_calculated)
     370        1966 :       CPASSERT(ASSOCIATED(tree_elem%parent))
     371             : 
     372             :       ! start the timing
     373        1966 :       CALL timeset(routineN, handle)
     374             : 
     375        1966 :       act_elem => tree_elem
     376             : 
     377             :       ! ------------------------------------------------------
     378             :       ! check this element
     379        1966 :       parent_elem => search_parent_element(act_elem)
     380        1966 :       CPASSERT(.NOT. ASSOCIATED(act_elem, parent_elem))
     381             : 
     382             :       ! check status of parent element
     383        1966 :       SELECT CASE (parent_elem%stat)
     384             :       CASE (status_created, status_calculate_energy, status_calculate_MD, &
     385             :             status_calculate_NMC_steps, status_canceled_ener, &
     386             :             status_canceled_nmc, status_cancel_nmc, status_cancel_ener)
     387             :          parent_ready = .FALSE.
     388             :       CASE (status_accepted_result, status_rejected_result, &
     389             :             status_accepted, status_rejected, status_calculated)
     390           0 :          parent_ready = .TRUE.
     391             :       CASE DEFAULT
     392        1966 :          CPABORT("unknown parent element status"//cp_to_string(parent_elem%stat))
     393             :       END SELECT
     394             : 
     395             :       ! ready ?
     396           0 :       IF (parent_ready) THEN
     397             :          ! acceptance check
     398             :          CALL check_and_change_status_of_subtree_elem(act_elem=act_elem, &
     399        1966 :                                                       parent_elem=parent_elem, tmc_env=tmc_env)
     400             :       END IF
     401             :       !------------------------------------------------------
     402             :       ! check acc child
     403        1966 :       parent_elem => tree_elem
     404        1966 :       IF (ASSOCIATED(tree_elem%acc)) THEN
     405           0 :          act_elem => tree_elem%acc
     406           0 :          IF (act_elem%stat .EQ. status_calculated) THEN
     407             :             CALL check_and_change_status_of_subtree_elem(act_elem=act_elem, &
     408           0 :                                                          parent_elem=parent_elem, tmc_env=tmc_env)
     409             :          END IF
     410             :          !------------------------------------------------------
     411             :          ! check all nacc childs
     412             :          nacc_loop: DO
     413           0 :             IF (.NOT. ASSOCIATED(act_elem%nacc)) EXIT nacc_loop
     414           0 :             act_elem => act_elem%nacc
     415           0 :             IF (act_elem%stat .EQ. status_calculated) THEN
     416             :                CALL check_and_change_status_of_subtree_elem(act_elem=act_elem, &
     417           0 :                                                             parent_elem=parent_elem, tmc_env=tmc_env)
     418             :             END IF
     419             :          END DO nacc_loop
     420             :       END IF
     421             : 
     422             :       ! end the timing
     423        1966 :       CALL timestop(handle)
     424        1966 :    END SUBROUTINE check_acceptance_of_depending_subtree_nodes
     425             : 
     426             : ! **************************************************************************************************
     427             : !> \brief checking the elements and change the status
     428             : !> \param act_elem actual tree element (new configuration)
     429             : !> \param parent_elem parent tree element (old configuration)
     430             : !> \param tmc_env TMC environment parameters
     431             : !> \author Mandes 12.2012
     432             : ! **************************************************************************************************
     433        1966 :    SUBROUTINE check_and_change_status_of_subtree_elem(act_elem, parent_elem, tmc_env)
     434             :       TYPE(tree_type), POINTER                           :: act_elem, parent_elem
     435             :       TYPE(tmc_env_type), POINTER                        :: tmc_env
     436             : 
     437             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'check_and_change_status_of_subtree_elem'
     438             : 
     439             :       INTEGER                                            :: handle
     440             :       LOGICAL                                            :: flag, result_acc
     441             :       TYPE(global_tree_type), POINTER                    :: tmp_gt_ptr
     442             :       TYPE(gt_elem_list_type), POINTER                   :: tmp_gt_list_ptr
     443             : 
     444        1966 :       NULLIFY (tmp_gt_list_ptr, tmp_gt_ptr)
     445             : 
     446             :       ! start the timing
     447        1966 :       CALL timeset(routineN, handle)
     448             : 
     449        1966 :       CPASSERT(ASSOCIATED(tmc_env))
     450        1966 :       CPASSERT(ASSOCIATED(tmc_env%params))
     451        1966 :       CPASSERT(ASSOCIATED(act_elem))
     452        1966 :       CPASSERT(ASSOCIATED(parent_elem))
     453        1966 :       CPASSERT(act_elem%stat .EQ. status_calculated)
     454             :       MARK_USED(parent_elem)
     455             : 
     456        1966 :       flag = .FALSE.
     457             : 
     458        1966 :       tmp_gt_list_ptr => act_elem%gt_nodes_references
     459             :       ! check for all global tree elements referring to this subtree element
     460             :       ! same subtree element can be accepted at a certain temperature and rejected at another one
     461        3932 :       DO WHILE (ASSOCIATED(tmp_gt_list_ptr))
     462        1966 :          CPASSERT(tmp_gt_list_ptr%gt_elem%stat .NE. status_accepted_result)
     463        1966 :          CPASSERT(tmp_gt_list_ptr%gt_elem%stat .NE. status_rejected_result)
     464             : 
     465             :          CALL check_elements(gt_act_elem=tmp_gt_list_ptr%gt_elem, tmc_env=tmc_env, &
     466        1966 :                              check_done=flag, result_acc=result_acc)
     467        1966 :          IF (flag) THEN
     468             :             ! probability update
     469             :             CALL prob_update(move_types=tmc_env%params%move_types, elem=act_elem, &
     470        1966 :                              acc=result_acc, prob_opt=tmc_env%params%esimate_acc_prob)
     471             : 
     472             :             ! change status
     473             :             ! accepted case: delete (and cancel) not accepted branch
     474        1966 :             IF (result_acc .EQV. .TRUE.) THEN
     475         579 :                tmp_gt_list_ptr%gt_elem%stat = status_accepted
     476         579 :                IF (ASSOCIATED(tmp_gt_list_ptr%gt_elem%nacc)) THEN
     477           0 :                   tmp_gt_ptr => tmp_gt_list_ptr%gt_elem%nacc
     478             :                   CALL remove_unused_g_tree(begin_ptr=tmp_gt_ptr, &
     479             :                                             end_ptr=tmp_gt_list_ptr%gt_elem, removed=flag, &
     480           0 :                                             tmc_env=tmc_env)
     481             :                END IF
     482             :             ELSE
     483             :                ! not accepted case: delete (and cancel) accepted branch
     484        1387 :                tmp_gt_list_ptr%gt_elem%stat = status_rejected
     485        1387 :                IF (ASSOCIATED(tmp_gt_list_ptr%gt_elem%acc)) THEN
     486           0 :                   tmp_gt_ptr => tmp_gt_list_ptr%gt_elem%acc
     487             :                   CALL remove_unused_g_tree(begin_ptr=tmp_gt_ptr, &
     488             :                                             end_ptr=tmp_gt_list_ptr%gt_elem, removed=flag, &
     489           0 :                                             tmc_env=tmc_env)
     490             :                END IF
     491             :             END IF
     492             : 
     493        1966 :             IF (tmc_env%params%DRAW_TREE) &
     494             :                CALL create_global_tree_dot_color(gt_tree_element=tmp_gt_list_ptr%gt_elem, &
     495          36 :                                                  tmc_params=tmc_env%params)
     496             :          END IF
     497        1966 :          tmp_gt_list_ptr => tmp_gt_list_ptr%next
     498             :       END DO
     499             :       ! end the timing
     500        1966 :       CALL timestop(handle)
     501        1966 :    END SUBROUTINE check_and_change_status_of_subtree_elem
     502             : 
     503             : ! **************************************************************************************************
     504             : !> \brief change status flag of all subtree elements
     505             : !>        most important to draw the right subtrees
     506             : !> \param gt_ptr global tree pointer, the related/changed sub tree element
     507             : !>        should change status
     508             : !> \param stat the status of the global tree element
     509             : !> \param tmc_params TMC emvironment parameters for drawing
     510             : !> \author Mandes 12.2012
     511             : ! **************************************************************************************************
     512        3932 :    SUBROUTINE subtree_configuration_stat_change(gt_ptr, stat, tmc_params)
     513             :       TYPE(global_tree_type), POINTER                    :: gt_ptr
     514             :       INTEGER                                            :: stat
     515             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     516             : 
     517             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'subtree_configuration_stat_change'
     518             : 
     519             :       INTEGER                                            :: handle
     520             :       TYPE(tree_type), POINTER                           :: ptr
     521             : 
     522        1966 :       NULLIFY (ptr)
     523             : 
     524        1966 :       CPASSERT(ASSOCIATED(gt_ptr))
     525        1966 :       CPASSERT(ASSOCIATED(tmc_params))
     526             : 
     527             :       ! start the timing
     528        1966 :       CALL timeset(routineN, handle)
     529             : 
     530             :       ! don't respect swaped nodes (subtree nodes are already in the list,
     531             :       !   and we dont want them in marked in subtrees)
     532        1966 :       IF (.NOT. gt_ptr%swaped) THEN
     533        1966 :          ptr => gt_ptr%conf(gt_ptr%mv_conf)%elem
     534             :          ! dependent on the status (don't map each status 1:1,
     535             :          !   e.g. not the result ones)
     536        2545 :          SELECT CASE (stat)
     537             :          CASE (status_accepted_result)
     538         579 :             ptr%stat = status_accepted
     539             :          CASE (status_rejected_result)
     540        1387 :             ptr%stat = status_rejected
     541             :          CASE (status_accepted, status_rejected)
     542           0 :             ptr%stat = stat
     543             :          CASE DEFAULT
     544             :             CALL cp_abort(__LOCATION__, &
     545             :                           "unknown global tree status"// &
     546        1966 :                           cp_to_string(stat))
     547             :          END SELECT
     548             : 
     549        1966 :          IF (tmc_params%DRAW_TREE) &
     550          36 :             CALL create_dot_color(tree_element=ptr, tmc_params=tmc_params)
     551             :       END IF
     552             : 
     553             :       ! end the timing
     554        1966 :       CALL timestop(handle)
     555        1966 :    END SUBROUTINE subtree_configuration_stat_change
     556             : 
     557             : ! **************************************************************************************************
     558             : !> \brief checks if the element is ready for an acceptance check
     559             : !> \param elem sub tree element
     560             : !> \param ready return value
     561             : !> \author Mandes 12.2012
     562             : ! **************************************************************************************************
     563      520808 :    SUBROUTINE elem_ready_to_check(elem, ready)
     564             :       TYPE(tree_type), POINTER                           :: elem
     565             :       LOGICAL                                            :: ready
     566             : 
     567             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'elem_ready_to_check'
     568             : 
     569             :       INTEGER                                            :: handle
     570             : 
     571      260404 :       CPASSERT(ASSOCIATED(elem))
     572             : 
     573             :       ! start the timing
     574      260404 :       CALL timeset(routineN, handle)
     575             : 
     576      260404 :       ready = .FALSE.
     577             : 
     578      260404 :       SELECT CASE (elem%stat)
     579             :       CASE (status_created, status_calculate_energy, &
     580             :             status_calculate_MD, status_calculate_NMC_steps, status_calc_approx_ener, &
     581             :             status_cancel_nmc, status_cancel_ener, status_canceled_nmc, status_canceled_ener)
     582      123797 :          ready = .FALSE.
     583             :       CASE (status_calculated, status_accepted_result, status_accepted, status_rejected, status_rejected_result)
     584      123797 :          ready = .TRUE.
     585             :       CASE DEFAULT
     586             :          CALL cp_abort(__LOCATION__, &
     587             :                        "status of actual tree node is unknown" &
     588      260404 :                        //cp_to_string(elem%stat))
     589             :       END SELECT
     590             :       ! end the timing
     591      260404 :       CALL timestop(handle)
     592      260404 :    END SUBROUTINE elem_ready_to_check
     593             : 
     594             : ! **************************************************************************************************
     595             : !> \brief checking the elements (find element and do acceptance check)
     596             : !> \param gt_act_elem actual global tree element
     597             : !> \param tmc_env TMC environment
     598             : !> \param check_done successful checked? (result)
     599             : !> \param result_acc checked configuration accepted? (result)
     600             : !> \author Mandes 12.2012
     601             : ! **************************************************************************************************
     602      264336 :    SUBROUTINE check_elements(gt_act_elem, tmc_env, check_done, result_acc)
     603             :       TYPE(global_tree_type), POINTER                    :: gt_act_elem
     604             :       TYPE(tmc_env_type), POINTER                        :: tmc_env
     605             :       LOGICAL                                            :: check_done, result_acc
     606             : 
     607             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'check_elements'
     608             : 
     609             :       INTEGER                                            :: handle
     610             :       LOGICAL                                            :: act_ready, tmp_ready
     611             :       TYPE(tree_type), POINTER                           :: act_element, tmp_element
     612             : 
     613      132168 :       NULLIFY (act_element, tmp_element)
     614      132168 :       check_done = .FALSE.
     615             : 
     616      132168 :       CPASSERT(ASSOCIATED(gt_act_elem))
     617      132168 :       CPASSERT(ASSOCIATED(tmc_env))
     618      132168 :       CPASSERT(ASSOCIATED(tmc_env%m_env))
     619             : 
     620             :       ! start the timing
     621      132168 :       CALL timeset(routineN, handle)
     622             : 
     623             :       !====================================================================================
     624             :       ! check if global tree element is already checked
     625      134134 :       SELECT CASE (gt_act_elem%stat)
     626             :       CASE (status_accepted, status_rejected)
     627        1966 :          check_done = .TRUE.
     628      132168 :          IF (gt_act_elem%stat .EQ. status_accepted) THEN
     629         579 :             result_acc = .TRUE.
     630        1387 :          ELSE IF (gt_act_elem%stat .EQ. status_rejected) THEN
     631        1387 :             result_acc = .FALSE.
     632             :          ELSE
     633             :             CALL cp_abort(__LOCATION__, &
     634             :                           "undefinite status of already checked elements:" &
     635           0 :                           //cp_to_string(gt_act_elem%stat))
     636             :          END IF
     637             :       CASE DEFAULT
     638             :       END SELECT
     639             : 
     640      132168 :       IF (.NOT. check_done) THEN
     641             :          !====================================================================================
     642             :          ! get elements related to global tree element to check
     643             :          CALL get_subtree_elements_to_check(gt_act_elem=gt_act_elem, elem1=act_element, &
     644      130202 :                                             elem2=tmp_element)
     645      130202 :          CPASSERT(ASSOCIATED(act_element))
     646      130202 :          CPASSERT(ASSOCIATED(tmp_element))
     647             : 
     648             :          ! check status of both elements (if they are already calculated and hence ready to check)
     649      130202 :          CALL elem_ready_to_check(elem=act_element, ready=act_ready)
     650      130202 :          CALL elem_ready_to_check(elem=tmp_element, ready=tmp_ready)
     651             : 
     652             :          ! both ready ? check
     653      130202 :          IF (act_ready .AND. tmp_ready) THEN
     654             :             ! 2 temperature (swap) check
     655        2134 :             IF (gt_act_elem%swaped) THEN
     656             :                CALL swap_acceptance_check(tree_elem=gt_act_elem, conf1=act_element, &
     657             :                                           conf2=tmp_element, tmc_params=tmc_env%params, &
     658         168 :                                           accept=result_acc)
     659             :                ! volume move check
     660        1966 :             ELSE IF (act_element%move_type .EQ. mv_type_volume_move) THEN
     661             :                CALL volume_acceptance_check(parent_elem=tmp_element, new_elem=act_element, &
     662             :                                             temperature=tmc_env%params%Temp(gt_act_elem%mv_conf), &
     663             :                                             rnd_nr=gt_act_elem%rnd_nr, &
     664         170 :                                             tmc_params=tmc_env%params, accept=result_acc)
     665             :             ELSE
     666        1796 :                IF (tmc_env%m_env%temp_decrease .NE. 1.0_dp) THEN
     667             :                   CALL acceptance_check(tree_element=act_element, parent_element=tmp_element, &
     668             :                                         tmc_params=tmc_env%params, temperature=gt_act_elem%Temp, &
     669             :                                         diff_pot_check=.TRUE., accept=result_acc, &
     670           0 :                                         rnd_nr=gt_act_elem%rnd_nr)
     671             :                ELSE
     672             :                   CALL acceptance_check(tree_element=act_element, parent_element=tmp_element, &
     673             :                                         tmc_params=tmc_env%params, &
     674             :                                         temperature=tmc_env%params%Temp(gt_act_elem%mv_conf), &
     675             :                                         diff_pot_check=.TRUE., accept=result_acc, &
     676        1796 :                                         rnd_nr=gt_act_elem%rnd_nr)
     677             :                END IF
     678             :             END IF
     679        2134 :             check_done = .TRUE.
     680             :          ELSE
     681             :             ! if element is not ready, update status of global tree element
     682             :             !TODO maybe update this part (how much status we need in global tree from sub tree??)
     683      256136 :             SELECT CASE (gt_act_elem%stat)
     684             :                !-- check if at least the MD move is already done
     685             :                !-- hence new configurations can be created on basis of this configuration
     686             :             CASE (status_calculate_MD, status_calculate_NMC_steps, status_calculate_energy, &
     687             :                   status_created, status_calc_approx_ener)
     688      128068 :                IF (gt_act_elem%conf(gt_act_elem%mv_conf)%elem%stat .NE. gt_act_elem%stat) &
     689        1959 :                   gt_act_elem%stat = gt_act_elem%conf(gt_act_elem%mv_conf)%elem%stat
     690             :             CASE (status_calculated)
     691             :             CASE DEFAULT
     692             :                CALL cp_abort(__LOCATION__, &
     693             :                              "status of actual checked node is unknown" &
     694      128068 :                              //cp_to_string(gt_act_elem%stat))
     695             :             END SELECT
     696             :          END IF
     697             :       END IF
     698             :       ! end the timing
     699      132168 :       CALL timestop(handle)
     700      132168 :    END SUBROUTINE check_elements
     701             : 
     702             : ! **************************************************************************************************
     703             : !> \brief searching tree nodes to check for Markov Chain, elements are marked
     704             : !>        and stored in lists ... (main entry point)
     705             : !> \param tmc_env TMC environment
     706             : !> \param result_acc checked configuration accepted? (result)
     707             : !> \param something_updated ...
     708             : !> \param
     709             : !> \param
     710             : !> \author Mandes 12.2012
     711             : ! **************************************************************************************************
     712      260432 :    SUBROUTINE tree_update(tmc_env, result_acc, something_updated)
     713             :       TYPE(tmc_env_type), POINTER                        :: tmc_env
     714             :       LOGICAL                                            :: result_acc, something_updated
     715             : 
     716             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tree_update'
     717             : 
     718             :       INTEGER                                            :: handle, itmp
     719             :       LOGICAL                                            :: found
     720             :       TYPE(global_tree_type), POINTER                    :: gt_act_elem
     721             :       TYPE(tree_type), POINTER                           :: act_element, tmp_element
     722             : 
     723      130216 :       NULLIFY (gt_act_elem, tmp_element, act_element)
     724             : 
     725      130216 :       CPASSERT(ASSOCIATED(tmc_env))
     726             : 
     727             :       ! start the timing
     728      130216 :       CALL timeset(routineN, handle)
     729             : 
     730      130216 :       result_acc = .FALSE.
     731      130216 :       something_updated = .FALSE.
     732             : 
     733      130216 :       gt_act_elem => tmc_env%m_env%gt_act
     734             :       search_calculated_element_loop: DO
     735      132350 :          NULLIFY (act_element, tmp_element)
     736             :          ! search for element to check
     737      132350 :          CALL search_next_gt_element_to_check(ptr=gt_act_elem, found=found)
     738      132350 :          IF (.NOT. found) EXIT search_calculated_element_loop
     739             : 
     740             :          ! check the elements status end, if possible, do acceptance check
     741             :          CALL check_elements(gt_act_elem=gt_act_elem, tmc_env=tmc_env, &
     742      130202 :                              check_done=found, result_acc=result_acc)
     743             :          ! if check was not possible, update the status of the global tree element and return
     744      130202 :          IF (.NOT. found) EXIT search_calculated_element_loop
     745             : 
     746             :          ! get elements related to global tree element, which were checked
     747             :          CALL get_subtree_elements_to_check(gt_act_elem=gt_act_elem, &
     748        2134 :                                             elem1=act_element, elem2=tmp_element)
     749             : 
     750             :          !========================================================================
     751             :          !-- set result counters
     752             :          ! counter for certain temperature
     753             :          tmc_env%m_env%result_count(gt_act_elem%mv_conf) = &
     754        2134 :             tmc_env%m_env%result_count(gt_act_elem%mv_conf) + 1
     755             :          ! in case of swapped also count the result for
     756             :          !  the other swapped temperature
     757        2134 :          IF (gt_act_elem%swaped) &
     758             :             tmc_env%m_env%result_count(gt_act_elem%mv_conf + 1) = &
     759         168 :             tmc_env%m_env%result_count(gt_act_elem%mv_conf + 1) + 1
     760             :          ! count also for global tree Markov Chain
     761        2134 :          tmc_env%m_env%result_count(0) = tmc_env%m_env%result_count(0) + 1
     762             : 
     763             :          ! flag for doing tree cleaning with canceling and deallocation
     764        2134 :          something_updated = .TRUE.
     765             : 
     766             :          !IF((.NOT.gt_act_elem%swaped) .AND. act_element%temp_created.NE.0 .AND. &
     767             :          !   act_element%temp_created.NE.gt_act_elem%mv_conf)&
     768             :          !  count_wrong_temp_move(gt_act_elem%mv_conf) = &
     769             :          !       count_wrong_temp_move(gt_act_elem%mv_conf)+1
     770             : 
     771             :          !========================================================================
     772             :          !-- sort element in lists
     773             :          !-- case NOT ACCEPTED
     774        2134 :          IF (.NOT. result_acc) THEN !result NOT accepted
     775        1425 :             IF (gt_act_elem%prob_acc .GT. 0.0_dp) THEN
     776        1397 :                IF (gt_act_elem%prob_acc .GE. 0.5_dp) THEN
     777             :                   ! wrong estimate (estimated accepted)
     778         197 :                   tmc_env%m_env%estim_corr_wrong(4) = tmc_env%m_env%estim_corr_wrong(4) + 1
     779             :                   IF (DEBUG .GT. 0) WRITE (tmc_env%m_env%io_unit, *) &
     780             :                      "Wrong guess for NACC (elem/estim acc prob)", gt_act_elem%nr, gt_act_elem%prob_acc
     781             :                ELSE
     782        1200 :                   tmc_env%m_env%estim_corr_wrong(3) = tmc_env%m_env%estim_corr_wrong(3) + 1
     783             :                END IF
     784             :             END IF
     785        1425 :             gt_act_elem%stat = status_rejected_result
     786        1425 :             gt_act_elem%prob_acc = 0.0_dp
     787        1425 :             itmp = gt_act_elem%conf(gt_act_elem%mv_conf)%elem%sub_tree_nr
     788             : 
     789             :             !-- case ACCEPTED
     790             :          ELSE
     791         709 :             IF (gt_act_elem%prob_acc .GT. 0.0_dp) THEN
     792         694 :                IF (gt_act_elem%prob_acc .LE. 0.5_dp) THEN
     793             :                   ! wrong estimate (estimated NOT accepted)
     794         242 :                   tmc_env%m_env%estim_corr_wrong(2) = tmc_env%m_env%estim_corr_wrong(2) + 1
     795             :                   IF (DEBUG .GT. 0) WRITE (tmc_env%m_env%io_unit, *) &
     796             :                      "wrong guess for ACC  (elem/estim acc prob)", gt_act_elem%nr, gt_act_elem%prob_acc
     797             :                ELSE
     798         452 :                   tmc_env%m_env%estim_corr_wrong(1) = tmc_env%m_env%estim_corr_wrong(1) + 1
     799             :                END IF
     800             :             END IF
     801         709 :             gt_act_elem%stat = status_accepted_result
     802         709 :             gt_act_elem%prob_acc = 1.0_dp
     803             : 
     804             :             ! save result in the result list
     805         709 :             IF (.NOT. gt_act_elem%swaped) THEN
     806             :                !-- saving moved/changed configuration of result node
     807             :                tmc_env%m_env%result_list(gt_act_elem%mv_conf)%elem => &
     808         579 :                   gt_act_elem%conf(gt_act_elem%mv_conf)%elem
     809             :             ELSE
     810             :                ! ATTENTION: act_element != gt_act_elem%conf(gt_act_elem%mv_conf),
     811             :                !   because we take the last accepted conf
     812         130 :                tmc_env%m_env%result_list(gt_act_elem%mv_conf)%elem => act_element
     813         130 :                tmc_env%m_env%result_list(gt_act_elem%mv_conf + 1)%elem => tmp_element
     814             :             END IF
     815         709 :             tmc_env%m_env%gt_act => gt_act_elem
     816             :          END IF ! result acceptance check
     817             :          !======================================================================
     818             :          !-- set status accepted or rejected of certain subtree elements
     819        2134 :          IF (.NOT. gt_act_elem%swaped) &
     820             :             CALL subtree_configuration_stat_change(gt_ptr=gt_act_elem, &
     821             :                                                    stat=gt_act_elem%stat, &
     822        1966 :                                                    tmc_params=tmc_env%params)
     823             : 
     824        2134 :          IF (tmc_env%params%DRAW_TREE) &
     825             :             CALL create_global_tree_dot_color(gt_tree_element=gt_act_elem, &
     826          74 :                                               tmc_params=tmc_env%params)
     827             : 
     828             :          ! probability update
     829             :          CALL prob_update(move_types=tmc_env%params%move_types, &
     830             :                           pt_el=gt_act_elem, &
     831        2134 :                           prob_opt=tmc_env%params%esimate_acc_prob)
     832             : 
     833             :          !writes only different configurations with repetition file if possible
     834             :          CALL write_result_list_element(result_list=tmc_env%m_env%result_list, &
     835             :                                         result_count=tmc_env%m_env%result_count, &
     836             :                                         conf_updated=gt_act_elem%mv_conf, accepted=result_acc, &
     837        2134 :                                         tmc_params=tmc_env%params)
     838        2134 :          IF (gt_act_elem%swaped) &
     839             :             CALL write_result_list_element(result_list=tmc_env%m_env%result_list, &
     840             :                                            result_count=tmc_env%m_env%result_count, &
     841             :                                            conf_updated=gt_act_elem%mv_conf + 1, accepted=result_acc, &
     842         168 :                                            tmc_params=tmc_env%params)
     843             : 
     844             :          ! save for analysis
     845      132350 :          IF (tmc_env%tmc_comp_set%para_env_m_ana%num_pe .GT. 1 .AND. result_acc) THEN
     846             :             CALL add_to_list(elem=tmc_env%m_env%result_list(gt_act_elem%mv_conf)%elem, &
     847             :                              list=tmc_env%m_env%analysis_list, &
     848             :                              temp_ind=gt_act_elem%mv_conf, &
     849           0 :                              nr=tmc_env%m_env%result_count(gt_act_elem%mv_conf))
     850           0 :             IF (gt_act_elem%swaped) THEN
     851             :                CALL add_to_list(elem=tmc_env%m_env%result_list(gt_act_elem%mv_conf + 1)%elem, &
     852             :                                 list=tmc_env%m_env%analysis_list, &
     853             :                                 temp_ind=gt_act_elem%mv_conf + 1, &
     854           0 :                                 nr=tmc_env%m_env%result_count(gt_act_elem%mv_conf + 1))
     855             :             END IF
     856             :          END IF
     857             :       END DO search_calculated_element_loop
     858             : 
     859             :       ! end the timing
     860      130216 :       CALL timestop(handle)
     861             : 
     862      130216 :    END SUBROUTINE tree_update
     863             : 
     864             :    !============================================================================
     865             :    ! tree node acceptance probability
     866             :    !============================================================================
     867             : ! **************************************************************************************************
     868             : !> \brief checks if element is ready for accaptance probability update
     869             : !>        checks status and the amount of already received intermediate energies
     870             : !> \param elem sub tree element to update
     871             : !> \return return value
     872             : !> \author Mandes 12.2012
     873             : ! **************************************************************************************************
     874           0 :    FUNCTION ready_for_update_acc_prob(elem) RESULT(ready)
     875             :       TYPE(tree_type), POINTER                           :: elem
     876             :       LOGICAL                                            :: ready
     877             : 
     878           0 :       CPASSERT(ASSOCIATED(elem))
     879           0 :       ready = .FALSE.
     880             :       IF ((elem%scf_energies_count .GE. 4) &
     881             :           .AND. (elem%stat .NE. status_deleted) .AND. (elem%stat .NE. status_deleted_result) &
     882           0 :           .AND. (elem%stat .NE. status_canceled_ener)) &
     883           0 :          ready = .TRUE.
     884           0 :    END FUNCTION ready_for_update_acc_prob
     885             : 
     886             : ! **************************************************************************************************
     887             : !> \brief updates the subtree acceptance probability
     888             : !>        the swap probabilities are handled within the certain checks of the
     889             : !>        global tree elements (pt references)
     890             : !> \param tree_elem sub tree element to update
     891             : !> \param tmc_env TMC environment
     892             : !> \author Mandes 12.2012
     893             : ! **************************************************************************************************
     894           0 :    SUBROUTINE check_elements_for_acc_prob_update(tree_elem, tmc_env)
     895             :       TYPE(tree_type), POINTER                           :: tree_elem
     896             :       TYPE(tmc_env_type), POINTER                        :: tmc_env
     897             : 
     898             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'check_elements_for_acc_prob_update'
     899             : 
     900             :       INTEGER                                            :: handle
     901             :       TYPE(tree_type), POINTER                           :: act_elem, parent_elem
     902             : 
     903           0 :       NULLIFY (parent_elem, act_elem)
     904             : 
     905           0 :       CPASSERT(ASSOCIATED(tree_elem))
     906           0 :       CPASSERT(ASSOCIATED(tmc_env))
     907             : 
     908             :       ! start the timing
     909           0 :       CALL timeset(routineN, handle)
     910             : 
     911           0 :       act_elem => tree_elem
     912             :       !-- nothing to do if option is disabled or element not ready
     913           0 :       IF (tmc_env%params%esimate_acc_prob .AND. &
     914             :           ready_for_update_acc_prob(act_elem)) THEN
     915             :          ! ------------------------------------------------------
     916             :          ! check this element
     917             :          !   for usual moves and swapping
     918             :          ! get parent subtree elment for case of not swapped configurations
     919           0 :          parent_elem => search_parent_element(act_elem)
     920             :          ! update the prob of accaptance
     921             :          CALL update_prob_gt_node_list(reference_list=act_elem%gt_nodes_references, &
     922             :                                        act_elem=act_elem, parent_elem=parent_elem, act_parent=.TRUE., &
     923           0 :                                        tmc_env=tmc_env)
     924             : 
     925             :          !------------------------------------------------------
     926             :          ! check the childs of this element
     927           0 :          parent_elem => tree_elem
     928             : 
     929             :          ! check ACC child (parent element is the entered/updated element)
     930           0 :          IF (ASSOCIATED(tree_elem%acc)) THEN
     931           0 :             act_elem => tree_elem%acc
     932             :             ! update the prob of accaptance
     933             :             CALL update_prob_gt_node_list(reference_list=act_elem%gt_nodes_references, &
     934             :                                           act_elem=act_elem, parent_elem=parent_elem, act_parent=.FALSE., &
     935           0 :                                           tmc_env=tmc_env)
     936             :          END IF
     937             : 
     938             :          ! check all NACC childs of next accepted one
     939           0 :          nacc_loop: DO
     940           0 :             IF (.NOT. ASSOCIATED(act_elem%nacc)) EXIT nacc_loop
     941           0 :             act_elem => act_elem%nacc
     942             :             CALL update_prob_gt_node_list(reference_list=act_elem%gt_nodes_references, &
     943             :                                           act_elem=act_elem, parent_elem=parent_elem, act_parent=.FALSE., &
     944           0 :                                           tmc_env=tmc_env)
     945             :          END DO nacc_loop
     946             :       END IF
     947             :       ! end the timing
     948           0 :       CALL timestop(handle)
     949           0 :    END SUBROUTINE check_elements_for_acc_prob_update
     950             : 
     951             : ! **************************************************************************************************
     952             : !> \brief search all pt nodes in reference list and update the estimated
     953             : !>        acceptance probabilities
     954             : !> \param reference_list list of global tree element referring to the updated
     955             : !>        subtree element
     956             : !> \param act_elem actual sub tree element updated or child of the updated one
     957             : !> \param parent_elem parent or the actual updated subtree element
     958             : !> \param act_parent flag if updated element is the actual (TRUE) or
     959             : !>        the parent (FALSE) element
     960             : !> \param tmc_env TMC environment
     961             : !> \author Mandes 12.2012
     962             : ! **************************************************************************************************
     963           0 :    SUBROUTINE update_prob_gt_node_list(reference_list, act_elem, parent_elem, act_parent, tmc_env)
     964             :       TYPE(gt_elem_list_type), POINTER                   :: reference_list
     965             :       TYPE(tree_type), POINTER                           :: act_elem, parent_elem
     966             :       LOGICAL                                            :: act_parent
     967             :       TYPE(tmc_env_type), POINTER                        :: tmc_env
     968             : 
     969             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'update_prob_gt_node_list'
     970             : 
     971             :       INTEGER                                            :: handle, integration_steps
     972             :       REAL(KIND=dp)                                      :: kB, tmp_prob
     973             :       TYPE(gt_elem_list_type), POINTER                   :: tmp_pt_ptr
     974             :       TYPE(tree_type), POINTER                           :: tmp_elem, tmp_parent_elem
     975             : 
     976           0 :       kB = boltzmann/joule
     977             : 
     978           0 :       NULLIFY (tmp_pt_ptr, tmp_elem, tmp_parent_elem)
     979             : 
     980           0 :       IF (.NOT. ASSOCIATED(reference_list)) RETURN ! could be canceled already
     981             : 
     982           0 :       CPASSERT(ASSOCIATED(reference_list%gt_elem))
     983           0 :       CPASSERT(ASSOCIATED(act_elem))
     984           0 :       CPASSERT(ASSOCIATED(parent_elem))
     985             : 
     986             :       ! start the timing
     987           0 :       CALL timeset(routineN, handle)
     988             : 
     989             :       ! check if element is ready
     990           0 :       IF (ready_for_update_acc_prob(act_elem)) THEN
     991             :          ! set the number of integration steps used for 3 point approximation
     992             :          ! of the exact energy, using the sub step energies
     993             :          ! still fixed value
     994           0 :          integration_steps = 100
     995             : 
     996           0 :          tmp_pt_ptr => reference_list
     997             :          ! do for all references using the updated subtree node
     998           0 :          reference_loop: DO WHILE (ASSOCIATED(tmp_pt_ptr))
     999           0 :             tmp_prob = -1.0_dp
    1000           0 :             IF (tmp_pt_ptr%gt_elem%swaped) THEN
    1001           0 :                NULLIFY (tmp_elem, tmp_parent_elem)
    1002             :                ! in case of swap use the other swaped configuration as related one
    1003             :                CALL get_subtree_elements_to_check(gt_act_elem=tmp_pt_ptr%gt_elem, &
    1004           0 :                                                   elem1=tmp_elem, elem2=tmp_parent_elem)
    1005             :                ! NOT if parent is the updated one, and check for correct elements ready
    1006           0 :                IF (act_parent .AND. ASSOCIATED(act_elem, tmp_elem) .AND. &
    1007             :                    ready_for_update_acc_prob(elem=tmp_parent_elem)) THEN
    1008             :                   ! using ln(rnd)/-(beta_i-beta_j) < U_j-U_i)
    1009             :                   tmp_prob = compute_estimated_prob(elem_old=tmp_parent_elem, &
    1010             :                                                     elem_new=act_elem, &
    1011             :                                                     E_classical_diff=0.0_dp, &
    1012             :                                                     rnd_nr=tmp_pt_ptr%gt_elem%rnd_nr, &
    1013             :                                                     beta=1.0_dp/(kB*(tmc_env%params%Temp(tmp_pt_ptr%gt_elem%mv_conf) - &
    1014             :                                                                      tmc_env%params%Temp(tmp_pt_ptr%gt_elem%mv_conf + 1))), &
    1015           0 :                                                     tmc_params=tmc_env%params)
    1016             :                ELSE
    1017           0 :                   tmp_pt_ptr => tmp_pt_ptr%next
    1018           0 :                   CYCLE reference_loop
    1019             :                END IF
    1020             :             ELSE
    1021             :                ! if no swap, use the parent configuration as related one
    1022           0 :                IF (ready_for_update_acc_prob(parent_elem)) THEN
    1023             :                   tmp_prob = compute_estimated_prob(elem_old=parent_elem, &
    1024             :                                                     elem_new=act_elem, &
    1025             :                                                     E_classical_diff=act_elem%e_pot_approx - parent_elem%e_pot_approx, &
    1026             :                                                     rnd_nr=tmp_pt_ptr%gt_elem%rnd_nr, &
    1027             :                                                     beta=1.0_dp/(kB*tmc_env%params%Temp(tmp_pt_ptr%gt_elem%mv_conf)), &
    1028           0 :                                                     tmc_params=tmc_env%params)
    1029             :                END IF
    1030             :             END IF
    1031             :             !successful probability update
    1032           0 :             IF (tmp_prob .GE. 0.0_dp) THEN
    1033           0 :                tmp_pt_ptr%gt_elem%prob_acc = tmp_prob
    1034             :                !-- speculative canceling for the related direction
    1035           0 :                IF (tmc_env%params%SPECULATIVE_CANCELING) &
    1036             :                   CALL search_canceling_elements(pt_elem_in=tmp_pt_ptr%gt_elem, &
    1037           0 :                                                  prob=tmp_pt_ptr%gt_elem%prob_acc, tmc_env=tmc_env)
    1038             :             END IF
    1039             : 
    1040             :             ! get next related global tree pointer
    1041           0 :             tmp_pt_ptr => tmp_pt_ptr%next
    1042             :          END DO reference_loop ! global tree references of actual subtree element
    1043             :       END IF
    1044             :       ! end the timing
    1045           0 :       CALL timestop(handle)
    1046             :    END SUBROUTINE update_prob_gt_node_list
    1047             : 
    1048             : END MODULE tmc_tree_acceptance

Generated by: LCOV version 1.15