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 Module with utility to perform MD Nudged Elastic Band Calculation
      10             : !> \note
      11             : !>      Numerical accuracy for parallel runs:
      12             : !>       Each replica starts the SCF run from the one optimized
      13             : !>       in a previous run. It may happen then energies and derivatives
      14             : !>       of a serial run and a parallel run could be slightly different
      15             : !>       'cause of a different starting density matrix.
      16             : !>       Exact results are obtained using:
      17             : !>          EXTRAPOLATION USE_GUESS in QS section (Teo 09.2006)
      18             : !> \author Teodoro Laino 10.2006
      19             : ! **************************************************************************************************
      20             : MODULE neb_md_utils
      21             :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      22             :    USE global_types,                    ONLY: global_environment_type
      23             :    USE input_constants,                 ONLY: band_md_opt,&
      24             :                                               do_band_collective
      25             :    USE input_section_types,             ONLY: section_vals_get,&
      26             :                                               section_vals_get_subs_vals,&
      27             :                                               section_vals_type,&
      28             :                                               section_vals_val_get
      29             :    USE kinds,                           ONLY: dp
      30             :    USE neb_types,                       ONLY: neb_type,&
      31             :                                               neb_var_type
      32             :    USE particle_types,                  ONLY: get_particle_pos_or_vel,&
      33             :                                               particle_type,&
      34             :                                               update_particle_pos_or_vel
      35             :    USE physcon,                         ONLY: kelvin,&
      36             :                                               massunit
      37             : #include "../base/base_uses.f90"
      38             : 
      39             :    IMPLICIT NONE
      40             :    PRIVATE
      41             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'neb_md_utils'
      42             : 
      43             :    PUBLIC :: neb_initialize_velocity, &
      44             :              control_vels_a, &
      45             :              control_vels_b, &
      46             :              get_temperatures
      47             : 
      48             : CONTAINS
      49             : 
      50             : ! **************************************************************************************************
      51             : !> \brief Initialize velocities of replica in an MD optimized algorithm within
      52             : !>        NEB
      53             : !> \param vels ...
      54             : !> \param neb_section ...
      55             : !> \param particle_set ...
      56             : !> \param i_rep ...
      57             : !> \param iw ...
      58             : !> \param globenv ...
      59             : !> \param neb_env ...
      60             : !> \par History
      61             : !>      25.11.2010 Consider core-shell model (MK)
      62             : !> \author Teodoro Laino 09.2006
      63             : ! **************************************************************************************************
      64         184 :    SUBROUTINE neb_initialize_velocity(vels, neb_section, particle_set, i_rep, iw, &
      65             :                                       globenv, neb_env)
      66             : 
      67             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vels
      68             :       TYPE(section_vals_type), POINTER                   :: neb_section
      69             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      70             :       INTEGER, INTENT(IN)                                :: i_rep, iw
      71             :       TYPE(global_environment_type), POINTER             :: globenv
      72             :       TYPE(neb_type), POINTER                            :: neb_env
      73             : 
      74             :       INTEGER                                            :: iatom, ivar, natom, nparticle, nvar
      75             :       REAL(KIND=dp)                                      :: akin, mass, mass_tot, sc, temp, &
      76             :                                                             temp_ext, tmp_r1
      77             :       REAL(KIND=dp), DIMENSION(3)                        :: v, vcom
      78             :       TYPE(section_vals_type), POINTER                   :: md_section
      79             : 
      80         184 :       IF (neb_env%opt_type == band_md_opt) THEN
      81          70 :          md_section => section_vals_get_subs_vals(neb_section, "OPTIMIZE_BAND%MD")
      82          70 :          CALL section_vals_val_get(md_section, "TEMPERATURE", r_val=temp_ext)
      83             :          ! Initialize velocity according to external temperature
      84          70 :          nparticle = SIZE(vels, 1)
      85          70 :          natom = SIZE(particle_set)
      86          70 :          mass_tot = 0.0_dp
      87          70 :          vcom(1:3) = 0.0_dp
      88          70 :          CALL globenv%gaussian_rng_stream%fill(vels(:, i_rep))
      89             :          ! Check always if BAND is working in Cartesian or in internal coordinates
      90             :          ! If working in cartesian coordinates let's get rid of the COM
      91             :          ! Compute also the total mass (both in Cartesian and internal)
      92          70 :          IF (neb_env%use_colvar) THEN
      93          10 :             nvar = nparticle
      94          10 :             mass_tot = REAL(nvar, KIND=dp)*massunit
      95             :          ELSE
      96        1080 :             DO iatom = 1, natom
      97        1020 :                mass = particle_set(iatom)%atomic_kind%mass
      98        1020 :                mass_tot = mass_tot + mass
      99        1020 :                v(1:3) = get_particle_pos_or_vel(iatom, particle_set, vels(:, i_rep))
     100        4140 :                vcom(1:3) = vcom(1:3) + mass*v(1:3)
     101             :             END DO
     102         240 :             vcom(1:3) = vcom(1:3)/mass_tot
     103             :          END IF
     104             :          ! Compute kinetic energy and temperature
     105          70 :          akin = 0.0_dp
     106          70 :          IF (neb_env%use_colvar) THEN
     107          20 :             DO ivar = 1, nvar
     108          20 :                akin = akin + 0.5_dp*massunit*vels(ivar, i_rep)*vels(ivar, i_rep)
     109             :             END DO
     110             :          ELSE
     111        1080 :             DO iatom = 1, natom
     112        1020 :                mass = particle_set(iatom)%atomic_kind%mass
     113        4080 :                v(1:3) = -vcom(1:3)
     114        1020 :                CALL update_particle_pos_or_vel(iatom, particle_set, v(1:3), vels(:, i_rep))
     115        4140 :                akin = akin + 0.5_dp*mass*DOT_PRODUCT(v(1:3), v(1:3))
     116             :             END DO
     117          60 :             nvar = 3*natom
     118             :          END IF
     119          70 :          temp = 2.0_dp*akin/REAL(nvar, KIND=dp)
     120             :          ! Scale velocities to get the correct initial temperature and
     121          70 :          sc = SQRT(temp_ext/temp)
     122        3140 :          vels(:, i_rep) = vels(:, i_rep)*sc
     123             :          ! Re-compute kinetic energya and temperature and velocity of COM
     124          70 :          akin = 0.0_dp
     125          70 :          vcom = 0.0_dp
     126          70 :          IF (neb_env%use_colvar) THEN
     127          20 :             DO ivar = 1, nvar
     128          20 :                akin = akin + 0.5_dp*massunit*vels(ivar, i_rep)*vels(ivar, i_rep)
     129             :             END DO
     130             :          ELSE
     131        1080 :             DO iatom = 1, natom
     132        1020 :                mass = particle_set(iatom)%atomic_kind%mass
     133        1020 :                v(1:3) = get_particle_pos_or_vel(iatom, particle_set, vels(:, i_rep))
     134        4080 :                vcom(1:3) = vcom(1:3) + mass*v(1:3)
     135        4140 :                akin = akin + 0.5_dp*mass*DOT_PRODUCT(v(1:3), v(1:3))
     136             :             END DO
     137             :          END IF
     138         280 :          vcom(1:3) = vcom(1:3)/mass_tot
     139             :          ! Dump information
     140          70 :          IF (iw > 0) THEN
     141           5 :             temp = 2.0_dp*akin/REAL(nvar, KIND=dp)
     142           5 :             tmp_r1 = cp_unit_from_cp2k(temp, "K")
     143             :             WRITE (iw, '(A,T61,F18.2,A2)') &
     144           5 :                ' NEB| Initial Temperature ', tmp_r1, " K"
     145             :             WRITE (iw, '(A,T61,F20.12)') &
     146           5 :                ' NEB| Centre of mass velocity in direction x:', vcom(1), &
     147           5 :                ' NEB| Centre of mass velocity in direction y:', vcom(2), &
     148          10 :                ' NEB| Centre of mass velocity in direction z:', vcom(3)
     149           5 :             WRITE (iw, '(T2,"NEB|",75("*"))')
     150             :          END IF
     151             :       ELSE
     152       32490 :          vels(:, i_rep) = 0.0_dp
     153             :       END IF
     154             : 
     155         184 :    END SUBROUTINE neb_initialize_velocity
     156             : 
     157             : ! **************************************************************************************************
     158             : !> \brief Control on  velocities - I part
     159             : !> \param vels ...
     160             : !> \param particle_set ...
     161             : !> \param tc_section ...
     162             : !> \param vc_section ...
     163             : !> \param output_unit ...
     164             : !> \param istep ...
     165             : !> \author Teodoro Laino 09.2006
     166             : ! **************************************************************************************************
     167         160 :    SUBROUTINE control_vels_a(vels, particle_set, tc_section, vc_section, &
     168             :                              output_unit, istep)
     169             :       TYPE(neb_var_type), POINTER                        :: vels
     170             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     171             :       TYPE(section_vals_type), POINTER                   :: tc_section, vc_section
     172             :       INTEGER, INTENT(IN)                                :: output_unit, istep
     173             : 
     174             :       INTEGER                                            :: i, temp_tol_steps
     175             :       LOGICAL                                            :: explicit
     176             :       REAL(KIND=dp)                                      :: ext_temp, f_annealing, scale, temp_tol, &
     177             :                                                             temploc, tmp_r1, tmp_r2
     178         160 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: temperatures
     179             : 
     180             : ! Temperature control
     181             : 
     182         160 :       CALL section_vals_get(tc_section, explicit=explicit)
     183         160 :       IF (explicit) THEN
     184          60 :          CALL section_vals_val_get(tc_section, "TEMP_TOL_STEPS", i_val=temp_tol_steps)
     185          60 :          CALL section_vals_val_get(tc_section, "TEMPERATURE", r_val=ext_temp)
     186          60 :          CALL section_vals_val_get(tc_section, "TEMP_TOL", r_val=temp_tol)
     187         180 :          ALLOCATE (temperatures(SIZE(vels%wrk, 2)))
     188             :          ! Computes temperatures
     189          60 :          CALL get_temperatures(vels, particle_set, temperatures, factor=1.0_dp)
     190             :          ! Possibly rescale
     191          60 :          IF (istep <= temp_tol_steps) THEN
     192         260 :             DO i = 2, SIZE(vels%wrk, 2) - 1
     193         220 :                temploc = temperatures(i)
     194         260 :                IF (ABS(temploc - ext_temp) > temp_tol) THEN
     195         158 :                   IF (output_unit > 0) THEN
     196          79 :                      tmp_r1 = cp_unit_from_cp2k(temploc, "K")
     197          79 :                      tmp_r2 = cp_unit_from_cp2k(ext_temp, "K")
     198             :                      WRITE (output_unit, '(T2,"NEB| Replica Nr.",I5,'// &
     199             :                             '"  - Velocity rescaled from: ",F12.6," to: ",F12.6,".")') &
     200          79 :                         i, tmp_r1, tmp_r2
     201             : 
     202             :                   END IF
     203         158 :                   scale = SQRT(ext_temp/temploc)
     204        8216 :                   vels%wrk(:, i) = scale*vels%wrk(:, i)
     205             :                END IF
     206             :             END DO
     207             :          END IF
     208         120 :          DEALLOCATE (temperatures)
     209             :       END IF
     210             :       ! Annealing
     211         160 :       CALL section_vals_get(vc_section, explicit=explicit)
     212         160 :       IF (explicit) THEN
     213         160 :          CALL section_vals_val_get(vc_section, "ANNEALING", r_val=f_annealing)
     214         740 :          DO i = 2, SIZE(vels%wrk, 2) - 1
     215       27320 :             vels%wrk(:, i) = f_annealing*vels%wrk(:, i)
     216             :          END DO
     217             :       END IF
     218         160 :    END SUBROUTINE control_vels_a
     219             : 
     220             : ! **************************************************************************************************
     221             : !> \brief Control on velocities - II part
     222             : !> \param vels ...
     223             : !> \param forces ...
     224             : !> \param vc_section ...
     225             : !> \author Teodoro Laino 09.2006
     226             : ! **************************************************************************************************
     227         160 :    SUBROUTINE control_vels_b(vels, forces, vc_section)
     228             :       TYPE(neb_var_type), POINTER                        :: vels, forces
     229             :       TYPE(section_vals_type), POINTER                   :: vc_section
     230             : 
     231             :       INTEGER                                            :: i
     232             :       LOGICAL                                            :: explicit, lval
     233             :       REAL(KIND=dp)                                      :: factor, norm
     234             : 
     235             : ! Check the sign of V.dot.F
     236             : 
     237         160 :       CALL section_vals_get(vc_section, explicit=explicit)
     238         160 :       IF (explicit) THEN
     239         160 :          CALL section_vals_val_get(vc_section, "PROJ_VELOCITY_VERLET", l_val=lval)
     240         160 :          IF (lval) THEN
     241         560 :             DO i = 2, SIZE(vels%wrk, 2) - 1
     242       18840 :                norm = DOT_PRODUCT(forces%wrk(:, i), forces%wrk(:, i))
     243       18840 :                factor = DOT_PRODUCT(vels%wrk(:, i), forces%wrk(:, i))
     244         560 :                IF (factor > 0 .AND. (norm >= EPSILON(0.0_dp))) THEN
     245       36608 :                   vels%wrk(:, i) = factor/norm*forces%wrk(:, i)
     246             :                ELSE
     247         536 :                   vels%wrk(:, i) = 0.0_dp
     248             :                END IF
     249             :             END DO
     250             :          END IF
     251         160 :          CALL section_vals_val_get(vc_section, "SD_LIKE", l_val=lval)
     252         160 :          IF (lval) THEN
     253           0 :             DO i = 2, SIZE(vels%wrk, 2) - 1
     254           0 :                vels%wrk(:, i) = 0.0_dp
     255             :             END DO
     256             :          END IF
     257             :       END IF
     258         160 :    END SUBROUTINE control_vels_b
     259             : 
     260             : ! **************************************************************************************************
     261             : !> \brief Computes temperatures
     262             : !> \param vels ...
     263             : !> \param particle_set ...
     264             : !> \param temperatures ...
     265             : !> \param ekin ...
     266             : !> \param factor ...
     267             : !> \par History
     268             : !>      24.11.2010 rewritten to include core-shell model (MK)
     269             : !> \author Teodoro Laino 09.2006
     270             : ! **************************************************************************************************
     271         349 :    SUBROUTINE get_temperatures(vels, particle_set, temperatures, ekin, factor)
     272             : 
     273             :       TYPE(neb_var_type), POINTER                        :: vels
     274             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     275             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: temperatures
     276             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: ekin
     277             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: factor
     278             : 
     279             :       INTEGER                                            :: i_rep, iatom, ivar, n_rep, natom, nvar
     280             :       REAL(KIND=dp)                                      :: akin, mass, myfactor, temploc
     281             :       REAL(KIND=dp), DIMENSION(3)                        :: v
     282             : 
     283         349 :       myfactor = kelvin
     284         349 :       IF (PRESENT(factor)) myfactor = factor
     285        2425 :       IF (PRESENT(ekin)) ekin(:) = 0.0_dp
     286        2536 :       temperatures(:) = 0.0_dp
     287         349 :       nvar = SIZE(vels%wrk, 1)
     288         349 :       n_rep = SIZE(vels%wrk, 2)
     289         349 :       natom = SIZE(particle_set)
     290        1838 :       DO i_rep = 2, n_rep - 1
     291        1489 :          akin = 0.0_dp
     292        1489 :          IF (vels%in_use == do_band_collective) THEN
     293          72 :             DO ivar = 1, nvar
     294          72 :                akin = akin + 0.5_dp*massunit*vels%wrk(ivar, i_rep)*vels%wrk(ivar, i_rep)
     295             :             END DO
     296             :          ELSE
     297       44832 :             DO iatom = 1, natom
     298       43379 :                mass = particle_set(iatom)%atomic_kind%mass
     299       43379 :                v(1:3) = get_particle_pos_or_vel(iatom, particle_set, vels%wrk(:, i_rep))
     300      174969 :                akin = akin + 0.5_dp*mass*DOT_PRODUCT(v(1:3), v(1:3))
     301             :             END DO
     302        1453 :             nvar = 3*natom
     303             :          END IF
     304        1489 :          IF (PRESENT(ekin)) ekin(i_rep) = akin
     305        1489 :          temploc = 2.0_dp*akin/REAL(nvar, KIND=dp)
     306        1838 :          temperatures(i_rep) = temploc*myfactor
     307             :       END DO
     308             : 
     309         349 :    END SUBROUTINE get_temperatures
     310             : 
     311             : END MODULE neb_md_utils