LCOV - code coverage report
Current view: top level - src/motion - integrator_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 552 636 86.8 %
Date: 2024-04-20 06:29:22 Functions: 20 22 90.9 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Provides integrator utility routines for the integrators
      10             : !> \par History
      11             : !>      Teodoro Laino [tlaino] 02.2008 - Splitting from integrator and creation
      12             : !>      Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
      13             : !>                                       (patch by Marcel Baer)
      14             : ! **************************************************************************************************
      15             : MODULE integrator_utils
      16             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      17             :                                               get_atomic_kind
      18             :    USE barostat_types,                  ONLY: do_clv_x,&
      19             :                                               do_clv_xy,&
      20             :                                               do_clv_xyz,&
      21             :                                               do_clv_xz,&
      22             :                                               do_clv_y,&
      23             :                                               do_clv_yz,&
      24             :                                               do_clv_z
      25             :    USE cell_types,                      ONLY: cell_type
      26             :    USE constraint,                      ONLY: rattle_roll_control
      27             :    USE constraint_util,                 ONLY: pv_constraint
      28             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      29             :    USE extended_system_types,           ONLY: debug_isotropic_limit,&
      30             :                                               debug_uniaxial_limit,&
      31             :                                               npt_info_type
      32             :    USE input_constants,                 ONLY: &
      33             :         isokin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, &
      34             :         nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, npt_ia_ensemble, nve_ensemble, &
      35             :         nvt_ensemble
      36             :    USE kinds,                           ONLY: dp
      37             :    USE md_environment_types,            ONLY: get_md_env,&
      38             :                                               md_environment_type
      39             :    USE message_passing,                 ONLY: mp_comm_type,&
      40             :                                               mp_para_env_type
      41             :    USE molecule_kind_types,             ONLY: local_fixd_constraint_type,&
      42             :                                               molecule_kind_type
      43             :    USE molecule_types,                  ONLY: get_molecule,&
      44             :                                               global_constraint_type,&
      45             :                                               molecule_type
      46             :    USE particle_types,                  ONLY: particle_type,&
      47             :                                               update_particle_set
      48             :    USE shell_potential_types,           ONLY: shell_kind_type
      49             :    USE simpar_types,                    ONLY: simpar_type
      50             :    USE thermostat_types,                ONLY: set_thermostats,&
      51             :                                               thermostats_type
      52             :    USE virial_types,                    ONLY: virial_type
      53             : #include "../base/base_uses.f90"
      54             : 
      55             :    IMPLICIT NONE
      56             : 
      57             :    PRIVATE
      58             : 
      59             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'integrator_utils'
      60             : 
      61             : ! **************************************************************************************************
      62             :    TYPE old_variables_type
      63             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: v
      64             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: r
      65             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: eps
      66             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: veps
      67             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: h
      68             :    END TYPE old_variables_type
      69             : 
      70             : ! **************************************************************************************************
      71             :    TYPE tmp_variables_type
      72             :       INTEGER, POINTER :: itimes
      73             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: pos
      74             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: vel
      75             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: shell_pos
      76             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: shell_vel
      77             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: core_pos
      78             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: core_vel
      79             :       REAL(KIND=dp) :: max_vel, max_vel_sc
      80             :       REAL(KIND=dp) :: max_dvel, max_dvel_sc
      81             :       REAL(KIND=dp) :: max_dr, max_dsc
      82             :       REAL(KIND=dp) :: arg_r(3), arg_v(3), u(3, 3), e_val(3), s, ds
      83             :       REAL(KIND=dp), DIMENSION(3) :: poly_r, &
      84             :                                      poly_v, scale_r, scale_v
      85             :    END TYPE tmp_variables_type
      86             : 
      87             :    INTERFACE set
      88             :       MODULE PROCEDURE set_particle_set, set_vel
      89             :    END INTERFACE
      90             : 
      91             :    INTERFACE update_pv
      92             :       MODULE PROCEDURE update_pv_particle_set, update_pv_velocity
      93             :    END INTERFACE
      94             : 
      95             :    INTERFACE damp_v
      96             :       MODULE PROCEDURE damp_v_particle_set, damp_v_velocity
      97             :    END INTERFACE
      98             : 
      99             :    PUBLIC :: old_variables_type, tmp_variables_type, allocate_old, deallocate_old, &
     100             :              allocate_tmp, update_dealloc_tmp, damp_v, set, update_pv, get_s_ds, &
     101             :              rattle_roll_setup, damp_veps, update_veps, vv_first, &
     102             :              vv_second, variable_timestep
     103             : 
     104             : CONTAINS
     105             : 
     106             : ! **************************************************************************************************
     107             : !> \brief ...
     108             : !> \param old ...
     109             : !> \param particle_set ...
     110             : !> \param npt ...
     111             : !> \par History
     112             : !>      none
     113             : !> \author CJM
     114             : ! **************************************************************************************************
     115        2520 :    SUBROUTINE allocate_old(old, particle_set, npt)
     116             :       TYPE(old_variables_type), POINTER                  :: old
     117             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     118             :       TYPE(npt_info_type), POINTER                       :: npt(:, :)
     119             : 
     120             :       INTEGER                                            :: idim, jdim, natoms
     121             : 
     122        2520 :       natoms = SIZE(particle_set)
     123        2520 :       idim = SIZE(npt, 1)
     124        2520 :       jdim = SIZE(npt, 2)
     125        2520 :       CPASSERT(.NOT. ASSOCIATED(old))
     126        2520 :       ALLOCATE (old)
     127             : 
     128        7560 :       ALLOCATE (old%v(natoms, 3))
     129     2163108 :       old%v = 0.0_dp
     130        7560 :       ALLOCATE (old%r(natoms, 3))
     131     2163108 :       old%r = 0.0_dp
     132       10080 :       ALLOCATE (old%eps(idim, jdim))
     133       16520 :       old%eps = 0.0_dp
     134       10080 :       ALLOCATE (old%veps(idim, jdim))
     135       16520 :       old%veps = 0.0_dp
     136        2520 :       ALLOCATE (old%h(3, 3))
     137       32760 :       old%h = 0.0_dp
     138             : 
     139        2520 :    END SUBROUTINE allocate_old
     140             : 
     141             : ! **************************************************************************************************
     142             : !> \brief ...
     143             : !> \param old ...
     144             : !> \par History
     145             : !>      none
     146             : !> \author CJM
     147             : ! **************************************************************************************************
     148        2520 :    SUBROUTINE deallocate_old(old)
     149             :       TYPE(old_variables_type), POINTER                  :: old
     150             : 
     151        2520 :       IF (ASSOCIATED(old)) THEN
     152        2520 :          IF (ASSOCIATED(old%v)) THEN
     153        2520 :             DEALLOCATE (old%v)
     154             :          END IF
     155        2520 :          IF (ASSOCIATED(old%r)) THEN
     156        2520 :             DEALLOCATE (old%r)
     157             :          END IF
     158        2520 :          IF (ASSOCIATED(old%eps)) THEN
     159        2520 :             DEALLOCATE (old%eps)
     160             :          END IF
     161        2520 :          IF (ASSOCIATED(old%veps)) THEN
     162        2520 :             DEALLOCATE (old%veps)
     163             :          END IF
     164        2520 :          IF (ASSOCIATED(old%h)) THEN
     165        2520 :             DEALLOCATE (old%h)
     166             :          END IF
     167        2520 :          DEALLOCATE (old)
     168             :       END IF
     169             : 
     170        2520 :    END SUBROUTINE deallocate_old
     171             : 
     172             : ! **************************************************************************************************
     173             : !> \brief allocate temporary variables to store positions and velocities
     174             : !>        used by the velocity-verlet integrator
     175             : !> \param md_env ...
     176             : !> \param tmp ...
     177             : !> \param nparticle ...
     178             : !> \param nshell ...
     179             : !> \param shell_adiabatic ...
     180             : !> \par History
     181             : !>      none
     182             : !> \author MI (February 2008)
     183             : ! **************************************************************************************************
     184       42611 :    SUBROUTINE allocate_tmp(md_env, tmp, nparticle, nshell, shell_adiabatic)
     185             :       TYPE(md_environment_type), POINTER                 :: md_env
     186             :       TYPE(tmp_variables_type), POINTER                  :: tmp
     187             :       INTEGER, INTENT(IN)                                :: nparticle, nshell
     188             :       LOGICAL, INTENT(IN)                                :: shell_adiabatic
     189             : 
     190       42611 :       CPASSERT(.NOT. ASSOCIATED(tmp))
     191       42611 :       ALLOCATE (tmp)
     192             : 
     193             :       ! Nullify Components
     194       42611 :       NULLIFY (tmp%pos)
     195       42611 :       NULLIFY (tmp%vel)
     196       42611 :       NULLIFY (tmp%shell_pos)
     197       42611 :       NULLIFY (tmp%shell_vel)
     198       42611 :       NULLIFY (tmp%core_pos)
     199       42611 :       NULLIFY (tmp%core_vel)
     200       42611 :       NULLIFY (tmp%itimes)
     201             : 
     202             :       !     *** Allocate work storage for positions and velocities ***
     203      127833 :       ALLOCATE (tmp%pos(3, nparticle))
     204             : 
     205      127833 :       ALLOCATE (tmp%vel(3, nparticle))
     206             : 
     207    23778595 :       tmp%pos(:, :) = 0.0_dp
     208    23778595 :       tmp%vel(:, :) = 0.0_dp
     209             : 
     210       42611 :       IF (shell_adiabatic) THEN
     211             :          !     *** Allocate work storage for positions and velocities ***
     212        6870 :          ALLOCATE (tmp%shell_pos(3, nshell))
     213             : 
     214        6870 :          ALLOCATE (tmp%core_pos(3, nshell))
     215             : 
     216        6870 :          ALLOCATE (tmp%shell_vel(3, nshell))
     217             : 
     218        6870 :          ALLOCATE (tmp%core_vel(3, nshell))
     219             : 
     220      906050 :          tmp%shell_pos(:, :) = 0.0_dp
     221      906050 :          tmp%shell_vel(:, :) = 0.0_dp
     222      906050 :          tmp%core_pos(:, :) = 0.0_dp
     223      906050 :          tmp%core_vel(:, :) = 0.0_dp
     224             :       END IF
     225             : 
     226      170444 :       tmp%arg_r = 0.0_dp
     227      170444 :       tmp%arg_v = 0.0_dp
     228      553943 :       tmp%u = 0.0_dp
     229      170444 :       tmp%e_val = 0.0_dp
     230       42611 :       tmp%max_vel = 0.0_dp
     231       42611 :       tmp%max_vel_sc = 0.0_dp
     232       42611 :       tmp%max_dvel = 0.0_dp
     233       42611 :       tmp%max_dvel_sc = 0.0_dp
     234      170444 :       tmp%scale_r = 1.0_dp
     235      170444 :       tmp%scale_v = 1.0_dp
     236      170444 :       tmp%poly_r = 1.0_dp
     237      170444 :       tmp%poly_v = 1.0_dp
     238             : 
     239       42611 :       CALL get_md_env(md_env=md_env, itimes=tmp%itimes)
     240             : 
     241       42611 :    END SUBROUTINE allocate_tmp
     242             : 
     243             : ! **************************************************************************************************
     244             : !> \brief update positions and deallocate temporary variable
     245             : !> \param tmp ...
     246             : !> \param particle_set ...
     247             : !> \param shell_particle_set ...
     248             : !> \param core_particle_set ...
     249             : !> \param para_env ...
     250             : !> \param shell_adiabatic ...
     251             : !> \param pos ...
     252             : !> \param vel ...
     253             : !> \param should_deall_vel ...
     254             : !> \par History
     255             : !>      none
     256             : !> \author MI (February 2008)
     257             : ! **************************************************************************************************
     258       86372 :    SUBROUTINE update_dealloc_tmp(tmp, particle_set, shell_particle_set, &
     259             :                                  core_particle_set, para_env, shell_adiabatic, pos, vel, &
     260             :                                  should_deall_vel)
     261             : 
     262             :       TYPE(tmp_variables_type), POINTER                  :: tmp
     263             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set, shell_particle_set, &
     264             :                                                             core_particle_set
     265             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     266             :       LOGICAL, INTENT(IN)                                :: shell_adiabatic
     267             :       LOGICAL, INTENT(IN), OPTIONAL                      :: pos, vel, should_deall_vel
     268             : 
     269             :       LOGICAL                                            :: my_deall, my_pos, my_vel
     270             : 
     271       86372 :       CPASSERT(ASSOCIATED(tmp))
     272       86372 :       my_pos = .FALSE.
     273       86372 :       my_vel = .FALSE.
     274       86372 :       my_deall = .TRUE.
     275       86372 :       IF (PRESENT(pos)) my_pos = pos
     276       86372 :       IF (PRESENT(vel)) my_vel = vel
     277       86372 :       IF (PRESENT(should_deall_vel)) my_deall = should_deall_vel
     278             : 
     279             :       !      *** Broadcast the new particle positions ***
     280       86372 :       IF (my_pos) THEN
     281       42611 :          CALL update_particle_set(particle_set, para_env, pos=tmp%pos)
     282       42611 :          DEALLOCATE (tmp%pos)
     283             : 
     284       42611 :          IF (shell_adiabatic) THEN
     285        2290 :             CALL update_particle_set(shell_particle_set, para_env, pos=tmp%shell_pos)
     286        2290 :             CALL update_particle_set(core_particle_set, para_env, pos=tmp%core_pos)
     287        2290 :             DEALLOCATE (tmp%shell_pos)
     288        2290 :             DEALLOCATE (tmp%core_pos)
     289             :          END IF
     290             :       END IF
     291             : 
     292             :       !     *** Broadcast the new particle velocities ***
     293       86372 :       IF (my_vel) THEN
     294       43761 :          CALL update_particle_set(particle_set, para_env, vel=tmp%vel)
     295       43761 :          IF (shell_adiabatic) THEN
     296        2290 :             CALL update_particle_set(shell_particle_set, para_env, vel=tmp%shell_vel)
     297        2290 :             CALL update_particle_set(core_particle_set, para_env, vel=tmp%core_vel)
     298             :          END IF
     299       43761 :          IF (my_deall) THEN
     300       42611 :             DEALLOCATE (tmp%vel)
     301       42611 :             IF (shell_adiabatic) THEN
     302        2290 :                DEALLOCATE (tmp%shell_vel)
     303        2290 :                DEALLOCATE (tmp%core_vel)
     304             :             END IF
     305       42611 :             CPASSERT(.NOT. ASSOCIATED(tmp%pos))
     306       42611 :             CPASSERT(.NOT. ASSOCIATED(tmp%shell_pos))
     307       42611 :             CPASSERT(.NOT. ASSOCIATED(tmp%core_pos))
     308       42611 :             DEALLOCATE (tmp)
     309             :          END IF
     310             :       END IF
     311             : 
     312       86372 :    END SUBROUTINE update_dealloc_tmp
     313             : 
     314             : ! **************************************************************************************************
     315             : !> \brief ...
     316             : !> \param tmp ...
     317             : !> \param nparticle_kind ...
     318             : !> \param atomic_kind_set ...
     319             : !> \param local_particles ...
     320             : !> \param particle_set ...
     321             : !> \param dt ...
     322             : !> \param para_env ...
     323             : !> \param tmpv ...
     324             : !> \par History
     325             : !>      - Created [2004-07]
     326             : !> \author Joost VandeVondele
     327             : ! **************************************************************************************************
     328          12 :    SUBROUTINE get_s_ds(tmp, nparticle_kind, atomic_kind_set, local_particles, particle_set, &
     329             :                        dt, para_env, tmpv)
     330             :       TYPE(tmp_variables_type), POINTER                  :: tmp
     331             :       INTEGER                                            :: nparticle_kind
     332             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     333             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     334             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     335             :       REAL(KIND=dp)                                      :: dt
     336             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     337             :       LOGICAL, INTENT(IN), OPTIONAL                      :: tmpv
     338             : 
     339             :       INTEGER                                            :: iparticle, iparticle_kind, &
     340             :                                                             iparticle_local, nparticle_local
     341             :       LOGICAL                                            :: my_tmpv
     342             :       REAL(KIND=dp)                                      :: a, b, K, mass, rb
     343             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     344             : 
     345          12 :       my_tmpv = .FALSE.
     346          12 :       IF (PRESENT(tmpv)) my_tmpv = tmpv
     347             : 
     348          12 :       K = 0.0_dp
     349          12 :       a = 0.0_dp
     350          12 :       b = 0.0_dp
     351          36 :       DO iparticle_kind = 1, nparticle_kind
     352          24 :          atomic_kind => atomic_kind_set(iparticle_kind)
     353          24 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     354          36 :          IF (mass /= 0.0_dp) THEN
     355          24 :             nparticle_local = local_particles%n_el(iparticle_kind)
     356          24 :             IF (my_tmpv) THEN
     357          21 :                DO iparticle_local = 1, nparticle_local
     358           9 :                   iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     359          36 :                   K = K + 0.5_dp*mass*DOT_PRODUCT(tmp%vel(:, iparticle), tmp%vel(:, iparticle))
     360          36 :                   a = a + DOT_PRODUCT(tmp%vel(:, iparticle), particle_set(iparticle)%f(:))
     361          48 :                   b = b + (1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:), particle_set(iparticle)%f(:))
     362             :                END DO
     363             :             ELSE
     364          21 :                DO iparticle_local = 1, nparticle_local
     365           9 :                   iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     366          36 :                   K = K + 0.5_dp*mass*DOT_PRODUCT(particle_set(iparticle)%v(:), particle_set(iparticle)%v(:))
     367          36 :                   a = a + DOT_PRODUCT(particle_set(iparticle)%v(:), particle_set(iparticle)%f(:))
     368          48 :                   b = b + (1.0_dp/mass)*DOT_PRODUCT(particle_set(iparticle)%f(:), particle_set(iparticle)%f(:))
     369             :                END DO
     370             :             END IF
     371             :          END IF
     372             :       END DO
     373          12 :       CALL para_env%sum(K)
     374          12 :       CALL para_env%sum(a)
     375          12 :       CALL para_env%sum(b)
     376          12 :       a = a/(2.0_dp*K)
     377          12 :       b = b/(2.0_dp*K)
     378          12 :       rb = SQRT(b)
     379          12 :       tmp%s = (a/b)*(COSH(dt*rb/2.0_dp) - 1) + SINH(dt*rb/2.0_dp)/rb
     380          12 :       tmp%ds = (a/b)*(SINH(dt*rb/2.0_dp)*rb) + COSH(dt*rb/2.0_dp)
     381             : 
     382          12 :    END SUBROUTINE get_s_ds
     383             : 
     384             : ! **************************************************************************************************
     385             : !> \brief ...
     386             : !> \param old ...
     387             : !> \param atomic_kind_set ...
     388             : !> \param particle_set ...
     389             : !> \param local_particles ...
     390             : !> \param cell ...
     391             : !> \param npt ...
     392             : !> \param char ...
     393             : !> \par History
     394             : !>      none
     395             : !> \author CJM
     396             : ! **************************************************************************************************
     397        2504 :    SUBROUTINE set_particle_set(old, atomic_kind_set, particle_set, local_particles, cell, npt, char)
     398             :       TYPE(old_variables_type), POINTER                  :: old
     399             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
     400             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     401             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     402             :       TYPE(cell_type), POINTER                           :: cell
     403             :       TYPE(npt_info_type), DIMENSION(:, :), POINTER      :: npt
     404             :       CHARACTER(LEN=*), INTENT(IN)                       :: char
     405             : 
     406             :       INTEGER                                            :: idim, iparticle, iparticle_kind, &
     407             :                                                             iparticle_local, nparticle_kind, &
     408             :                                                             nparticle_local
     409             : 
     410        2504 :       nparticle_kind = SIZE(atomic_kind_set)
     411             :       SELECT CASE (char)
     412             :       CASE ('F') ! forward assigning the old
     413        2040 :          DO iparticle_kind = 1, nparticle_kind
     414        1362 :             nparticle_local = local_particles%n_el(iparticle_kind)
     415      124267 :             DO iparticle_local = 1, nparticle_local
     416      122227 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     417      490270 :                DO idim = 1, 3
     418      366681 :                   old%v(iparticle, idim) = particle_set(iparticle)%v(idim)
     419      488908 :                   old%r(iparticle, idim) = particle_set(iparticle)%r(idim)
     420             :                END DO
     421             :             END DO
     422             :          END DO
     423        2134 :          old%eps(:, :) = npt(:, :)%eps
     424        2134 :          old%veps(:, :) = npt(:, :)%v
     425       16950 :          old%h(:, :) = cell%hmat(:, :)
     426             :       CASE ('B') ! back assigning the original variables
     427        5498 :          DO iparticle_kind = 1, nparticle_kind
     428        3672 :             nparticle_local = local_particles%n_el(iparticle_kind)
     429      361572 :             DO iparticle_local = 1, nparticle_local
     430      356074 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     431     1427968 :                DO idim = 1, 3
     432     1068222 :                   particle_set(iparticle)%v(idim) = old%v(iparticle, idim)
     433     1424296 :                   particle_set(iparticle)%r(idim) = old%r(iparticle, idim)
     434             :                END DO
     435             :             END DO
     436             :          END DO
     437        5678 :          npt(:, :)%eps = old%eps(:, :)
     438        5678 :          npt(:, :)%v = old%veps(:, :)
     439       48154 :          cell%hmat(:, :) = old%h(:, :)
     440             :       END SELECT
     441             : 
     442        2504 :    END SUBROUTINE set_particle_set
     443             : 
     444             : ! **************************************************************************************************
     445             : !> \brief ...
     446             : !> \param old ...
     447             : !> \param atomic_kind_set ...
     448             : !> \param particle_set ...
     449             : !> \param vel ...
     450             : !> \param local_particles ...
     451             : !> \param cell ...
     452             : !> \param npt ...
     453             : !> \param char ...
     454             : !> \par History
     455             : !>      none
     456             : !> \author CJM
     457             : ! **************************************************************************************************
     458        2472 :    SUBROUTINE set_vel(old, atomic_kind_set, particle_set, vel, local_particles, cell, npt, char)
     459             :       TYPE(old_variables_type), POINTER                  :: old
     460             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
     461             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     462             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     463             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     464             :       TYPE(cell_type), POINTER                           :: cell
     465             :       TYPE(npt_info_type), DIMENSION(:, :), POINTER      :: npt
     466             :       CHARACTER(LEN=*), INTENT(IN)                       :: char
     467             : 
     468             :       INTEGER                                            :: idim, iparticle, iparticle_kind, &
     469             :                                                             iparticle_local, nparticle_kind, &
     470             :                                                             nparticle_local
     471             : 
     472        2472 :       nparticle_kind = SIZE(atomic_kind_set)
     473             :       SELECT CASE (char)
     474             :       CASE ('F') ! forward assigning the old
     475        2040 :          DO iparticle_kind = 1, nparticle_kind
     476        1362 :             nparticle_local = local_particles%n_el(iparticle_kind)
     477      124267 :             DO iparticle_local = 1, nparticle_local
     478      122227 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     479      490270 :                DO idim = 1, 3
     480      366681 :                   old%v(iparticle, idim) = vel(idim, iparticle)
     481      488908 :                   old%r(iparticle, idim) = particle_set(iparticle)%r(idim)
     482             :                END DO
     483             :             END DO
     484             :          END DO
     485        2134 :          old%eps(:, :) = npt(:, :)%eps
     486        2134 :          old%veps(:, :) = npt(:, :)%v
     487       16950 :          old%h(:, :) = cell%hmat(:, :)
     488             :       CASE ('B') ! back assigning the original variables
     489        5400 :          DO iparticle_kind = 1, nparticle_kind
     490        3606 :             nparticle_local = local_particles%n_el(iparticle_kind)
     491      317027 :             DO iparticle_local = 1, nparticle_local
     492      311627 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     493     1250114 :                DO idim = 1, 3
     494      934881 :                   vel(idim, iparticle) = old%v(iparticle, idim)
     495     1246508 :                   particle_set(iparticle)%r(idim) = old%r(iparticle, idim)
     496             :                END DO
     497             :             END DO
     498             :          END DO
     499        5582 :          npt(:, :)%eps = old%eps(:, :)
     500        5582 :          npt(:, :)%v = old%veps(:, :)
     501       47322 :          cell%hmat(:, :) = old%h(:, :)
     502             :       END SELECT
     503             : 
     504        2472 :    END SUBROUTINE set_vel
     505             : 
     506             : ! **************************************************************************************************
     507             : !> \brief overloaded routine provides damping for particles via nph_uniaxial_damped dynamics
     508             : !> \param molecule_kind_set ...
     509             : !> \param molecule_set ...
     510             : !> \param particle_set ...
     511             : !> \param local_molecules ...
     512             : !> \param gamma1 ...
     513             : !> \param npt ...
     514             : !> \param dt ...
     515             : !> \param group ...
     516             : !> \par History
     517             : !>      none
     518             : !> \author CJM
     519             : ! **************************************************************************************************
     520          20 :    SUBROUTINE damp_v_particle_set(molecule_kind_set, molecule_set, &
     521             :                                   particle_set, local_molecules, gamma1, npt, dt, group)
     522             : 
     523             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     524             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     525             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     526             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     527             :       REAL(KIND=dp), INTENT(IN)                          :: gamma1
     528             :       TYPE(npt_info_type), INTENT(IN)                    :: npt
     529             :       REAL(KIND=dp), INTENT(IN)                          :: dt
     530             : 
     531             :       CLASS(mp_comm_type), INTENT(IN)                     :: group
     532             : 
     533             :       INTEGER                                            :: first_atom, ikind, imol, imol_local, &
     534             :                                                             ipart, last_atom, nmol_local
     535             :       REAL(KIND=dp)                                      :: alpha, ikin, kin, mass, scale
     536             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     537             :       TYPE(molecule_type), POINTER                       :: molecule
     538             : 
     539          20 :       kin = 0.0_dp
     540        1420 :       DO ikind = 1, SIZE(molecule_kind_set)
     541             :          ! Compute the total kinetic energy
     542        1400 :          nmol_local = local_molecules%n_el(ikind)
     543        2120 :          DO imol_local = 1, nmol_local
     544         700 :             imol = local_molecules%list(ikind)%array(imol_local)
     545         700 :             molecule => molecule_set(imol)
     546             :             CALL get_molecule(molecule, first_atom=first_atom, &
     547         700 :                               last_atom=last_atom)
     548        2800 :             DO ipart = first_atom, last_atom
     549         700 :                atomic_kind => particle_set(ipart)%atomic_kind
     550         700 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     551             :                kin = kin + mass*particle_set(ipart)%v(1)* &
     552         700 :                      particle_set(ipart)%v(1)
     553             :                kin = kin + mass*particle_set(ipart)%v(2)* &
     554         700 :                      particle_set(ipart)%v(2)
     555             :                kin = kin + mass*particle_set(ipart)%v(3)* &
     556        1400 :                      particle_set(ipart)%v(3)
     557             :             END DO
     558             :          END DO
     559             :       END DO
     560             :       !
     561          20 :       CALL group%sum(kin)
     562             :       !
     563          20 :       ikin = 1.0_dp/kin
     564          20 :       scale = 1.0_dp
     565          20 :       alpha = 2.0_dp*npt%mass*npt%v*npt%v*gamma1*ikin
     566          20 :       scale = scale*SQRT(1.0_dp + alpha*0.5_dp*dt)
     567             :       ! Scale
     568        1420 :       DO ikind = 1, SIZE(molecule_kind_set)
     569        1400 :          nmol_local = local_molecules%n_el(ikind)
     570        2120 :          DO imol_local = 1, nmol_local
     571         700 :             imol = local_molecules%list(ikind)%array(imol_local)
     572         700 :             molecule => molecule_set(imol)
     573             :             CALL get_molecule(molecule, first_atom=first_atom, &
     574         700 :                               last_atom=last_atom)
     575        2800 :             DO ipart = first_atom, last_atom
     576         700 :                particle_set(ipart)%v(1) = particle_set(ipart)%v(1)*scale
     577         700 :                particle_set(ipart)%v(2) = particle_set(ipart)%v(2)*scale
     578        1400 :                particle_set(ipart)%v(3) = particle_set(ipart)%v(3)*scale
     579             :             END DO
     580             :          END DO
     581             :       END DO
     582             : 
     583          20 :    END SUBROUTINE damp_v_particle_set
     584             : 
     585             : ! **************************************************************************************************
     586             : !> \brief overloaded provides damping for particles via nph_uniaxial_damped dynamics
     587             : !> \param molecule_kind_set ...
     588             : !> \param molecule_set ...
     589             : !> \param particle_set ...
     590             : !> \param local_molecules ...
     591             : !> \param vel ...
     592             : !> \param gamma1 ...
     593             : !> \param npt ...
     594             : !> \param dt ...
     595             : !> \param group ...
     596             : !> \par History
     597             : !>      none
     598             : !> \author CJM
     599             : ! **************************************************************************************************
     600          20 :    SUBROUTINE damp_v_velocity(molecule_kind_set, molecule_set, &
     601          20 :                               particle_set, local_molecules, vel, gamma1, npt, dt, group)
     602             : 
     603             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     604             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     605             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     606             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     607             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     608             :       REAL(KIND=dp), INTENT(IN)                          :: gamma1
     609             :       TYPE(npt_info_type), INTENT(IN)                    :: npt
     610             :       REAL(KIND=dp), INTENT(IN)                          :: dt
     611             : 
     612             :       CLASS(mp_comm_type), INTENT(IN)                     :: group
     613             : 
     614             :       INTEGER                                            :: first_atom, ikind, imol, imol_local, &
     615             :                                                             ipart, last_atom, nmol_local
     616             :       REAL(KIND=dp)                                      :: alpha, ikin, kin, mass, scale
     617             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     618             :       TYPE(molecule_type), POINTER                       :: molecule
     619             : 
     620          20 :       kin = 0.0_dp
     621             :       ! Compute the total kinetic energy
     622        1420 :       DO ikind = 1, SIZE(molecule_kind_set)
     623        1400 :          nmol_local = local_molecules%n_el(ikind)
     624        2120 :          DO imol_local = 1, nmol_local
     625         700 :             imol = local_molecules%list(ikind)%array(imol_local)
     626         700 :             molecule => molecule_set(imol)
     627             :             CALL get_molecule(molecule, first_atom=first_atom, &
     628         700 :                               last_atom=last_atom)
     629        2800 :             DO ipart = first_atom, last_atom
     630         700 :                atomic_kind => particle_set(ipart)%atomic_kind
     631         700 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     632         700 :                kin = kin + mass*vel(1, ipart)*vel(1, ipart)
     633         700 :                kin = kin + mass*vel(2, ipart)*vel(2, ipart)
     634        1400 :                kin = kin + mass*vel(3, ipart)*vel(3, ipart)
     635             :             END DO
     636             :          END DO
     637             :       END DO
     638             :       !
     639          20 :       CALL group%sum(kin)
     640             :       !
     641          20 :       ikin = 1.0_dp/kin
     642          20 :       scale = 1.0_dp
     643          20 :       alpha = 2.0_dp*npt%mass*npt%v*npt%v*gamma1*ikin
     644          20 :       scale = scale*SQRT(1.0_dp + alpha*0.5_dp*dt)
     645             :       ! Scale
     646        1420 :       DO ikind = 1, SIZE(molecule_kind_set)
     647        1400 :          nmol_local = local_molecules%n_el(ikind)
     648        2120 :          DO imol_local = 1, nmol_local
     649         700 :             imol = local_molecules%list(ikind)%array(imol_local)
     650         700 :             molecule => molecule_set(imol)
     651             :             CALL get_molecule(molecule, first_atom=first_atom, &
     652         700 :                               last_atom=last_atom)
     653        2800 :             DO ipart = first_atom, last_atom
     654         700 :                vel(1, ipart) = vel(1, ipart)*scale
     655         700 :                vel(2, ipart) = vel(2, ipart)*scale
     656        1400 :                vel(3, ipart) = vel(3, ipart)*scale
     657             :             END DO
     658             :          END DO
     659             :       END DO
     660             : 
     661          20 :    END SUBROUTINE damp_v_velocity
     662             : 
     663             : ! **************************************************************************************************
     664             : !> \brief provides damping for barostat via nph_uniaxial_damped dynamics
     665             : !> \param npt ...
     666             : !> \param gamma1 ...
     667             : !> \param dt ...
     668             : !> \par History
     669             : !>      none
     670             : !> \author CJM
     671             : ! **************************************************************************************************
     672          80 :    ELEMENTAL SUBROUTINE damp_veps(npt, gamma1, dt)
     673             : 
     674             :       TYPE(npt_info_type), INTENT(INOUT)                 :: npt
     675             :       REAL(KIND=dp), INTENT(IN)                          :: gamma1, dt
     676             : 
     677             :       REAL(KIND=dp)                                      :: scale
     678             : 
     679          80 :       scale = 1.0_dp
     680          80 :       scale = scale*EXP(-gamma1*0.25_dp*dt)
     681             :       ! Scale
     682          80 :       npt%v = npt%v*scale
     683             : 
     684          80 :    END SUBROUTINE damp_veps
     685             : 
     686             : ! **************************************************************************************************
     687             : !> \brief update veps using multiplier obtained from SHAKE
     688             : !> \param old ...
     689             : !> \param gci ...
     690             : !> \param atomic_kind_set ...
     691             : !> \param particle_set ...
     692             : !> \param local_particles ...
     693             : !> \param molecule_kind_set ...
     694             : !> \param molecule_set ...
     695             : !> \param local_molecules ...
     696             : !> \param vel ...
     697             : !> \param dt ...
     698             : !> \param cell ...
     699             : !> \param npt ...
     700             : !> \param simpar ...
     701             : !> \param virial ...
     702             : !> \param vector_v ...
     703             : !> \param roll_tol ...
     704             : !> \param iroll ...
     705             : !> \param infree ...
     706             : !> \param first ...
     707             : !> \param para_env ...
     708             : !> \param u ...
     709             : !> \par History
     710             : !>      none
     711             : !> \author CJM
     712             : ! **************************************************************************************************
     713        1794 :    SUBROUTINE rattle_roll_setup(old, gci, atomic_kind_set, particle_set, local_particles, &
     714        1794 :                                 molecule_kind_set, molecule_set, local_molecules, vel, dt, cell, npt, simpar, virial, &
     715        1794 :                                 vector_v, roll_tol, iroll, infree, first, para_env, u)
     716             : 
     717             :       TYPE(old_variables_type), POINTER                  :: old
     718             :       TYPE(global_constraint_type), POINTER              :: gci
     719             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
     720             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     721             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     722             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     723             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     724             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     725             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     726             :       REAL(KIND=dp), INTENT(IN)                          :: dt
     727             :       TYPE(cell_type), POINTER                           :: cell
     728             :       TYPE(npt_info_type), DIMENSION(:, :), POINTER      :: npt
     729             :       TYPE(simpar_type), INTENT(IN)                      :: simpar
     730             :       TYPE(virial_type), POINTER                         :: virial
     731             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: vector_v
     732             :       REAL(KIND=dp), INTENT(OUT)                         :: roll_tol
     733             :       INTEGER, INTENT(INOUT)                             :: iroll
     734             :       REAL(KIND=dp), INTENT(IN)                          :: infree
     735             :       LOGICAL, INTENT(INOUT)                             :: first
     736             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     737             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: u(:, :)
     738             : 
     739             :       REAL(KIND=dp)                                      :: kin
     740             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_kin
     741             :       TYPE(npt_info_type), DIMENSION(3, 3)               :: npt_loc
     742             : 
     743        1794 :       IF (first) THEN
     744             :          CALL update_pv(gci, simpar, atomic_kind_set, vel, particle_set, &
     745             :                         local_molecules, molecule_set, molecule_kind_set, &
     746         678 :                         local_particles, kin, pv_kin, virial, para_env)
     747         678 :          CALL update_veps(cell, npt, simpar, pv_kin, kin, virial, infree)
     748             : 
     749             :       END IF
     750        1794 :       first = .FALSE.
     751             : 
     752             :       ! assigning local variable
     753        1794 :       SELECT CASE (simpar%ensemble)
     754             :       CASE (npt_i_ensemble, npt_ia_ensemble)
     755       23062 :          npt_loc(:, :)%v = 0.0_dp
     756       23062 :          npt_loc(:, :)%mass = 0.0_dp
     757        1774 :          npt_loc(1, 1)%v = npt(1, 1)%v
     758        1774 :          npt_loc(2, 2)%v = npt(1, 1)%v
     759        1774 :          npt_loc(3, 3)%v = npt(1, 1)%v
     760             :          npt_loc(1, 1)%mass = npt(1, 1)%mass
     761             :          npt_loc(2, 2)%mass = npt(1, 1)%mass
     762        1774 :          npt_loc(3, 3)%mass = npt(1, 1)%mass
     763             :       CASE (npt_f_ensemble)
     764        2034 :          npt_loc = npt
     765             :       END SELECT
     766             : 
     767             :       ! resetting
     768        1794 :       CALL set(old, atomic_kind_set, particle_set, vel, local_particles, cell, npt, 'B')
     769             :       CALL rattle_roll_control(gci, local_molecules, molecule_set, molecule_kind_set, &
     770             :                                particle_set, vel, dt, simpar, vector_v, npt_loc%v, roll_tol, iroll, &
     771       46624 :                                para_env, u, cell, local_particles)
     772             : 
     773        1794 :    END SUBROUTINE rattle_roll_setup
     774             : 
     775             : ! **************************************************************************************************
     776             : !> \brief Overloaded routine to compute veps given the particles
     777             : !>      structure or a local copy of the velocity array
     778             : !> \param gci ...
     779             : !> \param simpar ...
     780             : !> \param atomic_kind_set ...
     781             : !> \param particle_set ...
     782             : !> \param local_molecules ...
     783             : !> \param molecule_set ...
     784             : !> \param molecule_kind_set ...
     785             : !> \param local_particles ...
     786             : !> \param kin ...
     787             : !> \param pv_kin ...
     788             : !> \param virial ...
     789             : !> \param int_group ...
     790             : !> \par History
     791             : !>      none
     792             : !> \author CJM
     793             : ! **************************************************************************************************
     794        3668 :    SUBROUTINE update_pv_particle_set(gci, simpar, atomic_kind_set, particle_set, &
     795             :                                      local_molecules, molecule_set, molecule_kind_set, &
     796             :                                      local_particles, kin, pv_kin, virial, int_group)
     797             :       TYPE(global_constraint_type), POINTER              :: gci
     798             :       TYPE(simpar_type), INTENT(IN)                      :: simpar
     799             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
     800             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     801             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     802             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     803             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     804             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     805             :       REAL(KIND=dp), INTENT(OUT)                         :: kin
     806             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT)        :: pv_kin
     807             :       TYPE(virial_type), INTENT(INOUT)                   :: virial
     808             : 
     809             :       CLASS(mp_comm_type), INTENT(IN)                     :: int_group
     810             : 
     811             :       INTEGER                                            :: i, iparticle, iparticle_kind, &
     812             :                                                             iparticle_local, j, nparticle_local
     813             :       REAL(KIND=dp)                                      :: mass
     814             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     815             : 
     816        3668 :       pv_kin = 0.0_dp
     817        3668 :       kin = 0.0_dp
     818       11584 :       DO iparticle_kind = 1, SIZE(atomic_kind_set)
     819        7916 :          atomic_kind => atomic_kind_set(iparticle_kind)
     820        7916 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     821        7916 :          nparticle_local = local_particles%n_el(iparticle_kind)
     822      604269 :          DO iparticle_local = 1, nparticle_local
     823      592685 :             iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     824     2378656 :             DO i = 1, 3
     825     7112220 :                DO j = 1, 3
     826             :                   pv_kin(i, j) = pv_kin(i, j) + &
     827             :                                  mass*particle_set(iparticle)%v(i)* &
     828     7112220 :                                  particle_set(iparticle)%v(j)
     829             :                END DO
     830             :                kin = kin + mass*particle_set(iparticle)%v(i)* &
     831     2370740 :                      particle_set(iparticle)%v(i)
     832             :             END DO
     833             :          END DO
     834             :       END DO
     835             : 
     836        3668 :       CALL int_group%sum(pv_kin)
     837        3668 :       CALL int_group%sum(kin)
     838             : 
     839             :       ! updating the constraint virial
     840        3668 :       IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
     841             :                                                 molecule_set, &
     842             :                                                 molecule_kind_set, &
     843             :                                                 particle_set, virial, &
     844        1826 :                                                 int_group)
     845             : 
     846        3668 :    END SUBROUTINE update_pv_particle_set
     847             : 
     848             : ! **************************************************************************************************
     849             : !> \brief Overloaded routine to compute kinetic virials given the particles
     850             : !>      structure or a local copy of the velocity array
     851             : !> \param gci ...
     852             : !> \param simpar ...
     853             : !> \param atomic_kind_set ...
     854             : !> \param vel ...
     855             : !> \param particle_set ...
     856             : !> \param local_molecules ...
     857             : !> \param molecule_set ...
     858             : !> \param molecule_kind_set ...
     859             : !> \param local_particles ...
     860             : !> \param kin ...
     861             : !> \param pv_kin ...
     862             : !> \param virial ...
     863             : !> \param int_group ...
     864             : !> \par History
     865             : !>      none
     866             : !> \author CJM
     867             : ! **************************************************************************************************
     868        4314 :    SUBROUTINE update_pv_velocity(gci, simpar, atomic_kind_set, vel, particle_set, &
     869             :                                  local_molecules, molecule_set, molecule_kind_set, &
     870        4314 :                                  local_particles, kin, pv_kin, virial, int_group)
     871             :       TYPE(global_constraint_type), POINTER              :: gci
     872             :       TYPE(simpar_type), INTENT(IN)                      :: simpar
     873             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
     874             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     875             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     876             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     877             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     878             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     879             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     880             :       REAL(KIND=dp), INTENT(OUT)                         :: kin
     881             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: pv_kin
     882             :       TYPE(virial_type), INTENT(INOUT)                   :: virial
     883             : 
     884             :       CLASS(mp_comm_type), INTENT(IN)                     :: int_group
     885             : 
     886             :       INTEGER                                            :: i, iparticle, iparticle_kind, &
     887             :                                                             iparticle_local, j, nparticle_local
     888             :       REAL(KIND=dp)                                      :: mass
     889             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     890             : 
     891       56082 :       pv_kin = 0.0_dp
     892        4314 :       kin = 0._dp
     893       13526 :       DO iparticle_kind = 1, SIZE(atomic_kind_set)
     894        9212 :          atomic_kind => atomic_kind_set(iparticle_kind)
     895        9212 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     896        9212 :          nparticle_local = local_particles%n_el(iparticle_kind)
     897      683991 :          DO iparticle_local = 1, nparticle_local
     898      670465 :             iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     899     2691072 :             DO i = 1, 3
     900     8045580 :                DO j = 1, 3
     901             :                   pv_kin(i, j) = pv_kin(i, j) + &
     902     8045580 :                                  mass*vel(i, iparticle)*vel(j, iparticle)
     903             :                END DO
     904     2681860 :                kin = kin + mass*vel(i, iparticle)*vel(i, iparticle)
     905             :             END DO
     906             :          END DO
     907             :       END DO
     908             : 
     909      107850 :       CALL int_group%sum(pv_kin)
     910        4314 :       CALL int_group%sum(kin)
     911             : 
     912             :       ! updating the constraint virial
     913        4314 :       IF (simpar%constraint) CALL pv_constraint(gci, local_molecules, &
     914             :                                                 molecule_set, &
     915             :                                                 molecule_kind_set, &
     916             :                                                 particle_set, virial, &
     917        2472 :                                                 int_group)
     918             : 
     919        4314 :    END SUBROUTINE update_pv_velocity
     920             : 
     921             : ! **************************************************************************************************
     922             : !> \brief Routine to compute veps
     923             : !> \param box ...
     924             : !> \param npt ...
     925             : !> \param simpar ...
     926             : !> \param pv_kin ...
     927             : !> \param kin ...
     928             : !> \param virial ...
     929             : !> \param infree ...
     930             : !> \param virial_components ...
     931             : !> \par History
     932             : !>      none
     933             : !> \author CJM
     934             : ! **************************************************************************************************
     935        7982 :    SUBROUTINE update_veps(box, npt, simpar, pv_kin, kin, virial, infree, virial_components)
     936             : 
     937             :       TYPE(cell_type), POINTER                           :: box
     938             :       TYPE(npt_info_type), DIMENSION(:, :), &
     939             :          INTENT(INOUT)                                   :: npt
     940             :       TYPE(simpar_type), INTENT(IN)                      :: simpar
     941             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: pv_kin
     942             :       REAL(KIND=dp), INTENT(IN)                          :: kin
     943             :       TYPE(virial_type), INTENT(INOUT)                   :: virial
     944             :       REAL(KIND=dp), INTENT(IN)                          :: infree
     945             :       INTEGER, INTENT(IN), OPTIONAL                      :: virial_components
     946             : 
     947             :       INTEGER                                            :: ii
     948             :       REAL(KIND=dp)                                      :: fdotr, trace, v, v0, v0i, vi
     949             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: unit
     950             : 
     951        7982 :       CPASSERT(ASSOCIATED(box))
     952             : 
     953             :       ! Initialize unit
     954             : 
     955        7982 :       unit = 0.0_dp
     956        7982 :       unit(1, 1) = 1.0_dp
     957        7982 :       unit(2, 2) = 1.0_dp
     958        7982 :       unit(3, 3) = 1.0_dp
     959             : 
     960             :       IF (simpar%ensemble == npt_i_ensemble .OR. &
     961        7982 :           simpar%ensemble == npt_ia_ensemble .OR. &
     962             :           simpar%ensemble == npe_i_ensemble) THEN
     963             :          ! get force on barostat
     964             :          fdotr = 0.0_dp
     965       24160 :          DO ii = 1, 3
     966             :             fdotr = fdotr + virial%pv_virial(ii, ii) + &
     967       24160 :                     virial%pv_constraint(ii, ii)
     968             :          END DO
     969             : 
     970             :          npt(:, :)%f = (1.0_dp + (3.0_dp*infree))*kin + fdotr - &
     971       18120 :                        3.0_dp*simpar%p_ext*box%deth
     972        1942 :       ELSEIF (simpar%ensemble == npt_f_ensemble .OR. &
     973             :               simpar%ensemble == npe_f_ensemble) THEN
     974             :          npt(:, :)%f = virial%pv_virial(:, :) + &
     975             :                        pv_kin(:, :) + virial%pv_constraint(:, :) - &
     976             :                        unit(:, :)*simpar%p_ext*box%deth + &
     977       23686 :                        infree*kin*unit(:, :)
     978             :          IF (debug_isotropic_limit) THEN
     979             :             trace = npt(1, 1)%f + npt(2, 2)%f + npt(3, 3)%f
     980             :             trace = trace/3.0_dp
     981             :             npt(:, :)%f = trace*unit(:, :)
     982             :          END IF
     983         120 :       ELSEIF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
     984             :               simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
     985         120 :          v = box%deth
     986         120 :          vi = 1._dp/v
     987         120 :          v0 = simpar%v0
     988         120 :          v0i = 1._dp/v0
     989             :          ! orthorhombic box ONLY
     990             :          ! Chooses only the compressive solution
     991         120 :          IF (v < v0) THEN
     992             :             npt(1, 1)%f = virial%pv_virial(1, 1) + &
     993             :                           pv_kin(1, 1) + virial%pv_constraint(1, 1) - &
     994             :                           simpar%p0*v - simpar%v_shock*simpar%v_shock* &
     995           0 :                           v*v0i*(1._dp - v*v0i) + infree*kin
     996             :          ELSE
     997             :             npt(1, 1)%f = virial%pv_virial(1, 1) + &
     998             :                           pv_kin(1, 1) + virial%pv_constraint(1, 1) - &
     999         120 :                           simpar%p0*v + infree*kin
    1000             :          END IF
    1001             :          IF (debug_uniaxial_limit) THEN
    1002             :             ! orthorhombic box ONLY
    1003             :             npt(1, 1)%f = virial%pv_virial(1, 1) + &
    1004             :                           pv_kin(1, 1) + virial%pv_constraint(1, 1) - &
    1005             :                           simpar%p0*box%deth + infree*kin
    1006             :          END IF
    1007             :       END IF
    1008             : 
    1009             :       ! update barostat velocities
    1010             :       npt(:, :)%v = npt(:, :)%v + &
    1011       42166 :                     0.5_dp*simpar%dt*npt(:, :)%f/npt(:, :)%mass
    1012             : 
    1013             :       ! Screen the dynamics of the barostat according user request
    1014        7982 :       IF (PRESENT(virial_components)) THEN
    1015        1812 :          SELECT CASE (virial_components)
    1016             :          CASE (do_clv_xyz)
    1017             :             ! Do nothing..
    1018             :          CASE (do_clv_x)
    1019           0 :             npt(1, 2)%v = 0.0_dp
    1020           0 :             npt(1, 3)%v = 0.0_dp
    1021           0 :             npt(1, 2)%f = 0.0_dp
    1022           0 :             npt(1, 3)%f = 0.0_dp
    1023           0 :             npt(2, 1)%v = 0.0_dp
    1024           0 :             npt(2, 2)%v = 0.0_dp
    1025           0 :             npt(2, 3)%v = 0.0_dp
    1026           0 :             npt(2, 1)%f = 0.0_dp
    1027           0 :             npt(2, 2)%f = 0.0_dp
    1028           0 :             npt(2, 3)%f = 0.0_dp
    1029           0 :             npt(3, 1)%v = 0.0_dp
    1030           0 :             npt(3, 2)%v = 0.0_dp
    1031           0 :             npt(3, 3)%v = 0.0_dp
    1032           0 :             npt(3, 1)%f = 0.0_dp
    1033           0 :             npt(3, 2)%f = 0.0_dp
    1034           0 :             npt(3, 3)%f = 0.0_dp
    1035             :          CASE (do_clv_y)
    1036           0 :             npt(1, 1)%v = 0.0_dp
    1037           0 :             npt(1, 2)%v = 0.0_dp
    1038           0 :             npt(1, 3)%v = 0.0_dp
    1039           0 :             npt(1, 1)%f = 0.0_dp
    1040           0 :             npt(1, 2)%f = 0.0_dp
    1041           0 :             npt(1, 3)%f = 0.0_dp
    1042           0 :             npt(2, 1)%v = 0.0_dp
    1043           0 :             npt(2, 3)%v = 0.0_dp
    1044           0 :             npt(2, 1)%f = 0.0_dp
    1045           0 :             npt(2, 3)%f = 0.0_dp
    1046           0 :             npt(3, 1)%v = 0.0_dp
    1047           0 :             npt(3, 2)%v = 0.0_dp
    1048           0 :             npt(3, 3)%v = 0.0_dp
    1049           0 :             npt(3, 1)%f = 0.0_dp
    1050           0 :             npt(3, 2)%f = 0.0_dp
    1051           0 :             npt(3, 3)%f = 0.0_dp
    1052             :          CASE (do_clv_z)
    1053          80 :             npt(1, 1)%v = 0.0_dp
    1054          80 :             npt(1, 2)%v = 0.0_dp
    1055          80 :             npt(1, 3)%v = 0.0_dp
    1056          80 :             npt(1, 1)%f = 0.0_dp
    1057          80 :             npt(1, 2)%f = 0.0_dp
    1058          80 :             npt(1, 3)%f = 0.0_dp
    1059          80 :             npt(2, 1)%v = 0.0_dp
    1060          80 :             npt(2, 2)%v = 0.0_dp
    1061          80 :             npt(2, 3)%v = 0.0_dp
    1062          80 :             npt(2, 1)%f = 0.0_dp
    1063          80 :             npt(2, 2)%f = 0.0_dp
    1064          80 :             npt(2, 3)%f = 0.0_dp
    1065          80 :             npt(3, 1)%v = 0.0_dp
    1066          80 :             npt(3, 2)%v = 0.0_dp
    1067          80 :             npt(3, 1)%f = 0.0_dp
    1068          80 :             npt(3, 2)%f = 0.0_dp
    1069             :          CASE (do_clv_xy)
    1070           0 :             npt(1, 3)%v = 0.0_dp
    1071           0 :             npt(1, 3)%f = 0.0_dp
    1072           0 :             npt(2, 3)%v = 0.0_dp
    1073           0 :             npt(2, 3)%f = 0.0_dp
    1074           0 :             npt(3, 1)%v = 0.0_dp
    1075           0 :             npt(3, 2)%v = 0.0_dp
    1076           0 :             npt(3, 3)%v = 0.0_dp
    1077           0 :             npt(3, 1)%f = 0.0_dp
    1078           0 :             npt(3, 2)%f = 0.0_dp
    1079           0 :             npt(3, 3)%f = 0.0_dp
    1080             :          CASE (do_clv_xz)
    1081           0 :             npt(1, 2)%v = 0.0_dp
    1082           0 :             npt(1, 2)%f = 0.0_dp
    1083           0 :             npt(2, 1)%v = 0.0_dp
    1084           0 :             npt(2, 2)%v = 0.0_dp
    1085           0 :             npt(2, 3)%v = 0.0_dp
    1086           0 :             npt(2, 1)%f = 0.0_dp
    1087           0 :             npt(2, 2)%f = 0.0_dp
    1088           0 :             npt(2, 3)%f = 0.0_dp
    1089           0 :             npt(3, 2)%v = 0.0_dp
    1090           0 :             npt(3, 2)%f = 0.0_dp
    1091             :          CASE (do_clv_yz)
    1092           0 :             npt(1, 1)%v = 0.0_dp
    1093           0 :             npt(1, 2)%v = 0.0_dp
    1094           0 :             npt(1, 3)%v = 0.0_dp
    1095           0 :             npt(1, 1)%f = 0.0_dp
    1096           0 :             npt(1, 2)%f = 0.0_dp
    1097           0 :             npt(1, 3)%f = 0.0_dp
    1098           0 :             npt(2, 1)%v = 0.0_dp
    1099           0 :             npt(2, 1)%f = 0.0_dp
    1100           0 :             npt(3, 1)%v = 0.0_dp
    1101        1812 :             npt(3, 1)%f = 0.0_dp
    1102             :          END SELECT
    1103             :       END IF
    1104             : 
    1105        7982 :    END SUBROUTINE update_veps
    1106             : 
    1107             : ! **************************************************************************************************
    1108             : !> \brief First half of the velocity-verlet algorithm : update velocity by half
    1109             : !>        step and positions by full step
    1110             : !> \param tmp ...
    1111             : !> \param atomic_kind_set ...
    1112             : !> \param local_particles ...
    1113             : !> \param particle_set ...
    1114             : !> \param core_particle_set ...
    1115             : !> \param shell_particle_set ...
    1116             : !> \param nparticle_kind ...
    1117             : !> \param shell_adiabatic ...
    1118             : !> \param dt ...
    1119             : !> \param u ...
    1120             : !> \param lfixd_list ...
    1121             : !> \par History
    1122             : !>      none
    1123             : !> \author MI (February 2008)
    1124             : ! **************************************************************************************************
    1125       44235 :    SUBROUTINE vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1126             :                        core_particle_set, shell_particle_set, nparticle_kind, &
    1127       44235 :                        shell_adiabatic, dt, u, lfixd_list)
    1128             : 
    1129             :       TYPE(tmp_variables_type), POINTER                  :: tmp
    1130             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1131             :       TYPE(distribution_1d_type), POINTER                :: local_particles
    1132             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set, core_particle_set, &
    1133             :                                                             shell_particle_set
    1134             :       INTEGER, INTENT(IN)                                :: nparticle_kind
    1135             :       LOGICAL, INTENT(IN)                                :: shell_adiabatic
    1136             :       REAL(KIND=dp)                                      :: dt
    1137             :       REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL           :: u
    1138             :       TYPE(local_fixd_constraint_type), DIMENSION(:), &
    1139             :          OPTIONAL                                        :: lfixd_list
    1140             : 
    1141             :       INTEGER                                            :: iparticle, iparticle_kind, &
    1142             :                                                             iparticle_local, nparticle_local, &
    1143             :                                                             shell_index
    1144             :       LOGICAL                                            :: is_shell
    1145             :       REAL(KIND=dp)                                      :: dm, dmc, dms, dsc(3), dvsc(3), &
    1146             :                                                             fac_massc, fac_masss, mass
    1147             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1148             :       TYPE(shell_kind_type), POINTER                     :: shell
    1149             : 
    1150       44235 :       NULLIFY (atomic_kind, shell)
    1151       44235 :       tmp%max_vel = 0.0_dp
    1152       44235 :       tmp%max_vel_sc = 0.0_dp
    1153       44235 :       tmp%max_dvel = 0.0_dp
    1154       44235 :       tmp%max_dvel_sc = 0.0_dp
    1155       44235 :       tmp%max_dr = 0.0_dp
    1156       44235 :       tmp%max_dsc = 0.0_dp
    1157       44235 :       dsc = 0.0_dp
    1158       44235 :       dvsc = 0.0_dp
    1159             : 
    1160             :       !     *** Velocity Verlet (first part) ***
    1161      154510 :       DO iparticle_kind = 1, nparticle_kind
    1162      110275 :          atomic_kind => atomic_kind_set(iparticle_kind)
    1163      110275 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell_active=is_shell)
    1164      154510 :          IF (mass /= 0.0_dp) THEN
    1165      109869 :             dm = 0.5_dp*dt/mass
    1166      109869 :             IF (is_shell .AND. shell_adiabatic) THEN
    1167        5412 :                CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
    1168        5412 :                dms = 0.5_dp*dt/shell%mass_shell
    1169        5412 :                dmc = 0.5_dp*dt/shell%mass_core
    1170        5412 :                fac_masss = shell%mass_shell/mass
    1171        5412 :                fac_massc = shell%mass_core/mass
    1172        5412 :                nparticle_local = local_particles%n_el(iparticle_kind)
    1173        5412 :                IF (PRESENT(u)) THEN
    1174       54620 :                   DO iparticle_local = 1, nparticle_local
    1175       52704 :                      iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1176       52704 :                      shell_index = particle_set(iparticle)%shell_index
    1177             :                      ! Transform positions and velocities and forces of the shell
    1178             :                      CALL transform_first(shell_particle_set, tmp%shell_pos, tmp%shell_vel, &
    1179       52704 :                                           shell_index, u, dms, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
    1180             : 
    1181             :                      ! Transform positions and velocities and forces of the core
    1182             :                      CALL transform_first(core_particle_set, tmp%core_pos, tmp%core_vel, &
    1183       52704 :                                           shell_index, u, dmc, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
    1184             : 
    1185             :                      ! Derive velocities and positions of the COM
    1186             :                      tmp%vel(:, iparticle) = fac_masss*tmp%shell_vel(:, shell_index) + &
    1187      210816 :                                              fac_massc*tmp%core_vel(:, shell_index)
    1188             :                      tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
    1189      210816 :                                              fac_massc*tmp%core_pos(:, shell_index)
    1190             : 
    1191             :                      tmp%max_vel = MAX(tmp%max_vel, ABS(tmp%vel(1, iparticle)), &
    1192       52704 :                                        ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
    1193             :                      tmp%max_vel_sc = MAX(tmp%max_vel_sc, &
    1194             :                                           ABS(tmp%shell_vel(1, shell_index) - tmp%core_vel(1, shell_index)), &
    1195             :                                           ABS(tmp%shell_vel(2, shell_index) - tmp%core_vel(2, shell_index)), &
    1196       52704 :                                           ABS(tmp%shell_vel(3, shell_index) - tmp%core_vel(3, shell_index)))
    1197             :                      tmp%max_dr = MAX(tmp%max_dr, &
    1198             :                                       ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
    1199             :                                       ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
    1200       52704 :                                       ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
    1201             :                      tmp%max_dvel = MAX(tmp%max_dvel, &
    1202             :                                         ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
    1203             :                                         ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
    1204       52704 :                                         ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
    1205             : 
    1206             :                      dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
    1207      210816 :                               shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
    1208       52704 :                      tmp%max_dsc = MAX(tmp%max_dsc, ABS(dsc(1)), ABS(dsc(2)), ABS(dsc(3)))
    1209             :                      dvsc(:) = tmp%shell_vel(:, shell_index) - tmp%core_vel(:, shell_index) - &
    1210      210816 :                                shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
    1211       54620 :                      tmp%max_dvel_sc = MAX(tmp%max_dvel_sc, ABS(dvsc(1)), ABS(dvsc(2)), ABS(dvsc(3)))
    1212             :                   END DO ! iparticle_local
    1213             :                ELSE
    1214       88530 :                   DO iparticle_local = 1, nparticle_local
    1215       85034 :                      iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1216       85034 :                      shell_index = particle_set(iparticle)%shell_index
    1217             :                      tmp%shell_vel(:, shell_index) = &
    1218             :                         shell_particle_set(shell_index)%v(:)*tmp%scale_v(:)*tmp%scale_v(:) + &
    1219      680272 :                         tmp%scale_v(:)*tmp%poly_v(:)*dms*shell_particle_set(shell_index)%f(:)
    1220             :                      tmp%shell_pos(:, shell_index) = &
    1221             :                         shell_particle_set(shell_index)%r(:)*tmp%scale_r(:)*tmp%scale_r(:) + &
    1222      680272 :                         tmp%scale_r(:)*tmp%poly_r(:)*dt*tmp%shell_vel(:, shell_index)
    1223             :                      tmp%core_vel(:, shell_index) = &
    1224             :                         core_particle_set(shell_index)%v(:)*tmp%scale_v(:)*tmp%scale_v(:) + &
    1225      680272 :                         tmp%scale_v(:)*tmp%poly_v(:)*dmc*core_particle_set(shell_index)%f(:)
    1226             :                      tmp%core_pos(:, shell_index) = &
    1227             :                         core_particle_set(shell_index)%r(:)*tmp%scale_r(:)*tmp%scale_r(:) + &
    1228      680272 :                         tmp%scale_r(:)*tmp%poly_r(:)*dt*tmp%core_vel(:, shell_index)
    1229             : 
    1230             :                      tmp%vel(:, iparticle) = fac_masss*tmp%shell_vel(:, shell_index) + &
    1231      340136 :                                              fac_massc*tmp%core_vel(:, shell_index)
    1232             :                      tmp%pos(:, iparticle) = fac_masss*tmp%shell_pos(:, shell_index) + &
    1233      340136 :                                              fac_massc*tmp%core_pos(:, shell_index)
    1234             :                      tmp%max_vel = MAX(tmp%max_vel, &
    1235       85034 :                                        ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
    1236             :                      tmp%max_vel_sc = MAX(tmp%max_vel_sc, &
    1237             :                                           ABS(tmp%shell_vel(1, shell_index) - tmp%core_vel(1, shell_index)), &
    1238             :                                           ABS(tmp%shell_vel(2, shell_index) - tmp%core_vel(2, shell_index)), &
    1239       85034 :                                           ABS(tmp%shell_vel(3, shell_index) - tmp%core_vel(3, shell_index)))
    1240             :                      tmp%max_dr = MAX(tmp%max_dr, &
    1241             :                                       ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
    1242             :                                       ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
    1243       85034 :                                       ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
    1244             :                      tmp%max_dvel = MAX(tmp%max_dvel, &
    1245             :                                         ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
    1246             :                                         ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
    1247       85034 :                                         ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
    1248             :                      dsc(:) = tmp%shell_pos(:, shell_index) - tmp%core_pos(:, shell_index) - &
    1249      340136 :                               shell_particle_set(shell_index)%r(:) + core_particle_set(shell_index)%r(:)
    1250       85034 :                      tmp%max_dsc = MAX(tmp%max_dsc, ABS(dsc(1)), ABS(dsc(2)), ABS(dsc(3)))
    1251             :                      dvsc(:) = tmp%shell_vel(:, shell_index) - tmp%core_vel(:, shell_index) - &
    1252      340136 :                                shell_particle_set(shell_index)%v(:) + core_particle_set(shell_index)%v(:)
    1253       88530 :                      tmp%max_dvel_sc = MAX(tmp%max_dvel_sc, ABS(dvsc(1)), ABS(dvsc(2)), ABS(dvsc(3)))
    1254             :                   END DO ! iparticle_local
    1255             :                END IF
    1256             :             ELSE
    1257      104457 :                nparticle_local = local_particles%n_el(iparticle_kind)
    1258      104457 :                IF (PRESENT(u)) THEN
    1259       44038 :                   DO iparticle_local = 1, nparticle_local
    1260       43412 :                      iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1261             :                      ! Transform positions and velocities and forces
    1262             :                      CALL transform_first(particle_set, tmp%pos, tmp%vel, &
    1263       43412 :                                           iparticle, u, dm, dt, tmp%poly_v, tmp%poly_r, tmp%scale_v, tmp%scale_r)
    1264             : 
    1265             :                      tmp%max_vel = MAX(tmp%max_vel, &
    1266       43412 :                                        ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
    1267             :                      tmp%max_dr = MAX(tmp%max_dr, ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
    1268             :                                       ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
    1269       43412 :                                       ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
    1270             :                      tmp%max_dvel = MAX(tmp%max_dvel, &
    1271             :                                         ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
    1272             :                                         ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
    1273       44038 :                                         ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
    1274             :                   END DO ! iparticle_local
    1275             :                ELSE
    1276     3174720 :                   DO iparticle_local = 1, nparticle_local
    1277     3070889 :                      iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1278             :                      tmp%vel(1, iparticle) = &
    1279             :                         particle_set(iparticle)%v(1)*tmp%scale_v(1)*tmp%scale_v(1) + &
    1280     3070889 :                         tmp%scale_v(1)*tmp%poly_v(1)*dm*particle_set(iparticle)%f(1)
    1281             :                      tmp%vel(2, iparticle) = &
    1282             :                         particle_set(iparticle)%v(2)*tmp%scale_v(2)*tmp%scale_v(2) + &
    1283     3070889 :                         tmp%scale_v(2)*tmp%poly_v(2)*dm*particle_set(iparticle)%f(2)
    1284             :                      tmp%vel(3, iparticle) = &
    1285             :                         particle_set(iparticle)%v(3)*tmp%scale_v(3)*tmp%scale_v(3) + &
    1286     3070889 :                         tmp%scale_v(3)*tmp%poly_v(3)*dm*particle_set(iparticle)%f(3)
    1287             :                      tmp%pos(1, iparticle) = &
    1288             :                         particle_set(iparticle)%r(1)*tmp%scale_r(1)*tmp%scale_r(1) + &
    1289     3070889 :                         tmp%scale_r(1)*tmp%poly_r(1)*dt*tmp%vel(1, iparticle)
    1290             :                      tmp%pos(2, iparticle) = &
    1291             :                         particle_set(iparticle)%r(2)*tmp%scale_r(2)*tmp%scale_r(2) + &
    1292     3070889 :                         tmp%scale_r(2)*tmp%poly_r(2)*dt*tmp%vel(2, iparticle)
    1293             :                      tmp%pos(3, iparticle) = &
    1294             :                         particle_set(iparticle)%r(3)*tmp%scale_r(3)*tmp%scale_r(3) + &
    1295     3070889 :                         tmp%scale_r(3)*tmp%poly_r(3)*dt*tmp%vel(3, iparticle)
    1296             : 
    1297             :                      ! overwrite positions of frozen particles
    1298     3070889 :                      IF (PRESENT(lfixd_list)) THEN
    1299      251570 :                         IF (ANY(lfixd_list(:)%ifixd_index == iparticle)) THEN
    1300         200 :                            tmp%pos(1, iparticle) = particle_set(iparticle)%r(1)
    1301         200 :                            tmp%pos(2, iparticle) = particle_set(iparticle)%r(2)
    1302         200 :                            tmp%pos(3, iparticle) = particle_set(iparticle)%r(3)
    1303             :                         END IF
    1304             :                      END IF
    1305             : 
    1306             :                      tmp%max_vel = MAX(tmp%max_vel, &
    1307     3070889 :                                        ABS(tmp%vel(1, iparticle)), ABS(tmp%vel(2, iparticle)), ABS(tmp%vel(3, iparticle)))
    1308             :                      tmp%max_dr = MAX(tmp%max_dr, &
    1309             :                                       ABS(particle_set(iparticle)%r(1) - tmp%pos(1, iparticle)), &
    1310             :                                       ABS(particle_set(iparticle)%r(2) - tmp%pos(2, iparticle)), &
    1311     3070889 :                                       ABS(particle_set(iparticle)%r(3) - tmp%pos(3, iparticle)))
    1312             :                      tmp%max_dvel = MAX(tmp%max_dvel, &
    1313             :                                         ABS(particle_set(iparticle)%v(1) - tmp%vel(1, iparticle)), &
    1314             :                                         ABS(particle_set(iparticle)%v(2) - tmp%vel(2, iparticle)), &
    1315     3174720 :                                         ABS(particle_set(iparticle)%v(3) - tmp%vel(3, iparticle)))
    1316             :                   END DO
    1317             :                END IF
    1318             :             END IF
    1319             :          END IF
    1320             :       END DO
    1321       44235 :    END SUBROUTINE vv_first
    1322             : 
    1323             : ! **************************************************************************************************
    1324             : !> \brief Second half of the velocity-verlet algorithm : update velocity by half
    1325             : !>        step  using the new forces
    1326             : !> \param tmp ...
    1327             : !> \param atomic_kind_set ...
    1328             : !> \param local_particles ...
    1329             : !> \param particle_set ...
    1330             : !> \param core_particle_set ...
    1331             : !> \param shell_particle_set ...
    1332             : !> \param nparticle_kind ...
    1333             : !> \param shell_adiabatic ...
    1334             : !> \param dt ...
    1335             : !> \param u ...
    1336             : !> \par History
    1337             : !>      none
    1338             : !> \author MI (February 2008)
    1339             : ! **************************************************************************************************
    1340       43761 :    SUBROUTINE vv_second(tmp, atomic_kind_set, local_particles, particle_set, &
    1341             :                         core_particle_set, shell_particle_set, nparticle_kind, &
    1342             :                         shell_adiabatic, dt, u)
    1343             : 
    1344             :       TYPE(tmp_variables_type), POINTER                  :: tmp
    1345             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1346             :       TYPE(distribution_1d_type), POINTER                :: local_particles
    1347             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set, core_particle_set, &
    1348             :                                                             shell_particle_set
    1349             :       INTEGER, INTENT(IN)                                :: nparticle_kind
    1350             :       LOGICAL, INTENT(IN)                                :: shell_adiabatic
    1351             :       REAL(KIND=dp)                                      :: dt
    1352             :       REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL           :: u
    1353             : 
    1354             :       INTEGER                                            :: iparticle, iparticle_kind, &
    1355             :                                                             iparticle_local, nparticle_local, &
    1356             :                                                             shell_index
    1357             :       LOGICAL                                            :: is_shell
    1358             :       REAL(KIND=dp)                                      :: dm, dmc, dms, fac_massc, fac_masss, mass
    1359             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1360             :       TYPE(shell_kind_type), POINTER                     :: shell
    1361             : 
    1362       43761 :       NULLIFY (atomic_kind, shell)
    1363             : 
    1364             :       !     *** Velocity Verlet (second part) ***
    1365      152828 :       DO iparticle_kind = 1, nparticle_kind
    1366      109067 :          atomic_kind => atomic_kind_set(iparticle_kind)
    1367             :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
    1368      109067 :                               shell_active=is_shell)
    1369      152828 :          IF (mass /= 0.0_dp) THEN
    1370      108661 :             dm = 0.5_dp*dt/mass
    1371      108661 :             IF (is_shell .AND. shell_adiabatic) THEN
    1372        4520 :                CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
    1373        4520 :                dms = 0.5_dp*dt/shell%mass_shell
    1374        4520 :                dmc = 0.5_dp*dt/shell%mass_core
    1375        4520 :                fac_masss = shell%mass_shell/mass
    1376        4520 :                fac_massc = shell%mass_core/mass
    1377        4520 :                nparticle_local = local_particles%n_el(iparticle_kind)
    1378        4520 :                IF (PRESENT(u)) THEN
    1379       35860 :                   DO iparticle_local = 1, nparticle_local
    1380       34560 :                      iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1381       34560 :                      shell_index = particle_set(iparticle)%shell_index
    1382             :                      ! Transform velocities and forces shell
    1383             :                      CALL transform_second(shell_particle_set, tmp%shell_vel, shell_index, &
    1384       34560 :                                            u, dms, tmp%poly_v, tmp%scale_v)
    1385             : 
    1386             :                      ! Transform velocities and forces core
    1387             :                      CALL transform_second(core_particle_set, tmp%core_vel, shell_index, &
    1388       34560 :                                            u, dmc, tmp%poly_v, tmp%scale_v)
    1389             : 
    1390             :                      ! Derive velocties of the COM
    1391             :                      tmp%vel(1, iparticle) = fac_masss*tmp%shell_vel(1, shell_index) + &
    1392       34560 :                                              fac_massc*tmp%core_vel(1, shell_index)
    1393             :                      tmp%vel(2, iparticle) = fac_masss*tmp%shell_vel(2, shell_index) + &
    1394       34560 :                                              fac_massc*tmp%core_vel(2, shell_index)
    1395             :                      tmp%vel(3, iparticle) = fac_masss*tmp%shell_vel(3, shell_index) + &
    1396       35860 :                                              fac_massc*tmp%core_vel(3, shell_index)
    1397             :                   END DO ! iparticle_local
    1398             :                ELSE
    1399       81630 :                   DO iparticle_local = 1, nparticle_local
    1400       78410 :                      iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1401       78410 :                      shell_index = particle_set(iparticle)%shell_index
    1402             :                      tmp%shell_vel(1, shell_index) = &
    1403             :                         tmp%shell_vel(1, shell_index)*tmp%scale_v(1)*tmp%scale_v(1) + &
    1404       78410 :                         tmp%scale_v(1)*tmp%poly_v(1)*dms*shell_particle_set(shell_index)%f(1)
    1405             :                      tmp%shell_vel(2, shell_index) = &
    1406             :                         tmp%shell_vel(2, shell_index)*tmp%scale_v(2)*tmp%scale_v(2) + &
    1407       78410 :                         tmp%scale_v(2)*tmp%poly_v(2)*dms*shell_particle_set(shell_index)%f(2)
    1408             :                      tmp%shell_vel(3, shell_index) = &
    1409             :                         tmp%shell_vel(3, shell_index)*tmp%scale_v(3)*tmp%scale_v(3) + &
    1410       78410 :                         tmp%scale_v(3)*tmp%poly_v(3)*dms*shell_particle_set(shell_index)%f(3)
    1411             :                      tmp%core_vel(1, shell_index) = &
    1412             :                         tmp%core_vel(1, shell_index)*tmp%scale_v(1)*tmp%scale_v(1) + &
    1413       78410 :                         tmp%scale_v(1)*tmp%poly_v(1)*dmc*core_particle_set(shell_index)%f(1)
    1414             :                      tmp%core_vel(2, shell_index) = &
    1415             :                         tmp%core_vel(2, shell_index)*tmp%scale_v(2)*tmp%scale_v(2) + &
    1416       78410 :                         tmp%scale_v(2)*tmp%poly_v(2)*dmc*core_particle_set(shell_index)%f(2)
    1417             :                      tmp%core_vel(3, shell_index) = &
    1418             :                         tmp%core_vel(3, shell_index)*tmp%scale_v(3)*tmp%scale_v(3) + &
    1419       78410 :                         tmp%scale_v(3)*tmp%poly_v(3)*dmc*core_particle_set(shell_index)%f(3)
    1420             : 
    1421             :                      tmp%vel(1, iparticle) = fac_masss*tmp%shell_vel(1, shell_index) + &
    1422       78410 :                                              fac_massc*tmp%core_vel(1, shell_index)
    1423             :                      tmp%vel(2, iparticle) = fac_masss*tmp%shell_vel(2, shell_index) + &
    1424       78410 :                                              fac_massc*tmp%core_vel(2, shell_index)
    1425             :                      tmp%vel(3, iparticle) = fac_masss*tmp%shell_vel(3, shell_index) + &
    1426       81630 :                                              fac_massc*tmp%core_vel(3, shell_index)
    1427             :                   END DO ! iparticle_local
    1428             :                END IF
    1429             :             ELSE
    1430      104141 :                nparticle_local = local_particles%n_el(iparticle_kind)
    1431      104141 :                IF (PRESENT(u)) THEN
    1432       44038 :                   DO iparticle_local = 1, nparticle_local
    1433       43412 :                      iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1434             :                      CALL transform_second(particle_set, tmp%vel, iparticle, &
    1435       44038 :                                            u, dm, tmp%poly_v, tmp%scale_v)
    1436             :                   END DO ! iparticle_local
    1437             :                ELSE
    1438     2944176 :                   DO iparticle_local = 1, nparticle_local
    1439     2840661 :                      iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
    1440             :                      tmp%vel(1, iparticle) = &
    1441             :                         tmp%vel(1, iparticle)*tmp%scale_v(1)*tmp%scale_v(1) + &
    1442     2840661 :                         tmp%scale_v(1)*tmp%poly_v(1)*dm*particle_set(iparticle)%f(1)
    1443             :                      tmp%vel(2, iparticle) = &
    1444             :                         tmp%vel(2, iparticle)*tmp%scale_v(2)*tmp%scale_v(2) + &
    1445     2840661 :                         tmp%scale_v(2)*tmp%poly_v(2)*dm*particle_set(iparticle)%f(2)
    1446             :                      tmp%vel(3, iparticle) = &
    1447             :                         tmp%vel(3, iparticle)*tmp%scale_v(3)*tmp%scale_v(3) + &
    1448     2944176 :                         tmp%scale_v(3)*tmp%poly_v(3)*dm*particle_set(iparticle)%f(3)
    1449             :                   END DO
    1450             :                END IF
    1451             :             END IF
    1452             :          END IF
    1453             :       END DO
    1454             : 
    1455       43761 :    END SUBROUTINE vv_second
    1456             : 
    1457             : ! **************************************************************************************************
    1458             : !> \brief Transform position and velocities for the first half of the
    1459             : !>        Velocity-Verlet integration in the npt_f ensemble
    1460             : !> \param particle_set ...
    1461             : !> \param pos ...
    1462             : !> \param vel ...
    1463             : !> \param index ...
    1464             : !> \param u ...
    1465             : !> \param dm ...
    1466             : !> \param dt ...
    1467             : !> \param poly_v ...
    1468             : !> \param poly_r ...
    1469             : !> \param scale_v ...
    1470             : !> \param scale_r ...
    1471             : !> \par History
    1472             : !>      none
    1473             : !> \author MI (February 2008)
    1474             : ! **************************************************************************************************
    1475      148820 :    SUBROUTINE transform_first(particle_set, pos, vel, index, u, dm, dt, poly_v, &
    1476             :                               poly_r, scale_v, scale_r)
    1477             : 
    1478             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1479             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pos, vel
    1480             :       INTEGER, INTENT(IN)                                :: index
    1481             :       REAL(KIND=dp), INTENT(IN)                          :: u(3, 3), dm, dt, poly_v(3), poly_r(3), &
    1482             :                                                             scale_v(3), scale_r(3)
    1483             : 
    1484             :       REAL(KIND=dp), DIMENSION(3)                        :: uf, ur, uv
    1485             : 
    1486      148820 :       ur = MATMUL(TRANSPOSE(u), particle_set(index)%r(:))
    1487      148820 :       uv = MATMUL(TRANSPOSE(u), particle_set(index)%v(:))
    1488      148820 :       uf = MATMUL(TRANSPOSE(u), particle_set(index)%f(:))
    1489             : 
    1490      148820 :       uv(1) = uv(1)*scale_v(1)*scale_v(1) + uf(1)*scale_v(1)*poly_v(1)*dm
    1491      148820 :       uv(2) = uv(2)*scale_v(2)*scale_v(2) + uf(2)*scale_v(2)*poly_v(2)*dm
    1492      148820 :       uv(3) = uv(3)*scale_v(3)*scale_v(3) + uf(3)*scale_v(3)*poly_v(3)*dm
    1493             : 
    1494      148820 :       ur(1) = ur(1)*scale_r(1)*scale_r(1) + uv(1)*scale_r(1)*poly_r(1)*dt
    1495      148820 :       ur(2) = ur(2)*scale_r(2)*scale_r(2) + uv(2)*scale_r(2)*poly_r(2)*dt
    1496      148820 :       ur(3) = ur(3)*scale_r(3)*scale_r(3) + uv(3)*scale_r(3)*poly_r(3)*dt
    1497             : 
    1498     2529940 :       pos(:, index) = MATMUL(u, ur)
    1499     2529940 :       vel(:, index) = MATMUL(u, uv)
    1500             : 
    1501      148820 :    END SUBROUTINE transform_first
    1502             : 
    1503             : ! **************************************************************************************************
    1504             : !> \brief Transform position and velocities for the second half of the
    1505             : !>        Velocity-Verlet integration in the npt_f ensemble
    1506             : !> \param particle_set ...
    1507             : !> \param vel ...
    1508             : !> \param index ...
    1509             : !> \param u ...
    1510             : !> \param dm ...
    1511             : !> \param poly_v ...
    1512             : !> \param scale_v ...
    1513             : !> \par History
    1514             : !>      none
    1515             : !> \author MI (February 2008)
    1516             : ! **************************************************************************************************
    1517      112532 :    SUBROUTINE transform_second(particle_set, vel, index, u, dm, poly_v, scale_v)
    1518             : 
    1519             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1520             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vel
    1521             :       INTEGER, INTENT(IN)                                :: index
    1522             :       REAL(KIND=dp), INTENT(IN)                          :: u(3, 3), dm, poly_v(3), scale_v(3)
    1523             : 
    1524             :       REAL(KIND=dp), DIMENSION(3)                        :: uf, uv
    1525             : 
    1526      112532 :       uv = MATMUL(TRANSPOSE(u), vel(:, index))
    1527      112532 :       uf = MATMUL(TRANSPOSE(u), particle_set(index)%f(:))
    1528             : 
    1529      112532 :       uv(1) = uv(1)*scale_v(1)*scale_v(1) + uf(1)*scale_v(1)*poly_v(1)*dm
    1530      112532 :       uv(2) = uv(2)*scale_v(2)*scale_v(2) + uf(2)*scale_v(2)*poly_v(2)*dm
    1531      112532 :       uv(3) = uv(3)*scale_v(3)*scale_v(3) + uf(3)*scale_v(3)*poly_v(3)*dm
    1532             : 
    1533     1913044 :       vel(:, index) = MATMUL(u, uv)
    1534             : 
    1535      112532 :    END SUBROUTINE transform_second
    1536             : 
    1537             : ! **************************************************************************************************
    1538             : !> \brief  Compute the timestep rescaling factor
    1539             : !> \param md_env ...
    1540             : !> \param tmp ...
    1541             : !> \param dt ...
    1542             : !> \param simpar ...
    1543             : !> \param para_env ...
    1544             : !> \param atomic_kind_set ...
    1545             : !> \param local_particles ...
    1546             : !> \param particle_set ...
    1547             : !> \param core_particle_set ...
    1548             : !> \param shell_particle_set ...
    1549             : !> \param nparticle_kind ...
    1550             : !> \param shell_adiabatic ...
    1551             : !> \param npt ...
    1552             : !> \par History
    1553             : !>      none
    1554             : !> \author MI (October 2008)
    1555             : ! **************************************************************************************************
    1556         480 :    SUBROUTINE variable_timestep(md_env, tmp, dt, simpar, para_env, atomic_kind_set, local_particles, &
    1557             :                                 particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
    1558             : 
    1559             :       TYPE(md_environment_type), POINTER                 :: md_env
    1560             :       TYPE(tmp_variables_type), POINTER                  :: tmp
    1561             :       REAL(KIND=dp), INTENT(INOUT)                       :: dt
    1562             :       TYPE(simpar_type), POINTER                         :: simpar
    1563             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1564             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1565             :       TYPE(distribution_1d_type), POINTER                :: local_particles
    1566             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set, core_particle_set, &
    1567             :                                                             shell_particle_set
    1568             :       INTEGER, INTENT(IN)                                :: nparticle_kind
    1569             :       LOGICAL, INTENT(IN)                                :: shell_adiabatic
    1570             :       TYPE(npt_info_type), OPTIONAL, POINTER             :: npt(:, :)
    1571             : 
    1572             :       INTEGER, SAVE                                      :: itime = 0
    1573             :       REAL(KIND=dp)                                      :: dt_fact, dt_fact_f, dt_fact_old, &
    1574             :                                                             dt_fact_sc, dt_fact_v, dt_old
    1575             :       REAL(KIND=dp), POINTER                             :: time
    1576             :       TYPE(thermostats_type), POINTER                    :: thermostats
    1577             : 
    1578         480 :       dt_fact = 1.0_dp
    1579         480 :       dt_fact_sc = 1.0_dp
    1580         480 :       dt_fact_f = 1.0_dp
    1581         480 :       dt_fact_v = 1.0_dp
    1582         480 :       dt_old = dt
    1583         480 :       dt_fact_old = simpar%dt_fact
    1584         480 :       simpar%dt_fact = 1.0_dp
    1585         480 :       NULLIFY (thermostats)
    1586             : 
    1587         480 :       itime = itime + 1
    1588         480 :       CALL para_env%max(tmp%max_dr)
    1589         480 :       IF (tmp%max_dr > simpar%dr_tol) THEN
    1590         278 :          CALL para_env%max(tmp%max_dvel)
    1591         278 :          dt_fact = SQRT(simpar%dr_tol*dt/tmp%max_dvel)/dt
    1592             :       END IF
    1593             : 
    1594         480 :       IF (shell_adiabatic) THEN
    1595         440 :          CALL para_env%max(tmp%max_dsc)
    1596         440 :          IF (tmp%max_dsc > simpar%dsc_tol) THEN
    1597         118 :             CALL para_env%max(tmp%max_dvel_sc)
    1598         118 :             dt_fact_sc = SQRT(simpar%dsc_tol*dt/tmp%max_dvel_sc)/dt
    1599             :          END IF
    1600             :       END IF
    1601             : 
    1602         480 :       dt_fact_f = MIN(dt_fact, dt_fact_sc, 1.0_dp)
    1603         480 :       IF (dt_fact_f < 1.0_dp) THEN
    1604         100 :          IF (dt_fact_f < 0.1_dp) dt_fact_f = 0.1_dp
    1605             :          ! repeat the first vv half-step with the rescaled time step
    1606         100 :          dt = dt*dt_fact_f
    1607         100 :          simpar%dt_fact = dt_fact_f
    1608             :          CALL rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
    1609         100 :                                 particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
    1610             :       END IF
    1611             : 
    1612         480 :       dt_fact = 1.0_dp
    1613         480 :       dt_fact_sc = 1.0_dp
    1614         480 :       CALL para_env%max(tmp%max_dr)
    1615         480 :       IF (tmp%max_dr > simpar%dr_tol) THEN
    1616         278 :          CALL para_env%max(tmp%max_vel)
    1617         278 :          dt_fact = simpar%dr_tol/tmp%max_vel/dt
    1618             :       END IF
    1619             : 
    1620         480 :       IF (shell_adiabatic) THEN
    1621         440 :          CALL para_env%max(tmp%max_dsc)
    1622         440 :          IF (tmp%max_dsc > simpar%dsc_tol) THEN
    1623         116 :             CALL para_env%max(tmp%max_vel_sc)
    1624         116 :             dt_fact_sc = simpar%dsc_tol/tmp%max_vel_sc/dt
    1625             :          END IF
    1626             :       END IF
    1627         480 :       dt_fact_v = MIN(dt_fact, dt_fact_sc, 1.0_dp)
    1628             : 
    1629         480 :       IF (dt_fact_v < 1.0_dp) THEN
    1630         376 :          IF (dt_fact_v < 0.1_dp) dt_fact_v = 0.1_dp
    1631             :          ! repeat the first vv half-step with the rescaled time step
    1632         376 :          dt = dt*dt_fact_v
    1633         376 :          simpar%dt_fact = dt_fact_f*dt_fact_v
    1634             :          CALL rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
    1635         554 :                                 particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
    1636             :       END IF
    1637             : 
    1638         480 :       simpar%dt_fact = dt_fact_f*dt_fact_v
    1639         480 :       IF (simpar%dt_fact < 1.0_dp) THEN
    1640         378 :          CALL get_md_env(md_env, t=time, thermostats=thermostats)
    1641         378 :          time = time - dt_old + dt_old*simpar%dt_fact
    1642         378 :          IF (ASSOCIATED(thermostats)) THEN
    1643         240 :             CALL set_thermostats(thermostats, dt_fact=simpar%dt_fact)
    1644             :          END IF
    1645             :       END IF
    1646             : 
    1647         480 :    END SUBROUTINE variable_timestep
    1648             : 
    1649             : ! **************************************************************************************************
    1650             : !> \brief Repeat the first step of velocity verlet with the rescaled time step
    1651             : !> \param tmp ...
    1652             : !> \param dt ...
    1653             : !> \param simpar ...
    1654             : !> \param atomic_kind_set ...
    1655             : !> \param local_particles ...
    1656             : !> \param particle_set ...
    1657             : !> \param core_particle_set ...
    1658             : !> \param shell_particle_set ...
    1659             : !> \param nparticle_kind ...
    1660             : !> \param shell_adiabatic ...
    1661             : !> \param npt ...
    1662             : !> \par History
    1663             : !>      none
    1664             : !> \author MI (October 2008)
    1665             : !>     I soliti ignoti
    1666             : ! **************************************************************************************************
    1667         476 :    SUBROUTINE rescaled_vv_first(tmp, dt, simpar, atomic_kind_set, local_particles, &
    1668             :                                 particle_set, core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, npt)
    1669             : 
    1670             :       TYPE(tmp_variables_type), POINTER                  :: tmp
    1671             :       REAL(KIND=dp), INTENT(IN)                          :: dt
    1672             :       TYPE(simpar_type), POINTER                         :: simpar
    1673             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1674             :       TYPE(distribution_1d_type), POINTER                :: local_particles
    1675             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set, core_particle_set, &
    1676             :                                                             shell_particle_set
    1677             :       INTEGER, INTENT(IN)                                :: nparticle_kind
    1678             :       LOGICAL, INTENT(IN)                                :: shell_adiabatic
    1679             :       TYPE(npt_info_type), OPTIONAL, POINTER             :: npt(:, :)
    1680             : 
    1681             :       REAL(KIND=dp), PARAMETER                           :: e2 = 1.0_dp/6.0_dp, e4 = e2/20.0_dp, &
    1682             :                                                             e6 = e4/42.0_dp, e8 = e6/72.0_dp
    1683             : 
    1684             :       REAL(KIND=dp)                                      :: arg_r(3), arg_v(3), e_val(3), infree, &
    1685             :                                                             trvg, u(3, 3)
    1686             : 
    1687         476 :       infree = 1.0_dp/REAL(simpar%nfree, dp)
    1688        1904 :       arg_r = tmp%arg_r
    1689        1904 :       arg_v = tmp%arg_v
    1690        6188 :       u = tmp%u
    1691        1904 :       e_val = tmp%e_val
    1692             : 
    1693         654 :       SELECT CASE (simpar%ensemble)
    1694             :       CASE (nve_ensemble, nvt_ensemble)
    1695             :          CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1696         178 :                        core_particle_set, shell_particle_set, nparticle_kind, shell_adiabatic, dt)
    1697             :       CASE (isokin_ensemble)
    1698           0 :          tmp%scale_v(1:3) = SQRT(1.0_dp/tmp%ds)
    1699           0 :          tmp%poly_v(1:3) = 2.0_dp*tmp%s/SQRT(tmp%ds)/dt
    1700             :          CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1701             :                        core_particle_set, shell_particle_set, nparticle_kind, &
    1702           0 :                        shell_adiabatic, dt)
    1703             : 
    1704             :       CASE (npt_i_ensemble, npt_ia_ensemble, npe_i_ensemble)
    1705           0 :          arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
    1706           0 :          tmp%poly_r(1:3) = 1.0_dp + e2*arg_r(1) + e4*arg_r(1)*arg_r(1) + e6*arg_r(1)**3 + e8*arg_r(1)**4
    1707           0 :          arg_v = arg_v*simpar%dt_fact*simpar%dt_fact
    1708           0 :          tmp%poly_v(1:3) = 1.0_dp + e2*arg_v(1) + e4*arg_v(1)*arg_v(1) + e6*arg_v(1)**3 + e8*arg_v(1)**4
    1709           0 :          tmp%scale_r(1:3) = EXP(0.5_dp*dt*npt(1, 1)%v)
    1710             :          tmp%scale_v(1:3) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
    1711           0 :                                 (1.0_dp + 3.0_dp*infree))
    1712             :          CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1713             :                        core_particle_set, shell_particle_set, nparticle_kind, &
    1714           0 :                        shell_adiabatic, dt)
    1715             : 
    1716             :       CASE (npt_f_ensemble, npe_f_ensemble)
    1717         298 :          trvg = npt(1, 1)%v + npt(2, 2)%v + npt(3, 3)%v
    1718        1192 :          arg_r(:) = arg_r(:)*simpar%dt_fact*simpar%dt_fact
    1719        1192 :          tmp%poly_r = 1._dp + e2*arg_r + e4*arg_r*arg_r + e6*arg_r**3 + e8*arg_r**4
    1720        1192 :          tmp%scale_r(:) = EXP(0.5_dp*dt*e_val(:))
    1721        1192 :          arg_v(:) = arg_v(:)*simpar%dt_fact*simpar%dt_fact
    1722        1192 :          tmp%poly_v = 1.0_dp + e2*arg_v + e4*arg_v*arg_v + e6*arg_v**3 + e8*arg_v**4
    1723             :          tmp%scale_v(:) = EXP(-0.25_dp*dt*( &
    1724        1192 :                               e_val(:) + trvg*infree))
    1725             : 
    1726             :          CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1727             :                        core_particle_set, shell_particle_set, nparticle_kind, &
    1728         298 :                        shell_adiabatic, dt, u)
    1729             : 
    1730             :       CASE (nph_uniaxial_ensemble, nph_uniaxial_damped_ensemble)
    1731           0 :          arg_r = arg_r*simpar%dt_fact*simpar%dt_fact
    1732           0 :          tmp%poly_r(1) = 1._dp + e2*arg_r(1) + e4*arg_r(1)*arg_r(1) + e6*arg_r(1)**3 + e8*arg_r(1)**4
    1733           0 :          arg_v(2) = arg_v(2)*simpar%dt_fact*simpar%dt_fact
    1734           0 :          arg_v(1) = arg_v(1)*simpar%dt_fact*simpar%dt_fact
    1735           0 :          tmp%poly_v(1) = 1._dp + e2*arg_v(1) + e4*arg_v(1)*arg_v(1) + e6*arg_v(1)**3 + e8*arg_v(1)**4
    1736           0 :          tmp%poly_v(2) = 1._dp + e2*arg_v(2) + e4*arg_v(2)*arg_v(2) + e6*arg_v(2)**3 + e8*arg_v(2)**4
    1737           0 :          tmp%poly_v(3) = 1._dp + e2*arg_v(2) + e4*arg_v(2)*arg_v(2) + e6*arg_v(2)**3 + e8*arg_v(2)**4
    1738           0 :          tmp%scale_r(1) = EXP(0.5_dp*dt*npt(1, 1)%v)
    1739             :          tmp%scale_v(1) = EXP(-0.25_dp*dt*npt(1, 1)%v* &
    1740           0 :                               (1._dp + infree))
    1741           0 :          tmp%scale_v(2) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree)
    1742           0 :          tmp%scale_v(3) = EXP(-0.25_dp*dt*npt(1, 1)%v*infree)
    1743             :          CALL vv_first(tmp, atomic_kind_set, local_particles, particle_set, &
    1744             :                        core_particle_set, shell_particle_set, nparticle_kind, &
    1745         476 :                        shell_adiabatic, dt)
    1746             : 
    1747             :       END SELECT
    1748             : 
    1749         476 :    END SUBROUTINE rescaled_vv_first
    1750             : 
    1751           0 : END MODULE integrator_utils

Generated by: LCOV version 1.15