LCOV - code coverage report
Current view: top level - src/motion - md_vel_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 843 935 90.2 %
Date: 2024-04-20 06:29:22 Functions: 28 28 100.0 %

          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  Collection of utilities for setting-up and handle velocities in MD
      10             : !>         runs
      11             : !> \author CJM
      12             : !> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008
      13             : !>         reorganization of the original routines/modules
      14             : !>      Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
      15             : !>                                       (patch by Marcel Baer)
      16             : ! **************************************************************************************************
      17             : MODULE md_vel_utils
      18             :    USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
      19             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      20             :                                               get_atomic_kind,&
      21             :                                               get_atomic_kind_set
      22             :    USE bibliography,                    ONLY: West2006,&
      23             :                                               cite_reference
      24             :    USE cell_types,                      ONLY: &
      25             :         cell_type, use_perd_none, use_perd_x, use_perd_xy, use_perd_xyz, use_perd_xz, use_perd_y, &
      26             :         use_perd_yz, use_perd_z
      27             :    USE cp_linked_list_input,            ONLY: cp_sll_val_next,&
      28             :                                               cp_sll_val_type
      29             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      30             :                                               cp_logger_type
      31             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      32             :                                               cp_print_key_unit_nr
      33             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      34             :                                               cp_subsys_type
      35             :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      36             :    USE extended_system_types,           ONLY: npt_info_type
      37             :    USE force_env_methods,               ONLY: force_env_calc_energy_force
      38             :    USE force_env_types,                 ONLY: force_env_get,&
      39             :                                               force_env_type
      40             :    USE force_env_utils,                 ONLY: force_env_rattle,&
      41             :                                               force_env_shake
      42             :    USE global_types,                    ONLY: global_environment_type
      43             :    USE input_constants,                 ONLY: &
      44             :         md_init_default, md_init_vib, npe_f_ensemble, npe_i_ensemble, &
      45             :         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
      46             :         npt_ia_ensemble, reftraj_ensemble
      47             :    USE input_cp2k_binary_restarts,      ONLY: read_binary_velocities
      48             :    USE input_restart_force_eval,        ONLY: update_subsys
      49             :    USE input_section_types,             ONLY: section_vals_get,&
      50             :                                               section_vals_get_subs_vals,&
      51             :                                               section_vals_list_get,&
      52             :                                               section_vals_type,&
      53             :                                               section_vals_val_get
      54             :    USE input_val_types,                 ONLY: val_get,&
      55             :                                               val_type
      56             :    USE kinds,                           ONLY: default_string_length,&
      57             :                                               dp
      58             :    USE mathconstants,                   ONLY: pi
      59             :    USE mathlib,                         ONLY: diamat_all
      60             :    USE md_ener_types,                   ONLY: md_ener_type
      61             :    USE md_environment_types,            ONLY: get_md_env,&
      62             :                                               md_environment_type
      63             :    USE md_util,                         ONLY: read_vib_eigs_unformatted
      64             :    USE message_passing,                 ONLY: mp_para_env_type
      65             :    USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
      66             :    USE molecule_kind_types,             ONLY: fixd_constraint_type,&
      67             :                                               get_molecule_kind,&
      68             :                                               get_molecule_kind_set,&
      69             :                                               molecule_kind_type
      70             :    USE parallel_rng_types,              ONLY: UNIFORM,&
      71             :                                               rng_stream_type
      72             :    USE particle_list_types,             ONLY: particle_list_type
      73             :    USE particle_types,                  ONLY: particle_type
      74             :    USE physcon,                         ONLY: kelvin
      75             :    USE shell_opt,                       ONLY: optimize_shell_core
      76             :    USE shell_potential_types,           ONLY: shell_kind_type
      77             :    USE simpar_types,                    ONLY: simpar_type
      78             :    USE thermal_region_types,            ONLY: thermal_region_type,&
      79             :                                               thermal_regions_type
      80             : #include "../base/base_uses.f90"
      81             : 
      82             :    IMPLICIT NONE
      83             : 
      84             :    PRIVATE
      85             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_vel_utils'
      86             : 
      87             :    PUBLIC :: temperature_control, &
      88             :              comvel_control, &
      89             :              angvel_control, &
      90             :              setup_velocities
      91             : 
      92             : CONTAINS
      93             : 
      94             : ! **************************************************************************************************
      95             : !> \brief compute center of mass position
      96             : !>      *** is only used by initialize_velocities below ***
      97             : !> \param part ...
      98             : !> \param is_fixed ...
      99             : !> \param rcom ...
     100             : !> \par History
     101             : !>      2007-11-6: created
     102             : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
     103             : ! **************************************************************************************************
     104         104 :    SUBROUTINE compute_rcom(part, is_fixed, rcom)
     105             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     106             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: is_fixed
     107             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: rcom
     108             : 
     109             :       INTEGER                                            :: i
     110             :       REAL(KIND=dp)                                      :: denom, mass
     111             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     112             : 
     113         104 :       rcom(:) = 0.0_dp
     114         104 :       denom = 0.0_dp
     115        1238 :       DO i = 1, SIZE(part)
     116        1134 :          atomic_kind => part(i)%atomic_kind
     117        1134 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     118        1238 :          SELECT CASE (is_fixed(i))
     119             :          CASE (use_perd_x, use_perd_y, use_perd_z, use_perd_xy, use_perd_xz, use_perd_yz, use_perd_none)
     120        1134 :             rcom(1) = rcom(1) + part(i)%r(1)*mass
     121        1134 :             rcom(2) = rcom(2) + part(i)%r(2)*mass
     122        1134 :             rcom(3) = rcom(3) + part(i)%r(3)*mass
     123        1134 :             denom = denom + mass
     124             :          END SELECT
     125             :       END DO
     126         416 :       rcom = rcom/denom
     127             : 
     128         104 :    END SUBROUTINE compute_rcom
     129             : 
     130             : ! **************************************************************************************************
     131             : !> \brief compute center of mass velocity
     132             : !>      *** is only used by initialize_velocities below ***
     133             : !> \param part ...
     134             : !> \param is_fixed ...
     135             : !> \param vcom ...
     136             : !> \param ecom ...
     137             : !> \par History
     138             : !>      2007-11-6: created
     139             : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
     140             : ! **************************************************************************************************
     141        3839 :    SUBROUTINE compute_vcom(part, is_fixed, vcom, ecom)
     142             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     143             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: is_fixed
     144             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: vcom
     145             :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: ecom
     146             : 
     147             :       INTEGER                                            :: i
     148             :       REAL(KIND=dp)                                      :: denom, mass
     149             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     150             : 
     151        3839 :       vcom = 0.0_dp
     152        3839 :       denom = 0.0_dp
     153      747856 :       DO i = 1, SIZE(part)
     154      744017 :          atomic_kind => part(i)%atomic_kind
     155      744017 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     156      747856 :          IF (mass .NE. 0.0) THEN
     157     1480666 :             SELECT CASE (is_fixed(i))
     158             :             CASE (use_perd_x, use_perd_y, use_perd_z, use_perd_xy, use_perd_xz, use_perd_yz, use_perd_none)
     159      737975 :                vcom(1) = vcom(1) + part(i)%v(1)*mass
     160      737975 :                vcom(2) = vcom(2) + part(i)%v(2)*mass
     161      737975 :                vcom(3) = vcom(3) + part(i)%v(3)*mass
     162      742691 :                denom = denom + mass
     163             :             END SELECT
     164             :          END IF
     165             :       END DO
     166       15356 :       vcom = vcom/denom
     167        3839 :       IF (PRESENT(ecom)) THEN
     168        4356 :          ecom = 0.5_dp*denom*SUM(vcom*vcom)
     169             :       END IF
     170             : 
     171        3839 :    END SUBROUTINE compute_vcom
     172             : 
     173             : ! **************************************************************************************************
     174             : !> \brief Copy atom velocities into core and shell velocities
     175             : !>      *** is only used by initialize_velocities below ***
     176             : !> \param part ...
     177             : !> \param shell_part ...
     178             : !> \param core_part ...
     179             : !> \par History
     180             : !>      2007-11-6: created
     181             : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
     182             : ! **************************************************************************************************
     183           8 :    SUBROUTINE clone_core_shell_vel(part, shell_part, core_part)
     184             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part, shell_part, core_part
     185             : 
     186             :       INTEGER                                            :: i
     187             :       LOGICAL                                            :: is_shell
     188             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     189             : 
     190         776 :       DO i = 1, SIZE(part)
     191         768 :          atomic_kind => part(i)%atomic_kind
     192         768 :          CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_shell)
     193         776 :          IF (is_shell) THEN
     194        6144 :             shell_part(part(i)%shell_index)%v(:) = part(i)%v(:)
     195        6144 :             core_part(part(i)%shell_index)%v(:) = part(i)%v(:)
     196             :          END IF
     197             :       END DO
     198             : 
     199           8 :    END SUBROUTINE clone_core_shell_vel
     200             : 
     201             : ! **************************************************************************************************
     202             : !> \brief Compute the kinetic energy. Does not subtract the center of mass kinetic
     203             : !>      energy.
     204             : !>      *** is only used by initialize_velocities below ***
     205             : !> \param part ...
     206             : !> \param ireg ...
     207             : !> \return ...
     208             : !> \par History
     209             : !>      2007-11-6: created
     210             : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
     211             : ! **************************************************************************************************
     212        4105 :    FUNCTION compute_ekin(part, ireg) RESULT(ekin)
     213             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     214             :       INTEGER, INTENT(IN), OPTIONAL                      :: ireg
     215             :       REAL(KIND=dp)                                      :: ekin
     216             : 
     217             :       INTEGER                                            :: i
     218             :       REAL(KIND=dp)                                      :: mass
     219             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     220             : 
     221        4105 :       NULLIFY (atomic_kind)
     222        4105 :       ekin = 0.0_dp
     223        4105 :       IF (PRESENT(ireg)) THEN
     224       13756 :          DO i = 1, SIZE(part)
     225       13756 :             IF (part(i)%t_region_index == ireg) THEN
     226        4560 :                atomic_kind => part(i)%atomic_kind
     227        4560 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     228       18240 :                ekin = ekin + 0.5_dp*mass*SUM(part(i)%v(:)*part(i)%v(:))
     229             :             END IF
     230             :          END DO
     231             :       ELSE
     232      747662 :          DO i = 1, SIZE(part)
     233      743825 :             atomic_kind => part(i)%atomic_kind
     234      743825 :             CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     235     2979137 :             ekin = ekin + 0.5_dp*mass*SUM(part(i)%v(:)*part(i)%v(:))
     236             :          END DO
     237             :       END IF
     238             : 
     239        4105 :    END FUNCTION compute_ekin
     240             : 
     241             : ! **************************************************************************************************
     242             : !> \brief Rescale the velocity to mimic the given external kinetic temperature.
     243             : !>      Optionally also rescale vcom.
     244             : !>      *** is only used by initialize_velocities below ***
     245             : !> \param part ...
     246             : !> \param simpar ...
     247             : !> \param ekin ...
     248             : !> \param vcom ...
     249             : !> \param ireg ...
     250             : !> \param nfree ...
     251             : !> \param temp ...
     252             : !> \par History
     253             : !>      2007-11-6: created
     254             : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
     255             : ! **************************************************************************************************
     256        2778 :    SUBROUTINE rescale_vel(part, simpar, ekin, vcom, ireg, nfree, temp)
     257             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     258             :       TYPE(simpar_type), POINTER                         :: simpar
     259             :       REAL(KIND=dp), INTENT(INOUT)                       :: ekin
     260             :       REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
     261             :          OPTIONAL                                        :: vcom
     262             :       INTEGER, INTENT(IN), OPTIONAL                      :: ireg, nfree
     263             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: temp
     264             : 
     265             :       INTEGER                                            :: i, my_ireg, my_nfree
     266             :       REAL(KIND=dp)                                      :: factor, my_temp
     267             : 
     268        2778 :       IF (PRESENT(ireg) .AND. PRESENT(nfree) .AND. PRESENT(temp)) THEN
     269           6 :          my_ireg = ireg
     270           6 :          my_nfree = nfree
     271           6 :          my_temp = temp
     272        2772 :       ELSEIF (PRESENT(nfree)) THEN
     273           0 :          my_ireg = 0
     274           0 :          my_nfree = nfree
     275           0 :          my_temp = simpar%temp_ext
     276             :       ELSE
     277        2772 :          my_ireg = 0
     278        2772 :          my_nfree = simpar%nfree
     279        2772 :          my_temp = simpar%temp_ext
     280             :       END IF
     281        2778 :       IF (my_nfree /= 0) THEN
     282        2766 :          factor = my_temp/(2.0_dp*ekin)*REAL(my_nfree, KIND=dp)
     283             :       ELSE
     284             :          factor = 0.0_dp
     285             :       END IF
     286             :       ! Note:
     287             :       ! this rescaling is still wrong, it should take the masses into account
     288             :       ! rescaling is generally not correct, so needs fixing
     289        2778 :       ekin = ekin*factor
     290        2778 :       factor = SQRT(factor)
     291        2778 :       IF (PRESENT(ireg)) THEN
     292         582 :          DO i = 1, SIZE(part)
     293        1158 :             IF (part(i)%t_region_index == my_ireg) part(i)%v(:) = factor*part(i)%v(:)
     294             :          END DO
     295             :       ELSE
     296      447570 :          DO i = 1, SIZE(part)
     297     1781964 :             part(i)%v(:) = factor*part(i)%v(:)
     298             :          END DO
     299        2772 :          IF (PRESENT(vcom)) THEN
     300          96 :             vcom = factor*vcom
     301             :          END IF
     302             :       END IF
     303             : 
     304        2778 :    END SUBROUTINE rescale_vel
     305             : 
     306             : ! **************************************************************************************************
     307             : !> \brief Rescale the velocity of separated regions independently
     308             : !> \param part ...
     309             : !> \param md_env ...
     310             : !> \param simpar ...
     311             : !> \par History
     312             : !>      2008-11
     313             : !> \author  MI
     314             : ! **************************************************************************************************
     315           2 :    SUBROUTINE rescale_vel_region(part, md_env, simpar)
     316             : 
     317             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     318             :       TYPE(md_environment_type), POINTER                 :: md_env
     319             :       TYPE(simpar_type), POINTER                         :: simpar
     320             : 
     321             :       INTEGER                                            :: ireg, nfree, nfree0, nfree_done
     322             :       REAL(KIND=dp)                                      :: ekin, temp
     323             :       TYPE(thermal_region_type), POINTER                 :: t_region
     324             :       TYPE(thermal_regions_type), POINTER                :: thermal_regions
     325             : 
     326           2 :       NULLIFY (thermal_regions, t_region)
     327             : 
     328           2 :       CALL get_md_env(md_env, thermal_regions=thermal_regions)
     329           2 :       nfree_done = 0
     330           6 :       DO ireg = 1, thermal_regions%nregions
     331           4 :          NULLIFY (t_region)
     332           4 :          t_region => thermal_regions%thermal_region(ireg)
     333           4 :          nfree = t_region%npart*3
     334           4 :          ekin = compute_ekin(part, ireg)
     335           4 :          temp = t_region%temp_expected
     336           4 :          CALL rescale_vel(part, simpar, ekin, ireg=ireg, nfree=nfree, temp=temp)
     337           4 :          nfree_done = nfree_done + nfree
     338           4 :          ekin = compute_ekin(part, ireg)
     339           4 :          temp = 2.0_dp*ekin/REAL(nfree, dp)*kelvin
     340           6 :          t_region%temperature = temp
     341             :       END DO
     342           2 :       nfree0 = simpar%nfree - nfree_done
     343           2 :       IF (nfree0 > 0) THEN
     344           2 :          ekin = compute_ekin(part, 0)
     345           2 :          CALL rescale_vel(part, simpar, ekin, ireg=0, nfree=nfree0, temp=simpar%temp_ext)
     346           2 :          ekin = compute_ekin(part, 0)
     347           2 :          temp = 2.0_dp*ekin/REAL(nfree0, dp)*kelvin
     348           2 :          thermal_regions%temp_reg0 = temp
     349             :       END IF
     350           2 :    END SUBROUTINE rescale_vel_region
     351             : 
     352             : ! **************************************************************************************************
     353             : !> \brief subtract center of mass velocity
     354             : !>      *** is only used by initialize_velocities below ***
     355             : !> \param part ...
     356             : !> \param is_fixed ...
     357             : !> \param vcom ...
     358             : !> \par History
     359             : !>      2007-11-6: created
     360             : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
     361             : ! **************************************************************************************************
     362        2750 :    SUBROUTINE subtract_vcom(part, is_fixed, vcom)
     363             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     364             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: is_fixed
     365             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: vcom
     366             : 
     367             :       INTEGER                                            :: i
     368             : 
     369      445870 :       DO i = 1, SIZE(part)
     370        2750 :          SELECT CASE (is_fixed(i))
     371             :          CASE (use_perd_x)
     372         512 :             part(i)%v(2) = part(i)%v(2) - vcom(2)
     373         512 :             part(i)%v(3) = part(i)%v(3) - vcom(3)
     374             :          CASE (use_perd_y)
     375         512 :             part(i)%v(1) = part(i)%v(1) - vcom(1)
     376         512 :             part(i)%v(3) = part(i)%v(3) - vcom(3)
     377             :          CASE (use_perd_z)
     378         512 :             part(i)%v(1) = part(i)%v(1) - vcom(1)
     379         512 :             part(i)%v(2) = part(i)%v(2) - vcom(2)
     380             :          CASE (use_perd_xy)
     381         512 :             part(i)%v(3) = part(i)%v(3) - vcom(3)
     382             :          CASE (use_perd_xz)
     383           0 :             part(i)%v(2) = part(i)%v(2) - vcom(2)
     384             :          CASE (use_perd_yz)
     385           0 :             part(i)%v(1) = part(i)%v(1) - vcom(1)
     386             :          CASE (use_perd_none)
     387     1756982 :             part(i)%v(:) = part(i)%v(:) - vcom(:)
     388             :          END SELECT
     389             :       END DO
     390        2750 :    END SUBROUTINE subtract_vcom
     391             : 
     392             : ! **************************************************************************************************
     393             : !> \brief compute the angular velocity
     394             : !>      *** is only used by initialize_velocities below ***
     395             : !> \param part ...
     396             : !> \param is_fixed ...
     397             : !> \param rcom ...
     398             : !> \param vang ...
     399             : !> \par History
     400             : !>      2007-11-9: created
     401             : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
     402             : ! **************************************************************************************************
     403         106 :    SUBROUTINE compute_vang(part, is_fixed, rcom, vang)
     404             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     405             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: is_fixed
     406             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rcom
     407             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: vang
     408             : 
     409             :       INTEGER                                            :: i
     410             :       REAL(KIND=dp)                                      :: mass, proj
     411             :       REAL(KIND=dp), DIMENSION(3)                        :: evals, mang, r
     412             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: iner
     413             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     414             : 
     415         106 :       NULLIFY (atomic_kind)
     416         106 :       mang(:) = 0.0_dp
     417         106 :       iner(:, :) = 0.0_dp
     418        1268 :       DO i = 1, SIZE(part)
     419             :          ! compute angular momentum and inertia tensor
     420         106 :          SELECT CASE (is_fixed(i))
     421             :          CASE (use_perd_x, use_perd_y, use_perd_z, use_perd_xy, use_perd_xz, use_perd_yz, use_perd_none)
     422        4648 :             r(:) = part(i)%r(:) - rcom(:)
     423        1162 :             atomic_kind => part(i)%atomic_kind
     424        1162 :             CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     425        1162 :             mang(1) = mang(1) + mass*(r(2)*part(i)%v(3) - r(3)*part(i)%v(2))
     426        1162 :             mang(2) = mang(2) + mass*(r(3)*part(i)%v(1) - r(1)*part(i)%v(3))
     427        1162 :             mang(3) = mang(3) + mass*(r(1)*part(i)%v(2) - r(2)*part(i)%v(1))
     428             : 
     429        1162 :             iner(1, 1) = iner(1, 1) + mass*(r(2)*r(2) + r(3)*r(3))
     430        1162 :             iner(2, 2) = iner(2, 2) + mass*(r(3)*r(3) + r(1)*r(1))
     431        1162 :             iner(3, 3) = iner(3, 3) + mass*(r(1)*r(1) + r(2)*r(2))
     432             : 
     433        1162 :             iner(1, 2) = iner(1, 2) - mass*r(1)*r(2)
     434        1162 :             iner(2, 3) = iner(2, 3) - mass*r(2)*r(3)
     435        2324 :             iner(3, 1) = iner(3, 1) - mass*r(3)*r(1)
     436             :          END SELECT
     437             :       END DO
     438         106 :       iner(2, 1) = iner(1, 2)
     439         106 :       iner(3, 2) = iner(2, 3)
     440         106 :       iner(1, 3) = iner(3, 1)
     441             : 
     442             :       ! Take the safest route, i.e. diagonalize the inertia tensor and solve
     443             :       ! the angular velocity only with the non-zero eigenvalues. A plain inversion
     444             :       ! would fail for linear molecules.
     445         106 :       CALL diamat_all(iner, evals)
     446             : 
     447         106 :       vang(:) = 0.0_dp
     448         424 :       DO i = 1, 3
     449         424 :          IF (evals(i) > 0.0_dp) THEN
     450        1228 :             proj = SUM(iner(:, i)*mang)/evals(i)
     451         307 :             vang(1) = vang(1) + proj*iner(1, i)
     452         307 :             vang(2) = vang(2) + proj*iner(2, i)
     453         307 :             vang(3) = vang(3) + proj*iner(3, i)
     454             :          END IF
     455             :       END DO
     456             : 
     457         106 :    END SUBROUTINE compute_vang
     458             : 
     459             : ! **************************************************************************************************
     460             : !> \brief subtract the angular velocity
     461             : !>      *** is only used by initialize_velocities below ***
     462             : !> \param part ...
     463             : !> \param is_fixed ...
     464             : !> \param rcom ...
     465             : !> \param vang ...
     466             : !> \par History
     467             : !>      2007-11-9: created
     468             : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
     469             : ! **************************************************************************************************
     470           6 :    SUBROUTINE subtract_vang(part, is_fixed, rcom, vang)
     471             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     472             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: is_fixed
     473             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rcom, vang
     474             : 
     475             :       INTEGER                                            :: i
     476             :       REAL(KIND=dp), DIMENSION(3)                        :: r
     477             : 
     478          52 :       DO i = 1, SIZE(part)
     479         184 :          r(:) = part(i)%r(:) - rcom(:)
     480           6 :          SELECT CASE (is_fixed(i))
     481             :          CASE (use_perd_x)
     482           0 :             part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
     483           0 :             part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
     484             :          CASE (use_perd_y)
     485           0 :             part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
     486           0 :             part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
     487             :          CASE (use_perd_z)
     488           0 :             part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
     489           0 :             part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
     490             :          CASE (use_perd_xy)
     491           0 :             part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
     492             :          CASE (use_perd_xz)
     493           0 :             part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
     494             :          CASE (use_perd_yz)
     495           0 :             part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
     496             :          CASE (use_perd_none)
     497          46 :             part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
     498          46 :             part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
     499          46 :             part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
     500             :          END SELECT
     501             :       END DO
     502             : 
     503           6 :    END SUBROUTINE subtract_vang
     504             : 
     505             : ! **************************************************************************************************
     506             : !> \brief Initializes the velocities to the Maxwellian distribution
     507             : !> \param simpar ...
     508             : !> \param part ...
     509             : !> \param force_env ...
     510             : !> \param globenv ...
     511             : !> \param md_env ...
     512             : !> \param molecule_kinds ...
     513             : !> \param label ...
     514             : !> \param print_section ...
     515             : !> \param subsys_section ...
     516             : !> \param shell_present ...
     517             : !> \param shell_part ...
     518             : !> \param core_part ...
     519             : !> \param force_rescaling ...
     520             : !> \param para_env ...
     521             : !> \param write_binary_restart_file ...
     522             : !> \par History
     523             : !>      - is_fixed removed from particle_type
     524             : !>      - 2007-11-07: Cleanup (TV)
     525             : !>      - 2007-11-09: Added angvel_zero feature
     526             : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>
     527             : ! **************************************************************************************************
     528        1818 :    SUBROUTINE initialize_velocities(simpar, &
     529             :                                     part, &
     530             :                                     force_env, &
     531             :                                     globenv, &
     532             :                                     md_env, &
     533             :                                     molecule_kinds, &
     534             :                                     label, &
     535             :                                     print_section, &
     536             :                                     subsys_section, &
     537             :                                     shell_present, &
     538             :                                     shell_part, &
     539             :                                     core_part, &
     540             :                                     force_rescaling, &
     541             :                                     para_env, &
     542             :                                     write_binary_restart_file)
     543             : 
     544             :       TYPE(simpar_type), POINTER                         :: simpar
     545             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     546             :       TYPE(force_env_type), POINTER                      :: force_env
     547             :       TYPE(global_environment_type), POINTER             :: globenv
     548             :       TYPE(md_environment_type), POINTER                 :: md_env
     549             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
     550             :       CHARACTER(LEN=*), INTENT(IN)                       :: label
     551             :       TYPE(section_vals_type), POINTER                   :: print_section, subsys_section
     552             :       LOGICAL, INTENT(IN)                                :: shell_present
     553             :       TYPE(particle_type), DIMENSION(:), POINTER         :: shell_part, core_part
     554             :       LOGICAL, INTENT(IN)                                :: force_rescaling
     555             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     556             :       LOGICAL, INTENT(IN)                                :: write_binary_restart_file
     557             : 
     558             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_velocities'
     559             : 
     560             :       INTEGER                                            :: handle, i, ifixd, imolecule_kind, iw, &
     561             :                                                             natoms
     562             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: is_fixed
     563             :       LOGICAL                                            :: success
     564             :       REAL(KIND=dp)                                      :: ecom, ekin, mass, mass_tot, temp, tmp_r1
     565             :       REAL(KIND=dp), DIMENSION(3)                        :: rcom, vang, vcom
     566             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     567             :       TYPE(cell_type), POINTER                           :: cell
     568             :       TYPE(cp_logger_type), POINTER                      :: logger
     569        1818 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     570        1818 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     571             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     572             :       TYPE(section_vals_type), POINTER                   :: md_section, root_section, vib_section
     573             : 
     574        1818 :       CALL timeset(routineN, handle)
     575             : 
     576             :       ! Initializing parameters
     577        1818 :       natoms = SIZE(part)
     578        1818 :       NULLIFY (atomic_kind, fixd_list, logger, molecule_kind)
     579        1818 :       NULLIFY (molecule_kind_set)
     580             : 
     581             :       ! Logging
     582        1818 :       logger => cp_get_default_logger()
     583        1818 :       iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
     584             : 
     585             :       ! Build a list of all fixed atoms (if any)
     586        5454 :       ALLOCATE (is_fixed(natoms))
     587             : 
     588      490090 :       is_fixed = use_perd_none
     589        1818 :       molecule_kind_set => molecule_kinds%els
     590       62516 :       DO imolecule_kind = 1, molecule_kinds%n_els
     591       60698 :          molecule_kind => molecule_kind_set(imolecule_kind)
     592       60698 :          CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
     593       62516 :          IF (ASSOCIATED(fixd_list)) THEN
     594        5432 :             DO ifixd = 1, SIZE(fixd_list)
     595        5432 :                IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
     596             :             END DO
     597             :          END IF
     598             :       END DO
     599             : 
     600             :       ! Compute the total mass when needed
     601        1818 :       IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
     602             :           simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
     603             :          mass_tot = 0.0_dp
     604        1006 :          DO i = 1, natoms
     605        1000 :             atomic_kind => part(i)%atomic_kind
     606        1000 :             CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     607        1006 :             mass_tot = mass_tot + mass
     608             :          END DO
     609           6 :          simpar%v_shock = simpar%v_shock*SQRT(mass_tot)
     610             :       END IF
     611             : 
     612             :       CALL read_input_velocities(simpar, part, force_env, md_env, subsys_section, &
     613        1818 :                                  shell_present, shell_part, core_part, force_rescaling, para_env, is_fixed, success)
     614        1818 :       IF (.NOT. success) THEN
     615        3102 :          SELECT CASE (simpar%initialization_method)
     616             :          CASE (md_init_default)
     617             :             CALL generate_velocities(simpar, part, force_env, globenv, md_env, shell_present, &
     618        1550 :                                      shell_part, core_part, is_fixed, iw)
     619             :          CASE (md_init_vib)
     620           2 :             CALL force_env_get(force_env=force_env, root_section=root_section)
     621           2 :             md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
     622           2 :             vib_section => section_vals_get_subs_vals(root_section, "VIBRATIONAL_ANALYSIS")
     623             :             CALL generate_coords_vels_vib(simpar, &
     624             :                                           part, &
     625             :                                           md_section, &
     626             :                                           vib_section, &
     627             :                                           force_env, &
     628             :                                           globenv, &
     629             :                                           shell_present, &
     630             :                                           shell_part, &
     631             :                                           core_part, &
     632           2 :                                           is_fixed)
     633             :             ! update restart file for the modified coordinates and velocities
     634        1554 :             CALL update_subsys(subsys_section, force_env, .FALSE., write_binary_restart_file)
     635             :          END SELECT
     636             :       END IF
     637             : 
     638        1818 :       IF (iw > 0) THEN
     639             :          WRITE (iw, '(/,T2,A)') &
     640         822 :             'MD_VEL| '//TRIM(ADJUSTL(label))
     641             :          ! Recompute vcom, ecom and ekin for IO
     642         822 :          CALL compute_vcom(part, is_fixed, vcom, ecom)
     643         822 :          ekin = compute_ekin(part) - ecom
     644         822 :          IF (simpar%nfree == 0) THEN
     645           6 :             CPASSERT(ekin == 0.0_dp)
     646           6 :             temp = 0.0_dp
     647             :          ELSE
     648         816 :             temp = 2.0_dp*ekin/REAL(simpar%nfree, KIND=dp)
     649             :          END IF
     650         822 :          tmp_r1 = cp_unit_from_cp2k(temp, "K")
     651             :          WRITE (iw, '(T2,A,T61,F20.6)') &
     652         822 :             'MD_VEL| Initial temperature [K]', tmp_r1
     653             :          WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
     654         822 :             'MD_VEL| COM velocity', vcom(1:3)
     655             : 
     656             :          ! compute and log rcom and vang if not periodic
     657         822 :          CALL force_env_get(force_env, cell=cell)
     658        3288 :          IF (SUM(cell%perd(1:3)) == 0) THEN
     659          60 :             CALL compute_rcom(part, is_fixed, rcom)
     660          60 :             CALL compute_vang(part, is_fixed, rcom, vang)
     661             :             WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
     662          60 :                'MD_VEL| COM position', rcom(1:3)
     663             :             WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
     664          60 :                'MD_VEL| Angular velocity', vang(1:3)
     665             :          END IF
     666             :       END IF
     667             : 
     668        1818 :       DEALLOCATE (is_fixed)
     669        1818 :       CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
     670        1818 :       CALL timestop(handle)
     671             : 
     672        3636 :    END SUBROUTINE initialize_velocities
     673             : 
     674             : ! **************************************************************************************************
     675             : !> \brief Read velocities from binary restart file if available
     676             : !> \param simpar ...
     677             : !> \param part ...
     678             : !> \param force_env ...
     679             : !> \param md_env ...
     680             : !> \param subsys_section ...
     681             : !> \param shell_present ...
     682             : !> \param shell_part ...
     683             : !> \param core_part ...
     684             : !> \param force_rescaling ...
     685             : !> \param para_env ...
     686             : !> \param is_fixed ...
     687             : !> \param success ...
     688             : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>
     689             : ! **************************************************************************************************
     690       12726 :    SUBROUTINE read_input_velocities(simpar, part, force_env, md_env, subsys_section, &
     691        1818 :                                     shell_present, shell_part, core_part, force_rescaling, para_env, is_fixed, success)
     692             :       TYPE(simpar_type), POINTER                         :: simpar
     693             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
     694             :       TYPE(force_env_type), POINTER                      :: force_env
     695             :       TYPE(md_environment_type), POINTER                 :: md_env
     696             :       TYPE(section_vals_type), POINTER                   :: subsys_section
     697             :       LOGICAL, INTENT(IN)                                :: shell_present
     698             :       TYPE(particle_type), DIMENSION(:), POINTER         :: shell_part, core_part
     699             :       LOGICAL, INTENT(IN)                                :: force_rescaling
     700             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     701             :       INTEGER, DIMENSION(:), INTENT(INOUT)               :: is_fixed
     702             :       LOGICAL, INTENT(OUT)                               :: success
     703             : 
     704             :       INTEGER                                            :: i, natoms, nshell, shell_index
     705             :       LOGICAL :: atomvel_explicit, atomvel_read, corevel_explicit, corevel_read, is_ok, &
     706             :          rescale_regions, shellvel_explicit, shellvel_read
     707             :       REAL(KIND=dp)                                      :: ecom, ekin, fac_massc, fac_masss, mass
     708             :       REAL(KIND=dp), DIMENSION(3)                        :: vc, vcom, vs
     709        1818 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: vel
     710             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     711             :       TYPE(cp_sll_val_type), POINTER                     :: atom_list, core_list, shell_list
     712             :       TYPE(section_vals_type), POINTER                   :: atomvel_section, corevel_section, &
     713             :                                                             shellvel_section
     714             :       TYPE(shell_kind_type), POINTER                     :: shell
     715             :       TYPE(thermal_regions_type), POINTER                :: thermal_regions
     716             :       TYPE(val_type), POINTER                            :: val
     717             : 
     718             : ! Initializing parameters
     719             : 
     720        1818 :       success = .FALSE.
     721        1818 :       natoms = SIZE(part)
     722        1818 :       atomvel_read = .FALSE.
     723        1818 :       corevel_read = .FALSE.
     724        1818 :       shellvel_read = .FALSE.
     725        1818 :       NULLIFY (vel, atomic_kind, atom_list, core_list, shell_list)
     726        1818 :       NULLIFY (atomvel_section, shellvel_section, corevel_section)
     727        1818 :       NULLIFY (shell, thermal_regions, val)
     728             : 
     729             :       ! Core-Shell Model
     730        1818 :       nshell = 0
     731        1818 :       IF (shell_present) THEN
     732         132 :          CPASSERT(ASSOCIATED(core_part))
     733         132 :          CPASSERT(ASSOCIATED(shell_part))
     734         132 :          nshell = SIZE(shell_part)
     735             :       END IF
     736             : 
     737        1818 :       atomvel_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
     738        1818 :       shellvel_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
     739        1818 :       corevel_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
     740             : 
     741             :       ! Read or initialize the particle velocities
     742        1818 :       CALL section_vals_get(atomvel_section, explicit=atomvel_explicit)
     743        1818 :       CALL section_vals_get(shellvel_section, explicit=shellvel_explicit)
     744        1818 :       CALL section_vals_get(corevel_section, explicit=corevel_explicit)
     745        1818 :       CPASSERT(shellvel_explicit .EQV. corevel_explicit)
     746             : 
     747             :       CALL read_binary_velocities("", part, force_env%root_section, para_env, &
     748        1818 :                                   subsys_section, atomvel_read)
     749             :       CALL read_binary_velocities("SHELL", shell_part, force_env%root_section, para_env, &
     750        1818 :                                   subsys_section, shellvel_read)
     751             :       CALL read_binary_velocities("CORE", core_part, force_env%root_section, para_env, &
     752        1818 :                                   subsys_section, corevel_read)
     753             : 
     754        1818 :       IF (.NOT. (atomvel_explicit .OR. atomvel_read)) RETURN
     755         266 :       success = .TRUE.
     756             : 
     757         266 :       IF (.NOT. atomvel_read) THEN
     758             :          ! Read the atom velocities if explicitly given in the input file
     759         220 :          CALL section_vals_list_get(atomvel_section, "_DEFAULT_KEYWORD_", list=atom_list)
     760       52950 :          DO i = 1, natoms
     761       52730 :             is_ok = cp_sll_val_next(atom_list, val)
     762       52730 :             CALL val_get(val, r_vals=vel)
     763      369330 :             part(i)%v = vel
     764             :          END DO
     765             :       END IF
     766       57412 :       DO i = 1, natoms
     767         266 :          SELECT CASE (is_fixed(i))
     768             :          CASE (use_perd_x)
     769           0 :             part(i)%v(1) = 0.0_dp
     770             :          CASE (use_perd_y)
     771           0 :             part(i)%v(2) = 0.0_dp
     772             :          CASE (use_perd_z)
     773           0 :             part(i)%v(3) = 0.0_dp
     774             :          CASE (use_perd_xy)
     775           0 :             part(i)%v(1) = 0.0_dp
     776           0 :             part(i)%v(2) = 0.0_dp
     777             :          CASE (use_perd_xz)
     778           0 :             part(i)%v(1) = 0.0_dp
     779           0 :             part(i)%v(3) = 0.0_dp
     780             :          CASE (use_perd_yz)
     781           0 :             part(i)%v(2) = 0.0_dp
     782           0 :             part(i)%v(3) = 0.0_dp
     783             :          CASE (use_perd_xyz)
     784       57224 :             part(i)%v = 0.0_dp
     785             :          END SELECT
     786             :       END DO
     787         266 :       IF (shell_present) THEN
     788          48 :          IF (shellvel_explicit) THEN
     789             :             ! If the atoms positions are given (?) and core and shell velocities are
     790             :             ! present in the input, read the latter.
     791          24 :             CALL section_vals_list_get(shellvel_section, "_DEFAULT_KEYWORD_", list=shell_list)
     792          24 :             CALL section_vals_list_get(corevel_section, "_DEFAULT_KEYWORD_", list=core_list)
     793        2328 :             DO i = 1, nshell
     794        2304 :                is_ok = cp_sll_val_next(shell_list, val)
     795        2304 :                CALL val_get(val, r_vals=vel)
     796       16128 :                shell_part(i)%v = vel
     797        2304 :                is_ok = cp_sll_val_next(core_list, val)
     798        2304 :                CALL val_get(val, r_vals=vel)
     799       16152 :                core_part(i)%v = vel
     800             :             END DO
     801             :          ELSE
     802          24 :             IF (.NOT. (shellvel_read .AND. corevel_read)) THEN
     803             :                ! Otherwise, just copy atom velocties into shell and core velocities.
     804           8 :                CALL clone_core_shell_vel(part, shell_part, core_part)
     805             :             END IF
     806             :          END IF
     807             :       END IF
     808             : 
     809             :       ! compute vcom, ecom and ekin
     810         266 :       CALL compute_vcom(part, is_fixed, vcom, ecom)
     811         266 :       ekin = compute_ekin(part) - ecom
     812             : 
     813         266 :       IF (simpar%do_thermal_region) THEN
     814          12 :          CALL get_md_env(md_env, thermal_regions=thermal_regions)
     815          12 :          IF (ASSOCIATED(thermal_regions)) THEN
     816          12 :             rescale_regions = thermal_regions%force_rescaling
     817             :          END IF
     818             :       ELSE
     819             :          rescale_regions = .FALSE.
     820             :       END IF
     821         266 :       IF (simpar%nfree /= 0 .AND. (force_rescaling .OR. rescale_regions)) THEN
     822          24 :          IF (simpar%do_thermal_region) THEN
     823           0 :             CALL rescale_vel_region(part, md_env, simpar)
     824             :          ELSE
     825          24 :             CALL rescale_vel(part, simpar, ekin, vcom=vcom)
     826             :          END IF
     827             : 
     828             :          ! After rescaling, the core and shell velocities must also adapt.
     829        1894 :          DO i = 1, natoms
     830        1870 :             shell_index = part(i)%shell_index
     831        2136 :             IF (shell_present .AND. shell_index /= 0) THEN
     832           0 :                atomic_kind => part(i)%atomic_kind
     833           0 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell=shell)
     834           0 :                fac_masss = shell%mass_shell/mass
     835           0 :                fac_massc = shell%mass_core/mass
     836           0 :                vs = shell_part(shell_index)%v
     837           0 :                vc = core_part(shell_index)%v
     838             : 
     839           0 :                shell_part(shell_index)%v(1) = part(i)%v(1) + fac_massc*(vs(1) - vc(1))
     840           0 :                shell_part(shell_index)%v(2) = part(i)%v(2) + fac_massc*(vs(2) - vc(2))
     841           0 :                shell_part(shell_index)%v(3) = part(i)%v(3) + fac_massc*(vs(3) - vc(3))
     842           0 :                core_part(shell_index)%v(1) = part(i)%v(1) + fac_masss*(vc(1) - vs(1))
     843           0 :                core_part(shell_index)%v(2) = part(i)%v(2) + fac_masss*(vc(2) - vs(2))
     844           0 :                core_part(shell_index)%v(3) = part(i)%v(3) + fac_masss*(vc(3) - vs(3))
     845             :             END IF
     846             :          END DO
     847             :       END IF
     848        1818 :    END SUBROUTINE read_input_velocities
     849             : 
     850             : ! **************************************************************************************************
     851             : !> \brief Initializing velocities AND positions randomly on all processors, based on vibrational
     852             : !>        modes of the system, so that the starting coordinates are already very close to
     853             : !>        canonical ensumble corresponding to temperature of a head bath.
     854             : !> \param simpar          : MD simulation parameters
     855             : !> \param particles       : global array of all particles
     856             : !> \param md_section      : MD input subsection
     857             : !> \param vib_section     : vibrational analysis input section
     858             : !> \param force_env       : force environment container
     859             : !> \param global_env      : global environment container
     860             : !> \param shell_present   : if core-shell model is used
     861             : !> \param shell_particles : global array of all shell particles in shell model
     862             : !> \param core_particles  : global array of all core particles in shell model
     863             : !> \param is_fixed        : array of size of total number of atoms, that determines if any
     864             : !>                          cartesian components are fixed
     865             : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
     866             : ! **************************************************************************************************
     867           2 :    SUBROUTINE generate_coords_vels_vib(simpar, &
     868             :                                        particles, &
     869             :                                        md_section, &
     870             :                                        vib_section, &
     871             :                                        force_env, &
     872             :                                        global_env, &
     873             :                                        shell_present, &
     874             :                                        shell_particles, &
     875             :                                        core_particles, &
     876           2 :                                        is_fixed)
     877             :       TYPE(simpar_type), POINTER                         :: simpar
     878             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles
     879             :       TYPE(section_vals_type), POINTER                   :: md_section, vib_section
     880             :       TYPE(force_env_type), POINTER                      :: force_env
     881             :       TYPE(global_environment_type), POINTER             :: global_env
     882             :       LOGICAL, INTENT(IN)                                :: shell_present
     883             :       TYPE(particle_type), DIMENSION(:), POINTER         :: shell_particles, core_particles
     884             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: is_fixed
     885             : 
     886             :       INTEGER                                            :: dof, fixed_dof, iatom, ii, imode, &
     887             :                                                             my_dof, natoms, shell_index
     888             :       REAL(KIND=dp)                                      :: Erand, mass, my_phase, ratio, temperature
     889           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues, phase, random
     890           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dr, eigenvectors
     891             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     892             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     893           2 :       TYPE(rng_stream_type), ALLOCATABLE                 :: random_stream
     894             : 
     895           2 :       CALL cite_reference(West2006)
     896           2 :       natoms = SIZE(particles)
     897           2 :       temperature = simpar%temp_ext
     898           2 :       my_dof = 3*natoms
     899           6 :       ALLOCATE (eigenvalues(my_dof))
     900           8 :       ALLOCATE (eigenvectors(my_dof, my_dof))
     901           6 :       ALLOCATE (phase(my_dof))
     902           6 :       ALLOCATE (random(my_dof))
     903           6 :       ALLOCATE (dr(3, natoms))
     904           2 :       CALL force_env_get(force_env=force_env, para_env=para_env)
     905             :       ! read vibration modes
     906             :       CALL read_vib_eigs_unformatted(md_section, &
     907             :                                      vib_section, &
     908             :                                      para_env, &
     909             :                                      dof, &
     910             :                                      eigenvalues, &
     911           2 :                                      eigenvectors)
     912           2 :       IF (my_dof .NE. dof) THEN
     913             :          CALL cp_abort(__LOCATION__, &
     914             :                        "number of degrees of freedom in vibrational analysis data "// &
     915           0 :                        "do not match total number of cartesian degrees of freedom")
     916             :       END IF
     917             :       ! read phases
     918           2 :       CALL section_vals_val_get(md_section, "INITIAL_VIBRATION%PHASE", r_val=my_phase)
     919           2 :       my_phase = MIN(1.0_dp, my_phase)
     920             :       ! generate random numbers
     921           2 :       random_stream = rng_stream_type(name="MD_INIT_VIB", distribution_type=UNIFORM)
     922           2 :       CALL random_stream%fill(random)
     923           2 :       IF (my_phase .LT. 0.0_dp) THEN
     924           0 :          CALL random_stream%fill(phase)
     925             :       ELSE
     926          20 :          phase = my_phase
     927             :       END IF
     928           2 :       DEALLOCATE (random_stream)
     929             : 
     930             :       ! the first three modes are acoustic with zero frequencies,
     931             :       ! exclude these from considerations
     932           2 :       my_dof = dof - 3
     933             :       ! randomly selects energy from distribution about kT, all
     934             :       ! energies are scaled so that the sum over vibration modes gives
     935             :       ! exactly my_dof*kT. Note that k = 1.0 in atomic units
     936           2 :       Erand = 0.0_dp
     937          14 :       DO imode = 4, dof
     938          14 :          Erand = Erand - temperature*LOG(1.0_dp - random(imode))
     939             :       END DO
     940             :       ! need to take into account of fixed constraints too
     941           2 :       fixed_dof = 0
     942           8 :       DO iatom = 1, natoms
     943           2 :          SELECT CASE (is_fixed(iatom))
     944             :          CASE (use_perd_x, use_perd_y, use_perd_z)
     945           0 :             fixed_dof = fixed_dof + 1
     946             :          CASE (use_perd_xy, use_perd_xz, use_perd_yz)
     947           0 :             fixed_dof = fixed_dof + 2
     948             :          CASE (use_perd_xyz)
     949           6 :             fixed_dof = fixed_dof + 3
     950             :          END SELECT
     951             :       END DO
     952           2 :       my_dof = my_dof - fixed_dof
     953           2 :       ratio = REAL(my_dof, KIND=dp)*temperature/Erand
     954             :       ! update  velocities AND positions
     955           8 :       DO iatom = 1, natoms
     956           6 :          atomic_kind => particles(iatom)%atomic_kind
     957           6 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     958           8 :          SELECT CASE (is_fixed(iatom))
     959             :          CASE (use_perd_x)
     960           0 :             DO ii = 2, 3
     961             :                dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
     962           0 :                                                 eigenvectors, random, phase, dof, ratio)
     963             :                particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
     964             :                                                          eigenvectors, random, phase, dof, &
     965           0 :                                                          ratio)
     966             :             END DO
     967             :          CASE (use_perd_y)
     968           0 :             DO ii = 1, 3, 2
     969             :                dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
     970           0 :                                                 eigenvectors, random, phase, dof, ratio)
     971             :                particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
     972             :                                                          eigenvectors, random, phase, dof, &
     973           0 :                                                          ratio)
     974             :             END DO
     975             :          CASE (use_perd_z)
     976           0 :             DO ii = 1, 2
     977             :                dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
     978           0 :                                                 eigenvectors, random, phase, dof, ratio)
     979             :                particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
     980             :                                                          eigenvectors, random, phase, dof, &
     981           0 :                                                          ratio)
     982             :             END DO
     983             :          CASE (use_perd_xy)
     984             :             dr(3, iatom) = dr_from_vib_data(iatom, 3, mass, temperature, eigenvalues, &
     985           0 :                                             eigenvectors, random, phase, dof, ratio)
     986             :             particles(iatom)%v(3) = dv_from_vib_data(iatom, 3, mass, temperature, &
     987             :                                                      eigenvectors, random, phase, dof, &
     988           0 :                                                      ratio)
     989             :          CASE (use_perd_xz)
     990             :             dr(2, iatom) = dr_from_vib_data(iatom, 2, mass, temperature, eigenvalues, &
     991           0 :                                             eigenvectors, random, phase, dof, ratio)
     992             :             particles(iatom)%v(2) = dv_from_vib_data(iatom, 2, mass, temperature, &
     993             :                                                      eigenvectors, random, phase, dof, &
     994           0 :                                                      ratio)
     995             :          CASE (use_perd_yz)
     996             :             dr(1, iatom) = dr_from_vib_data(iatom, 1, mass, temperature, eigenvalues, &
     997           0 :                                             eigenvectors, random, phase, dof, ratio)
     998             :             particles(iatom)%v(1) = dv_from_vib_data(iatom, 1, mass, temperature, &
     999             :                                                      eigenvectors, random, phase, dof, &
    1000           0 :                                                      ratio)
    1001             :          CASE (use_perd_none)
    1002          24 :             DO ii = 1, 3
    1003             :                dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
    1004          18 :                                                 eigenvectors, random, phase, dof, ratio)
    1005             :                particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
    1006             :                                                          eigenvectors, random, phase, dof, &
    1007          24 :                                                          ratio)
    1008             :             END DO
    1009             :          END SELECT
    1010             :       END DO ! iatom
    1011             :       ! free memory
    1012           2 :       DEALLOCATE (eigenvalues)
    1013           2 :       DEALLOCATE (eigenvectors)
    1014           2 :       DEALLOCATE (phase)
    1015           2 :       DEALLOCATE (random)
    1016             :       ! update particle coordinates
    1017           8 :       DO iatom = 1, natoms
    1018          26 :          particles(iatom)%r(:) = particles(iatom)%r(:) + dr(:, iatom)
    1019             :       END DO
    1020             :       ! update core-shell model coordinates
    1021           2 :       IF (shell_present) THEN
    1022             :          ! particles have moved, and for core-shell model this means
    1023             :          ! the cores and shells must also move by the same amount. The
    1024             :          ! shell positions will then be optimised if needed
    1025           0 :          shell_index = particles(iatom)%shell_index
    1026           0 :          IF (shell_index .NE. 0) THEN
    1027             :             core_particles(shell_index)%r(:) = core_particles(shell_index)%r(:) + &
    1028           0 :                                                dr(:, iatom)
    1029             :             shell_particles(shell_index)%r(:) = shell_particles(shell_index)%r(:) + &
    1030           0 :                                                 dr(:, iatom)
    1031             :          END IF
    1032             :          CALL optimize_shell_core(force_env, &
    1033             :                                   particles, &
    1034             :                                   shell_particles, &
    1035             :                                   core_particles, &
    1036           0 :                                   global_env)
    1037             :       END IF
    1038             :       ! cleanup
    1039           2 :       DEALLOCATE (dr)
    1040           6 :    END SUBROUTINE generate_coords_vels_vib
    1041             : 
    1042             : ! **************************************************************************************************
    1043             : !> \brief calculates componbent of initial velocity of an atom from vibreational modes
    1044             : !> \param iatom        : global atomic index
    1045             : !> \param icart        : cartesian index (1, 2 or 3)
    1046             : !> \param mass         : atomic mass
    1047             : !> \param temperature  : target temperature of canonical ensemble
    1048             : !> \param eigenvalues  : array containing all cartesian vibrational mode eigenvalues (frequencies)
    1049             : !> \param eigenvectors : array containing all corresponding vibrational mode eigenvectors
    1050             : !>                       (displacements)
    1051             : !> \param random       : array containing uniform distributed random numbers, must be the size
    1052             : !>                       of 3*natoms. Numbers must be between 0 and 1
    1053             : !> \param phase        : array containing numbers between 0 and 1 that determines for each
    1054             : !>                       vibration mode the ratio of potential energy vs kinetic energy contribution
    1055             : !>                       to total energy
    1056             : !> \param dof          : total number of degrees of freedom, = 3*natoms
    1057             : !> \param scale        : scale to make sure the sum of vibrational modes give the correct energy
    1058             : !> \return : outputs icart-th cartesian component of initial position of atom iatom
    1059             : !> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
    1060             : ! **************************************************************************************************
    1061          18 :    PURE FUNCTION dr_from_vib_data(iatom, &
    1062             :                                   icart, &
    1063             :                                   mass, &
    1064             :                                   temperature, &
    1065          36 :                                   eigenvalues, &
    1066          18 :                                   eigenvectors, &
    1067          18 :                                   random, &
    1068          18 :                                   phase, &
    1069             :                                   dof, &
    1070             :                                   scale) &
    1071             :       RESULT(res)
    1072             :       INTEGER, INTENT(IN)                                :: iatom, icart
    1073             :       REAL(KIND=dp), INTENT(IN)                          :: mass, temperature
    1074             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: eigenvalues
    1075             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: eigenvectors
    1076             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: random, phase
    1077             :       INTEGER, INTENT(IN)                                :: dof
    1078             :       REAL(KIND=dp), INTENT(IN)                          :: scale
    1079             :       REAL(KIND=dp)                                      :: res
    1080             : 
    1081             :       INTEGER                                            :: imode, ind
    1082             : 
    1083          18 :       res = 0.0_dp
    1084             :       ! assuming the eigenvalues are sorted in ascending order, the
    1085             :       ! first three modes are acoustic with zero frequencies. These are
    1086             :       ! excluded from considerations, and should have been reflected in
    1087             :       ! the calculation of scale outside this function
    1088          18 :       IF (mass .GT. 0.0_dp) THEN
    1089             :          ! eigenvector rows assumed to be grouped in atomic blocks
    1090          18 :          ind = (iatom - 1)*3 + icart
    1091         126 :          DO imode = 4, dof
    1092             :             res = res + &
    1093             :                   SQRT(-2.0_dp*scale*temperature*LOG(1 - random(imode))/mass)/ &
    1094             :                   eigenvalues(imode)* &
    1095             :                   eigenvectors(ind, imode)* &
    1096         126 :                   COS(2.0_dp*pi*phase(imode))
    1097             :          END DO
    1098             :       END IF
    1099          18 :    END FUNCTION dr_from_vib_data
    1100             : 
    1101             : ! **************************************************************************************************
    1102             : !> \brief calculates componbent of initial velocity of an atom from vibreational modes
    1103             : !> \param iatom        : global atomic index
    1104             : !> \param icart        : cartesian index (1, 2 or 3)
    1105             : !> \param mass         : atomic mass
    1106             : !> \param temperature  : target temperature of canonical ensemble
    1107             : !> \param eigenvectors : array containing all corresponding vibrational mode eigenvectors
    1108             : !>                       (displacements)
    1109             : !> \param random       : array containing uniform distributed random numbers, must be the size
    1110             : !>                       of 3*natoms. Numbers must be between 0 and 1
    1111             : !> \param phase        : array containing numbers between 0 and 1 that determines for each
    1112             : !>                       vibration mode the ratio of potential energy vs kinetic energy contribution
    1113             : !>                       to total energy
    1114             : !> \param dof          : total number of degrees of freedom, = 3*natoms
    1115             : !> \param scale        : scale to make sure the sum of vibrational modes give the correct energy
    1116             : !> \return : outputs icart-th cartesian component of initial velocity of atom iatom
    1117             : !> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
    1118             : ! **************************************************************************************************
    1119          18 :    PURE FUNCTION dv_from_vib_data(iatom, &
    1120             :                                   icart, &
    1121             :                                   mass, &
    1122             :                                   temperature, &
    1123          36 :                                   eigenvectors, &
    1124          18 :                                   random, &
    1125          18 :                                   phase, &
    1126             :                                   dof, &
    1127             :                                   scale) &
    1128             :       RESULT(res)
    1129             :       INTEGER, INTENT(IN)                                :: iatom, icart
    1130             :       REAL(KIND=dp), INTENT(IN)                          :: mass, temperature
    1131             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: eigenvectors
    1132             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: random, phase
    1133             :       INTEGER, INTENT(IN)                                :: dof
    1134             :       REAL(KIND=dp), INTENT(IN)                          :: scale
    1135             :       REAL(KIND=dp)                                      :: res
    1136             : 
    1137             :       INTEGER                                            :: imode, ind
    1138             : 
    1139          18 :       res = 0.0_dp
    1140             :       ! assuming the eigenvalues are sorted in ascending order, the
    1141             :       ! first three modes are acoustic with zero frequencies. These are
    1142             :       ! excluded from considerations, and should have been reflected in
    1143             :       ! the calculation of scale outside this function
    1144          18 :       IF (mass .GT. 0.0_dp) THEN
    1145             :          ! eigenvector rows assumed to be grouped in atomic blocks
    1146          18 :          ind = (iatom - 1)*3 + icart
    1147         126 :          DO imode = 4, dof
    1148             :             res = res - &
    1149             :                   SQRT(-2.0_dp*scale*temperature*LOG(1 - random(imode))/mass)* &
    1150             :                   eigenvectors(ind, imode)* &
    1151         126 :                   SIN(2.0_dp*pi*phase(imode))
    1152             :          END DO
    1153             :       END IF
    1154          18 :    END FUNCTION dv_from_vib_data
    1155             : 
    1156             : ! **************************************************************************************************
    1157             : !> \brief Initializing velocities deterministically on all processors, if not given in input
    1158             : !> \param simpar ...
    1159             : !> \param part ...
    1160             : !> \param force_env ...
    1161             : !> \param globenv ...
    1162             : !> \param md_env ...
    1163             : !> \param shell_present ...
    1164             : !> \param shell_part ...
    1165             : !> \param core_part ...
    1166             : !> \param is_fixed ...
    1167             : !> \param iw ...
    1168             : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
    1169             : ! **************************************************************************************************
    1170        1550 :    SUBROUTINE generate_velocities(simpar, part, force_env, globenv, md_env, &
    1171        1550 :                                   shell_present, shell_part, core_part, is_fixed, iw)
    1172             :       TYPE(simpar_type), POINTER                         :: simpar
    1173             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
    1174             :       TYPE(force_env_type), POINTER                      :: force_env
    1175             :       TYPE(global_environment_type), POINTER             :: globenv
    1176             :       TYPE(md_environment_type), POINTER                 :: md_env
    1177             :       LOGICAL, INTENT(IN)                                :: shell_present
    1178             :       TYPE(particle_type), DIMENSION(:), POINTER         :: shell_part, core_part
    1179             :       INTEGER, DIMENSION(:), INTENT(INOUT)               :: is_fixed
    1180             :       INTEGER, INTENT(IN)                                :: iw
    1181             : 
    1182             :       INTEGER                                            :: i, natoms
    1183             :       REAL(KIND=dp)                                      :: mass
    1184             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1185             : 
    1186        1550 :       NULLIFY (atomic_kind)
    1187        1550 :       natoms = SIZE(part)
    1188             : 
    1189      432670 :       DO i = 1, natoms
    1190      431120 :          atomic_kind => part(i)%atomic_kind
    1191      431120 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    1192      431120 :          part(i)%v(1) = 0.0_dp
    1193      431120 :          part(i)%v(2) = 0.0_dp
    1194      431120 :          part(i)%v(3) = 0.0_dp
    1195      432670 :          IF (mass .NE. 0.0) THEN
    1196      431186 :             SELECT CASE (is_fixed(i))
    1197             :             CASE (use_perd_x)
    1198         512 :                part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1199         512 :                part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1200             :             CASE (use_perd_y)
    1201         512 :                part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1202         512 :                part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1203             :             CASE (use_perd_z)
    1204         512 :                part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1205         512 :                part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1206             :             CASE (use_perd_xy)
    1207         512 :                part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1208             :             CASE (use_perd_xz)
    1209           0 :                part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1210             :             CASE (use_perd_yz)
    1211           0 :                part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1212             :             CASE (use_perd_none)
    1213      425508 :                part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1214      425508 :                part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1215      856182 :                part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
    1216             :             END SELECT
    1217             :          END IF
    1218             :       END DO
    1219             : 
    1220        1550 :       CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
    1221        1550 :       CALL soften_velocities(simpar, part, force_env, md_env, is_fixed, iw)
    1222             : 
    1223             :       ! Initialize the core and the shell velocity. Atom velocities are just
    1224             :       ! copied so that the initial relative core-shell velocity is zero.
    1225        1550 :       IF (shell_present) THEN
    1226          84 :          CALL optimize_shell_core(force_env, part, shell_part, core_part, globenv)
    1227             :       END IF
    1228        1550 :    END SUBROUTINE generate_velocities
    1229             : 
    1230             : ! **************************************************************************************************
    1231             : !> \brief Direct velocities along a low-curvature direction in order to
    1232             : !>        favors MD trajectories to cross rapidly over small energy barriers
    1233             : !>        into neighboring basins.
    1234             : !> \param simpar ...
    1235             : !> \param part ...
    1236             : !> \param force_env ...
    1237             : !> \param md_env ...
    1238             : !> \param is_fixed ...
    1239             : !> \param iw ...
    1240             : !> \author Ole Schuett
    1241             : ! **************************************************************************************************
    1242        1550 :    SUBROUTINE soften_velocities(simpar, part, force_env, md_env, is_fixed, iw)
    1243             :       TYPE(simpar_type), POINTER                         :: simpar
    1244             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
    1245             :       TYPE(force_env_type), POINTER                      :: force_env
    1246             :       TYPE(md_environment_type), POINTER                 :: md_env
    1247             :       INTEGER, DIMENSION(:), INTENT(INOUT)               :: is_fixed
    1248             :       INTEGER, INTENT(IN)                                :: iw
    1249             : 
    1250             :       INTEGER                                            :: i, k
    1251        1490 :       REAL(KIND=dp), DIMENSION(SIZE(part), 3)            :: F, F_t, N, x0
    1252             : 
    1253        1550 :       IF (simpar%soften_nsteps <= 0) RETURN !nothing todo
    1254             : 
    1255         660 :       IF (ANY(is_fixed /= use_perd_none)) &
    1256           0 :          CPABORT("Velocitiy softening with constraints is not supported.")
    1257             : 
    1258             :       !backup positions
    1259         660 :       DO i = 1, SIZE(part)
    1260        2460 :          x0(i, :) = part(i)%r
    1261             :       END DO
    1262             : 
    1263        1260 :       DO k = 1, simpar%soften_nsteps
    1264             : 
    1265             :          !use normalized velocities as displace direction
    1266       13200 :          DO i = 1, SIZE(part)
    1267       49200 :             N(i, :) = part(i)%v
    1268             :          END DO
    1269       80400 :          N = N/SQRT(SUM(N**2))
    1270             : 
    1271             :          ! displace system temporarly to calculate forces
    1272       13200 :          DO i = 1, SIZE(part)
    1273       49200 :             part(i)%r = part(i)%r + simpar%soften_delta*N(i, :)
    1274             :          END DO
    1275        1200 :          CALL force_env_calc_energy_force(force_env)
    1276             : 
    1277             :          ! calculate velocity update direction F_t
    1278       13200 :          DO i = 1, SIZE(part)
    1279       49200 :             F(i, :) = part(i)%f
    1280             :          END DO
    1281       80400 :          F_t = F - N*SUM(N*F)
    1282             : 
    1283             :          ! restore positions and update velocities
    1284       13200 :          DO i = 1, SIZE(part)
    1285       48000 :             part(i)%r = x0(i, :)
    1286       49200 :             part(i)%v = part(i)%v + simpar%soften_alpha*F_t(i, :)
    1287             :          END DO
    1288             : 
    1289        1260 :          CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
    1290             :       END DO
    1291             : 
    1292          60 :       IF (iw > 0) THEN
    1293           0 :          WRITE (iw, "(A,T71, I10)") " Velocities softening Steps: ", simpar%soften_nsteps
    1294           0 :          WRITE (iw, "(A,T71, E10.3)") " Velocities softening NORM(F_t): ", SQRT(SUM(F_t**2))
    1295             :       END IF
    1296          60 :    END SUBROUTINE soften_velocities
    1297             : 
    1298             : ! **************************************************************************************************
    1299             : !> \brief Scale velocities according to temperature and remove rigid body motion.
    1300             : !> \param simpar ...
    1301             : !> \param part ...
    1302             : !> \param force_env ...
    1303             : !> \param md_env ...
    1304             : !> \param is_fixed ...
    1305             : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
    1306             : ! **************************************************************************************************
    1307        2750 :    SUBROUTINE normalize_velocities(simpar, part, force_env, md_env, is_fixed)
    1308             :       TYPE(simpar_type), POINTER                         :: simpar
    1309             :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
    1310             :       TYPE(force_env_type), POINTER                      :: force_env
    1311             :       TYPE(md_environment_type), POINTER                 :: md_env
    1312             :       INTEGER, DIMENSION(:), INTENT(INOUT)               :: is_fixed
    1313             : 
    1314             :       REAL(KIND=dp)                                      :: ekin
    1315             :       REAL(KIND=dp), DIMENSION(3)                        :: rcom, vang, vcom
    1316             :       TYPE(cell_type), POINTER                           :: cell
    1317             : 
    1318        2750 :       NULLIFY (cell)
    1319             : 
    1320             :       ! Subtract the vcom
    1321        2750 :       CALL compute_vcom(part, is_fixed, vcom)
    1322        2750 :       CALL subtract_vcom(part, is_fixed, vcom)
    1323             :       ! If requested and the system is not periodic, subtract the angular velocity
    1324        2750 :       CALL force_env_get(force_env, cell=cell)
    1325       11000 :       IF (SUM(cell%perd(1:3)) == 0 .AND. simpar%angvel_zero) THEN
    1326           4 :          CALL compute_rcom(part, is_fixed, rcom)
    1327           4 :          CALL compute_vang(part, is_fixed, rcom, vang)
    1328           4 :          CALL subtract_vang(part, is_fixed, rcom, vang)
    1329             :       END IF
    1330             :       ! Rescale the velocities
    1331        2750 :       IF (simpar%do_thermal_region) THEN
    1332           2 :          CALL rescale_vel_region(part, md_env, simpar)
    1333             :       ELSE
    1334        2748 :          ekin = compute_ekin(part)
    1335        2748 :          CALL rescale_vel(part, simpar, ekin)
    1336             :       END IF
    1337        2750 :    END SUBROUTINE normalize_velocities
    1338             : 
    1339             : ! **************************************************************************************************
    1340             : !> \brief Computes Ekin, VCOM and Temp for particles
    1341             : !> \param subsys ...
    1342             : !> \param md_ener ...
    1343             : !> \param vsubtract ...
    1344             : !> \par History
    1345             : !>     Teodoro Laino - University of Zurich - 09.2007 [tlaino]
    1346             : ! **************************************************************************************************
    1347          42 :    SUBROUTINE reset_vcom(subsys, md_ener, vsubtract)
    1348             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1349             :       TYPE(md_ener_type), POINTER                        :: md_ener
    1350             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: vsubtract
    1351             : 
    1352             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'reset_vcom'
    1353             : 
    1354             :       INTEGER                                            :: atom, handle, iatom, ikind, natom, &
    1355             :                                                             shell_index
    1356          42 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list
    1357             :       LOGICAL                                            :: is_shell
    1358             :       REAL(KIND=dp)                                      :: ekin_old, imass_c, imass_s, mass, v2
    1359             :       REAL(KIND=dp), DIMENSION(3)                        :: tmp, v
    1360             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    1361             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1362             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
    1363             :                                                             shell_particles
    1364             :       TYPE(shell_kind_type), POINTER                     :: shell
    1365             : 
    1366          42 :       NULLIFY (particles, atomic_kind, atomic_kinds, atom_list, shell)
    1367          42 :       CALL timeset(routineN, handle)
    1368             : 
    1369             :       CALL cp_subsys_get(subsys, &
    1370             :                          atomic_kinds=atomic_kinds, &
    1371             :                          particles=particles, &
    1372             :                          shell_particles=shell_particles, &
    1373          42 :                          core_particles=core_particles)
    1374             : 
    1375          42 :       ekin_old = md_ener%ekin
    1376             :       ! Possibly subtract a quantity from all velocities
    1377         126 :       DO ikind = 1, atomic_kinds%n_els
    1378          84 :          atomic_kind => atomic_kinds%els(ikind)
    1379             :          CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, &
    1380          84 :                               natom=natom, mass=mass, shell_active=is_shell, shell=shell)
    1381         126 :          IF (is_shell) THEN
    1382         336 :             tmp = 0.5_dp*vsubtract*mass
    1383          84 :             imass_s = 1.0_dp/shell%mass_shell
    1384          84 :             imass_c = 1.0_dp/shell%mass_core
    1385        3780 :             DO iatom = 1, natom
    1386        3696 :                atom = atom_list(iatom)
    1387        3696 :                shell_index = particles%els(atom)%shell_index
    1388       14784 :                shell_particles%els(shell_index)%v = shell_particles%els(shell_index)%v - tmp*imass_s
    1389       14784 :                core_particles%els(shell_index)%v = core_particles%els(shell_index)%v - tmp*imass_c
    1390       14868 :                particles%els(atom)%v = particles%els(atom)%v - vsubtract
    1391             :             END DO
    1392             :          ELSE
    1393           0 :             DO iatom = 1, natom
    1394           0 :                atom = atom_list(iatom)
    1395           0 :                particles%els(atom)%v = particles%els(atom)%v - vsubtract
    1396             :             END DO
    1397             :          END IF
    1398             :       END DO
    1399             :       ! Compute Kinetic Energy and COM Velocity
    1400         168 :       md_ener%vcom = 0.0_dp
    1401          42 :       md_ener%total_mass = 0.0_dp
    1402          42 :       md_ener%ekin = 0.0_dp
    1403         126 :       DO ikind = 1, atomic_kinds%n_els
    1404          84 :          atomic_kind => atomic_kinds%els(ikind)
    1405          84 :          CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom)
    1406          84 :          v2 = 0.0_dp
    1407          84 :          v = 0.0_dp
    1408        3780 :          DO iatom = 1, natom
    1409        3696 :             atom = atom_list(iatom)
    1410       14784 :             v2 = v2 + SUM(particles%els(atom)%v**2)
    1411        3696 :             v(1) = v(1) + particles%els(atom)%v(1)
    1412        3696 :             v(2) = v(2) + particles%els(atom)%v(2)
    1413        3780 :             v(3) = v(3) + particles%els(atom)%v(3)
    1414             :          END DO
    1415          84 :          md_ener%ekin = md_ener%ekin + 0.5_dp*mass*v2
    1416          84 :          md_ener%vcom(1) = md_ener%vcom(1) + mass*v(1)
    1417          84 :          md_ener%vcom(2) = md_ener%vcom(2) + mass*v(2)
    1418          84 :          md_ener%vcom(3) = md_ener%vcom(3) + mass*v(3)
    1419         210 :          md_ener%total_mass = md_ener%total_mass + REAL(natom, KIND=dp)*mass
    1420             :       END DO
    1421         168 :       md_ener%vcom = md_ener%vcom/md_ener%total_mass
    1422          42 :       md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
    1423          42 :       IF (md_ener%nfree /= 0) THEN
    1424          42 :          md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
    1425             :       END IF
    1426          42 :       CALL timestop(handle)
    1427             : 
    1428          42 :    END SUBROUTINE reset_vcom
    1429             : 
    1430             : ! **************************************************************************************************
    1431             : !> \brief Scale velocities to get the correct temperature
    1432             : !> \param subsys ...
    1433             : !> \param md_ener ...
    1434             : !> \param temp_expected ...
    1435             : !> \param temp_tol ...
    1436             : !> \param iw ...
    1437             : !> \par History
    1438             : !>     Teodoro Laino - University of Zurich - 09.2007 [tlaino]
    1439             : ! **************************************************************************************************
    1440       12814 :    SUBROUTINE scale_velocity(subsys, md_ener, temp_expected, temp_tol, iw)
    1441             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1442             :       TYPE(md_ener_type), POINTER                        :: md_ener
    1443             :       REAL(KIND=dp), INTENT(IN)                          :: temp_expected, temp_tol
    1444             :       INTEGER, INTENT(IN)                                :: iw
    1445             : 
    1446             :       REAL(KIND=dp)                                      :: ekin_old, scale, temp_old
    1447             : 
    1448       12814 :       IF (ABS(temp_expected - md_ener%temp_part/kelvin) > temp_tol) THEN
    1449        2562 :          scale = 0.0_dp
    1450        2562 :          IF (md_ener%temp_part > 0.0_dp) scale = SQRT((temp_expected/md_ener%temp_part)*kelvin)
    1451        2562 :          ekin_old = md_ener%ekin
    1452        2562 :          temp_old = md_ener%temp_part
    1453        2562 :          md_ener%ekin = 0.0_dp
    1454        2562 :          md_ener%temp_part = 0.0_dp
    1455       10248 :          md_ener%vcom = 0.0_dp
    1456        2562 :          md_ener%total_mass = 0.0_dp
    1457             : 
    1458        2562 :          CALL scale_velocity_low(subsys, scale, ireg=0, ekin=md_ener%ekin, vcom=md_ener%vcom)
    1459        2562 :          IF (md_ener%nfree /= 0) THEN
    1460        2562 :             md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
    1461             :          END IF
    1462        2562 :          md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
    1463        2562 :          IF (iw > 0) THEN
    1464             :             WRITE (UNIT=iw, FMT='(/,T2,A)') &
    1465        1281 :                'MD_VEL| Temperature scaled to requested temperature'
    1466             :             WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
    1467        1281 :                'MD_VEL| Old temperature [K]', temp_old, &
    1468        2562 :                'MD_VEL| New temperature [K]', md_ener%temp_part
    1469             :          END IF
    1470             :       END IF
    1471             : 
    1472       12814 :    END SUBROUTINE scale_velocity
    1473             : 
    1474             : ! **************************************************************************************************
    1475             : !> \brief Scale velocities of set of regions
    1476             : !> \param md_env ...
    1477             : !> \param subsys ...
    1478             : !> \param md_ener ...
    1479             : !> \param simpar ...
    1480             : !> \param iw ...
    1481             : !> \par author MI
    1482             : ! **************************************************************************************************
    1483          96 :    SUBROUTINE scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
    1484             : 
    1485             :       TYPE(md_environment_type), POINTER                 :: md_env
    1486             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1487             :       TYPE(md_ener_type), POINTER                        :: md_ener
    1488             :       TYPE(simpar_type), POINTER                         :: simpar
    1489             :       INTEGER, INTENT(IN)                                :: iw
    1490             : 
    1491             :       INTEGER                                            :: ireg, nfree, nfree_done, nregions
    1492             :       REAL(KIND=dp)                                      :: ekin, ekin_old, ekin_total_new, fscale, &
    1493             :                                                             vcom(3), vcom_total(3)
    1494          96 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: temp_new, temp_old
    1495             :       TYPE(particle_list_type), POINTER                  :: particles
    1496          96 :       TYPE(particle_type), DIMENSION(:), POINTER         :: part
    1497             :       TYPE(thermal_region_type), POINTER                 :: t_region
    1498             :       TYPE(thermal_regions_type), POINTER                :: thermal_regions
    1499             : 
    1500          96 :       NULLIFY (particles, part, thermal_regions, t_region)
    1501          96 :       CALL cp_subsys_get(subsys, particles=particles)
    1502          96 :       part => particles%els
    1503          96 :       CALL get_md_env(md_env, thermal_regions=thermal_regions)
    1504             : 
    1505          96 :       nregions = thermal_regions%nregions
    1506          96 :       nfree_done = 0
    1507          96 :       ekin_total_new = 0.0_dp
    1508          96 :       ekin_old = md_ener%ekin
    1509             :       vcom_total = 0.0_dp
    1510         480 :       ALLOCATE (temp_new(0:nregions), temp_old(0:nregions))
    1511         352 :       temp_new = 0.0_dp
    1512         352 :       temp_old = 0.0_dp
    1513             :       !loop regions
    1514         256 :       DO ireg = 1, nregions
    1515         160 :          NULLIFY (t_region)
    1516         160 :          t_region => thermal_regions%thermal_region(ireg)
    1517         160 :          nfree = 3*t_region%npart
    1518         160 :          ekin = compute_ekin(part, ireg)
    1519         160 :          IF (nfree > 0) t_region%temperature = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
    1520         160 :          temp_old(ireg) = t_region%temperature
    1521         160 :          IF (t_region%temp_tol > 0.0_dp .AND. &
    1522             :              ABS(t_region%temp_expected - t_region%temperature/kelvin) > t_region%temp_tol) THEN
    1523           2 :             fscale = SQRT((t_region%temp_expected/t_region%temperature)*kelvin)
    1524           2 :             CALL scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
    1525           2 :             t_region%temperature = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
    1526           2 :             temp_new(ireg) = t_region%temperature
    1527             :          END IF
    1528         160 :          nfree_done = nfree_done + nfree
    1529         256 :          ekin_total_new = ekin_total_new + ekin
    1530             :       END DO
    1531          96 :       nfree = simpar%nfree - nfree_done
    1532          96 :       ekin = compute_ekin(part, ireg=0)
    1533          96 :       IF (nfree > 0) thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
    1534          96 :       temp_old(0) = thermal_regions%temp_reg0
    1535          96 :       IF (simpar%temp_tol > 0.0_dp .AND. nfree > 0) THEN
    1536           0 :          IF (ABS(simpar%temp_ext - thermal_regions%temp_reg0/kelvin) > simpar%temp_tol) THEN
    1537           0 :             fscale = SQRT((simpar%temp_ext/thermal_regions%temp_reg0)*kelvin)
    1538           0 :             CALL scale_velocity_low(subsys, fscale, 0, ekin, vcom)
    1539           0 :             thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
    1540           0 :             temp_new(0) = thermal_regions%temp_reg0
    1541             :          END IF
    1542             :       END IF
    1543          96 :       ekin_total_new = ekin_total_new + ekin
    1544             : 
    1545          96 :       md_ener%ekin = ekin_total_new
    1546          96 :       IF (md_ener%nfree /= 0) THEN
    1547          96 :          md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
    1548             :       END IF
    1549          96 :       md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
    1550          96 :       IF (iw > 0) THEN
    1551         176 :          DO ireg = 0, nregions
    1552         176 :             IF (temp_new(ireg) > 0.0_dp) THEN
    1553             :                WRITE (UNIT=iw, FMT='(/,T2,A,I0,A)') &
    1554           1 :                   'MD_VEL| Temperature region ', ireg, ' scaled to requested temperature'
    1555             :                WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
    1556           1 :                   'MD_VEL| Old temperature [K]', temp_old(ireg), &
    1557           2 :                   'MD_VEL| New temperature [K]', temp_new(ireg)
    1558             :             END IF
    1559             :          END DO
    1560             :       END IF
    1561          96 :       DEALLOCATE (temp_new, temp_old)
    1562             : 
    1563          96 :    END SUBROUTINE scale_velocity_region
    1564             : 
    1565             : ! **************************************************************************************************
    1566             : !> \brief Scale velocities  for a specific region
    1567             : !> \param subsys ...
    1568             : !> \param fscale ...
    1569             : !> \param ireg ...
    1570             : !> \param ekin ...
    1571             : !> \param vcom ...
    1572             : !> \par author MI
    1573             : ! **************************************************************************************************
    1574        2564 :    SUBROUTINE scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
    1575             : 
    1576             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1577             :       REAL(KIND=dp), INTENT(IN)                          :: fscale
    1578             :       INTEGER, INTENT(IN)                                :: ireg
    1579             :       REAL(KIND=dp), INTENT(OUT)                         :: ekin, vcom(3)
    1580             : 
    1581             :       INTEGER                                            :: atom, iatom, ikind, my_ireg, natom, &
    1582             :                                                             shell_index
    1583        2564 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list
    1584             :       LOGICAL                                            :: is_shell
    1585             :       REAL(KIND=dp)                                      :: imass, mass, tmass, v2
    1586             :       REAL(KIND=dp), DIMENSION(3)                        :: tmp, v, vc, vs
    1587             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    1588             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1589             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
    1590             :                                                             shell_particles
    1591             :       TYPE(shell_kind_type), POINTER                     :: shell
    1592             : 
    1593        2564 :       NULLIFY (atomic_kinds, particles, shell_particles, core_particles, shell, atom_list)
    1594             : 
    1595        2564 :       my_ireg = ireg
    1596        2564 :       ekin = 0.0_dp
    1597        2564 :       tmass = 0.0_dp
    1598        2564 :       vcom = 0.0_dp
    1599             : 
    1600             :       CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
    1601        2564 :                          shell_particles=shell_particles, core_particles=core_particles)
    1602             : 
    1603        8848 :       DO ikind = 1, atomic_kinds%n_els
    1604        6284 :          atomic_kind => atomic_kinds%els(ikind)
    1605             :          CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, &
    1606        6284 :                               natom=natom, shell_active=is_shell, shell=shell)
    1607        6284 :          IF (is_shell) THEN
    1608         124 :             imass = 1.0_dp/mass
    1609         124 :             v2 = 0.0_dp
    1610         124 :             v = 0.0_dp
    1611        5740 :             DO iatom = 1, natom
    1612        5616 :                atom = atom_list(iatom)
    1613             :                !check region
    1614        5616 :                IF (particles%els(atom)%t_region_index /= my_ireg) CYCLE
    1615             : 
    1616       21080 :                particles%els(atom)%v(:) = fscale*particles%els(atom)%v
    1617        5270 :                shell_index = particles%els(atom)%shell_index
    1618       21080 :                vs = shell_particles%els(shell_index)%v
    1619       21080 :                vc = core_particles%els(shell_index)%v
    1620        5270 :                tmp(1) = imass*(vs(1) - vc(1))
    1621        5270 :                tmp(2) = imass*(vs(2) - vc(2))
    1622        5270 :                tmp(3) = imass*(vs(3) - vc(3))
    1623             : 
    1624        5270 :                shell_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) + tmp(1)*shell%mass_core
    1625        5270 :                shell_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) + tmp(2)*shell%mass_core
    1626        5270 :                shell_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) + tmp(3)*shell%mass_core
    1627             : 
    1628        5270 :                core_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) - tmp(1)*shell%mass_shell
    1629        5270 :                core_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) - tmp(2)*shell%mass_shell
    1630        5270 :                core_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) - tmp(3)*shell%mass_shell
    1631             : 
    1632             :                ! kinetic energy and velocity of COM
    1633       21080 :                v2 = v2 + SUM(particles%els(atom)%v**2)
    1634        5270 :                v(1) = v(1) + particles%els(atom)%v(1)
    1635        5270 :                v(2) = v(2) + particles%els(atom)%v(2)
    1636        5270 :                v(3) = v(3) + particles%els(atom)%v(3)
    1637        5740 :                tmass = tmass + mass
    1638             :             END DO
    1639             :          ELSE
    1640        6160 :             v2 = 0.0_dp
    1641        6160 :             v = 0.0_dp
    1642       22382 :             DO iatom = 1, natom
    1643       16222 :                atom = atom_list(iatom)
    1644             :                !check region
    1645       16222 :                IF (particles%els(atom)%t_region_index /= my_ireg) CYCLE
    1646             : 
    1647       64888 :                particles%els(atom)%v(:) = fscale*particles%els(atom)%v
    1648             :                ! kinetic energy and velocity of COM
    1649       64888 :                v2 = v2 + SUM(particles%els(atom)%v**2)
    1650       16222 :                v(1) = v(1) + particles%els(atom)%v(1)
    1651       16222 :                v(2) = v(2) + particles%els(atom)%v(2)
    1652       16222 :                v(3) = v(3) + particles%els(atom)%v(3)
    1653       22382 :                tmass = tmass + mass
    1654             :             END DO
    1655             :          END IF
    1656        6284 :          ekin = ekin + 0.5_dp*mass*v2
    1657        6284 :          vcom(1) = vcom(1) + mass*v(1)
    1658        6284 :          vcom(2) = vcom(2) + mass*v(2)
    1659       15132 :          vcom(3) = vcom(3) + mass*v(3)
    1660             : 
    1661             :       END DO
    1662       10256 :       vcom = vcom/tmass
    1663             : 
    1664        2564 :    END SUBROUTINE scale_velocity_low
    1665             : 
    1666             : ! **************************************************************************************************
    1667             : !> \brief Scale internal motion of CORE-SHELL model to the correct temperature
    1668             : !> \param subsys ...
    1669             : !> \param md_ener ...
    1670             : !> \param temp_expected ...
    1671             : !> \param temp_tol ...
    1672             : !> \param iw ...
    1673             : !> \par History
    1674             : !>     Teodoro Laino - University of Zurich - 09.2007 [tlaino]
    1675             : ! **************************************************************************************************
    1676        1060 :    SUBROUTINE scale_velocity_internal(subsys, md_ener, temp_expected, temp_tol, iw)
    1677             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1678             :       TYPE(md_ener_type), POINTER                        :: md_ener
    1679             :       REAL(KIND=dp), INTENT(IN)                          :: temp_expected, temp_tol
    1680             :       INTEGER, INTENT(IN)                                :: iw
    1681             : 
    1682             :       INTEGER                                            :: atom, iatom, ikind, natom, shell_index
    1683        1060 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list
    1684             :       LOGICAL                                            :: is_shell
    1685             :       REAL(KIND=dp)                                      :: ekin_shell_old, fac_mass, mass, scale, &
    1686             :                                                             temp_shell_old, v2
    1687             :       REAL(KIND=dp), DIMENSION(3)                        :: tmp, v, vc, vs
    1688             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    1689             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1690             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
    1691             :                                                             shell_particles
    1692             :       TYPE(shell_kind_type), POINTER                     :: shell
    1693             : 
    1694        1060 :       NULLIFY (atom_list, atomic_kinds, atomic_kind, core_particles, particles, shell_particles, shell)
    1695        1060 :       IF (ABS(temp_expected - md_ener%temp_shell/kelvin) > temp_tol) THEN
    1696          80 :          scale = 0.0_dp
    1697          80 :          IF (md_ener%temp_shell > EPSILON(0.0_dp)) scale = SQRT((temp_expected/md_ener%temp_shell)*kelvin)
    1698          80 :          ekin_shell_old = md_ener%ekin_shell
    1699          80 :          temp_shell_old = md_ener%temp_shell
    1700          80 :          md_ener%ekin_shell = 0.0_dp
    1701          80 :          md_ener%temp_shell = 0.0_dp
    1702             : 
    1703             :          CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, shell_particles=shell_particles, &
    1704          80 :                             core_particles=core_particles)
    1705             : 
    1706         240 :          DO ikind = 1, atomic_kinds%n_els
    1707         160 :             atomic_kind => atomic_kinds%els(ikind)
    1708             :             CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom, &
    1709         160 :                                  shell_active=is_shell, shell=shell)
    1710         240 :             IF (is_shell) THEN
    1711         160 :                fac_mass = 1.0_dp/mass
    1712         160 :                v2 = 0.0_dp
    1713         776 :                DO iatom = 1, natom
    1714         616 :                   atom = atom_list(iatom)
    1715         616 :                   shell_index = particles%els(atom)%shell_index
    1716        2464 :                   vs = shell_particles%els(shell_index)%v
    1717        2464 :                   vc = core_particles%els(shell_index)%v
    1718        2464 :                   v = particles%els(atom)%v
    1719         616 :                   tmp(1) = fac_mass*(vc(1) - vs(1))
    1720         616 :                   tmp(2) = fac_mass*(vc(2) - vs(2))
    1721         616 :                   tmp(3) = fac_mass*(vc(3) - vs(3))
    1722             : 
    1723         616 :                   shell_particles%els(shell_index)%v(1) = v(1) - shell%mass_core*scale*tmp(1)
    1724         616 :                   shell_particles%els(shell_index)%v(2) = v(2) - shell%mass_core*scale*tmp(2)
    1725         616 :                   shell_particles%els(shell_index)%v(3) = v(3) - shell%mass_core*scale*tmp(3)
    1726             : 
    1727         616 :                   core_particles%els(shell_index)%v(1) = v(1) + shell%mass_shell*scale*tmp(1)
    1728         616 :                   core_particles%els(shell_index)%v(2) = v(2) + shell%mass_shell*scale*tmp(2)
    1729         616 :                   core_particles%els(shell_index)%v(3) = v(3) + shell%mass_shell*scale*tmp(3)
    1730             : 
    1731        2464 :                   vs = shell_particles%els(shell_index)%v
    1732        2464 :                   vc = core_particles%els(shell_index)%v
    1733         616 :                   tmp(1) = vc(1) - vs(1)
    1734         616 :                   tmp(2) = vc(2) - vs(2)
    1735         616 :                   tmp(3) = vc(3) - vs(3)
    1736        2624 :                   v2 = v2 + SUM(tmp**2)
    1737             :                END DO
    1738         160 :                md_ener%ekin_shell = md_ener%ekin_shell + 0.5_dp*shell%mass_core*shell%mass_shell*fac_mass*v2
    1739             :             END IF
    1740             :          END DO
    1741          80 :          IF (md_ener%nfree_shell > 0) THEN
    1742          80 :             md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/REAL(md_ener%nfree_shell, KIND=dp)*kelvin
    1743             :          END IF
    1744          80 :          md_ener%constant = md_ener%constant - ekin_shell_old + md_ener%ekin_shell
    1745          80 :          IF (iw > 0) THEN
    1746             :             WRITE (UNIT=iw, FMT='(/,T2,A)') &
    1747          40 :                'MD_VEL| Temperature of shell internal motion scaled to requested temperature'
    1748             :             WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
    1749          40 :                'MD_VEL| Old temperature [K]', temp_shell_old, &
    1750          80 :                'MD_VEL| New temperature [K]', md_ener%temp_shell
    1751             :          END IF
    1752             :       END IF
    1753             : 
    1754        1060 :    END SUBROUTINE scale_velocity_internal
    1755             : 
    1756             : ! **************************************************************************************************
    1757             : !> \brief Scale barostat velocities to get the desired temperature
    1758             : !> \param md_env ...
    1759             : !> \param md_ener ...
    1760             : !> \param temp_expected ...
    1761             : !> \param temp_tol ...
    1762             : !> \param iw ...
    1763             : !> \par History
    1764             : !>     MI 02.2008
    1765             : ! **************************************************************************************************
    1766          40 :    SUBROUTINE scale_velocity_baro(md_env, md_ener, temp_expected, temp_tol, iw)
    1767             :       TYPE(md_environment_type), POINTER                 :: md_env
    1768             :       TYPE(md_ener_type), POINTER                        :: md_ener
    1769             :       REAL(KIND=dp), INTENT(IN)                          :: temp_expected, temp_tol
    1770             :       INTEGER, INTENT(IN)                                :: iw
    1771             : 
    1772             :       INTEGER                                            :: i, j, nfree
    1773             :       REAL(KIND=dp)                                      :: ekin_old, scale, temp_old
    1774          40 :       TYPE(npt_info_type), POINTER                       :: npt(:, :)
    1775             :       TYPE(simpar_type), POINTER                         :: simpar
    1776             : 
    1777          40 :       NULLIFY (npt, simpar)
    1778          40 :       CALL get_md_env(md_env, simpar=simpar, npt=npt)
    1779          40 :       IF (ABS(temp_expected - md_ener%temp_baro/kelvin) > temp_tol) THEN
    1780           2 :          scale = 0.0_dp
    1781           2 :          IF (md_ener%temp_baro > 0.0_dp) scale = SQRT((temp_expected/md_ener%temp_baro)*kelvin)
    1782           2 :          ekin_old = md_ener%baro_kin
    1783           2 :          temp_old = md_ener%temp_baro
    1784           2 :          md_ener%baro_kin = 0.0_dp
    1785           2 :          md_ener%temp_baro = 0.0_dp
    1786             :          IF (simpar%ensemble == npt_i_ensemble .OR. simpar%ensemble == npe_i_ensemble &
    1787           2 :              .OR. simpar%ensemble == npt_ia_ensemble) THEN
    1788           0 :             npt(1, 1)%v = npt(1, 1)%v*scale
    1789           0 :             md_ener%baro_kin = 0.5_dp*npt(1, 1)%v**2*npt(1, 1)%mass
    1790           2 :          ELSE IF (simpar%ensemble == npt_f_ensemble .OR. simpar%ensemble == npe_f_ensemble) THEN
    1791           2 :             md_ener%baro_kin = 0.0_dp
    1792           8 :             DO i = 1, 3
    1793          26 :                DO j = 1, 3
    1794          18 :                   npt(i, j)%v = npt(i, j)%v*scale
    1795          24 :                   md_ener%baro_kin = md_ener%baro_kin + 0.5_dp*npt(i, j)%v**2*npt(i, j)%mass
    1796             :                END DO
    1797             :             END DO
    1798             :          END IF
    1799           2 :          nfree = SIZE(npt, 1)*SIZE(npt, 2)
    1800           2 :          md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/REAL(nfree, dp)*kelvin
    1801           2 :          IF (iw > 0) THEN
    1802             :             WRITE (UNIT=iw, FMT='(/,T2,A)') &
    1803           1 :                'MD_VEL| Temperature of barostat motion scaled to requested temperature'
    1804             :             WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
    1805           1 :                'MD_VEL| Old temperature [K]', temp_old, &
    1806           2 :                'MD_VEL| New temperature [K]', md_ener%temp_baro
    1807             :          END IF
    1808             :       END IF
    1809             : 
    1810          40 :    END SUBROUTINE scale_velocity_baro
    1811             : 
    1812             : ! **************************************************************************************************
    1813             : !> \brief Perform all temperature manipulations during a QS MD run.
    1814             : !> \param simpar ...
    1815             : !> \param md_env ...
    1816             : !> \param md_ener ...
    1817             : !> \param force_env ...
    1818             : !> \param logger ...
    1819             : !> \par History
    1820             : !>     Creation (15.09.2003,MK)
    1821             : !>     adapted to force_env (05.10.2003,fawzi)
    1822             : !>     Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
    1823             : ! **************************************************************************************************
    1824       42781 :    SUBROUTINE temperature_control(simpar, md_env, md_ener, force_env, logger)
    1825             : 
    1826             :       TYPE(simpar_type), POINTER                         :: simpar
    1827             :       TYPE(md_environment_type), POINTER                 :: md_env
    1828             :       TYPE(md_ener_type), POINTER                        :: md_ener
    1829             :       TYPE(force_env_type), POINTER                      :: force_env
    1830             :       TYPE(cp_logger_type), POINTER                      :: logger
    1831             : 
    1832             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'temperature_control'
    1833             : 
    1834             :       INTEGER                                            :: handle, iw
    1835             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1836             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1837             : 
    1838       42781 :       CALL timeset(routineN, handle)
    1839       42781 :       NULLIFY (subsys, para_env)
    1840       42781 :       CPASSERT(ASSOCIATED(simpar))
    1841       42781 :       CPASSERT(ASSOCIATED(md_ener))
    1842       42781 :       CPASSERT(ASSOCIATED(force_env))
    1843       42781 :       CALL force_env_get(force_env, subsys=subsys, para_env=para_env)
    1844             :       iw = cp_print_key_unit_nr(logger, force_env%root_section, "MOTION%MD%PRINT%PROGRAM_RUN_INFO", &
    1845       42781 :                                 extension=".mdLog")
    1846             : 
    1847             :       ! Control the particle motion
    1848       42781 :       IF (simpar%do_thermal_region) THEN
    1849          96 :          CALL scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
    1850             :       ELSE
    1851       42685 :          IF (simpar%temp_tol > 0.0_dp) THEN
    1852       12770 :             CALL scale_velocity(subsys, md_ener, simpar%temp_ext, simpar%temp_tol, iw)
    1853             :          END IF
    1854             :       END IF
    1855             :       ! Control the internal core-shell motion
    1856       42781 :       IF (simpar%temp_sh_tol > 0.0_dp) THEN
    1857        1060 :          CALL scale_velocity_internal(subsys, md_ener, simpar%temp_sh_ext, simpar%temp_sh_tol, iw)
    1858             :       END IF
    1859             :       ! Control cell motion
    1860       45301 :       SELECT CASE (simpar%ensemble)
    1861             :       CASE (nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
    1862             :             npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble, npt_ia_ensemble)
    1863       42781 :          IF (simpar%temp_baro_tol > 0.0_dp) THEN
    1864          40 :             CALL scale_velocity_baro(md_env, md_ener, simpar%temp_baro_ext, simpar%temp_baro_tol, iw)
    1865             :          END IF
    1866             :       END SELECT
    1867             : 
    1868             :       CALL cp_print_key_finished_output(iw, logger, force_env%root_section, &
    1869       42781 :                                         "MOTION%MD%PRINT%PROGRAM_RUN_INFO")
    1870       42781 :       CALL timestop(handle)
    1871       42781 :    END SUBROUTINE temperature_control
    1872             : 
    1873             : ! **************************************************************************************************
    1874             : !> \brief Set to 0 the velocity of the COM along MD runs, if required.
    1875             : !> \param md_ener ...
    1876             : !> \param force_env ...
    1877             : !> \param md_section ...
    1878             : !> \param logger ...
    1879             : !> \par History
    1880             : !>      Creation (29.04.2007,MI)
    1881             : !>      Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
    1882             : ! **************************************************************************************************
    1883       85562 :    SUBROUTINE comvel_control(md_ener, force_env, md_section, logger)
    1884             : 
    1885             :       TYPE(md_ener_type), POINTER                        :: md_ener
    1886             :       TYPE(force_env_type), POINTER                      :: force_env
    1887             :       TYPE(section_vals_type), POINTER                   :: md_section
    1888             :       TYPE(cp_logger_type), POINTER                      :: logger
    1889             : 
    1890             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'comvel_control'
    1891             : 
    1892             :       INTEGER                                            :: handle, iw
    1893             :       LOGICAL                                            :: explicit
    1894             :       REAL(KIND=dp)                                      :: comvel_tol, temp_old, vel_com
    1895             :       REAL(KIND=dp), DIMENSION(3)                        :: vcom_old
    1896             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1897             : 
    1898       42781 :       CALL timeset(routineN, handle)
    1899       42781 :       NULLIFY (subsys)
    1900       42781 :       CPASSERT(ASSOCIATED(force_env))
    1901       42781 :       CALL force_env_get(force_env, subsys=subsys)
    1902             : 
    1903             :       ! Print COMVEL and COM Position
    1904       42781 :       iw = cp_print_key_unit_nr(logger, md_section, "PRINT%CENTER_OF_MASS", extension=".mdLog")
    1905       42781 :       IF (iw > 0) THEN
    1906             :          WRITE (UNIT=iw, FMT="(/,T2,A)") &
    1907        7742 :             "MD_VEL| Centre of mass motion (COM)"
    1908             :          WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
    1909        7742 :             "MD_VEL| VCOM [a.u.]", md_ener%vcom(1:3)
    1910             :       END IF
    1911       42781 :       CALL cp_print_key_finished_output(iw, logger, md_section, "PRINT%CENTER_OF_MASS")
    1912             : 
    1913             :       ! If requested rescale COMVEL
    1914       42781 :       CALL section_vals_val_get(md_section, "COMVEL_TOL", explicit=explicit)
    1915       42781 :       IF (explicit) THEN
    1916         826 :          CALL section_vals_val_get(md_section, "COMVEL_TOL", r_val=comvel_tol)
    1917             :          iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
    1918         826 :                                    extension=".mdLog")
    1919         826 :          vel_com = SQRT(md_ener%vcom(1)**2 + md_ener%vcom(2)**2 + md_ener%vcom(3)**2)
    1920             : 
    1921             :          ! Subtract the velocity of the COM, if requested
    1922         826 :          IF (vel_com > comvel_tol) THEN
    1923          42 :             temp_old = md_ener%temp_part/kelvin
    1924         168 :             vcom_old = md_ener%vcom
    1925          42 :             CALL reset_vcom(subsys, md_ener, vsubtract=vcom_old)
    1926          42 :             CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
    1927          42 :             IF (iw > 0) THEN
    1928             :                WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
    1929          21 :                   "MD_VEL| Old VCOM [a.u.]", vcom_old(1:3)
    1930             :                WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
    1931          21 :                   "MD_VEL| New VCOM [a.u.]", md_ener%vcom(1:3)
    1932             :             END IF
    1933             :          END IF
    1934             :          CALL cp_print_key_finished_output(iw, logger, md_section, &
    1935         826 :                                            "PRINT%PROGRAM_RUN_INFO")
    1936             :       END IF
    1937             : 
    1938       42781 :       CALL timestop(handle)
    1939       42781 :    END SUBROUTINE comvel_control
    1940             : 
    1941             : ! **************************************************************************************************
    1942             : !> \brief Set to 0 the angular velocity along MD runs, if required.
    1943             : !> \param md_ener ...
    1944             : !> \param force_env ...
    1945             : !> \param md_section ...
    1946             : !> \param logger ...
    1947             : !> \par History
    1948             : !>      Creation (10.2009) Teodoro Laino [tlaino]
    1949             : ! **************************************************************************************************
    1950       42781 :    SUBROUTINE angvel_control(md_ener, force_env, md_section, logger)
    1951             : 
    1952             :       TYPE(md_ener_type), POINTER                        :: md_ener
    1953             :       TYPE(force_env_type), POINTER                      :: force_env
    1954             :       TYPE(section_vals_type), POINTER                   :: md_section
    1955             :       TYPE(cp_logger_type), POINTER                      :: logger
    1956             : 
    1957             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'angvel_control'
    1958             : 
    1959             :       INTEGER                                            :: handle, ifixd, imolecule_kind, iw, natoms
    1960       42781 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: is_fixed
    1961             :       LOGICAL                                            :: explicit
    1962             :       REAL(KIND=dp)                                      :: angvel_tol, rcom(3), temp_old, vang(3), &
    1963             :                                                             vang_new(3)
    1964             :       TYPE(cell_type), POINTER                           :: cell
    1965             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    1966       42781 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
    1967             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    1968       42781 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
    1969             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
    1970             :       TYPE(particle_list_type), POINTER                  :: particles
    1971             : 
    1972       42781 :       CALL timeset(routineN, handle)
    1973             :       ! If requested rescale ANGVEL
    1974       42781 :       CALL section_vals_val_get(md_section, "ANGVEL_TOL", explicit=explicit)
    1975       42781 :       IF (explicit) THEN
    1976          40 :          NULLIFY (subsys, cell)
    1977          40 :          CPASSERT(ASSOCIATED(force_env))
    1978          40 :          CALL force_env_get(force_env, subsys=subsys, cell=cell)
    1979             : 
    1980         160 :          IF (SUM(cell%perd(1:3)) == 0) THEN
    1981          40 :             CALL section_vals_val_get(md_section, "ANGVEL_TOL", r_val=angvel_tol)
    1982             :             iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
    1983          40 :                                       extension=".mdLog")
    1984             : 
    1985             :             CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, &
    1986          40 :                                particles=particles)
    1987             : 
    1988          40 :             natoms = SIZE(particles%els)
    1989             :             ! Build a list of all fixed atoms (if any)
    1990         120 :             ALLOCATE (is_fixed(natoms))
    1991             : 
    1992         600 :             is_fixed = use_perd_none
    1993          40 :             molecule_kind_set => molecule_kinds%els
    1994         600 :             DO imolecule_kind = 1, molecule_kinds%n_els
    1995         560 :                molecule_kind => molecule_kind_set(imolecule_kind)
    1996         560 :                CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
    1997         600 :                IF (ASSOCIATED(fixd_list)) THEN
    1998           0 :                   DO ifixd = 1, SIZE(fixd_list)
    1999           0 :                      IF (.NOT. fixd_list(ifixd)%restraint%active) &
    2000           0 :                         is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
    2001             :                   END DO
    2002             :                END IF
    2003             :             END DO
    2004             : 
    2005             :             ! If requested and the system is not periodic, subtract the angular velocity
    2006          40 :             CALL compute_rcom(particles%els, is_fixed, rcom)
    2007          40 :             CALL compute_vang(particles%els, is_fixed, rcom, vang)
    2008             :             ! SQRT(DOT_PRODUCT(vang,vang))>angvel_tol
    2009         160 :             IF (DOT_PRODUCT(vang, vang) > (angvel_tol*angvel_tol)) THEN
    2010           2 :                CALL subtract_vang(particles%els, is_fixed, rcom, vang)
    2011             : 
    2012             :                ! Rescale velocities after removal
    2013           2 :                temp_old = md_ener%temp_part/kelvin
    2014           2 :                CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
    2015           2 :                CALL compute_vang(particles%els, is_fixed, rcom, vang_new)
    2016           2 :                IF (iw > 0) THEN
    2017             :                   WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
    2018           1 :                      'MD_VEL| Old VANG [a.u.]', vang(1:3)
    2019             :                   WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
    2020           1 :                      'MD_VEL| New VANG [a.u.]', vang_new(1:3)
    2021             :                END IF
    2022             :             END IF
    2023             : 
    2024          40 :             DEALLOCATE (is_fixed)
    2025             : 
    2026             :             CALL cp_print_key_finished_output(iw, logger, md_section, &
    2027          80 :                                               "PRINT%PROGRAM_RUN_INFO")
    2028             :          END IF
    2029             :       END IF
    2030             : 
    2031       42781 :       CALL timestop(handle)
    2032       85562 :    END SUBROUTINE angvel_control
    2033             : 
    2034             : ! **************************************************************************************************
    2035             : !> \brief Initialize Velocities for MD runs
    2036             : !> \param force_env ...
    2037             : !> \param simpar ...
    2038             : !> \param globenv ...
    2039             : !> \param md_env ...
    2040             : !> \param md_section ...
    2041             : !> \param constraint_section ...
    2042             : !> \param write_binary_restart_file ...
    2043             : !> \par History
    2044             : !>     Teodoro Laino - University of Zurich - 09.2007 [tlaino]
    2045             : ! **************************************************************************************************
    2046        3588 :    SUBROUTINE setup_velocities(force_env, simpar, globenv, md_env, md_section, &
    2047             :                                constraint_section, write_binary_restart_file)
    2048             : 
    2049             :       TYPE(force_env_type), POINTER                      :: force_env
    2050             :       TYPE(simpar_type), POINTER                         :: simpar
    2051             :       TYPE(global_environment_type), POINTER             :: globenv
    2052             :       TYPE(md_environment_type), POINTER                 :: md_env
    2053             :       TYPE(section_vals_type), POINTER                   :: md_section, constraint_section
    2054             :       LOGICAL, INTENT(IN)                                :: write_binary_restart_file
    2055             : 
    2056             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'setup_velocities'
    2057             : 
    2058             :       INTEGER                                            :: handle, nconstraint, nconstraint_fixd
    2059             :       LOGICAL                                            :: apply_cns0, shell_adiabatic, &
    2060             :                                                             shell_present
    2061             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
    2062             :       TYPE(cell_type), POINTER                           :: cell
    2063             :       TYPE(cp_subsys_type), POINTER                      :: subsys
    2064             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    2065             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2066             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
    2067             :                                                             shell_particles
    2068        1794 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
    2069        1794 :                                                             shell_particle_set
    2070             :       TYPE(section_vals_type), POINTER                   :: force_env_section, print_section, &
    2071             :                                                             subsys_section
    2072             : 
    2073        1794 :       CALL timeset(routineN, handle)
    2074             : 
    2075        1794 :       NULLIFY (atomic_kinds, cell, para_env, subsys, molecule_kinds, core_particles, particles)
    2076        1794 :       NULLIFY (shell_particles, core_particle_set, particle_set, shell_particle_set)
    2077        1794 :       NULLIFY (force_env_section, print_section, subsys_section)
    2078             : 
    2079        1794 :       print_section => section_vals_get_subs_vals(md_section, "PRINT")
    2080        1794 :       apply_cns0 = .FALSE.
    2081        1794 :       IF (simpar%constraint) THEN
    2082         310 :          CALL section_vals_val_get(constraint_section, "CONSTRAINT_INIT", l_val=apply_cns0)
    2083             :       END IF
    2084             :       ! Always initialize velocities and possibly restart them
    2085             :       CALL force_env_get(force_env, subsys=subsys, cell=cell, para_env=para_env, &
    2086        1794 :                          force_env_section=force_env_section)
    2087        1794 :       subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
    2088             : 
    2089             :       CALL cp_subsys_get(subsys, &
    2090             :                          atomic_kinds=atomic_kinds, &
    2091             :                          core_particles=core_particles, &
    2092             :                          molecule_kinds=molecule_kinds, &
    2093             :                          particles=particles, &
    2094        1794 :                          shell_particles=shell_particles)
    2095             : 
    2096             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, &
    2097             :                                shell_present=shell_present, &
    2098        1794 :                                shell_adiabatic=shell_adiabatic)
    2099             : 
    2100        1794 :       NULLIFY (core_particle_set)
    2101             :       NULLIFY (particle_set)
    2102        1794 :       NULLIFY (shell_particle_set)
    2103        1794 :       particle_set => particles%els
    2104             : 
    2105        1794 :       IF (shell_present .AND. shell_adiabatic) THEN
    2106             :          ! Constraints are not yet implemented for core-shell models generally
    2107             :          CALL get_molecule_kind_set(molecule_kind_set=molecule_kinds%els, &
    2108             :                                     nconstraint=nconstraint, &
    2109         132 :                                     nconstraint_fixd=nconstraint_fixd)
    2110         132 :          IF (nconstraint - nconstraint_fixd /= 0) &
    2111           0 :             CPABORT("Only the fixed atom constraint is implemented for core-shell models")
    2112             : !MK    CPPostcondition(.NOT.simpar%constraint,cp_failure_level,routineP,failure)
    2113         132 :          CPASSERT(ASSOCIATED(shell_particles))
    2114         132 :          CPASSERT(ASSOCIATED(core_particles))
    2115         132 :          shell_particle_set => shell_particles%els
    2116         132 :          core_particle_set => core_particles%els
    2117             :       END IF
    2118             : 
    2119             :       CALL initialize_velocities(simpar, &
    2120             :                                  particle_set, &
    2121             :                                  molecule_kinds=molecule_kinds, &
    2122             :                                  force_env=force_env, &
    2123             :                                  globenv=globenv, &
    2124             :                                  md_env=md_env, &
    2125             :                                  label="Velocities initialization", &
    2126             :                                  print_section=print_section, &
    2127             :                                  subsys_section=subsys_section, &
    2128             :                                  shell_present=(shell_present .AND. shell_adiabatic), &
    2129             :                                  shell_part=shell_particle_set, &
    2130             :                                  core_part=core_particle_set, &
    2131             :                                  force_rescaling=.FALSE., &
    2132             :                                  para_env=para_env, &
    2133        3456 :                                  write_binary_restart_file=write_binary_restart_file)
    2134             : 
    2135             :       ! Apply constraints if required and rescale velocities..
    2136        1794 :       IF (simpar%ensemble /= reftraj_ensemble) THEN
    2137        1756 :          IF (apply_cns0) THEN
    2138          24 :             CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
    2139             :             CALL force_env_shake(force_env, &
    2140             :                                  shake_tol=simpar%shake_tol, &
    2141             :                                  log_unit=simpar%info_constraint, &
    2142             :                                  lagrange_mult=simpar%lagrange_multipliers, &
    2143             :                                  dump_lm=simpar%dump_lm, &
    2144          24 :                                  compold=.TRUE.)
    2145             :             CALL force_env_rattle(force_env, shake_tol=simpar%shake_tol, &
    2146             :                                   log_unit=simpar%info_constraint, lagrange_mult=simpar%lagrange_multipliers, &
    2147          24 :                                   dump_lm=simpar%dump_lm, reset=.TRUE.)
    2148          24 :             IF (simpar%do_respa) THEN
    2149             :                CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
    2150           0 :                                                 calc_force=.TRUE.)
    2151             :                CALL force_env_shake(force_env%sub_force_env(1)%force_env, &
    2152             :                                     shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
    2153           0 :                                     lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, compold=.TRUE.)
    2154             :                CALL force_env_rattle(force_env%sub_force_env(1)%force_env, &
    2155             :                                      shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
    2156           0 :                                      lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, reset=.TRUE.)
    2157             :             END IF
    2158             :             ! Reinitialize velocities rescaling properly after rattle
    2159          24 :             subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
    2160          24 :             CALL update_subsys(subsys_section, force_env, .FALSE., write_binary_restart_file)
    2161             :             CALL initialize_velocities(simpar, &
    2162             :                                        particle_set, &
    2163             :                                        molecule_kinds=molecule_kinds, &
    2164             :                                        force_env=force_env, &
    2165             :                                        globenv=globenv, &
    2166             :                                        md_env=md_env, &
    2167             :                                        label="Re-Initializing velocities after applying constraints", &
    2168             :                                        print_section=print_section, &
    2169             :                                        subsys_section=subsys_section, &
    2170             :                                        shell_present=(shell_present .AND. shell_adiabatic), &
    2171             :                                        shell_part=shell_particle_set, &
    2172             :                                        core_part=core_particle_set, &
    2173             :                                        force_rescaling=.TRUE., &
    2174             :                                        para_env=para_env, &
    2175          48 :                                        write_binary_restart_file=write_binary_restart_file)
    2176             :          END IF
    2177             :       END IF
    2178             : 
    2179             :       ! Perform setup for a cascade run
    2180        1794 :       CALL initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
    2181             : 
    2182        1794 :       CALL timestop(handle)
    2183             : 
    2184        1794 :    END SUBROUTINE setup_velocities
    2185             : 
    2186             : ! **************************************************************************************************
    2187             : !> \brief   Perform the initialization for a cascade run
    2188             : !> \param simpar ...
    2189             : !> \param particle_set ...
    2190             : !> \param molecule_kinds ...
    2191             : !> \param md_section ...
    2192             : !> \date    05.02.2012
    2193             : !> \author  Matthias Krack (MK)
    2194             : !> \version 1.0
    2195             : ! **************************************************************************************************
    2196        1794 :    SUBROUTINE initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
    2197             : 
    2198             :       TYPE(simpar_type), POINTER                         :: simpar
    2199             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2200             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
    2201             :       TYPE(section_vals_type), POINTER                   :: md_section
    2202             : 
    2203             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_cascade'
    2204             : 
    2205             :       CHARACTER(len=2*default_string_length)             :: line
    2206             :       INTEGER                                            :: handle, iatom, ifixd, imolecule_kind, &
    2207             :                                                             iparticle, iw, natom, nparticle
    2208        1794 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_index, is_fixed
    2209             :       LOGICAL                                            :: init_cascade, is_ok, no_read_error
    2210             :       REAL(KIND=dp)                                      :: ecom, ekin, energy, norm, temp, &
    2211             :                                                             temperature
    2212        1794 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: matom, weight
    2213        1794 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: vatom
    2214             :       REAL(KIND=dp), DIMENSION(3)                        :: vcom
    2215             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    2216             :       TYPE(cp_logger_type), POINTER                      :: logger
    2217             :       TYPE(cp_sll_val_type), POINTER                     :: atom_list
    2218        1794 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
    2219        1794 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
    2220             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
    2221             :       TYPE(section_vals_type), POINTER                   :: atom_list_section, cascade_section, &
    2222             :                                                             print_section
    2223             :       TYPE(val_type), POINTER                            :: val
    2224             : 
    2225        1794 :       CALL timeset(routineN, handle)
    2226             : 
    2227        1794 :       NULLIFY (atom_list)
    2228        1794 :       NULLIFY (atom_list_section)
    2229        1794 :       NULLIFY (atomic_kind)
    2230        1794 :       NULLIFY (cascade_section)
    2231        1794 :       NULLIFY (fixd_list)
    2232        1794 :       NULLIFY (molecule_kind)
    2233        1794 :       NULLIFY (molecule_kind_set)
    2234        1794 :       NULLIFY (logger)
    2235        1794 :       NULLIFY (val)
    2236             : 
    2237        1794 :       logger => cp_get_default_logger()
    2238        1794 :       print_section => section_vals_get_subs_vals(md_section, "PRINT")
    2239        1794 :       iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
    2240             : 
    2241        1794 :       cascade_section => section_vals_get_subs_vals(md_section, "CASCADE")
    2242        1794 :       CALL section_vals_val_get(cascade_section, "_SECTION_PARAMETERS_", l_val=init_cascade)
    2243             : 
    2244        1794 :       nparticle = SIZE(particle_set)
    2245             : 
    2246        1794 :       IF (init_cascade) THEN
    2247             : 
    2248           2 :          CALL section_vals_val_get(cascade_section, "ENERGY", r_val=energy)
    2249           2 :          IF (energy < 0.0_dp) &
    2250           0 :             CPABORT("Error occurred reading &CASCADE section: Negative energy found")
    2251             : 
    2252           2 :          IF (iw > 0) THEN
    2253           1 :             ekin = cp_unit_from_cp2k(energy, "keV")
    2254             :             WRITE (UNIT=iw, FMT="(/,T2,A,T61,F20.6)") &
    2255           1 :                "CASCADE| Energy [keV]", ekin
    2256             :             WRITE (UNIT=iw, FMT="(T2,A)") &
    2257           1 :                "CASCADE|"
    2258             :          END IF
    2259             : 
    2260             :          ! Read the atomic velocities given in the input file
    2261           2 :          atom_list_section => section_vals_get_subs_vals(cascade_section, "ATOM_LIST")
    2262           2 :          CALL section_vals_val_get(atom_list_section, "_DEFAULT_KEYWORD_", n_rep_val=natom)
    2263           2 :          CALL section_vals_list_get(atom_list_section, "_DEFAULT_KEYWORD_", list=atom_list)
    2264           2 :          IF (natom <= 0) &
    2265           0 :             CPABORT("Error occurred reading &CASCADE section: No atom list found")
    2266             : 
    2267           2 :          IF (iw > 0) THEN
    2268             :             WRITE (UNIT=iw, FMT="(T2,A,T11,A,3(11X,A),9X,A)") &
    2269           1 :                "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "weight"
    2270             :          END IF
    2271             : 
    2272           6 :          ALLOCATE (atom_index(natom))
    2273           6 :          ALLOCATE (matom(natom))
    2274           6 :          ALLOCATE (vatom(3, natom))
    2275           6 :          ALLOCATE (weight(natom))
    2276             : 
    2277           8 :          DO iatom = 1, natom
    2278           6 :             is_ok = cp_sll_val_next(atom_list, val)
    2279           6 :             CALL val_get(val, c_val=line)
    2280             :             ! Read atomic index, velocity vector, and weight
    2281           6 :             no_read_error = .FALSE.
    2282           6 :             READ (UNIT=line, FMT=*, ERR=999) atom_index(iatom), vatom(1:3, iatom), weight(iatom)
    2283             :             no_read_error = .TRUE.
    2284             : 999         IF (.NOT. no_read_error) &
    2285           0 :                CPABORT("Error occurred reading &CASCADE section. Last line read <"//TRIM(line)//">")
    2286           6 :             IF ((atom_index(iatom) <= 0) .OR. ((atom_index(iatom) > nparticle))) &
    2287           0 :                CPABORT("Error occurred reading &CASCADE section: Invalid atom index found")
    2288           6 :             IF (weight(iatom) < 0.0_dp) &
    2289           0 :                CPABORT("Error occurred reading &CASCADE section: Negative weight found")
    2290           8 :             IF (iw > 0) THEN
    2291             :                WRITE (UNIT=iw, FMT="(T2,A,I10,4(1X,F14.6))") &
    2292           3 :                   "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), weight(iatom)
    2293             :             END IF
    2294             :          END DO
    2295             : 
    2296             :          ! Normalise velocities and weights
    2297             :          norm = 0.0_dp
    2298           8 :          DO iatom = 1, natom
    2299           6 :             iparticle = atom_index(iatom)
    2300           6 :             IF (particle_set(iparticle)%shell_index /= 0) &
    2301           0 :                CPWARN("Warning: The primary knock-on atom is a core-shell atom")
    2302           6 :             atomic_kind => particle_set(iparticle)%atomic_kind
    2303           6 :             CALL get_atomic_kind(atomic_kind=atomic_kind, mass=matom(iatom))
    2304           8 :             norm = norm + matom(iatom)*weight(iatom)
    2305             :          END DO
    2306           8 :          weight(:) = matom(:)*weight(:)*energy/norm
    2307           8 :          DO iatom = 1, natom
    2308          24 :             norm = SQRT(DOT_PRODUCT(vatom(1:3, iatom), vatom(1:3, iatom)))
    2309          26 :             vatom(1:3, iatom) = vatom(1:3, iatom)/norm
    2310             :          END DO
    2311             : 
    2312           2 :          IF (iw > 0) THEN
    2313             :             WRITE (UNIT=iw, FMT="(T2,A)") &
    2314           1 :                "CASCADE|", &
    2315           1 :                "CASCADE| Normalised velocities and additional kinetic energy [keV]", &
    2316           2 :                "CASCADE|"
    2317             :             WRITE (UNIT=iw, FMT="(T2,A,T11,A,3(11X,A),9X,A)") &
    2318           1 :                "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "E(kin)"
    2319           4 :             DO iatom = 1, natom
    2320           3 :                ekin = cp_unit_from_cp2k(weight(iatom), "keV")
    2321             :                WRITE (UNIT=iw, FMT="(T2,A,I10,4(1X,F14.6))") &
    2322           4 :                   "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), ekin
    2323             :             END DO
    2324             :          END IF
    2325             : 
    2326             :          ! Apply velocity modifications
    2327           8 :          DO iatom = 1, natom
    2328           6 :             iparticle = atom_index(iatom)
    2329             :             particle_set(iparticle)%v(:) = particle_set(iparticle)%v(:) + &
    2330          26 :                                            SQRT(2.0_dp*weight(iatom)/matom(iatom))*vatom(1:3, iatom)
    2331             :          END DO
    2332             : 
    2333           2 :          DEALLOCATE (atom_index)
    2334           2 :          DEALLOCATE (matom)
    2335           2 :          DEALLOCATE (vatom)
    2336           2 :          DEALLOCATE (weight)
    2337             : 
    2338           6 :          IF (iw > 0) THEN
    2339             :             ! Build a list of all fixed atoms (if any)
    2340           3 :             ALLOCATE (is_fixed(nparticle))
    2341          97 :             is_fixed = use_perd_none
    2342           1 :             molecule_kind_set => molecule_kinds%els
    2343           2 :             DO imolecule_kind = 1, molecule_kinds%n_els
    2344           1 :                molecule_kind => molecule_kind_set(imolecule_kind)
    2345           1 :                CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
    2346           2 :                IF (ASSOCIATED(fixd_list)) THEN
    2347           0 :                   DO ifixd = 1, SIZE(fixd_list)
    2348           0 :                      IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
    2349             :                   END DO
    2350             :                END IF
    2351             :             END DO
    2352             :             ! Compute vcom, ecom and ekin for printout
    2353           1 :             CALL compute_vcom(particle_set, is_fixed, vcom, ecom)
    2354           1 :             ekin = compute_ekin(particle_set) - ecom
    2355           1 :             IF (simpar%nfree == 0) THEN
    2356           0 :                CPASSERT(ekin == 0.0_dp)
    2357           0 :                temp = 0.0_dp
    2358             :             ELSE
    2359           1 :                temp = 2.0_dp*ekin/REAL(simpar%nfree, KIND=dp)
    2360             :             END IF
    2361           1 :             temperature = cp_unit_from_cp2k(temp, "K")
    2362             :             WRITE (UNIT=iw, FMT="(T2,A)") &
    2363           1 :                "CASCADE|"
    2364             :             WRITE (UNIT=iw, FMT="(T2,A,T61,F20.6)") &
    2365           1 :                "CASCADE| Temperature after cascade initialization [K]", temperature
    2366             :             WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,ES16.8))") &
    2367           1 :                "CASCADE| COM velocity", vcom(1:3)
    2368             : !MK          ! compute and log rcom and vang if not periodic
    2369             : !MK          CALL force_env_get(force_env,cell=cell)
    2370             : !MK          IF (SUM(cell%perd(1:3)) == 0) THEN
    2371             : !MK             CALL compute_rcom(particle_set,is_fixed,rcom)
    2372             : !MK             CALL compute_vang(particle_set,is_fixed,rcom,vang)
    2373             : !MK             WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM position:',rcom(1:3)
    2374             : !MK             WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' Angular velocity:',vang(1:3)
    2375             : !MK          END IF
    2376           2 :             DEALLOCATE (is_fixed)
    2377             :          END IF
    2378             : 
    2379             :       END IF
    2380             : 
    2381        1794 :       CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
    2382             : 
    2383        1794 :       CALL timestop(handle)
    2384             : 
    2385        3588 :    END SUBROUTINE initialize_cascade
    2386             : 
    2387             : END MODULE md_vel_utils

Generated by: LCOV version 1.15