LCOV - code coverage report
Current view: top level - src/tmc - tmc_moves.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:9f1e18e) Lines: 370 418 88.5 %
Date: 2021-09-22 20:58:48 Functions: 11 12 91.7 %

          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 different move types are applied
      10             : !> \par History
      11             : !>      11.2012 created [Mandes Schoenherr]
      12             : !> \author Mandes 11/2012
      13             : ! **************************************************************************************************
      14             : 
      15             : MODULE tmc_moves
      16             :    USE cell_types,                      ONLY: cell_type,&
      17             :                                               get_cell,&
      18             :                                               pbc
      19             :    USE cp_log_handling,                 ONLY: cp_to_string
      20             :    USE kinds,                           ONLY: dp
      21             :    USE mathconstants,                   ONLY: pi
      22             :    USE mathlib,                         ONLY: dihedral_angle,&
      23             :                                               rotate_vector
      24             :    USE parallel_rng_types,              ONLY: rng_stream_type
      25             :    USE physcon,                         ONLY: boltzmann,&
      26             :                                               joule
      27             :    USE tmc_calculations,                ONLY: center_of_mass,&
      28             :                                               geometrical_center,&
      29             :                                               get_scaled_cell,&
      30             :                                               nearest_distance
      31             :    USE tmc_move_types,                  ONLY: &
      32             :         mv_type_MD, mv_type_atom_swap, mv_type_atom_trans, mv_type_gausian_adapt, mv_type_mol_rot, &
      33             :         mv_type_mol_trans, mv_type_proton_reorder, mv_type_volume_move, tmc_move_type
      34             :    USE tmc_tree_types,                  ONLY: status_frozen,&
      35             :                                               status_ok,&
      36             :                                               tree_type
      37             :    USE tmc_types,                       ONLY: tmc_atom_type,&
      38             :                                               tmc_param_type
      39             : #include "../base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             : 
      43             :    PRIVATE
      44             : 
      45             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_moves'
      46             : 
      47             :    PUBLIC :: change_pos
      48             :    PUBLIC :: elements_in_new_subbox
      49             : 
      50             :    INTEGER, PARAMETER :: not_selected = 0
      51             :    INTEGER, PARAMETER :: proton_donor = -1
      52             :    INTEGER, PARAMETER :: proton_acceptor = 1
      53             : 
      54             : CONTAINS
      55             : ! **************************************************************************************************
      56             : !> \brief applying the preselected move type
      57             : !> \param tmc_params TMC parameters with dimensions ...
      58             : !> \param move_types ...
      59             : !> \param rng_stream random number stream
      60             : !> \param elem configuration to change
      61             : !> \param mv_conf temperature index for determinig the move size
      62             : !> \param new_subbox flag if new sub box should be crated
      63             : !> \param move_rejected return flag if during configurational change
      64             : !>        configuration should still be accepted (not if e.g. atom/molecule
      65             : !>        leave the sub box
      66             : !> \author Mandes 12.2012
      67             : ! **************************************************************************************************
      68        2047 :    SUBROUTINE change_pos(tmc_params, move_types, rng_stream, elem, mv_conf, &
      69             :                          new_subbox, move_rejected)
      70             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
      71             :       TYPE(tmc_move_type), POINTER                       :: move_types
      72             :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      73             :       TYPE(tree_type), POINTER                           :: elem
      74             :       INTEGER                                            :: mv_conf
      75             :       LOGICAL                                            :: new_subbox, move_rejected
      76             : 
      77             :       INTEGER                                            :: act_nr_elem_mv, counter, d, i, ind, &
      78             :                                                             ind_e, m, nr_molec, nr_sub_box_elem
      79        2047 :       INTEGER, DIMENSION(:), POINTER                     :: mol_in_sb
      80             :       REAL(KIND=dp)                                      :: rnd
      81        2047 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: direction, elem_center
      82             : 
      83        2047 :       NULLIFY (direction, elem_center, mol_in_sb)
      84             : 
      85           0 :       CPASSERT(ASSOCIATED(tmc_params))
      86        2047 :       CPASSERT(ASSOCIATED(move_types))
      87        2047 :       CPASSERT(ASSOCIATED(elem))
      88             : 
      89        2047 :       move_rejected = .FALSE.
      90             : 
      91             :       CALL rng_stream%set(bg=elem%rng_seed(:, :, 1), &
      92        2047 :                           cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3))
      93             : 
      94        2047 :       IF (new_subbox) THEN
      95        1740 :          IF (ALL(tmc_params%sub_box_size .GT. 0.0_dp)) THEN
      96             :             CALL elements_in_new_subbox(tmc_params=tmc_params, &
      97             :                                         rng_stream=rng_stream, elem=elem, &
      98           0 :                                         nr_of_sub_box_elements=nr_sub_box_elem)
      99             :          ELSE
     100      115854 :             elem%elem_stat(:) = status_ok
     101             :          END IF
     102             :       END IF
     103             : 
     104             :       ! at least one atom should be in the sub box
     105       15385 :       CPASSERT(ANY(elem%elem_stat(:) .EQ. status_ok))
     106        2047 :       IF (tmc_params%nr_elem_mv .EQ. 0) THEN
     107             :          ! move all elements (could be all atoms or all molecules)
     108             :          act_nr_elem_mv = 0
     109             :       ELSE
     110         988 :          act_nr_elem_mv = tmc_params%nr_elem_mv
     111             :       END IF
     112             :       !-- select the type of move (looked up in list, using the move type index)
     113             :       !-- for each move type exist single moves of certain number of elements
     114             :       !-- or move of all elements
     115             :       !-- one element is a position or velocity of an atom.
     116             :       !-- Always all dimension are changed.
     117        2047 :       SELECT CASE (elem%move_type)
     118             :       CASE (mv_type_gausian_adapt)
     119             :          ! just for Gaussian Adaptation
     120           0 :          CPABORT("gaussian adaptation is not imlemented yet.")
     121             : !TODO       CALL new_pos_gauss_adapt(acc=ASSOCIATED(elem%parent%acc, elem), &
     122             : !                    pos=elem%pos, covari=elem%frc, pot=elem%potential, &
     123             : !                    step_size=elem%ekin, pos_aver=elem%vel, temp=elem%ekin_before_md, &
     124             : !                    rng_seed=elem%rng_seed, rng_seed_last_acc=last_acc_elem%rng_seed)
     125             :          !-- atom translation
     126             :       CASE (mv_type_atom_trans)
     127        1158 :          IF (act_nr_elem_mv .EQ. 0) &
     128         264 :             act_nr_elem_mv = SIZE(elem%pos)/tmc_params%dim_per_elem
     129        3474 :          ALLOCATE (elem_center(tmc_params%dim_per_elem))
     130        1158 :          i = 1
     131             :          move_elements_loop: DO
     132             :             ! select atom
     133       10043 :             IF (tmc_params%nr_elem_mv .EQ. 0) THEN
     134        9149 :                ind = (i - 1)*(tmc_params%dim_per_elem) + 1
     135             :             ELSE
     136         894 :                rnd = rng_stream%next()
     137             :                ind = tmc_params%dim_per_elem* &
     138         894 :                      INT(rnd*(SIZE(elem%pos)/tmc_params%dim_per_elem)) + 1
     139             :             END IF
     140             :             ! apply move
     141       10043 :             IF (elem%elem_stat(ind) .EQ. status_ok) THEN
     142             :                ! displace atom
     143       23676 :                DO d = 0, tmc_params%dim_per_elem - 1
     144       17757 :                   rnd = rng_stream%next()
     145             :                   elem%pos(ind + d) = elem%pos(ind + d) + (rnd - 0.5)*2.0* &
     146       23676 :                                       move_types%mv_size(mv_type_atom_trans, mv_conf)
     147             :                END DO
     148             :                ! check if new position is in subbox
     149       47352 :                elem_center = elem%pos(ind:ind + tmc_params%dim_per_elem - 1)
     150        5919 :                IF (.NOT. check_pos_in_subbox(pos=elem_center, &
     151             :                                              subbox_center=elem%subbox_center, &
     152             :                                              box_scale=elem%box_scale, tmc_params=tmc_params) &
     153             :                    ) THEN
     154           6 :                   move_rejected = .TRUE.
     155           6 :                   EXIT move_elements_loop
     156             :                END IF
     157             :             ELSE
     158             :                ! element was not in sub box, search new one instead
     159        4124 :                IF (tmc_params%nr_elem_mv .GT. 0) i = i - 1
     160             :             END IF
     161       10037 :             i = i + 1
     162       10037 :             IF (i .GT. act_nr_elem_mv) EXIT move_elements_loop
     163             :          END DO move_elements_loop
     164        1158 :          DEALLOCATE (elem_center)
     165             : 
     166             :          !-- molecule translation
     167             :       CASE (mv_type_mol_trans)
     168        8052 :          nr_molec = MAXVAL(elem%mol(:))
     169             :          ! if all particles should be displaced, set the amount of molecules
     170         207 :          IF (act_nr_elem_mv .EQ. 0) &
     171         201 :             act_nr_elem_mv = nr_molec
     172         621 :          ALLOCATE (mol_in_sb(nr_molec))
     173         621 :          ALLOCATE (elem_center(tmc_params%dim_per_elem))
     174        2812 :          mol_in_sb(:) = status_frozen
     175             :          ! check if any molecule is in sub_box
     176        2812 :          DO m = 1, nr_molec
     177             :             CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
     178        2605 :                                  start_ind=ind, end_ind=ind_e)
     179             :             CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
     180        2605 :                                     center=elem_center)
     181        2605 :             IF (check_pos_in_subbox(pos=elem_center, &
     182             :                                     subbox_center=elem%subbox_center, &
     183             :                                     box_scale=elem%box_scale, tmc_params=tmc_params) &
     184             :                 ) &
     185        4680 :                mol_in_sb(m) = status_ok
     186             :          END DO
     187             :          ! displace the selected amount of molecules
     188         558 :          IF (ANY(mol_in_sb(:) .EQ. status_ok)) THEN
     189         621 :             ALLOCATE (direction(tmc_params%dim_per_elem))
     190        1638 :             counter = 1
     191        1638 :             move_molecule_loop: DO
     192             :                ! select molecule
     193        1638 :                IF (tmc_params%nr_elem_mv .EQ. 0) THEN
     194        1632 :                   m = counter
     195             :                ELSE
     196           6 :                   rnd = rng_stream%next()
     197           6 :                   m = INT(rnd*nr_molec) + 1
     198             :                END IF
     199             :                CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
     200        1638 :                                     start_ind=ind, end_ind=ind_e)
     201             :                ! when "molecule" is single atom, search a new one
     202        1638 :                IF (ind .EQ. ind_e) CYCLE move_molecule_loop
     203             : 
     204             :                ! calculate displacement
     205             :                !  move only molecules, with geom. center in subbox
     206        1638 :                IF (mol_in_sb(m) .EQ. status_ok) THEN
     207             :                   ! calculate displacement
     208        5124 :                   DO d = 1, tmc_params%dim_per_elem
     209        3843 :                      rnd = rng_stream%next()
     210             :                      direction(d) = (rnd - 0.5)*2.0_dp*move_types%mv_size( &
     211        5124 :                                     mv_type_mol_trans, mv_conf)
     212             :                   END DO
     213             :                   ! check if displaced position is still in subbox
     214       10248 :                   elem_center(:) = elem_center(:) + direction(:)
     215        1281 :                   IF (check_pos_in_subbox(pos=elem_center, &
     216             :                                           subbox_center=elem%subbox_center, &
     217             :                                           box_scale=elem%box_scale, tmc_params=tmc_params) &
     218             :                       ) THEN
     219             :                      ! apply move
     220        5134 :                      atom_in_mol_loop: DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
     221       16708 :                         dim_loop: DO d = 0, tmc_params%dim_per_elem - 1
     222       15432 :                            elem%pos(i + d) = elem%pos(i + d) + direction(d + 1)
     223             :                         END DO dim_loop
     224             :                      END DO atom_in_mol_loop
     225             :                   ELSE
     226             :                      ! the whole move is rejected, because one element is outside the subbox
     227           5 :                      move_rejected = .TRUE.
     228           5 :                      EXIT move_molecule_loop
     229             :                   END IF
     230             :                ELSE
     231             :                   ! element was not in sub box, search new one instead
     232         357 :                   IF (tmc_params%nr_elem_mv .GT. 0) counter = counter - 1
     233             :                END IF
     234        1633 :                counter = counter + 1
     235        1633 :                IF (counter .GT. act_nr_elem_mv) EXIT move_molecule_loop
     236             :             END DO move_molecule_loop
     237         207 :             DEALLOCATE (direction)
     238             :          END IF
     239         207 :          DEALLOCATE (elem_center)
     240         207 :          DEALLOCATE (mol_in_sb)
     241             : 
     242             :          !-- molecule rotation
     243             :       CASE (mv_type_mol_rot)
     244       10810 :          nr_molec = MAXVAL(elem%mol(:))
     245         261 :          IF (act_nr_elem_mv .EQ. 0) &
     246         257 :             act_nr_elem_mv = nr_molec
     247         783 :          ALLOCATE (mol_in_sb(nr_molec))
     248         783 :          ALLOCATE (elem_center(tmc_params%dim_per_elem))
     249        3766 :          mol_in_sb(:) = status_frozen
     250             :          ! check if any molecule is in sub_box
     251        3766 :          DO m = 1, nr_molec
     252             :             CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
     253        3505 :                                  start_ind=ind, end_ind=ind_e)
     254             :             CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
     255        3505 :                                     center=elem_center)
     256        3505 :             IF (check_pos_in_subbox(pos=elem_center, &
     257             :                                     subbox_center=elem%subbox_center, &
     258             :                                     box_scale=elem%box_scale, tmc_params=tmc_params) &
     259             :                 ) &
     260        5796 :                mol_in_sb(m) = status_ok
     261             :          END DO
     262             :          ! rotate the selected amount of molecules
     263         961 :          IF (ANY(mol_in_sb(:) .EQ. status_ok)) THEN
     264             :             counter = 1
     265        3125 :             rot_molecule_loop: DO
     266             :                ! select molecule
     267        3125 :                IF (tmc_params%nr_elem_mv .EQ. 0) THEN
     268        3121 :                   m = counter
     269             :                ELSE
     270           4 :                   rnd = rng_stream%next()
     271           4 :                   m = INT(rnd*nr_molec) + 1
     272             :                END IF
     273             :                CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
     274        3125 :                                     start_ind=ind, end_ind=ind_e)
     275             :                ! when "molecule" is single atom, search a new one
     276        3125 :                IF (ind .EQ. ind_e) CYCLE rot_molecule_loop
     277             : 
     278             :                ! apply move
     279        3125 :                IF (mol_in_sb(m) .EQ. status_ok) THEN
     280             :                   CALL do_mol_rot(pos=elem%pos, ind_start=ind, ind_end=ind_e, &
     281             :                                   max_angle=move_types%mv_size( &
     282             :                                   mv_type_mol_rot, mv_conf), &
     283             :                                   move_types=move_types, rng_stream=rng_stream, &
     284        1650 :                                   dim_per_elem=tmc_params%dim_per_elem)
     285             :                   ! update sub box status of single atom
     286        6634 :                   DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
     287       39872 :                      elem_center = elem%pos(i:i + tmc_params%dim_per_elem - 1)
     288        4984 :                      IF (check_pos_in_subbox(pos=elem_center, &
     289             :                                              subbox_center=elem%subbox_center, &
     290             :                                              box_scale=elem%box_scale, tmc_params=tmc_params) &
     291        1650 :                          ) THEN
     292       19772 :                         elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_ok
     293             :                      ELSE
     294         164 :                         elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_frozen
     295             :                      END IF
     296             :                   END DO
     297             :                ELSE
     298             :                   ! element was not in sub box, search new one instead
     299        1475 :                   IF (tmc_params%nr_elem_mv .GT. 0) counter = counter - 1
     300             :                END IF
     301        3125 :                counter = counter + 1
     302        3125 :                IF (counter .GT. act_nr_elem_mv) EXIT rot_molecule_loop
     303             :             END DO rot_molecule_loop
     304             :          END IF
     305         261 :          DEALLOCATE (elem_center)
     306         261 :          DEALLOCATE (mol_in_sb)
     307             : 
     308             :          !-- velocity changes for MD
     309             :          !-- here all velocities are changed
     310             :       CASE (mv_type_MD)
     311           0 :          CPASSERT(ASSOCIATED(tmc_params%atoms))
     312           0 :          change_all_velocities_loop: DO i = 1, SIZE(elem%pos)
     313             :             !-- attention, move type size is in atomic units of velocity
     314           0 :             IF (elem%elem_stat(i) .NE. status_frozen) THEN
     315             :                CALL vel_change(vel=elem%vel(i), &
     316             :                                atom_kind=tmc_params%atoms(INT(i/REAL(tmc_params%dim_per_elem, KIND=dp)) + 1), &
     317             :                                phi=move_types%mv_size(mv_type_MD, 1), & ! TODO parallel tempering move sizes for vel_change
     318             :                                temp=tmc_params%Temp(mv_conf), &
     319             :                                rnd_sign_change=.TRUE., & ! MD_vel_invert, &
     320           0 :                                rng_stream=rng_stream)
     321             :             END IF
     322             :          END DO change_all_velocities_loop
     323             : 
     324             :          !-- proton order and disorder
     325             :          !   a loop of molecules is build an in this loop proton acceptors become proton donators
     326             :          !   Therefor the molecules are rotated along the not involved O-H bond
     327             :       CASE (mv_type_proton_reorder)
     328             :          CALL search_and_do_proton_displace_loop(elem=elem, &
     329             :                                                  short_loop=move_rejected, rng_stream=rng_stream, &
     330          12 :                                                  tmc_params=tmc_params)
     331             : 
     332             :          !-- volume move
     333             :          ! the box is increased or decreased and with it the coordinates
     334             :       CASE (mv_type_volume_move)
     335             :          CALL change_volume(conf=elem, T_ind=mv_conf, move_types=move_types, &
     336             :                             rng_stream=rng_stream, tmc_params=tmc_params, &
     337         224 :                             mv_cen_of_mass=tmc_params%mv_cen_of_mass)
     338             : 
     339             :          !-- atom swap
     340             :          ! two atoms of different types are swapped
     341             :       CASE (mv_type_atom_swap)
     342             :          CALL swap_atoms(conf=elem, move_types=move_types, rng_stream=rng_stream, &
     343         185 :                          tmc_params=tmc_params)
     344             : 
     345             :       CASE DEFAULT
     346             :          CALL cp_abort(__LOCATION__, &
     347             :                        "unknown move type "// &
     348        2047 :                        cp_to_string(elem%move_type))
     349             :       END SELECT
     350             : 
     351             :       CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), &
     352        2047 :                           cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3))
     353             : 
     354        2047 :    END SUBROUTINE change_pos
     355             : 
     356             : ! **************************************************************************************************
     357             : !> \brief gets the index of the first molecule element position and the size
     358             : !> \param tmc_params TMC parameters with dim_per_elem
     359             : !> \param mol_arr array with molecule information (which atom attend which mol)
     360             : !> \param mol the selected molecule number
     361             : !> \param start_ind start index of the first atom in molecule
     362             : !> \param end_ind index of the last atom in molecule
     363             : !> \author Mandes 10.2013
     364             : ! **************************************************************************************************
     365       27693 :    SUBROUTINE get_mol_indeces(tmc_params, mol_arr, mol, start_ind, end_ind)
     366             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     367             :       INTEGER, DIMENSION(:), INTENT(IN), POINTER         :: mol_arr
     368             :       INTEGER, INTENT(IN)                                :: mol
     369             :       INTEGER, INTENT(OUT)                               :: start_ind, end_ind
     370             : 
     371             :       INTEGER                                            :: i
     372             : 
     373       27693 :       start_ind = -1
     374       27693 :       end_ind = -1
     375             : 
     376       27693 :       CPASSERT(ASSOCIATED(mol_arr))
     377     6479417 :       CPASSERT(mol .LE. MAXVAL(mol_arr(:)))
     378             :       ! get start index
     379     3201629 :       loop_start: DO i = 1, SIZE(mol_arr)
     380     3201629 :          IF (mol_arr(i) .EQ. mol) THEN
     381       27693 :             start_ind = i
     382       27693 :             EXIT loop_start
     383             :          END IF
     384             :       END DO loop_start
     385             :       ! get end index
     386     3222274 :       loop_end: DO i = SIZE(mol_arr), i, -1
     387     3222274 :          IF (mol_arr(i) .EQ. mol) THEN
     388       27693 :             end_ind = i
     389       27693 :             EXIT loop_end
     390             :          END IF
     391             :       END DO loop_end
     392             :       ! check if all atoms inbetween attend to molecule
     393      110900 :       CPASSERT(ALL(mol_arr(start_ind:end_ind) .EQ. mol))
     394       27693 :       CPASSERT(start_ind .GT. 0)
     395       27693 :       CPASSERT(end_ind .GT. 0)
     396             :       ! convert to indeces mapped for the position array (multiple dim per atom)
     397       27693 :       start_ind = (start_ind - 1)*tmc_params%dim_per_elem + 1
     398       27693 :       end_ind = (end_ind - 1)*tmc_params%dim_per_elem + 1
     399       27693 :    END SUBROUTINE
     400             : 
     401             : ! **************************************************************************************************
     402             : !> \brief checks if a position is within the sub box
     403             : !>        returns true if position is inside
     404             : !> \param pos array with positions
     405             : !> \param subbox_center actual center of sub box
     406             : !> \param box_scale scaling factors for the cell
     407             : !> \param tmc_params TMC parameters with sub box size and cell
     408             : !> \return ...
     409             : !> \author Mandes 11.2012
     410             : ! **************************************************************************************************
     411       22326 :    FUNCTION check_pos_in_subbox(pos, subbox_center, box_scale, tmc_params) &
     412             :       RESULT(inside)
     413             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pos, subbox_center, box_scale
     414             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     415             :       LOGICAL                                            :: inside
     416             : 
     417             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'check_pos_in_subbox'
     418             : 
     419             :       INTEGER                                            :: handle
     420             :       LOGICAL                                            :: flag
     421       22326 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: pos_tmp
     422             : 
     423       22326 :       CPASSERT(ASSOCIATED(pos))
     424       22326 :       CPASSERT(ASSOCIATED(subbox_center))
     425       22326 :       CPASSERT(ASSOCIATED(box_scale))
     426             :       ! if pressure is defined, no scale should be 0
     427       89304 :       flag = .NOT. ((tmc_params%pressure .GT. 0.0_dp) .AND. (ANY(box_scale .EQ. 0.0_dp)))
     428           0 :       CPASSERT(flag)
     429       22326 :       CPASSERT(SIZE(pos) .EQ. 3)
     430       22326 :       CPASSERT(SIZE(pos) .EQ. SIZE(subbox_center))
     431             : 
     432             :       ! start the timing
     433       22326 :       CALL timeset(routineN, handle)
     434             : 
     435       66978 :       ALLOCATE (pos_tmp(SIZE(pos)))
     436             : 
     437       22326 :       inside = .TRUE.
     438             :       ! return if no subbox is defined
     439       42342 :       IF (.NOT. ANY(tmc_params%sub_box_size(:) .LE. 0.1_dp)) THEN
     440       26688 :          pos_tmp(:) = pos(:) - subbox_center(:)
     441             :          CALL get_scaled_cell(cell=tmc_params%cell, box_scale=box_scale, &
     442        6672 :                               vec=pos_tmp)
     443             :          ! check
     444       33520 :          IF (ANY(pos_tmp(:) .GE. tmc_params%sub_box_size(:)/2.0) .OR. &
     445             :              ANY(pos_tmp(:) .LE. -tmc_params%sub_box_size(:)/2.0)) THEN
     446        6142 :             inside = .FALSE.
     447             :          END IF
     448             :       END IF
     449       22326 :       DEALLOCATE (pos_tmp)
     450             :       ! end the timing
     451       22326 :       CALL timestop(handle)
     452       22326 :    END FUNCTION check_pos_in_subbox
     453             : 
     454             : ! **************************************************************************************************
     455             : !> \brief set a new random sub box center and counte the number of atoms in it
     456             : !> \param tmc_params ...
     457             : !> \param rng_stream ...
     458             : !> \param elem ...
     459             : !> \param nr_of_sub_box_elements ...
     460             : !> \param
     461             : !> \param
     462             : !> \author Mandes 11.2012
     463             : ! **************************************************************************************************
     464         114 :    SUBROUTINE elements_in_new_subbox(tmc_params, rng_stream, elem, &
     465             :                                      nr_of_sub_box_elements)
     466             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     467             :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     468             :       TYPE(tree_type), POINTER                           :: elem
     469             :       INTEGER, INTENT(OUT)                               :: nr_of_sub_box_elements
     470             : 
     471             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'elements_in_new_subbox'
     472             : 
     473             :       INTEGER                                            :: handle, i
     474             :       REAL(KIND=dp)                                      :: rnd
     475             :       REAL(KIND=dp), DIMENSION(3)                        :: box_size
     476          57 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: atom_tmp, center_of_sub_box
     477             : 
     478          57 :       NULLIFY (center_of_sub_box, atom_tmp)
     479             : 
     480           0 :       CPASSERT(ASSOCIATED(tmc_params))
     481          57 :       CPASSERT(ASSOCIATED(elem))
     482             : 
     483             :       ! start the timing
     484          57 :       CALL timeset(routineN, handle)
     485             : 
     486          99 :       IF (ANY(tmc_params%sub_box_size(:) .LE. 0.1_dp)) THEN
     487             :          !CPWARN("try to count elements in sub box without sub box.")
     488       37195 :          elem%elem_stat = status_ok
     489          43 :          nr_of_sub_box_elements = SIZE(elem%elem_stat)
     490             :       ELSE
     491          42 :          ALLOCATE (center_of_sub_box(tmc_params%dim_per_elem))
     492          42 :          ALLOCATE (atom_tmp(tmc_params%dim_per_elem))
     493          14 :          nr_of_sub_box_elements = 0
     494             :          ! -- define the center of the sub box
     495             :          CALL rng_stream%set(bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), &
     496          14 :                              ig=elem%rng_seed(:, :, 3))
     497             : 
     498          14 :          CALL get_cell(cell=tmc_params%cell, abc=box_size)
     499          56 :          DO i = 1, SIZE(tmc_params%sub_box_size)
     500          42 :             rnd = rng_stream%next()
     501          56 :             center_of_sub_box(i) = rnd*box_size(i)
     502             :          END DO
     503         112 :          elem%subbox_center(:) = center_of_sub_box(:)
     504             : 
     505             :          CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), &
     506          14 :                              ig=elem%rng_seed(:, :, 3))
     507             : 
     508             :          ! check all elements if they are in subbox
     509        4046 :          DO i = 1, SIZE(elem%pos), tmc_params%dim_per_elem
     510       32256 :             atom_tmp(:) = elem%pos(i:i + tmc_params%dim_per_elem - 1)
     511        4032 :             IF (check_pos_in_subbox(pos=atom_tmp, &
     512             :                                     subbox_center=center_of_sub_box, box_scale=elem%box_scale, &
     513          14 :                                     tmc_params=tmc_params)) THEN
     514         616 :                elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_ok
     515         154 :                nr_of_sub_box_elements = nr_of_sub_box_elements + 1
     516             :             ELSE
     517       15512 :                elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_frozen
     518             :             END IF
     519             :          END DO
     520          14 :          DEALLOCATE (atom_tmp)
     521          14 :          DEALLOCATE (center_of_sub_box)
     522             :       END IF
     523             :       ! end the timing
     524          57 :       CALL timestop(handle)
     525          57 :    END SUBROUTINE elements_in_new_subbox
     526             : 
     527             : ! **************************************************************************************************
     528             : !> \brief molecule rotation using quaternions
     529             : !> \param pos atom positions
     530             : !> \param ind_start starting index in the array
     531             : !> \param ind_end index of last atom in the array
     532             : !> \param max_angle maximal angle in each direction
     533             : !> \param move_types ...
     534             : !> \param rng_stream ramdon stream
     535             : !> \param dim_per_elem dimension per atom
     536             : !> \author Mandes 11.2012
     537             : ! **************************************************************************************************
     538        1650 :    SUBROUTINE do_mol_rot(pos, ind_start, ind_end, max_angle, move_types, &
     539             :                          rng_stream, dim_per_elem)
     540             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pos
     541             :       INTEGER                                            :: ind_start, ind_end
     542             :       REAL(KIND=dp)                                      :: max_angle
     543             :       TYPE(tmc_move_type), POINTER                       :: move_types
     544             :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     545             :       INTEGER                                            :: dim_per_elem
     546             : 
     547             :       INTEGER                                            :: i
     548             :       REAL(KIND=dp)                                      :: a1, a2, a3, q0, q1, q2, q3, rnd
     549             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: rot
     550        1650 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: elem_center
     551             : 
     552        1650 :       NULLIFY (elem_center)
     553             : 
     554           0 :       CPASSERT(ASSOCIATED(pos))
     555        1650 :       CPASSERT(dim_per_elem .EQ. 3)
     556        1650 :       CPASSERT(ind_start .GT. 0 .AND. ind_start .LT. SIZE(pos))
     557        1650 :       CPASSERT(ind_end .GT. 0 .AND. ind_end .LT. SIZE(pos))
     558        1650 :       CPASSERT(ASSOCIATED(move_types))
     559             :       MARK_USED(move_types)
     560             : 
     561             :       ! calculate rotation matrix (using quanternions)
     562        1650 :       rnd = rng_stream%next()
     563        1650 :       a1 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
     564        1650 :       rnd = rng_stream%next()
     565        1650 :       a2 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
     566        1650 :       rnd = rng_stream%next()
     567        1650 :       a3 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
     568        1650 :       q0 = COS(a2/2)*COS((a1 + a3)/2.0_dp)
     569        1650 :       q1 = SIN(a2/2)*COS((a1 - a3)/2.0_dp)
     570        1650 :       q2 = SIN(a2/2)*SIN((a1 - a3)/2.0_dp)
     571        1650 :       q3 = COS(a2/2)*SIN((a1 + a3)/2.0_dp)
     572             :       rot = RESHAPE((/q0*q0 + q1*q1 - q2*q2 - q3*q3, 2*(q1*q2 - q0*q3), 2*(q1*q3 + q0*q2), &
     573             :                       2*(q1*q2 + q0*q3), q0*q0 - q1*q1 + q2*q2 - q3*q3, 2*(q2*q3 - q0*q1), &
     574       16500 :                       2*(q1*q3 - q0*q2), 2*(q2*q3 + q0*q1), q0*q0 - q1*q1 - q2*q2 + q3*q3/), (/3, 3/))
     575             : 
     576        4950 :       ALLOCATE (elem_center(dim_per_elem))
     577             :       ! calculate geometrical center
     578             :       CALL geometrical_center(pos=pos(ind_start:ind_end + dim_per_elem - 1), &
     579        1650 :                               center=elem_center)
     580             : 
     581             :       ! proceed rotation
     582        1650 :       atom_loop: DO i = ind_start, ind_end + dim_per_elem - 1, dim_per_elem
     583      109648 :          pos(i:i + 2) = MATMUL(pos(i:i + 2) - elem_center(:), rot) + elem_center(:)
     584             :       END DO atom_loop
     585        1650 :       DEALLOCATE (elem_center)
     586        1650 :    END SUBROUTINE do_mol_rot
     587             : 
     588             : ! **************************************************************************************************
     589             : !> \brief velocity change should be gaussian distributed
     590             : !>        around the old velocity with respect to kB*T/m
     591             : !> \param vel velocity of atom (one direction)
     592             : !> \param atom_kind ...
     593             : !> \param phi angle for mixing old with random gaussian distributed velocity
     594             : !>        phi =90 degree -> only gaussian velocity around 0
     595             : !>        phi = 0 degree -> only old velocity (with sign change)
     596             : !> \param temp temperature for gaussian distributed velocity
     597             : !> \param rnd_sign_change if sign of old velocity should change randomly
     598             : !> \param rng_stream random number stream
     599             : !> \author Mandes 11.2012
     600             : ! **************************************************************************************************
     601           0 :    SUBROUTINE vel_change(vel, atom_kind, phi, temp, rnd_sign_change, rng_stream)
     602             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel
     603             :       TYPE(tmc_atom_type)                                :: atom_kind
     604             :       REAL(KIND=dp), INTENT(IN)                          :: phi, temp
     605             :       LOGICAL                                            :: rnd_sign_change
     606             :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     607             : 
     608             :       INTEGER                                            :: d
     609             :       REAL(KIND=dp)                                      :: delta_vel, kB, rnd1, rnd2, rnd3, rnd_g
     610             : 
     611           0 :       kB = boltzmann/joule
     612             : 
     613             :       !phi = move_types%mv_size(mv_type_MD,1) ! TODO parallel tempering move sizes for vel_change
     614             :       ! hence first producing a gaussian random number
     615           0 :       rnd1 = rng_stream%next()
     616           0 :       rnd2 = rng_stream%next()
     617             : 
     618           0 :       rnd_g = SQRT(-2.0_dp*LOG(rnd1))*COS(2.0_dp*PI*rnd2)
     619             :       !we can also produce a second one in the same step:
     620             :       !rnd_g2 = SQRT(-2.0_dp*LOG(rnd1))*SIN(2.0_dp*PI*rnd2)
     621             : 
     622             :       ! adapting the variance with respect to kB*T/m
     623           0 :       delta_vel = SQRT(kB*temp/atom_kind%mass)*rnd_g
     624             :       ! check if TODO random velocity sign change
     625             :       ! using detailed balance, velocity sign changes are necessary,
     626             :       ! which are done randomly and
     627             :       ! can be switched of using MD_vel_invert
     628             :       ! without still the balance condition should be fulfilled
     629             : 
     630           0 :       rnd3 = rng_stream%next()
     631           0 :       IF (rnd3 .GE. 0.5 .AND. rnd_sign_change) THEN
     632             :          d = -1
     633             :       ELSE
     634           0 :          d = 1
     635             :       END IF
     636           0 :       vel = SIN(phi)*delta_vel + COS(phi)*vel*d*1.0_dp
     637           0 :    END SUBROUTINE vel_change
     638             : 
     639             : ! **************************************************************************************************
     640             : !> \brief proton order and disorder (specialized move for ice Ih)
     641             : !>        a loop of molecules is build an
     642             : !>        in this loop proton acceptors become proton donators
     643             : !>        Therefor the molecules are rotated along the not involved O-H bond
     644             : !> \param elem sub tree element with actual positions
     645             : !> \param short_loop return if the a loop shorter than 6 molecules is found
     646             : !>        (should not be in ice structure)
     647             : !> \param rng_stream random number stream
     648             : !> \param tmc_params TMC parameters with numbers of dimensions per element
     649             : !>        number of atoms per molecule
     650             : !> \author Mandes 11.2012
     651             : ! **************************************************************************************************
     652          12 :    SUBROUTINE search_and_do_proton_displace_loop(elem, short_loop, rng_stream, &
     653             :                                                  tmc_params)
     654             :       TYPE(tree_type), POINTER                           :: elem
     655             :       LOGICAL                                            :: short_loop
     656             :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     657             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     658             : 
     659             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'search_and_do_proton_displace_loop'
     660             : 
     661             :       CHARACTER(LEN=1000)                                :: tmp_chr
     662             :       INTEGER                                            :: counter, donor_acceptor, handle, k, mol, &
     663             :                                                             nr_mol
     664          12 :       INTEGER, DIMENSION(:), POINTER                     :: mol_arr
     665             :       REAL(KIND=dp)                                      :: rnd
     666             : 
     667          12 :       NULLIFY (mol_arr)
     668             : 
     669           0 :       CPASSERT(ASSOCIATED(elem))
     670          12 :       CPASSERT(ASSOCIATED(tmc_params))
     671             : 
     672             :       ! start the timing
     673          12 :       CALL timeset(routineN, handle)
     674             : 
     675          12 :       short_loop = .FALSE.
     676          12 :       counter = 0
     677        3468 :       nr_mol = MAXVAL(elem%mol(:))
     678             :       ! ind_arr: one array element for each molecule
     679          36 :       ALLOCATE (mol_arr(nr_mol))
     680        1164 :       mol_arr(:) = -1
     681             :       donor_acceptor = not_selected
     682             :       ! select randomly if neighboring molecule is donor / acceptor
     683          12 :       IF (rng_stream%next() .LT. 0.5_dp) THEN
     684           7 :          donor_acceptor = proton_acceptor
     685             :       ELSE
     686           5 :          donor_acceptor = proton_donor
     687             :       END IF
     688             : 
     689             :       ! first step build loop
     690             :       !  select randomly one atom
     691          12 :       rnd = rng_stream%next()
     692             :       ! the randomly selected first atom
     693          12 :       mol = INT(rnd*nr_mol) + 1
     694          12 :       counter = counter + 1
     695          12 :       mol_arr(counter) = mol
     696             : 
     697             :       ! do until the loop is closed
     698             :       !  (until path connects back to any spot of the path)
     699         162 :       chain_completition_loop: DO
     700         174 :          counter = counter + 1
     701             :          ! find nearest neighbor
     702             :          !  (with same state, in the chain, proton donator or proton accptor)
     703             :          CALL find_nearest_proton_acceptor_donator(elem=elem, mol=mol, &
     704             :                                                    donor_acceptor=donor_acceptor, tmc_params=tmc_params, &
     705         174 :                                                    rng_stream=rng_stream)
     706       15784 :          IF (ANY(mol_arr(:) .EQ. mol)) &
     707             :             EXIT chain_completition_loop
     708         174 :          mol_arr(counter) = mol
     709             :       END DO chain_completition_loop
     710             :       counter = counter - 1 ! last searched element is equal to one other in list
     711             : 
     712             :       ! just take the loop of molecules out of the chain
     713          70 :       DO k = 1, counter
     714          70 :          IF (mol_arr(k) .EQ. mol) &
     715          12 :             EXIT
     716             :       END DO
     717         256 :       mol_arr(1:counter - k + 1) = mol_arr(k:counter)
     718          12 :       counter = counter - k + 1
     719             : 
     720             :       ! check if loop is minimum size of 6 molecules
     721          12 :       IF (counter .LT. 6) THEN
     722             :          CALL cp_warn(__LOCATION__, &
     723             :                       "short proton loop with"//cp_to_string(counter)// &
     724           0 :                       "molecules.")
     725           0 :          tmp_chr = ""
     726           0 :          WRITE (tmp_chr, *) mol_arr(1:counter)
     727           0 :          CPWARN("selected molecules:"//TRIM(tmp_chr))
     728           0 :          short_loop = .TRUE.
     729             :       ENDIF
     730             : 
     731             :       ! rotate the molecule along the not involved O-H bond
     732             :       !   (about the angle in of the neighboring chain elements)
     733             :       CALL rotate_molecules_in_chain(tmc_params=tmc_params, elem=elem, &
     734          12 :                                      mol_arr_in=mol_arr(1:counter), donor_acceptor=donor_acceptor)
     735          12 :       DEALLOCATE (mol_arr)
     736             : 
     737             :       ! end the timing
     738          12 :       CALL timestop(handle)
     739          24 :    END SUBROUTINE search_and_do_proton_displace_loop
     740             : 
     741             : ! **************************************************************************************************
     742             : !> \brief searches the next (first atom of) neighboring molecule
     743             : !>        which is proton donor / acceptor
     744             : !> \param elem sub tree element with actual positions
     745             : !> \param mol (in_out) actual regarded molecule, which neighbor is searched for
     746             : !> \param donor_acceptor type of searched neighbor
     747             : !>        (proton donor or proton acceptor)
     748             : !> \param tmc_params TMC parameters with numbers of dimensions per element
     749             : !>        number of atoms per molecule
     750             : !> \param rng_stream random number stream
     751             : !> \author Mandes 12.2012
     752             : ! **************************************************************************************************
     753         174 :    SUBROUTINE find_nearest_proton_acceptor_donator(elem, mol, donor_acceptor, &
     754             :                                                    tmc_params, rng_stream)
     755             :       TYPE(tree_type), POINTER                           :: elem
     756             :       INTEGER                                            :: mol, donor_acceptor
     757             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     758             :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
     759             : 
     760             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'find_nearest_proton_acceptor_donator'
     761             : 
     762             :       INTEGER                                            :: handle, ind, ind_e, ind_n, mol_tmp, &
     763             :                                                             nr_mol
     764             :       INTEGER, DIMENSION(2)                              :: neighbor_mol
     765             :       REAL(KIND=dp)                                      :: dist_tmp, rnd
     766         174 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: distH1, distH2, distO
     767             : 
     768         174 :       NULLIFY (distO, distH1, distH2)
     769           0 :       CPASSERT(ASSOCIATED(elem))
     770         174 :       CPASSERT(ASSOCIATED(tmc_params))
     771             : 
     772             :       ! start the timing
     773         174 :       CALL timeset(routineN, handle)
     774             : 
     775       50286 :       nr_mol = MAXVAL(elem%mol)
     776         522 :       ALLOCATE (distO(nr_mol))
     777         522 :       ALLOCATE (distH1(nr_mol))
     778         522 :       ALLOCATE (distH2(nr_mol))
     779             :       !-- initialize the distances to huge values
     780             :       ! distance of nearest proton of certain molecule to preselected O
     781       16878 :       distO(:) = HUGE(distO(1))
     782             :       ! distance of (first) proton of preselected molecule to certain molecule
     783       16878 :       distH1(:) = HUGE(distH1(1))
     784             :       ! distance of (second) proton of preselected molecule to certain molecule
     785       16878 :       distH2(:) = HUGE(distH2(1))
     786             : 
     787             :       ! get the indices of the old O atom (assuming the first atom of the molecule the first atom)
     788             :       CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=mol, &
     789         174 :                            start_ind=ind, end_ind=ind_e)
     790             : 
     791             :       ! calculate distances to all molecules
     792       16878 :       list_distances: DO mol_tmp = 1, nr_mol
     793       16704 :          IF (mol_tmp .EQ. mol) CYCLE list_distances
     794             :          ! index of the molecule (the O atom)
     795             :          ! assume the first atom of the molecule the first atom
     796             :          CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, &
     797       16530 :                               mol=mol_tmp, start_ind=ind_n, end_ind=ind_e)
     798             :          ! check if selected molecule is water respectively consists of 3 atoms
     799       16530 :          IF (MOD(ind_e - ind_n, 3) .GT. 0) THEN
     800             :             CALL cp_warn(__LOCATION__, &
     801             :                          "selected a molecule with more than 3 atoms, "// &
     802           0 :                          "the proton reordering does not support, skip molecule")
     803           0 :             CYCLE list_distances
     804             :          END IF
     805       16530 :          IF (donor_acceptor .EQ. proton_acceptor) THEN
     806        9785 :             IF (check_donor_acceptor(elem=elem, i_orig=ind, i_neighbor=ind_n, &
     807             :                                      tmc_params=tmc_params) .EQ. proton_acceptor) THEN
     808             :                !distance of fist proton to certain O
     809             :                distH1(mol_tmp) = nearest_distance( &
     810             :                                  x1=elem%pos(ind + tmc_params%dim_per_elem: &
     811             :                                              ind + 2*tmc_params%dim_per_elem - 1), &
     812             :                                  x2=elem%pos(ind_n:ind_n + tmc_params%dim_per_elem - 1), &
     813        4979 :                                  cell=tmc_params%cell, box_scale=elem%box_scale)
     814             :                !distance of second proton to certain O
     815             :                distH2(mol_tmp) = nearest_distance( &
     816             :                                  x1=elem%pos(ind + 2*tmc_params%dim_per_elem: &
     817             :                                              ind + 3*tmc_params%dim_per_elem - 1), &
     818             :                                  x2=elem%pos(ind_n:ind_n + tmc_params%dim_per_elem - 1), &
     819        4979 :                                  cell=tmc_params%cell, box_scale=elem%box_scale)
     820             :             END IF
     821             :          END IF
     822             :          !check for neighboring proton donors
     823       33234 :          IF (donor_acceptor .EQ. proton_donor) THEN
     824        6745 :             IF (check_donor_acceptor(elem=elem, i_orig=ind, i_neighbor=ind_n, &
     825             :                                      tmc_params=tmc_params) .EQ. proton_donor) THEN
     826             :                !distance of selected O to all first protons of other melecules
     827             :                distO(mol_tmp) = nearest_distance( &
     828             :                                 x1=elem%pos(ind:ind + tmc_params%dim_per_elem - 1), &
     829             :                                 x2=elem%pos(ind_n + tmc_params%dim_per_elem: &
     830             :                                             ind_n + 2*tmc_params%dim_per_elem - 1), &
     831        3315 :                                 cell=tmc_params%cell, box_scale=elem%box_scale)
     832             :                dist_tmp = nearest_distance( &
     833             :                           x1=elem%pos(ind:ind + tmc_params%dim_per_elem - 1), &
     834             :                           x2=elem%pos(ind_n + 2*tmc_params%dim_per_elem: &
     835             :                                       ind_n + 3*tmc_params%dim_per_elem - 1), &
     836        3315 :                           cell=tmc_params%cell, box_scale=elem%box_scale)
     837        3315 :                IF (dist_tmp .LT. distO(mol_tmp)) distO(mol_tmp) = dist_tmp
     838             :             END IF
     839             :          END IF
     840             :       END DO list_distances
     841             : 
     842         174 :       mol_tmp = 1
     843             :       ! select the nearest neighbors
     844             :       !check for neighboring proton acceptors
     845         174 :       IF (donor_acceptor .EQ. proton_acceptor) THEN
     846       10094 :          neighbor_mol(mol_tmp) = MINLOC(distH1(:), 1)
     847       10094 :          neighbor_mol(mol_tmp + 1) = MINLOC(distH2(:), 1)
     848             :          ! if both smallest distances points to the shortest molecule search also the second next shortest distance
     849         103 :          IF (neighbor_mol(mol_tmp) .EQ. neighbor_mol(mol_tmp + 1)) THEN
     850           0 :             distH1(neighbor_mol(mol_tmp)) = HUGE(distH1(1))
     851           0 :             distH2(neighbor_mol(mol_tmp + 1)) = HUGE(distH2(1))
     852           0 :             IF (MINVAL(distH1(:), 1) .LT. MINVAL(distH2(:), 1)) THEN
     853           0 :                neighbor_mol(mol_tmp) = MINLOC(distH1(:), 1)
     854             :             ELSE
     855           0 :                neighbor_mol(mol_tmp + 1) = MINLOC(distH2(:), 1)
     856             :             END IF
     857             :          END IF
     858         103 :          mol_tmp = mol_tmp + 2
     859             :       END IF
     860             : 
     861             :       !check for neighboring proton donors
     862         174 :       IF (donor_acceptor .EQ. proton_donor) THEN
     863        6958 :          neighbor_mol(mol_tmp) = MINLOC(distO(:), 1)
     864          71 :          distO(neighbor_mol(mol_tmp)) = HUGE(distO(1))
     865        6958 :          neighbor_mol(mol_tmp + 1) = MINLOC(distO(:), 1)
     866             :       END IF
     867             : 
     868             :       ! select randomly the next neighboring molecule
     869         174 :       rnd = rng_stream%next()
     870             :       ! the randomly selected atom: return value!
     871         174 :       mol_tmp = neighbor_mol(INT(rnd*SIZE(neighbor_mol(:))) + 1)
     872         174 :       mol = mol_tmp
     873             : 
     874         174 :       DEALLOCATE (distO)
     875         174 :       DEALLOCATE (distH1)
     876         174 :       DEALLOCATE (distH2)
     877             : 
     878             :       ! end the timing
     879         174 :       CALL timestop(handle)
     880         522 :    END SUBROUTINE find_nearest_proton_acceptor_donator
     881             : 
     882             : ! **************************************************************************************************
     883             : !> \brief checks if neighbor of the selected/orig element
     884             : !>        is a proron donator or acceptor
     885             : !> \param elem ...
     886             : !> \param i_orig ...
     887             : !> \param i_neighbor ...
     888             : !> \param tmc_params ...
     889             : !> \return ...
     890             : !> \author Mandes 11.2012
     891             : ! **************************************************************************************************
     892       16530 :    FUNCTION check_donor_acceptor(elem, i_orig, i_neighbor, tmc_params) &
     893             :       RESULT(donor_acceptor)
     894             :       TYPE(tree_type), POINTER                           :: elem
     895             :       INTEGER                                            :: i_orig, i_neighbor
     896             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     897             :       INTEGER                                            :: donor_acceptor
     898             : 
     899             :       REAL(KIND=dp), DIMENSION(4)                        :: distances
     900             : 
     901       16530 :       CPASSERT(ASSOCIATED(elem))
     902       16530 :       CPASSERT(i_orig .GE. 1 .AND. i_orig .LE. SIZE(elem%pos))
     903       16530 :       CPASSERT(i_neighbor .GE. 1 .AND. i_neighbor .LE. SIZE(elem%pos))
     904       16530 :       CPASSERT(ASSOCIATED(tmc_params))
     905             : 
     906             :       ! 1. proton of orig with neighbor O
     907             :       distances(1) = nearest_distance( &
     908             :                      x1=elem%pos(i_neighbor:i_neighbor + tmc_params%dim_per_elem - 1), &
     909             :                      x2=elem%pos(i_orig + tmc_params%dim_per_elem: &
     910             :                                  i_orig + 2*tmc_params%dim_per_elem - 1), &
     911       16530 :                      cell=tmc_params%cell, box_scale=elem%box_scale)
     912             :       ! 2. proton of orig with neighbor O
     913             :       distances(2) = nearest_distance( &
     914             :                      x1=elem%pos(i_neighbor:i_neighbor + tmc_params%dim_per_elem - 1), &
     915             :                      x2=elem%pos(i_orig + 2*tmc_params%dim_per_elem: &
     916             :                                  i_orig + 3*tmc_params%dim_per_elem - 1), &
     917       16530 :                      cell=tmc_params%cell, box_scale=elem%box_scale)
     918             :       ! 1. proton of neighbor with orig O
     919             :       distances(3) = nearest_distance( &
     920             :                      x1=elem%pos(i_orig:i_orig + tmc_params%dim_per_elem - 1), &
     921             :                      x2=elem%pos(i_neighbor + tmc_params%dim_per_elem: &
     922             :                                  i_neighbor + 2*tmc_params%dim_per_elem - 1), &
     923       16530 :                      cell=tmc_params%cell, box_scale=elem%box_scale)
     924             :       ! 2. proton of neigbor with orig O
     925             :       distances(4) = nearest_distance( &
     926             :                      x1=elem%pos(i_orig:i_orig + tmc_params%dim_per_elem - 1), &
     927             :                      x2=elem%pos(i_neighbor + 2*tmc_params%dim_per_elem: &
     928             :                                  i_neighbor + 3*tmc_params%dim_per_elem - 1), &
     929       16530 :                      cell=tmc_params%cell, box_scale=elem%box_scale)
     930             : 
     931       99180 :       IF (MINLOC(distances(:), 1) .LE. 2) THEN
     932             :          donor_acceptor = proton_acceptor
     933             :       ELSE
     934        8121 :          donor_acceptor = proton_donor
     935             :       END IF
     936       16530 :    END FUNCTION check_donor_acceptor
     937             : 
     938             : ! **************************************************************************************************
     939             : !> \brief rotates all the molecules in the chain
     940             : !>        the protons were flipped from the donor to the acceptor
     941             : !> \param tmc_params TMC environment parameters
     942             : !> \param elem sub tree element the pos of the molecules in chain should be
     943             : !>        changed by rotating
     944             : !> \param mol_arr_in array of indeces of molecules, should be rotated
     945             : !> \param donor_acceptor gives the direction of rotation
     946             : !> \author Mandes 11.2012
     947             : ! **************************************************************************************************
     948          12 :    SUBROUTINE rotate_molecules_in_chain(tmc_params, elem, mol_arr_in, &
     949             :                                         donor_acceptor)
     950             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     951             :       TYPE(tree_type), POINTER                           :: elem
     952             :       INTEGER, DIMENSION(:)                              :: mol_arr_in
     953             :       INTEGER                                            :: donor_acceptor
     954             : 
     955             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'rotate_molecules_in_chain'
     956             : 
     957             :       INTEGER                                            :: H_offset, handle, i, ind
     958          12 :       INTEGER, DIMENSION(:), POINTER                     :: ind_arr
     959             :       REAL(KIND=dp)                                      :: dihe_angle, dist_near, tmp
     960             :       REAL(KIND=dp), DIMENSION(3)                        :: rot_axis, tmp_1, tmp_2, vec_1O, &
     961             :                                                             vec_2H_f, vec_2H_m, vec_2O, vec_3O, &
     962             :                                                             vec_4O, vec_rotated
     963             :       TYPE(cell_type), POINTER                           :: tmp_cell
     964             : 
     965          12 :       NULLIFY (ind_arr, tmp_cell)
     966             : 
     967           0 :       CPASSERT(ASSOCIATED(tmc_params))
     968          12 :       CPASSERT(ASSOCIATED(elem))
     969             : 
     970             :       ! start the timing
     971          12 :       CALL timeset(routineN, handle)
     972             : 
     973          36 :       ALLOCATE (ind_arr(0:SIZE(mol_arr_in) + 1))
     974         128 :       DO i = 1, SIZE(mol_arr_in)
     975             :          CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, &
     976             :                               mol=mol_arr_in(i), &
     977         128 :                               start_ind=ind_arr(i), end_ind=ind)
     978             :       END DO
     979          12 :       ind_arr(0) = ind_arr(SIZE(ind_arr) - 2)
     980          12 :       ind_arr(SIZE(ind_arr) - 1) = ind_arr(1)
     981             : 
     982             :       ! get the scaled cell
     983          12 :       ALLOCATE (tmp_cell)
     984             :       CALL get_scaled_cell(cell=tmc_params%cell, box_scale=elem%box_scale, &
     985          12 :                            scaled_cell=tmp_cell)
     986             : 
     987             :       ! rotate single molecules
     988         128 :       DO i = 1, SIZE(ind_arr) - 2
     989             :          ! the 3 O atoms
     990         464 :          vec_1O(:) = elem%pos(ind_arr(i - 1):ind_arr(i - 1) + tmc_params%dim_per_elem - 1)
     991         464 :          vec_2O(:) = elem%pos(ind_arr(i):ind_arr(i) + tmc_params%dim_per_elem - 1)
     992         464 :          vec_3O(:) = elem%pos(ind_arr(i + 1):ind_arr(i + 1) + tmc_params%dim_per_elem - 1)
     993             :          ! the H atoms
     994             :          ! distinguished between the one fixed (rotation axis with 2 O)
     995             :          ! and the moved one
     996             :          ! if true the first H atom is between the O atoms
     997         116 :          IF (nearest_distance( &
     998             :              x1=elem%pos(ind_arr(i + donor_acceptor): &
     999             :                          ind_arr(i + donor_acceptor) + tmc_params%dim_per_elem - 1), &
    1000             :              x2=elem%pos(ind_arr(i) + tmc_params%dim_per_elem: &
    1001             :                          ind_arr(i) + 2*tmc_params%dim_per_elem - 1), &
    1002             :              cell=tmc_params%cell, box_scale=elem%box_scale) &
    1003             :              .LT. &
    1004             :              nearest_distance( &
    1005             :              x1=elem%pos(ind_arr(i + donor_acceptor): &
    1006             :                          ind_arr(i + donor_acceptor) + tmc_params%dim_per_elem - 1), &
    1007             :              x2=elem%pos(ind_arr(i) + 2*tmc_params%dim_per_elem: &
    1008             :                          ind_arr(i) + 3*tmc_params%dim_per_elem - 1), &
    1009             :              cell=tmc_params%cell, box_scale=elem%box_scale) &
    1010             :              ) THEN
    1011             :             vec_2H_m = elem%pos(ind_arr(i) + tmc_params%dim_per_elem: &
    1012         276 :                                 ind_arr(i) + 2*tmc_params%dim_per_elem - 1)
    1013             :             vec_2H_f = elem%pos(ind_arr(i) + 2*tmc_params%dim_per_elem: &
    1014         276 :                                 ind_arr(i) + 3*tmc_params%dim_per_elem - 1)
    1015             :             H_offset = 1
    1016             :          ELSE
    1017             :             vec_2H_f = elem%pos(ind_arr(i) + tmc_params%dim_per_elem: &
    1018         188 :                                 ind_arr(i) + 2*tmc_params%dim_per_elem - 1)
    1019             :             vec_2H_m = elem%pos(ind_arr(i) + 2*tmc_params%dim_per_elem: &
    1020         188 :                                 ind_arr(i) + 3*tmc_params%dim_per_elem - 1)
    1021             :             H_offset = 2
    1022             :          END IF
    1023             : 
    1024          12 :          IF (.TRUE.) THEN !TODO find a better switch for the pauling model
    1025             : 
    1026             :             ! do rotation (NOT pauling model)
    1027         464 :             tmp_1 = pbc(vec_2O - vec_1O, tmp_cell)
    1028         464 :             tmp_2 = pbc(vec_3O - vec_2H_f, tmp_cell)
    1029             : 
    1030         464 :             dihe_angle = donor_acceptor*dihedral_angle(tmp_1, vec_2H_f - vec_2O, tmp_2)
    1031         464 :             DO ind = ind_arr(i), ind_arr(i) + tmc_params%dim_per_elem*3 - 1, tmc_params%dim_per_elem
    1032             :                ! set rotation vector
    1033             :                !vec_rotated = rotate_vector(vec_2H_m-vec_2O, dihe_angle, vec_2H_f-vec_2O)
    1034             :                vec_rotated = rotate_vector(elem%pos(ind: &
    1035             :                                                     ind + tmc_params%dim_per_elem - 1) - vec_2O, &
    1036        2436 :                                            dihe_angle, vec_2H_f - vec_2O)
    1037             : 
    1038             :                ! set new position
    1039             :                !elem%pos(ind_arr(i)+H_offset*dim_per_elem:ind_arr(i)+(H_offset+1)*dim_per_elem-1) = vec_2O+vec_rotated
    1040        1508 :                elem%pos(ind:ind + tmc_params%dim_per_elem - 1) = vec_2O + vec_rotated
    1041             :             END DO
    1042             :          ELSE
    1043             :             ! using the pauling model
    1044             :             !  (see Aragones and Vega: Dielectric constant of ices...)
    1045             :             ! the rotation axis is defined using the 4th not involved O
    1046             :             !  (next to the not involved H)
    1047             :             ! O atom next to not involved proton for axis calculation
    1048             :             dist_near = HUGE(dist_near)
    1049             :             search_O_loop: DO ind = 1, SIZE(elem%pos), &
    1050             :                tmc_params%dim_per_elem*3
    1051             :                IF (ind .EQ. ind_arr(i)) CYCLE search_O_loop
    1052             :                tmp = nearest_distance(x1=vec_2H_f, &
    1053             :                                       x2=elem%pos(ind:ind + tmc_params%dim_per_elem - 1), &
    1054             :                                       cell=tmc_params%cell, box_scale=elem%box_scale)
    1055             :                IF (dist_near .GT. tmp) THEN
    1056             :                   dist_near = tmp
    1057             :                   vec_4O = elem%pos(ind:ind + tmc_params%dim_per_elem - 1)
    1058             :                END IF
    1059             :             END DO search_O_loop
    1060             :             rot_axis = pbc(-vec_2O(:) + vec_4O(:), tmp_cell)
    1061             :             tmp_1 = pbc(vec_2O - vec_1O, tmp_cell)
    1062             :             tmp_2 = pbc(vec_3O - vec_4O, tmp_cell)
    1063             :             dihe_angle = donor_acceptor*dihedral_angle(tmp_1, rot_axis, tmp_2)
    1064             :             vec_rotated = rotate_vector(vec_2H_m - vec_2O, dihe_angle, rot_axis)
    1065             :             ! set new position
    1066             :             elem%pos(ind_arr(i) + H_offset*tmc_params%dim_per_elem: &
    1067             :                      ind_arr(i) + (H_offset + 1)*tmc_params%dim_per_elem - 1) &
    1068             :                = vec_2O + vec_rotated
    1069             :             vec_rotated = rotate_vector(vec_2H_f - vec_2O, dihe_angle, rot_axis)
    1070             :             IF (H_offset .EQ. 1) THEN
    1071             :                H_offset = 2
    1072             :             ELSE
    1073             :                H_offset = 1
    1074             :             END IF
    1075             :             elem%pos(ind_arr(i) + H_offset*tmc_params%dim_per_elem: &
    1076             :                      ind_arr(i) + (H_offset + 1)*tmc_params%dim_per_elem - 1) &
    1077             :                = vec_2O + vec_rotated
    1078             :          END IF
    1079             :       END DO
    1080          12 :       DEALLOCATE (tmp_cell)
    1081          12 :       DEALLOCATE (ind_arr)
    1082             :       ! end the timing
    1083          12 :       CALL timestop(handle)
    1084          24 :    END SUBROUTINE rotate_molecules_in_chain
    1085             : 
    1086             : ! **************************************************************************************************
    1087             : !> \brief volume move, the box size is increased or decreased,
    1088             : !>        using the mv_size a the factor.
    1089             : !>        the coordinated are scaled moleculewise
    1090             : !>        (the is moved like the center of mass is moves)
    1091             : !> \param conf configuration to change with positions
    1092             : !> \param T_ind temperature index, to select the correct temperature
    1093             : !>        for move size
    1094             : !> \param move_types ...
    1095             : !> \param rng_stream random number generator stream
    1096             : !> \param tmc_params TMC parameters with e.g. dimensions of atoms and molecules
    1097             : !> \param mv_cen_of_mass ...
    1098             : !> \author Mandes 11.2012
    1099             : ! **************************************************************************************************
    1100         224 :    SUBROUTINE change_volume(conf, T_ind, move_types, rng_stream, tmc_params, &
    1101             :                             mv_cen_of_mass)
    1102             :       TYPE(tree_type), POINTER                           :: conf
    1103             :       INTEGER                                            :: T_ind
    1104             :       TYPE(tmc_move_type), POINTER                       :: move_types
    1105             :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    1106             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
    1107             :       LOGICAL                                            :: mv_cen_of_mass
    1108             : 
    1109             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'change_volume'
    1110             : 
    1111             :       INTEGER                                            :: atom, dir, handle, ind, ind_e, mol
    1112             :       REAL(KIND=dp)                                      :: rnd, vol
    1113             :       REAL(KIND=dp), DIMENSION(3)                        :: box_length_new, box_length_orig, &
    1114             :                                                             box_scale_old
    1115         224 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: disp, scaling
    1116             : 
    1117         224 :       NULLIFY (scaling, disp)
    1118             : 
    1119           0 :       CPASSERT(ASSOCIATED(conf))
    1120         224 :       CPASSERT(ASSOCIATED(move_types))
    1121         224 :       CPASSERT(ASSOCIATED(tmc_params))
    1122         224 :       CPASSERT(T_ind .GT. 0 .AND. T_ind .LE. tmc_params%nr_temp)
    1123         224 :       CPASSERT(tmc_params%dim_per_elem .EQ. 3)
    1124         224 :       CPASSERT(tmc_params%cell%orthorhombic)
    1125             : 
    1126             :       ! start the timing
    1127         224 :       CALL timeset(routineN, handle)
    1128             : 
    1129         672 :       ALLOCATE (scaling(tmc_params%dim_per_elem))
    1130         672 :       ALLOCATE (disp(tmc_params%dim_per_elem))
    1131             : 
    1132         896 :       box_scale_old(:) = conf%box_scale
    1133             :       ! get the cell vector length of the configuration (before move)
    1134             :       CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, &
    1135         224 :                            abc=box_length_new)
    1136             : 
    1137             :       IF (.FALSE.) THEN
    1138             :          ! the volume move in volume space (dV)
    1139             :          IF (tmc_params%v_isotropic) THEN
    1140             :             CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, &
    1141             :                                  abc=box_length_new, vol=vol)
    1142             :             rnd = rng_stream%next()
    1143             :             vol = vol + (rnd - 0.5_dp)*2.0_dp*move_types%mv_size(mv_type_volume_move, T_ind)
    1144             :             box_length_new(:) = vol**(1/REAL(3, KIND=dp))
    1145             :          ELSE
    1146             :             CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, &
    1147             :                                  abc=box_length_new, vol=vol)
    1148             :             rnd = rng_stream%next()
    1149             :             vol = vol + (rnd - 0.5_dp)*2.0_dp*move_types%mv_size(mv_type_volume_move, T_ind)
    1150             :             rnd = rng_stream%next()
    1151             :             dir = 1 + INT(rnd*3)
    1152             :             box_length_new(dir) = 1.0_dp
    1153             :             box_length_new(dir) = vol/PRODUCT(box_length_new(:))
    1154             :          END IF
    1155             :       ELSE
    1156             :          ! the volume move in box length space (dL)
    1157             :          ! increase / decrease box length in this direction
    1158             :          ! l_n = l_o +- rnd * mv_size
    1159         224 :          IF (tmc_params%v_isotropic) THEN
    1160         224 :             rnd = rng_stream%next()
    1161             :             box_length_new(:) = box_length_new(:) + &
    1162             :                                 (rnd - 0.5_dp)*2.0_dp* &
    1163         896 :                                 move_types%mv_size(mv_type_volume_move, T_ind)
    1164             :          ELSE
    1165             :             ! select a random direction
    1166           0 :             rnd = rng_stream%next()
    1167           0 :             dir = 1 + INT(rnd*3)
    1168           0 :             rnd = rng_stream%next()
    1169             :             box_length_new(dir) = box_length_new(dir) + &
    1170             :                                   (rnd - 0.5_dp)*2.0_dp* &
    1171           0 :                                   move_types%mv_size(mv_type_volume_move, T_ind)
    1172             :          END IF
    1173             :       END IF
    1174             : 
    1175             :       ! get the original box length
    1176         896 :       scaling(:) = 1.0_dp
    1177             :       CALL get_scaled_cell(cell=tmc_params%cell, &
    1178             :                            box_scale=scaling, &
    1179         224 :                            abc=box_length_orig)
    1180             :       ! get the new box scale
    1181         896 :       conf%box_scale(:) = box_length_new(:)/box_length_orig(:)
    1182             :       ! molecule scaling
    1183        1792 :       scaling(:) = conf%box_scale(:)/box_scale_old(:)
    1184             : 
    1185         224 :       IF (mv_cen_of_mass .EQV. .FALSE.) THEN
    1186             :          ! homogene scaling of atomic coordinates
    1187         224 :          DO atom = 1, SIZE(conf%pos), tmc_params%dim_per_elem
    1188             :             conf%pos(atom:atom + tmc_params%dim_per_elem - 1) = &
    1189      182880 :                conf%pos(atom:atom + tmc_params%dim_per_elem - 1)*scaling(:)
    1190             :          END DO
    1191             :       ELSE
    1192           0 :          DO mol = 1, MAXVAL(conf%mol(:))
    1193             :             ! move the molecule related to the molecule center of mass
    1194             :             ! get center of mass
    1195           0 :             CPASSERT(ASSOCIATED(tmc_params%atoms))
    1196             : 
    1197             :             CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=conf%mol, mol=mol, &
    1198           0 :                                  start_ind=ind, end_ind=ind_e)
    1199             :             CALL center_of_mass( &
    1200             :                pos=conf%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
    1201             :                atoms=tmc_params%atoms(INT(ind/REAL(tmc_params%dim_per_elem, KIND=dp)) + 1: &
    1202             :                                       INT(ind_e/REAL(tmc_params%dim_per_elem, KIND=dp)) + 1), &
    1203           0 :                center=disp)
    1204             :             ! calculate the center of mass DISPLACEMENT
    1205           0 :             disp(:) = disp(:)*(scaling(:) - 1.0_dp)
    1206             :             ! displace all atoms of the molecule
    1207           0 :             DO atom = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
    1208             :                conf%pos(atom:atom + tmc_params%dim_per_elem - 1) = &
    1209           0 :                   conf%pos(atom:atom + tmc_params%dim_per_elem - 1) + disp(:)
    1210             :             END DO
    1211             :          END DO
    1212             :       END IF
    1213             : 
    1214         224 :       DEALLOCATE (scaling)
    1215         224 :       DEALLOCATE (disp)
    1216             : 
    1217             :       ! end the timing
    1218         224 :       CALL timestop(handle)
    1219         448 :    END SUBROUTINE change_volume
    1220             : 
    1221             : ! **************************************************************************************************
    1222             : !> \brief volume move, two atoms of different types are swapped, both selected
    1223             : !>        randomly
    1224             : !> \param conf configuration to change with positions
    1225             : !> \param move_types ...
    1226             : !> \param rng_stream random number generator stream
    1227             : !> \param tmc_params TMC parameters with e.g. dimensions of atoms and molecules
    1228             : !> \author Mandes 11.2012
    1229             : ! **************************************************************************************************
    1230         185 :    SUBROUTINE swap_atoms(conf, move_types, rng_stream, tmc_params)
    1231             :       TYPE(tree_type), POINTER                           :: conf
    1232             :       TYPE(tmc_move_type), POINTER                       :: move_types
    1233             :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
    1234             :       TYPE(tmc_param_type), POINTER                      :: tmc_params
    1235             : 
    1236             :       INTEGER                                            :: a_1, a_2, ind_1, ind_2
    1237             :       LOGICAL                                            :: found
    1238         185 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: pos_tmp
    1239             : 
    1240         185 :       CPASSERT(ASSOCIATED(conf))
    1241         185 :       CPASSERT(ASSOCIATED(move_types))
    1242         185 :       CPASSERT(ASSOCIATED(tmc_params))
    1243         185 :       CPASSERT(ASSOCIATED(tmc_params%atoms))
    1244             : 
    1245             :       ! loop until two different atoms are found
    1246             :       atom_search_loop: DO
    1247             :          ! select one atom randomly
    1248             :          a_1 = INT(SIZE(conf%pos)/REAL(tmc_params%dim_per_elem, KIND=dp)* &
    1249         532 :                    rng_stream%next()) + 1
    1250             :          ! select the second atom randomly
    1251             :          a_2 = INT(SIZE(conf%pos)/REAL(tmc_params%dim_per_elem, KIND=dp)* &
    1252         532 :                    rng_stream%next()) + 1
    1253             :          ! check if they have different kinds
    1254         532 :          IF (tmc_params%atoms(a_1)%name .NE. tmc_params%atoms(a_2)%name) THEN
    1255             :             ! if present, check if atoms have different type related to the specified table
    1256         234 :             IF (ASSOCIATED(move_types%atom_lists)) THEN
    1257         338 :                DO ind_1 = 1, SIZE(move_types%atom_lists)
    1258             :                   IF (ANY(move_types%atom_lists(ind_1)%atoms(:) .EQ. &
    1259        1073 :                           tmc_params%atoms(a_1)%name) .AND. &
    1260             :                       ANY(move_types%atom_lists(ind_1)%atoms(:) .EQ. &
    1261          49 :                           tmc_params%atoms(a_2)%name)) THEN
    1262             :                      found = .TRUE.
    1263             :                      EXIT atom_search_loop
    1264             :                   END IF
    1265             :                END DO
    1266             :             ELSE
    1267             :                found = .TRUE.
    1268             :                EXIT atom_search_loop
    1269             :             END IF
    1270             :          END IF
    1271             :       END DO atom_search_loop
    1272             :       IF (found) THEN
    1273             :          ! perform coordinate exchange
    1274         555 :          ALLOCATE (pos_tmp(tmc_params%dim_per_elem))
    1275         185 :          ind_1 = (a_1 - 1)*tmc_params%dim_per_elem + 1
    1276         740 :          pos_tmp(:) = conf%pos(ind_1:ind_1 + tmc_params%dim_per_elem - 1)
    1277         185 :          ind_2 = (a_2 - 1)*tmc_params%dim_per_elem + 1
    1278             :          conf%pos(ind_1:ind_1 + tmc_params%dim_per_elem - 1) = &
    1279        1295 :             conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1)
    1280         740 :          conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1) = pos_tmp(:)
    1281         185 :          DEALLOCATE (pos_tmp)
    1282             :       END IF
    1283         185 :    END SUBROUTINE swap_atoms
    1284             : 
    1285             : END MODULE tmc_moves

Generated by: LCOV version 1.15