tmc_moves.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief different move types are applied
!> \par History
!>      11.2012 created [Mandes Schoenherr]
!> \author Mandes 11/2012
! **************************************************************************************************
 
MODULE tmc_moves
   USE cell_types,                      ONLY: cell_type,&
                                              get_cell,&
                                              pbc
   USE cp_log_handling,                 ONLY: cp_to_string
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: pi
   USE mathlib,                         ONLY: dihedral_angle,&
                                              rotate_vector
   USE parallel_rng_types,              ONLY: rng_stream_type
   USE physcon,                         ONLY: boltzmann,&
                                              joule
   USE tmc_calculations,                ONLY: center_of_mass,&
                                              geometrical_center,&
                                              get_scaled_cell,&
                                              nearest_distance
   USE tmc_move_types,                  ONLY: &
        mv_type_md, mv_type_atom_swap, mv_type_atom_trans, mv_type_gausian_adapt, mv_type_mol_rot, &
        mv_type_mol_trans, mv_type_proton_reorder, mv_type_volume_move, tmc_move_type
   USE tmc_tree_types,                  ONLY: status_frozen,&
                                              status_ok,&
                                              tree_type
   USE tmc_types,                       ONLY: tmc_atom_type,&
                                              tmc_param_type
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_moves'
 
   PUBLIC :: change_pos
   PUBLIC :: elements_in_new_subbox
 
   INTEGER, PARAMETER :: not_selected = 0
   INTEGER, PARAMETER :: proton_donor = -1
   INTEGER, PARAMETER :: proton_acceptor = 1
 
CONTAINS
! **************************************************************************************************
!> \brief applying the preselected move type
!> \param tmc_params TMC parameters with dimensions ...
!> \param move_types ...
!> \param rng_stream random number stream
!> \param elem configuration to change
!> \param mv_conf temperature index for determinig the move size
!> \param new_subbox flag if new sub box should be crated
!> \param move_rejected return flag if during configurational change
!>        configuration should still be accepted (not if e.g. atom/molecule
!>        leave the sub box
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE change_pos(tmc_params, move_types, rng_stream, elem, mv_conf, &
                         new_subbox, move_rejected)
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      TYPE(tmc_move_type), POINTER                       :: move_types
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      TYPE(tree_type), POINTER                           :: elem
      INTEGER                                            :: mv_conf
      LOGICAL                                            :: new_subbox, move_rejected
 
      INTEGER                                            :: act_nr_elem_mv, counter, d, i, ind, &
                                                            ind_e, m, nr_molec, nr_sub_box_elem
      INTEGER, DIMENSION(:), POINTER                     :: mol_in_sb
      REAL(kind=dp)                                      :: rnd
      REAL(kind=dp), DIMENSION(:), POINTER               :: direction, elem_center
 
      NULLIFY (direction, elem_center, mol_in_sb)
 
      cpassert(ASSOCIATED(tmc_params))
      cpassert(ASSOCIATED(move_types))
      cpassert(ASSOCIATED(elem))
 
      move_rejected = .false.
 
      CALL rng_stream%set(bg=elem%rng_seed(:, :, 1), &
                          cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3))
 
      IF (new_subbox) THEN
         IF (all(tmc_params%sub_box_size .GT. 0.0_dp)) THEN
            CALL elements_in_new_subbox(tmc_params=tmc_params, &
                                        rng_stream=rng_stream, elem=elem, &
                                        nr_of_sub_box_elements=nr_sub_box_elem)
         ELSE
            elem%elem_stat(:) = status_ok
         END IF
      END IF
 
      ! at least one atom should be in the sub box
      cpassert(any(elem%elem_stat(:) .EQ. status_ok))
      IF (tmc_params%nr_elem_mv .EQ. 0) THEN
         ! move all elements (could be all atoms or all molecules)
         act_nr_elem_mv = 0
      ELSE
         act_nr_elem_mv = tmc_params%nr_elem_mv
      END IF
      !-- select the type of move (looked up in list, using the move type index)
      !-- for each move type exist single moves of certain number of elements
      !-- or move of all elements
      !-- one element is a position or velocity of an atom.
      !-- Always all dimension are changed.
      SELECT CASE (elem%move_type)
      CASE (mv_type_gausian_adapt)
         ! just for Gaussian Adaptation
         cpabort("gaussian adaptation is not imlemented yet.")
!TODO       CALL new_pos_gauss_adapt(acc=ASSOCIATED(elem%parent%acc, elem), &
!                    pos=elem%pos, covari=elem%frc, pot=elem%potential, &
!                    step_size=elem%ekin, pos_aver=elem%vel, temp=elem%ekin_before_md, &
!                    rng_seed=elem%rng_seed, rng_seed_last_acc=last_acc_elem%rng_seed)
         !-- atom translation
      CASE (mv_type_atom_trans)
         IF (act_nr_elem_mv .EQ. 0) &
            act_nr_elem_mv = SIZE(elem%pos)/tmc_params%dim_per_elem
         ALLOCATE (elem_center(tmc_params%dim_per_elem))
         i = 1
         move_elements_loop: DO
            ! select atom
            IF (tmc_params%nr_elem_mv .EQ. 0) THEN
               ind = (i - 1)*(tmc_params%dim_per_elem) + 1
            ELSE
               rnd = rng_stream%next()
               ind = tmc_params%dim_per_elem* &
                     int(rnd*(SIZE(elem%pos)/tmc_params%dim_per_elem)) + 1
            END IF
            ! apply move
            IF (elem%elem_stat(ind) .EQ. status_ok) THEN
               ! displace atom
               DO d = 0, tmc_params%dim_per_elem - 1
                  rnd = rng_stream%next()
                  elem%pos(ind + d) = elem%pos(ind + d) + (rnd - 0.5)*2.0* &
                                      move_types%mv_size(mv_type_atom_trans, mv_conf)
               END DO
               ! check if new position is in subbox
               elem_center = elem%pos(ind:ind + tmc_params%dim_per_elem - 1)
               IF (.NOT. check_pos_in_subbox(pos=elem_center, &
                                             subbox_center=elem%subbox_center, &
                                             box_scale=elem%box_scale, tmc_params=tmc_params) &
                   ) THEN
                  move_rejected = .true.
                  EXIT move_elements_loop
               END IF
            ELSE
               ! element was not in sub box, search new one instead
               IF (tmc_params%nr_elem_mv .GT. 0) i = i - 1
            END IF
            i = i + 1
            IF (i .GT. act_nr_elem_mv) EXIT move_elements_loop
         END DO move_elements_loop
         DEALLOCATE (elem_center)
 
         !-- molecule translation
      CASE (mv_type_mol_trans)
         nr_molec = maxval(elem%mol(:))
         ! if all particles should be displaced, set the amount of molecules
         IF (act_nr_elem_mv .EQ. 0) &
            act_nr_elem_mv = nr_molec
         ALLOCATE (mol_in_sb(nr_molec))
         ALLOCATE (elem_center(tmc_params%dim_per_elem))
         mol_in_sb(:) = status_frozen
         ! check if any molecule is in sub_box
         DO m = 1, nr_molec
            CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
                                 start_ind=ind, end_ind=ind_e)
            CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
                                    center=elem_center)
            IF (check_pos_in_subbox(pos=elem_center, &
                                    subbox_center=elem%subbox_center, &
                                    box_scale=elem%box_scale, tmc_params=tmc_params) &
                ) &
               mol_in_sb(m) = status_ok
         END DO
         ! displace the selected amount of molecules
         IF (any(mol_in_sb(:) .EQ. status_ok)) THEN
            ALLOCATE (direction(tmc_params%dim_per_elem))
            counter = 1
            move_molecule_loop: DO
               ! select molecule
               IF (tmc_params%nr_elem_mv .EQ. 0) THEN
                  m = counter
               ELSE
                  rnd = rng_stream%next()
                  m = int(rnd*nr_molec) + 1
               END IF
               CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
                                    start_ind=ind, end_ind=ind_e)
               ! when "molecule" is single atom, search a new one
               IF (ind .EQ. ind_e) cycle move_molecule_loop
 
               ! calculate displacement
               !  move only molecules, with geom. center in subbox
               IF (mol_in_sb(m) .EQ. status_ok) THEN
                  ! calculate displacement
                  DO d = 1, tmc_params%dim_per_elem
                     rnd = rng_stream%next()
                     direction(d) = (rnd - 0.5)*2.0_dp*move_types%mv_size( &
                                    mv_type_mol_trans, mv_conf)
                  END DO
                  ! check if displaced position is still in subbox
                  elem_center(:) = elem_center(:) + direction(:)
                  IF (check_pos_in_subbox(pos=elem_center, &
                                          subbox_center=elem%subbox_center, &
                                          box_scale=elem%box_scale, tmc_params=tmc_params) &
                      ) THEN
                     ! apply move
                     atom_in_mol_loop: DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
                        dim_loop: DO d = 0, tmc_params%dim_per_elem - 1
                           elem%pos(i + d) = elem%pos(i + d) + direction(d + 1)
                        END DO dim_loop
                     END DO atom_in_mol_loop
                  ELSE
                     ! the whole move is rejected, because one element is outside the subbox
                     move_rejected = .true.
                     EXIT move_molecule_loop
                  END IF
               ELSE
                  ! element was not in sub box, search new one instead
                  IF (tmc_params%nr_elem_mv .GT. 0) counter = counter - 1
               END IF
               counter = counter + 1
               IF (counter .GT. act_nr_elem_mv) EXIT move_molecule_loop
            END DO move_molecule_loop
            DEALLOCATE (direction)
         END IF
         DEALLOCATE (elem_center)
         DEALLOCATE (mol_in_sb)
 
         !-- molecule rotation
      CASE (mv_type_mol_rot)
         nr_molec = maxval(elem%mol(:))
         IF (act_nr_elem_mv .EQ. 0) &
            act_nr_elem_mv = nr_molec
         ALLOCATE (mol_in_sb(nr_molec))
         ALLOCATE (elem_center(tmc_params%dim_per_elem))
         mol_in_sb(:) = status_frozen
         ! check if any molecule is in sub_box
         DO m = 1, nr_molec
            CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
                                 start_ind=ind, end_ind=ind_e)
            CALL geometrical_center(pos=elem%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
                                    center=elem_center)
            IF (check_pos_in_subbox(pos=elem_center, &
                                    subbox_center=elem%subbox_center, &
                                    box_scale=elem%box_scale, tmc_params=tmc_params) &
                ) &
               mol_in_sb(m) = status_ok
         END DO
         ! rotate the selected amount of molecules
         IF (any(mol_in_sb(:) .EQ. status_ok)) THEN
            counter = 1
            rot_molecule_loop: DO
               ! select molecule
               IF (tmc_params%nr_elem_mv .EQ. 0) THEN
                  m = counter
               ELSE
                  rnd = rng_stream%next()
                  m = int(rnd*nr_molec) + 1
               END IF
               CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=m, &
                                    start_ind=ind, end_ind=ind_e)
               ! when "molecule" is single atom, search a new one
               IF (ind .EQ. ind_e) cycle rot_molecule_loop
 
               ! apply move
               IF (mol_in_sb(m) .EQ. status_ok) THEN
                  CALL do_mol_rot(pos=elem%pos, ind_start=ind, ind_end=ind_e, &
                                  max_angle=move_types%mv_size( &
                                  mv_type_mol_rot, mv_conf), &
                                  move_types=move_types, rng_stream=rng_stream, &
                                  dim_per_elem=tmc_params%dim_per_elem)
                  ! update sub box status of single atom
                  DO i = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
                     elem_center = elem%pos(i:i + tmc_params%dim_per_elem - 1)
                     IF (check_pos_in_subbox(pos=elem_center, &
                                             subbox_center=elem%subbox_center, &
                                             box_scale=elem%box_scale, tmc_params=tmc_params) &
                         ) THEN
                        elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_ok
                     ELSE
                        elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_frozen
                     END IF
                  END DO
               ELSE
                  ! element was not in sub box, search new one instead
                  IF (tmc_params%nr_elem_mv .GT. 0) counter = counter - 1
               END IF
               counter = counter + 1
               IF (counter .GT. act_nr_elem_mv) EXIT rot_molecule_loop
            END DO rot_molecule_loop
         END IF
         DEALLOCATE (elem_center)
         DEALLOCATE (mol_in_sb)
 
         !-- velocity changes for MD
         !-- here all velocities are changed
      CASE (mv_type_md)
         cpassert(ASSOCIATED(tmc_params%atoms))
         change_all_velocities_loop: DO i = 1, SIZE(elem%pos)
            !-- attention, move type size is in atomic units of velocity
            IF (elem%elem_stat(i) .NE. status_frozen) THEN
               CALL vel_change(vel=elem%vel(i), &
                               atom_kind=tmc_params%atoms(int(i/real(tmc_params%dim_per_elem, kind=dp)) + 1), &
                               phi=move_types%mv_size(mv_type_md, 1), & ! TODO parallel tempering move sizes for vel_change
                               temp=tmc_params%Temp(mv_conf), &
                               rnd_sign_change=.true., & ! MD_vel_invert, &
                               rng_stream=rng_stream)
            END IF
         END DO change_all_velocities_loop
 
         !-- proton order and disorder
         !   a loop of molecules is build an in this loop proton acceptors become proton donators
         !   Therefor the molecules are rotated along the not involved O-H bond
      CASE (mv_type_proton_reorder)
         CALL search_and_do_proton_displace_loop(elem=elem, &
                                                 short_loop=move_rejected, rng_stream=rng_stream, &
                                                 tmc_params=tmc_params)
 
         !-- volume move
         ! the box is increased or decreased and with it the coordinates
      CASE (mv_type_volume_move)
         CALL change_volume(conf=elem, t_ind=mv_conf, move_types=move_types, &
                            rng_stream=rng_stream, tmc_params=tmc_params, &
                            mv_cen_of_mass=tmc_params%mv_cen_of_mass)
 
         !-- atom swap
         ! two atoms of different types are swapped
      CASE (mv_type_atom_swap)
         CALL swap_atoms(conf=elem, move_types=move_types, rng_stream=rng_stream, &
                         tmc_params=tmc_params)
 
      CASE DEFAULT
         CALL cp_abort(__location__, &
                       "unknown move type "// &
                       cp_to_string(elem%move_type))
      END SELECT
 
      CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), &
                          cg=elem%rng_seed(:, :, 2), ig=elem%rng_seed(:, :, 3))
 
   END SUBROUTINE change_pos
 
! **************************************************************************************************
!> \brief gets the index of the first molecule element position and the size
!> \param tmc_params TMC parameters with dim_per_elem
!> \param mol_arr array with molecule information (which atom attend which mol)
!> \param mol the selected molecule number
!> \param start_ind start index of the first atom in molecule
!> \param end_ind index of the last atom in molecule
!> \author Mandes 10.2013
! **************************************************************************************************
   SUBROUTINE get_mol_indeces(tmc_params, mol_arr, mol, start_ind, end_ind)
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      INTEGER, DIMENSION(:), INTENT(IN), POINTER         :: mol_arr
      INTEGER, INTENT(IN)                                :: mol
      INTEGER, INTENT(OUT)                               :: start_ind, end_ind
 
      INTEGER                                            :: i
 
      start_ind = -1
      end_ind = -1
 
      cpassert(ASSOCIATED(mol_arr))
      cpassert(mol .LE. maxval(mol_arr(:)))
      ! get start index
      loop_start: DO i = 1, SIZE(mol_arr)
         IF (mol_arr(i) .EQ. mol) THEN
            start_ind = i
            EXIT loop_start
         END IF
      END DO loop_start
      ! get end index
      loop_end: DO i = SIZE(mol_arr), i, -1
         IF (mol_arr(i) .EQ. mol) THEN
            end_ind = i
            EXIT loop_end
         END IF
      END DO loop_end
      ! check if all atoms inbetween attend to molecule
      cpassert(all(mol_arr(start_ind:end_ind) .EQ. mol))
      cpassert(start_ind .GT. 0)
      cpassert(end_ind .GT. 0)
      ! convert to indeces mapped for the position array (multiple dim per atom)
      start_ind = (start_ind - 1)*tmc_params%dim_per_elem + 1
      end_ind = (end_ind - 1)*tmc_params%dim_per_elem + 1
   END SUBROUTINE
 
! **************************************************************************************************
!> \brief checks if a position is within the sub box
!>        returns true if position is inside
!> \param pos array with positions
!> \param subbox_center actual center of sub box
!> \param box_scale scaling factors for the cell
!> \param tmc_params TMC parameters with sub box size and cell
!> \return ...
!> \author Mandes 11.2012
! **************************************************************************************************
   FUNCTION check_pos_in_subbox(pos, subbox_center, box_scale, tmc_params) &
      result(inside)
      REAL(kind=dp), DIMENSION(:), POINTER               :: pos, subbox_center, box_scale
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      LOGICAL                                            :: inside
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'check_pos_in_subbox'
 
      INTEGER                                            :: handle
      LOGICAL                                            :: flag
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: pos_tmp
 
      cpassert(ASSOCIATED(pos))
      cpassert(ASSOCIATED(subbox_center))
      cpassert(ASSOCIATED(box_scale))
      ! if pressure is defined, no scale should be 0
      flag = .NOT. ((tmc_params%pressure .GT. 0.0_dp) .AND. (any(box_scale .EQ. 0.0_dp)))
      cpassert(flag)
      cpassert(SIZE(pos) .EQ. 3)
      cpassert(SIZE(pos) .EQ. SIZE(subbox_center))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      ALLOCATE (pos_tmp(SIZE(pos)))
 
      inside = .true.
      ! return if no subbox is defined
      IF (.NOT. any(tmc_params%sub_box_size(:) .LE. 0.1_dp)) THEN
         pos_tmp(:) = pos(:) - subbox_center(:)
         CALL get_scaled_cell(cell=tmc_params%cell, box_scale=box_scale, &
                              vec=pos_tmp)
         ! check
         IF (any(pos_tmp(:) .GE. tmc_params%sub_box_size(:)/2.0) .OR. &
             any(pos_tmp(:) .LE. -tmc_params%sub_box_size(:)/2.0)) THEN
            inside = .false.
         END IF
      END IF
      DEALLOCATE (pos_tmp)
      ! end the timing
      CALL timestop(handle)
   END FUNCTION check_pos_in_subbox
 
! **************************************************************************************************
!> \brief set a new random sub box center and counte the number of atoms in it
!> \param tmc_params ...
!> \param rng_stream ...
!> \param elem ...
!> \param nr_of_sub_box_elements ...
!> \param
!> \param
!> \author Mandes 11.2012
! **************************************************************************************************
   SUBROUTINE elements_in_new_subbox(tmc_params, rng_stream, elem, &
                                     nr_of_sub_box_elements)
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      TYPE(tree_type), POINTER                           :: elem
      INTEGER, INTENT(OUT)                               :: nr_of_sub_box_elements
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'elements_in_new_subbox'
 
      INTEGER                                            :: handle, i
      REAL(kind=dp)                                      :: rnd
      REAL(kind=dp), DIMENSION(3)                        :: box_size
      REAL(kind=dp), DIMENSION(:), POINTER               :: atom_tmp, center_of_sub_box
 
      NULLIFY (center_of_sub_box, atom_tmp)
 
      cpassert(ASSOCIATED(tmc_params))
      cpassert(ASSOCIATED(elem))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      IF (any(tmc_params%sub_box_size(:) .LE. 0.1_dp)) THEN
         !CPWARN("try to count elements in sub box without sub box.")
         elem%elem_stat = status_ok
         nr_of_sub_box_elements = SIZE(elem%elem_stat)
      ELSE
         ALLOCATE (center_of_sub_box(tmc_params%dim_per_elem))
         ALLOCATE (atom_tmp(tmc_params%dim_per_elem))
         nr_of_sub_box_elements = 0
         ! -- define the center of the sub box
         CALL rng_stream%set(bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), &
                             ig=elem%rng_seed(:, :, 3))
 
         CALL get_cell(cell=tmc_params%cell, abc=box_size)
         DO i = 1, SIZE(tmc_params%sub_box_size)
            rnd = rng_stream%next()
            center_of_sub_box(i) = rnd*box_size(i)
         END DO
         elem%subbox_center(:) = center_of_sub_box(:)
 
         CALL rng_stream%get(bg=elem%rng_seed(:, :, 1), cg=elem%rng_seed(:, :, 2), &
                             ig=elem%rng_seed(:, :, 3))
 
         ! check all elements if they are in subbox
         DO i = 1, SIZE(elem%pos), tmc_params%dim_per_elem
            atom_tmp(:) = elem%pos(i:i + tmc_params%dim_per_elem - 1)
            IF (check_pos_in_subbox(pos=atom_tmp, &
                                    subbox_center=center_of_sub_box, box_scale=elem%box_scale, &
                                    tmc_params=tmc_params)) THEN
               elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_ok
               nr_of_sub_box_elements = nr_of_sub_box_elements + 1
            ELSE
               elem%elem_stat(i:i + tmc_params%dim_per_elem - 1) = status_frozen
            END IF
         END DO
         DEALLOCATE (atom_tmp)
         DEALLOCATE (center_of_sub_box)
      END IF
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE elements_in_new_subbox
 
! **************************************************************************************************
!> \brief molecule rotation using quaternions
!> \param pos atom positions
!> \param ind_start starting index in the array
!> \param ind_end index of last atom in the array
!> \param max_angle maximal angle in each direction
!> \param move_types ...
!> \param rng_stream ramdon stream
!> \param dim_per_elem dimension per atom
!> \author Mandes 11.2012
! **************************************************************************************************
   SUBROUTINE do_mol_rot(pos, ind_start, ind_end, max_angle, move_types, &
                         rng_stream, dim_per_elem)
      REAL(kind=dp), DIMENSION(:), POINTER               :: pos
      INTEGER                                            :: ind_start, ind_end
      REAL(kind=dp)                                      :: max_angle
      TYPE(tmc_move_type), POINTER                       :: move_types
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      INTEGER                                            :: dim_per_elem
 
      INTEGER                                            :: i
      REAL(kind=dp)                                      :: a1, a2, a3, q0, q1, q2, q3, rnd
      REAL(kind=dp), DIMENSION(3, 3)                     :: rot
      REAL(kind=dp), DIMENSION(:), POINTER               :: elem_center
 
      NULLIFY (elem_center)
 
      cpassert(ASSOCIATED(pos))
      cpassert(dim_per_elem .EQ. 3)
      cpassert(ind_start .GT. 0 .AND. ind_start .LT. SIZE(pos))
      cpassert(ind_end .GT. 0 .AND. ind_end .LT. SIZE(pos))
      cpassert(ASSOCIATED(move_types))
      mark_used(move_types)
 
      ! calculate rotation matrix (using quanternions)
      rnd = rng_stream%next()
      a1 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
      rnd = rng_stream%next()
      a2 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
      rnd = rng_stream%next()
      a3 = (rnd - 0.5)*2.0*max_angle !move_types%mv_size(mv_type_mol_rot,mv_conf)
      q0 = cos(a2/2)*cos((a1 + a3)/2.0_dp)
      q1 = sin(a2/2)*cos((a1 - a3)/2.0_dp)
      q2 = sin(a2/2)*sin((a1 - a3)/2.0_dp)
      q3 = cos(a2/2)*sin((a1 + a3)/2.0_dp)
      rot = reshape((/q0*q0 + q1*q1 - q2*q2 - q3*q3, 2*(q1*q2 - q0*q3), 2*(q1*q3 + q0*q2), &
                      2*(q1*q2 + q0*q3), q0*q0 - q1*q1 + q2*q2 - q3*q3, 2*(q2*q3 - q0*q1), &
                      2*(q1*q3 - q0*q2), 2*(q2*q3 + q0*q1), q0*q0 - q1*q1 - q2*q2 + q3*q3/), (/3, 3/))
 
      ALLOCATE (elem_center(dim_per_elem))
      ! calculate geometrical center
      CALL geometrical_center(pos=pos(ind_start:ind_end + dim_per_elem - 1), &
                              center=elem_center)
 
      ! proceed rotation
      atom_loop: DO i = ind_start, ind_end + dim_per_elem - 1, dim_per_elem
         pos(i:i + 2) = matmul(pos(i:i + 2) - elem_center(:), rot) + elem_center(:)
      END DO atom_loop
      DEALLOCATE (elem_center)
   END SUBROUTINE do_mol_rot
 
! **************************************************************************************************
!> \brief velocity change should be gaussian distributed
!>        around the old velocity with respect to kB*T/m
!> \param vel velocity of atom (one direction)
!> \param atom_kind ...
!> \param phi angle for mixing old with random gaussian distributed velocity
!>        phi =90 degree -> only gaussian velocity around 0
!>        phi = 0 degree -> only old velocity (with sign change)
!> \param temp temperature for gaussian distributed velocity
!> \param rnd_sign_change if sign of old velocity should change randomly
!> \param rng_stream random number stream
!> \author Mandes 11.2012
! **************************************************************************************************
   SUBROUTINE vel_change(vel, atom_kind, phi, temp, rnd_sign_change, rng_stream)
      REAL(kind=dp), INTENT(INOUT)                       :: vel
      TYPE(tmc_atom_type)                                :: atom_kind
      REAL(kind=dp), INTENT(IN)                          :: phi, temp
      LOGICAL                                            :: rnd_sign_change
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
      INTEGER                                            :: d
      REAL(kind=dp)                                      :: delta_vel, kb, rnd1, rnd2, rnd3, rnd_g
 
      kb = boltzmann/joule
 
      !phi = move_types%mv_size(mv_type_MD,1) ! TODO parallel tempering move sizes for vel_change
      ! hence first producing a gaussian random number
      rnd1 = rng_stream%next()
      rnd2 = rng_stream%next()
 
      rnd_g = sqrt(-2.0_dp*log(rnd1))*cos(2.0_dp*pi*rnd2)
      !we can also produce a second one in the same step:
      !rnd_g2 = SQRT(-2.0_dp*LOG(rnd1))*SIN(2.0_dp*PI*rnd2)
 
      ! adapting the variance with respect to kB*T/m
      delta_vel = sqrt(kb*temp/atom_kind%mass)*rnd_g
      ! check if TODO random velocity sign change
      ! using detailed balance, velocity sign changes are necessary,
      ! which are done randomly and
      ! can be switched of using MD_vel_invert
      ! without still the balance condition should be fulfilled
 
      rnd3 = rng_stream%next()
      IF (rnd3 .GE. 0.5 .AND. rnd_sign_change) THEN
         d = -1
      ELSE
         d = 1
      END IF
      vel = sin(phi)*delta_vel + cos(phi)*vel*d*1.0_dp
   END SUBROUTINE vel_change
 
! **************************************************************************************************
!> \brief proton order and disorder (specialized move for ice Ih)
!>        a loop of molecules is build an
!>        in this loop proton acceptors become proton donators
!>        Therefor the molecules are rotated along the not involved O-H bond
!> \param elem sub tree element with actual positions
!> \param short_loop return if the a loop shorter than 6 molecules is found
!>        (should not be in ice structure)
!> \param rng_stream random number stream
!> \param tmc_params TMC parameters with numbers of dimensions per element
!>        number of atoms per molecule
!> \author Mandes 11.2012
! **************************************************************************************************
   SUBROUTINE search_and_do_proton_displace_loop(elem, short_loop, rng_stream, &
                                                 tmc_params)
      TYPE(tree_type), POINTER                           :: elem
      LOGICAL                                            :: short_loop
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      TYPE(tmc_param_type), POINTER                      :: tmc_params
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'search_and_do_proton_displace_loop'
 
      CHARACTER(LEN=1000)                                :: tmp_chr
      INTEGER                                            :: counter, donor_acceptor, handle, k, mol, &
                                                            nr_mol
      INTEGER, DIMENSION(:), POINTER                     :: mol_arr
      REAL(kind=dp)                                      :: rnd
 
      NULLIFY (mol_arr)
 
      cpassert(ASSOCIATED(elem))
      cpassert(ASSOCIATED(tmc_params))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      short_loop = .false.
      counter = 0
      nr_mol = maxval(elem%mol(:))
      ! ind_arr: one array element for each molecule
      ALLOCATE (mol_arr(nr_mol))
      mol_arr(:) = -1
      donor_acceptor = not_selected
      ! select randomly if neighboring molecule is donor / acceptor
      IF (rng_stream%next() .LT. 0.5_dp) THEN
         donor_acceptor = proton_acceptor
      ELSE
         donor_acceptor = proton_donor
      END IF
 
      ! first step build loop
      !  select randomly one atom
      rnd = rng_stream%next()
      ! the randomly selected first atom
      mol = int(rnd*nr_mol) + 1
      counter = counter + 1
      mol_arr(counter) = mol
 
      ! do until the loop is closed
      !  (until path connects back to any spot of the path)
      chain_completition_loop: DO
         counter = counter + 1
         ! find nearest neighbor
         !  (with same state, in the chain, proton donator or proton accptor)
         CALL find_nearest_proton_acceptor_donator(elem=elem, mol=mol, &
                                                   donor_acceptor=donor_acceptor, tmc_params=tmc_params, &
                                                   rng_stream=rng_stream)
         IF (any(mol_arr(:) .EQ. mol)) &
            EXIT chain_completition_loop
         mol_arr(counter) = mol
      END DO chain_completition_loop
      counter = counter - 1 ! last searched element is equal to one other in list
 
      ! just take the loop of molecules out of the chain
      DO k = 1, counter
         IF (mol_arr(k) .EQ. mol) &
            EXIT
      END DO
      mol_arr(1:counter - k + 1) = mol_arr(k:counter)
      counter = counter - k + 1
 
      ! check if loop is minimum size of 6 molecules
      IF (counter .LT. 6) THEN
         CALL cp_warn(__location__, &
                      "short proton loop with"//cp_to_string(counter)// &
                      "molecules.")
         tmp_chr = ""
         WRITE (tmp_chr, *) mol_arr(1:counter)
         cpwarn("selected molecules:"//trim(tmp_chr))
         short_loop = .true.
      END IF
 
      ! rotate the molecule along the not involved O-H bond
      !   (about the angle in of the neighboring chain elements)
      CALL rotate_molecules_in_chain(tmc_params=tmc_params, elem=elem, &
                                     mol_arr_in=mol_arr(1:counter), donor_acceptor=donor_acceptor)
      DEALLOCATE (mol_arr)
 
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE search_and_do_proton_displace_loop
 
! **************************************************************************************************
!> \brief searches the next (first atom of) neighboring molecule
!>        which is proton donor / acceptor
!> \param elem sub tree element with actual positions
!> \param mol (in_out) actual regarded molecule, which neighbor is searched for
!> \param donor_acceptor type of searched neighbor
!>        (proton donor or proton acceptor)
!> \param tmc_params TMC parameters with numbers of dimensions per element
!>        number of atoms per molecule
!> \param rng_stream random number stream
!> \author Mandes 12.2012
! **************************************************************************************************
   SUBROUTINE find_nearest_proton_acceptor_donator(elem, mol, donor_acceptor, &
                                                   tmc_params, rng_stream)
      TYPE(tree_type), POINTER                           :: elem
      INTEGER                                            :: mol, donor_acceptor
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'find_nearest_proton_acceptor_donator'
 
      INTEGER                                            :: handle, ind, ind_e, ind_n, mol_tmp, &
                                                            nr_mol
      INTEGER, DIMENSION(2)                              :: neighbor_mol
      REAL(kind=dp)                                      :: dist_tmp, rnd
      REAL(kind=dp), DIMENSION(:), POINTER               :: disth1, disth2, disto
 
      NULLIFY (disto, disth1, disth2)
      cpassert(ASSOCIATED(elem))
      cpassert(ASSOCIATED(tmc_params))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      nr_mol = maxval(elem%mol)
      ALLOCATE (disto(nr_mol))
      ALLOCATE (disth1(nr_mol))
      ALLOCATE (disth2(nr_mol))
      !-- initialize the distances to huge values
      ! distance of nearest proton of certain molecule to preselected O
      disto(:) = huge(disto(1))
      ! distance of (first) proton of preselected molecule to certain molecule
      disth1(:) = huge(disth1(1))
      ! distance of (second) proton of preselected molecule to certain molecule
      disth2(:) = huge(disth2(1))
 
      ! get the indices of the old O atom (assuming the first atom of the molecule the first atom)
      CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, mol=mol, &
                           start_ind=ind, end_ind=ind_e)
 
      ! calculate distances to all molecules
      list_distances: DO mol_tmp = 1, nr_mol
         IF (mol_tmp .EQ. mol) cycle list_distances
         ! index of the molecule (the O atom)
         ! assume the first atom of the molecule the first atom
         CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, &
                              mol=mol_tmp, start_ind=ind_n, end_ind=ind_e)
         ! check if selected molecule is water respectively consists of 3 atoms
         IF (mod(ind_e - ind_n, 3) .GT. 0) THEN
            CALL cp_warn(__location__, &
                         "selected a molecule with more than 3 atoms, "// &
                         "the proton reordering does not support, skip molecule")
            cycle list_distances
         END IF
         IF (donor_acceptor .EQ. proton_acceptor) THEN
            IF (check_donor_acceptor(elem=elem, i_orig=ind, i_neighbor=ind_n, &
                                     tmc_params=tmc_params) .EQ. proton_acceptor) THEN
               !distance of fist proton to certain O
               disth1(mol_tmp) = nearest_distance( &
                                 x1=elem%pos(ind + tmc_params%dim_per_elem: &
                                             ind + 2*tmc_params%dim_per_elem - 1), &
                                 x2=elem%pos(ind_n:ind_n + tmc_params%dim_per_elem - 1), &
                                 cell=tmc_params%cell, box_scale=elem%box_scale)
               !distance of second proton to certain O
               disth2(mol_tmp) = nearest_distance( &
                                 x1=elem%pos(ind + 2*tmc_params%dim_per_elem: &
                                             ind + 3*tmc_params%dim_per_elem - 1), &
                                 x2=elem%pos(ind_n:ind_n + tmc_params%dim_per_elem - 1), &
                                 cell=tmc_params%cell, box_scale=elem%box_scale)
            END IF
         END IF
         !check for neighboring proton donors
         IF (donor_acceptor .EQ. proton_donor) THEN
            IF (check_donor_acceptor(elem=elem, i_orig=ind, i_neighbor=ind_n, &
                                     tmc_params=tmc_params) .EQ. proton_donor) THEN
               !distance of selected O to all first protons of other melecules
               disto(mol_tmp) = nearest_distance( &
                                x1=elem%pos(ind:ind + tmc_params%dim_per_elem - 1), &
                                x2=elem%pos(ind_n + tmc_params%dim_per_elem: &
                                            ind_n + 2*tmc_params%dim_per_elem - 1), &
                                cell=tmc_params%cell, box_scale=elem%box_scale)
               dist_tmp = nearest_distance( &
                          x1=elem%pos(ind:ind + tmc_params%dim_per_elem - 1), &
                          x2=elem%pos(ind_n + 2*tmc_params%dim_per_elem: &
                                      ind_n + 3*tmc_params%dim_per_elem - 1), &
                          cell=tmc_params%cell, box_scale=elem%box_scale)
               IF (dist_tmp .LT. disto(mol_tmp)) disto(mol_tmp) = dist_tmp
            END IF
         END IF
      END DO list_distances
 
      mol_tmp = 1
      ! select the nearest neighbors
      !check for neighboring proton acceptors
      IF (donor_acceptor .EQ. proton_acceptor) THEN
         neighbor_mol(mol_tmp) = minloc(disth1(:), 1)
         neighbor_mol(mol_tmp + 1) = minloc(disth2(:), 1)
         ! if both smallest distances points to the shortest molecule search also the second next shortest distance
         IF (neighbor_mol(mol_tmp) .EQ. neighbor_mol(mol_tmp + 1)) THEN
            disth1(neighbor_mol(mol_tmp)) = huge(disth1(1))
            disth2(neighbor_mol(mol_tmp + 1)) = huge(disth2(1))
            IF (minval(disth1(:), 1) .LT. minval(disth2(:), 1)) THEN
               neighbor_mol(mol_tmp) = minloc(disth1(:), 1)
            ELSE
               neighbor_mol(mol_tmp + 1) = minloc(disth2(:), 1)
            END IF
         END IF
         mol_tmp = mol_tmp + 2
      END IF
 
      !check for neighboring proton donors
      IF (donor_acceptor .EQ. proton_donor) THEN
         neighbor_mol(mol_tmp) = minloc(disto(:), 1)
         disto(neighbor_mol(mol_tmp)) = huge(disto(1))
         neighbor_mol(mol_tmp + 1) = minloc(disto(:), 1)
      END IF
 
      ! select randomly the next neighboring molecule
      rnd = rng_stream%next()
      ! the randomly selected atom: return value!
      mol_tmp = neighbor_mol(int(rnd*SIZE(neighbor_mol(:))) + 1)
      mol = mol_tmp
 
      DEALLOCATE (disto)
      DEALLOCATE (disth1)
      DEALLOCATE (disth2)
 
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE find_nearest_proton_acceptor_donator
 
! **************************************************************************************************
!> \brief checks if neighbor of the selected/orig element
!>        is a proron donator or acceptor
!> \param elem ...
!> \param i_orig ...
!> \param i_neighbor ...
!> \param tmc_params ...
!> \return ...
!> \author Mandes 11.2012
! **************************************************************************************************
   FUNCTION check_donor_acceptor(elem, i_orig, i_neighbor, tmc_params) &
      result(donor_acceptor)
      TYPE(tree_type), POINTER                           :: elem
      INTEGER                                            :: i_orig, i_neighbor
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      INTEGER                                            :: donor_acceptor
 
      REAL(kind=dp), DIMENSION(4)                        :: distances
 
      cpassert(ASSOCIATED(elem))
      cpassert(i_orig .GE. 1 .AND. i_orig .LE. SIZE(elem%pos))
      cpassert(i_neighbor .GE. 1 .AND. i_neighbor .LE. SIZE(elem%pos))
      cpassert(ASSOCIATED(tmc_params))
 
      ! 1. proton of orig with neighbor O
      distances(1) = nearest_distance( &
                     x1=elem%pos(i_neighbor:i_neighbor + tmc_params%dim_per_elem - 1), &
                     x2=elem%pos(i_orig + tmc_params%dim_per_elem: &
                                 i_orig + 2*tmc_params%dim_per_elem - 1), &
                     cell=tmc_params%cell, box_scale=elem%box_scale)
      ! 2. proton of orig with neighbor O
      distances(2) = nearest_distance( &
                     x1=elem%pos(i_neighbor:i_neighbor + tmc_params%dim_per_elem - 1), &
                     x2=elem%pos(i_orig + 2*tmc_params%dim_per_elem: &
                                 i_orig + 3*tmc_params%dim_per_elem - 1), &
                     cell=tmc_params%cell, box_scale=elem%box_scale)
      ! 1. proton of neighbor with orig O
      distances(3) = nearest_distance( &
                     x1=elem%pos(i_orig:i_orig + tmc_params%dim_per_elem - 1), &
                     x2=elem%pos(i_neighbor + tmc_params%dim_per_elem: &
                                 i_neighbor + 2*tmc_params%dim_per_elem - 1), &
                     cell=tmc_params%cell, box_scale=elem%box_scale)
      ! 2. proton of neigbor with orig O
      distances(4) = nearest_distance( &
                     x1=elem%pos(i_orig:i_orig + tmc_params%dim_per_elem - 1), &
                     x2=elem%pos(i_neighbor + 2*tmc_params%dim_per_elem: &
                                 i_neighbor + 3*tmc_params%dim_per_elem - 1), &
                     cell=tmc_params%cell, box_scale=elem%box_scale)
 
      IF (minloc(distances(:), 1) .LE. 2) THEN
         donor_acceptor = proton_acceptor
      ELSE
         donor_acceptor = proton_donor
      END IF
   END FUNCTION check_donor_acceptor
 
! **************************************************************************************************
!> \brief rotates all the molecules in the chain
!>        the protons were flipped from the donor to the acceptor
!> \param tmc_params TMC environment parameters
!> \param elem sub tree element the pos of the molecules in chain should be
!>        changed by rotating
!> \param mol_arr_in array of indeces of molecules, should be rotated
!> \param donor_acceptor gives the direction of rotation
!> \author Mandes 11.2012
! **************************************************************************************************
   SUBROUTINE rotate_molecules_in_chain(tmc_params, elem, mol_arr_in, &
                                        donor_acceptor)
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      TYPE(tree_type), POINTER                           :: elem
      INTEGER, DIMENSION(:)                              :: mol_arr_in
      INTEGER                                            :: donor_acceptor
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'rotate_molecules_in_chain'
 
      INTEGER                                            :: h_offset, handle, i, ind
      INTEGER, DIMENSION(:), POINTER                     :: ind_arr
      REAL(kind=dp)                                      :: dihe_angle, dist_near, tmp
      REAL(kind=dp), DIMENSION(3)                        :: rot_axis, tmp_1, tmp_2, vec_1o, &
                                                            vec_2h_f, vec_2h_m, vec_2o, vec_3o, &
                                                            vec_4o, vec_rotated
      TYPE(cell_type), POINTER                           :: tmp_cell
 
      NULLIFY (ind_arr, tmp_cell)
 
      cpassert(ASSOCIATED(tmc_params))
      cpassert(ASSOCIATED(elem))
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      ALLOCATE (ind_arr(0:SIZE(mol_arr_in) + 1))
      DO i = 1, SIZE(mol_arr_in)
         CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=elem%mol, &
                              mol=mol_arr_in(i), &
                              start_ind=ind_arr(i), end_ind=ind)
      END DO
      ind_arr(0) = ind_arr(SIZE(ind_arr) - 2)
      ind_arr(SIZE(ind_arr) - 1) = ind_arr(1)
 
      ! get the scaled cell
      ALLOCATE (tmp_cell)
      CALL get_scaled_cell(cell=tmc_params%cell, box_scale=elem%box_scale, &
                           scaled_cell=tmp_cell)
 
      ! rotate single molecules
      DO i = 1, SIZE(ind_arr) - 2
         ! the 3 O atoms
         vec_1o(:) = elem%pos(ind_arr(i - 1):ind_arr(i - 1) + tmc_params%dim_per_elem - 1)
         vec_2o(:) = elem%pos(ind_arr(i):ind_arr(i) + tmc_params%dim_per_elem - 1)
         vec_3o(:) = elem%pos(ind_arr(i + 1):ind_arr(i + 1) + tmc_params%dim_per_elem - 1)
         ! the H atoms
         ! distinguished between the one fixed (rotation axis with 2 O)
         ! and the moved one
         ! if true the first H atom is between the O atoms
         IF (nearest_distance( &
             x1=elem%pos(ind_arr(i + donor_acceptor): &
                         ind_arr(i + donor_acceptor) + tmc_params%dim_per_elem - 1), &
             x2=elem%pos(ind_arr(i) + tmc_params%dim_per_elem: &
                         ind_arr(i) + 2*tmc_params%dim_per_elem - 1), &
             cell=tmc_params%cell, box_scale=elem%box_scale) &
             .LT. &
             nearest_distance( &
             x1=elem%pos(ind_arr(i + donor_acceptor): &
                         ind_arr(i + donor_acceptor) + tmc_params%dim_per_elem - 1), &
             x2=elem%pos(ind_arr(i) + 2*tmc_params%dim_per_elem: &
                         ind_arr(i) + 3*tmc_params%dim_per_elem - 1), &
             cell=tmc_params%cell, box_scale=elem%box_scale) &
             ) THEN
            vec_2h_m = elem%pos(ind_arr(i) + tmc_params%dim_per_elem: &
                                ind_arr(i) + 2*tmc_params%dim_per_elem - 1)
            vec_2h_f = elem%pos(ind_arr(i) + 2*tmc_params%dim_per_elem: &
                                ind_arr(i) + 3*tmc_params%dim_per_elem - 1)
            h_offset = 1
         ELSE
            vec_2h_f = elem%pos(ind_arr(i) + tmc_params%dim_per_elem: &
                                ind_arr(i) + 2*tmc_params%dim_per_elem - 1)
            vec_2h_m = elem%pos(ind_arr(i) + 2*tmc_params%dim_per_elem: &
                                ind_arr(i) + 3*tmc_params%dim_per_elem - 1)
            h_offset = 2
         END IF
 
         IF (.true.) THEN !TODO find a better switch for the pauling model
 
            ! do rotation (NOT pauling model)
            tmp_1 = pbc(vec_2o - vec_1o, tmp_cell)
            tmp_2 = pbc(vec_3o - vec_2h_f, tmp_cell)
 
            dihe_angle = donor_acceptor*dihedral_angle(tmp_1, vec_2h_f - vec_2o, tmp_2)
            DO ind = ind_arr(i), ind_arr(i) + tmc_params%dim_per_elem*3 - 1, tmc_params%dim_per_elem
               ! set rotation vector
               !vec_rotated = rotate_vector(vec_2H_m-vec_2O, dihe_angle, vec_2H_f-vec_2O)
               vec_rotated = rotate_vector(elem%pos(ind: &
                                                    ind + tmc_params%dim_per_elem - 1) - vec_2o, &
                                           dihe_angle, vec_2h_f - vec_2o)
 
               ! set new position
               !elem%pos(ind_arr(i)+H_offset*dim_per_elem:ind_arr(i)+(H_offset+1)*dim_per_elem-1) = vec_2O+vec_rotated
               elem%pos(ind:ind + tmc_params%dim_per_elem - 1) = vec_2o + vec_rotated
            END DO
         ELSE
            ! using the pauling model
            !  (see Aragones and Vega: Dielectric constant of ices...)
            ! the rotation axis is defined using the 4th not involved O
            !  (next to the not involved H)
            ! O atom next to not involved proton for axis calculation
            dist_near = huge(dist_near)
            search_o_loop: DO ind = 1, SIZE(elem%pos), &
               tmc_params%dim_per_elem*3
               IF (ind .EQ. ind_arr(i)) cycle search_o_loop
               tmp = nearest_distance(x1=vec_2h_f, &
                                      x2=elem%pos(ind:ind + tmc_params%dim_per_elem - 1), &
                                      cell=tmc_params%cell, box_scale=elem%box_scale)
               IF (dist_near .GT. tmp) THEN
                  dist_near = tmp
                  vec_4o = elem%pos(ind:ind + tmc_params%dim_per_elem - 1)
               END IF
            END DO search_o_loop
            rot_axis = pbc(-vec_2o(:) + vec_4o(:), tmp_cell)
            tmp_1 = pbc(vec_2o - vec_1o, tmp_cell)
            tmp_2 = pbc(vec_3o - vec_4o, tmp_cell)
            dihe_angle = donor_acceptor*dihedral_angle(tmp_1, rot_axis, tmp_2)
            vec_rotated = rotate_vector(vec_2h_m - vec_2o, dihe_angle, rot_axis)
            ! set new position
            elem%pos(ind_arr(i) + h_offset*tmc_params%dim_per_elem: &
                     ind_arr(i) + (h_offset + 1)*tmc_params%dim_per_elem - 1) &
               = vec_2o + vec_rotated
            vec_rotated = rotate_vector(vec_2h_f - vec_2o, dihe_angle, rot_axis)
            IF (h_offset .EQ. 1) THEN
               h_offset = 2
            ELSE
               h_offset = 1
            END IF
            elem%pos(ind_arr(i) + h_offset*tmc_params%dim_per_elem: &
                     ind_arr(i) + (h_offset + 1)*tmc_params%dim_per_elem - 1) &
               = vec_2o + vec_rotated
         END IF
      END DO
      DEALLOCATE (tmp_cell)
      DEALLOCATE (ind_arr)
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE rotate_molecules_in_chain
 
! **************************************************************************************************
!> \brief volume move, the box size is increased or decreased,
!>        using the mv_size a the factor.
!>        the coordinated are scaled moleculewise
!>        (the is moved like the center of mass is moves)
!> \param conf configuration to change with positions
!> \param T_ind temperature index, to select the correct temperature
!>        for move size
!> \param move_types ...
!> \param rng_stream random number generator stream
!> \param tmc_params TMC parameters with e.g. dimensions of atoms and molecules
!> \param mv_cen_of_mass ...
!> \author Mandes 11.2012
! **************************************************************************************************
   SUBROUTINE change_volume(conf, T_ind, move_types, rng_stream, tmc_params, &
                            mv_cen_of_mass)
      TYPE(tree_type), POINTER                           :: conf
      INTEGER                                            :: t_ind
      TYPE(tmc_move_type), POINTER                       :: move_types
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      TYPE(tmc_param_type), POINTER                      :: tmc_params
      LOGICAL                                            :: mv_cen_of_mass
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'change_volume'
 
      INTEGER                                            :: atom, dir, handle, ind, ind_e, mol
      REAL(kind=dp)                                      :: rnd, vol
      REAL(kind=dp), DIMENSION(3)                        :: box_length_new, box_length_orig, &
                                                            box_scale_old
      REAL(kind=dp), DIMENSION(:), POINTER               :: disp, scaling
 
      NULLIFY (scaling, disp)
 
      cpassert(ASSOCIATED(conf))
      cpassert(ASSOCIATED(move_types))
      cpassert(ASSOCIATED(tmc_params))
      cpassert(t_ind .GT. 0 .AND. t_ind .LE. tmc_params%nr_temp)
      cpassert(tmc_params%dim_per_elem .EQ. 3)
      cpassert(tmc_params%cell%orthorhombic)
 
      ! start the timing
      CALL timeset(routinen, handle)
 
      ALLOCATE (scaling(tmc_params%dim_per_elem))
      ALLOCATE (disp(tmc_params%dim_per_elem))
 
      box_scale_old(:) = conf%box_scale
      ! get the cell vector length of the configuration (before move)
      CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, &
                           abc=box_length_new)
 
      IF (.false.) THEN
         ! the volume move in volume space (dV)
         IF (tmc_params%v_isotropic) THEN
            CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, &
                                 abc=box_length_new, vol=vol)
            rnd = rng_stream%next()
            vol = vol + (rnd - 0.5_dp)*2.0_dp*move_types%mv_size(mv_type_volume_move, t_ind)
            box_length_new(:) = vol**(1/real(3, kind=dp))
         ELSE
            CALL get_scaled_cell(cell=tmc_params%cell, box_scale=conf%box_scale, &
                                 abc=box_length_new, vol=vol)
            rnd = rng_stream%next()
            vol = vol + (rnd - 0.5_dp)*2.0_dp*move_types%mv_size(mv_type_volume_move, t_ind)
            rnd = rng_stream%next()
            dir = 1 + int(rnd*3)
            box_length_new(dir) = 1.0_dp
            box_length_new(dir) = vol/product(box_length_new(:))
         END IF
      ELSE
         ! the volume move in box length space (dL)
         ! increase / decrease box length in this direction
         ! l_n = l_o +- rnd * mv_size
         IF (tmc_params%v_isotropic) THEN
            rnd = rng_stream%next()
            box_length_new(:) = box_length_new(:) + &
                                (rnd - 0.5_dp)*2.0_dp* &
                                move_types%mv_size(mv_type_volume_move, t_ind)
         ELSE
            ! select a random direction
            rnd = rng_stream%next()
            dir = 1 + int(rnd*3)
            rnd = rng_stream%next()
            box_length_new(dir) = box_length_new(dir) + &
                                  (rnd - 0.5_dp)*2.0_dp* &
                                  move_types%mv_size(mv_type_volume_move, t_ind)
         END IF
      END IF
 
      ! get the original box length
      scaling(:) = 1.0_dp
      CALL get_scaled_cell(cell=tmc_params%cell, &
                           box_scale=scaling, &
                           abc=box_length_orig)
      ! get the new box scale
      conf%box_scale(:) = box_length_new(:)/box_length_orig(:)
      ! molecule scaling
      scaling(:) = conf%box_scale(:)/box_scale_old(:)
 
      IF (mv_cen_of_mass .EQV. .false.) THEN
         ! homogene scaling of atomic coordinates
         DO atom = 1, SIZE(conf%pos), tmc_params%dim_per_elem
            conf%pos(atom:atom + tmc_params%dim_per_elem - 1) = &
               conf%pos(atom:atom + tmc_params%dim_per_elem - 1)*scaling(:)
         END DO
      ELSE
         DO mol = 1, maxval(conf%mol(:))
            ! move the molecule related to the molecule center of mass
            ! get center of mass
            cpassert(ASSOCIATED(tmc_params%atoms))
 
            CALL get_mol_indeces(tmc_params=tmc_params, mol_arr=conf%mol, mol=mol, &
                                 start_ind=ind, end_ind=ind_e)
            CALL center_of_mass( &
               pos=conf%pos(ind:ind_e + tmc_params%dim_per_elem - 1), &
               atoms=tmc_params%atoms(int(ind/real(tmc_params%dim_per_elem, kind=dp)) + 1: &
                                      int(ind_e/real(tmc_params%dim_per_elem, kind=dp)) + 1), &
               center=disp)
            ! calculate the center of mass DISPLACEMENT
            disp(:) = disp(:)*(scaling(:) - 1.0_dp)
            ! displace all atoms of the molecule
            DO atom = ind, ind_e + tmc_params%dim_per_elem - 1, tmc_params%dim_per_elem
               conf%pos(atom:atom + tmc_params%dim_per_elem - 1) = &
                  conf%pos(atom:atom + tmc_params%dim_per_elem - 1) + disp(:)
            END DO
         END DO
      END IF
 
      DEALLOCATE (scaling)
      DEALLOCATE (disp)
 
      ! end the timing
      CALL timestop(handle)
   END SUBROUTINE change_volume
 
! **************************************************************************************************
!> \brief volume move, two atoms of different types are swapped, both selected
!>        randomly
!> \param conf configuration to change with positions
!> \param move_types ...
!> \param rng_stream random number generator stream
!> \param tmc_params TMC parameters with e.g. dimensions of atoms and molecules
!> \author Mandes 11.2012
! **************************************************************************************************
   SUBROUTINE swap_atoms(conf, move_types, rng_stream, tmc_params)
      TYPE(tree_type), POINTER                           :: conf
      TYPE(tmc_move_type), POINTER                       :: move_types
      TYPE(rng_stream_type), INTENT(INOUT)               :: rng_stream
      TYPE(tmc_param_type), POINTER                      :: tmc_params
 
      INTEGER                                            :: a_1, a_2, ind_1, ind_2
      LOGICAL                                            :: found
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: pos_tmp
 
      cpassert(ASSOCIATED(conf))
      cpassert(ASSOCIATED(move_types))
      cpassert(ASSOCIATED(tmc_params))
      cpassert(ASSOCIATED(tmc_params%atoms))
 
      ! loop until two different atoms are found
      atom_search_loop: DO
         ! select one atom randomly
         a_1 = int(SIZE(conf%pos)/real(tmc_params%dim_per_elem, kind=dp)* &
                   rng_stream%next()) + 1
         ! select the second atom randomly
         a_2 = int(SIZE(conf%pos)/real(tmc_params%dim_per_elem, kind=dp)* &
                   rng_stream%next()) + 1
         ! check if they have different kinds
         IF (tmc_params%atoms(a_1)%name .NE. tmc_params%atoms(a_2)%name) THEN
            ! if present, check if atoms have different type related to the specified table
            IF (ASSOCIATED(move_types%atom_lists)) THEN
               DO ind_1 = 1, SIZE(move_types%atom_lists)
                  IF (any(move_types%atom_lists(ind_1)%atoms(:) .EQ. &
                          tmc_params%atoms(a_1)%name) .AND. &
                      any(move_types%atom_lists(ind_1)%atoms(:) .EQ. &
                          tmc_params%atoms(a_2)%name)) THEN
                     found = .true.
                     EXIT atom_search_loop
                  END IF
               END DO
            ELSE
               found = .true.
               EXIT atom_search_loop
            END IF
         END IF
      END DO atom_search_loop
      IF (found) THEN
         ! perform coordinate exchange
         ALLOCATE (pos_tmp(tmc_params%dim_per_elem))
         ind_1 = (a_1 - 1)*tmc_params%dim_per_elem + 1
         pos_tmp(:) = conf%pos(ind_1:ind_1 + tmc_params%dim_per_elem - 1)
         ind_2 = (a_2 - 1)*tmc_params%dim_per_elem + 1
         conf%pos(ind_1:ind_1 + tmc_params%dim_per_elem - 1) = &
            conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1)
         conf%pos(ind_2:ind_2 + tmc_params%dim_per_elem - 1) = pos_tmp(:)
         DEALLOCATE (pos_tmp)
      END IF
   END SUBROUTINE swap_atoms
 
END MODULE tmc_moves