Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> 15.10.2007 Giovanni Bussi - Implementation validated.
11 : !> \author Teodoro Laino - 09.2007 - University of Zurich [tlaino]
12 : ! **************************************************************************************************
13 : MODULE csvr_system_utils
14 :
15 : USE kinds, ONLY: dp
16 : USE parallel_rng_types, ONLY: rng_stream_type
17 : #include "./base/base_uses.f90"
18 :
19 : IMPLICIT NONE
20 :
21 : PRIVATE
22 :
23 : LOGICAL, PARAMETER :: debug_this_module = .FALSE.
24 : PUBLIC :: rescaling_factor
25 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'csvr_system_utils'
26 :
27 : CONTAINS
28 :
29 : ! **************************************************************************************************
30 : !> \brief Stochastic velocity rescale, as described in
31 : !> Bussi, Donadio and Parrinello, J. Chem. Phys. 126, 014101 (2007)
32 : !>
33 : !> This subroutine implements Eq.(A7) and returns the new value for the kinetic energy,
34 : !> which can be used to rescale the velocities.
35 : !> The procedure can be applied to all atoms or to smaller groups.
36 : !> If it is applied to intersecting groups in sequence, the kinetic energy
37 : !> that is given as an input (kk) has to be up-to-date with respect to the previous
38 : !> rescalings.
39 : !>
40 : !> When applied to the entire system, and when performing standard molecular dynamics
41 : !> (fixed c.o.m. (center of mass))
42 : !> the degrees of freedom of the c.o.m. have to be discarded in the calculation of ndeg,
43 : !> and the c.o.m. momentum HAS TO BE SET TO ZERO.
44 : !> When applied to subgroups, one can chose to:
45 : !> (a) calculate the subgroup kinetic energy in the usual reference frame, and count
46 : !> the c.o.m. in ndeg
47 : !> (b) calculate the subgroup kinetic energy with respect to its c.o.m. motion, discard
48 : !> the c.o.m. in ndeg and apply the rescale factor with respect to the subgroup c.o.m.
49 : !> velocity.
50 : !> They should be almost equivalent.
51 : !> If the subgroups are expected to move one respect to the other, the choice (b)
52 : !> should be better.
53 : !>
54 : !> If a null relaxation time is required (taut=0.0), the procedure reduces to an istantaneous
55 : !> randomization of the kinetic energy, as described in paragraph IIA.
56 : !>
57 : !> HOW TO CALCULATE THE EFFECTIVE-ENERGY DRIFT
58 : !> The effective-energy (htilde) drift can be used to check the integrator against
59 : !> discretization errors.
60 : !> The easiest recipe is:
61 : !> htilde = h + conint
62 : !> where h is the total energy (kinetic + potential)
63 : !> and conint is a quantity accumulated along the trajectory as minus the sum of all
64 : !> the increments of kinetic energy due to the thermostat.
65 : !>
66 : !> Variables:
67 : !> kk ! present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
68 : !> sigma ! target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
69 : !> ndeg ! number of degrees of freedom of the atoms to be thermalized
70 : !> taut ! relaxation time of the thermostat, in units of 'how often this routine is called'
71 : !> \param kk ...
72 : !> \param sigma ...
73 : !> \param ndeg ...
74 : !> \param taut ...
75 : !> \param rng_stream ...
76 : !> \return ...
77 : !> \date 09.2007
78 : !> \author Giovanni Bussi - ETH Zurich, Lugano 10.2007
79 : ! **************************************************************************************************
80 114684 : FUNCTION rescaling_factor(kk, sigma, ndeg, taut, rng_stream) RESULT(my_res)
81 : REAL(KIND=dp), INTENT(IN) :: kk, sigma
82 : INTEGER, INTENT(IN) :: ndeg
83 : REAL(KIND=dp), INTENT(IN) :: taut
84 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
85 : REAL(KIND=dp) :: my_res
86 :
87 : REAL(KIND=dp) :: factor, resample, reverse, rr
88 :
89 114684 : my_res = 0.0_dp
90 114684 : IF (kk > 0.0_dp) THEN
91 114586 : IF (taut > 0.1_dp) THEN
92 114586 : factor = EXP(-1.0_dp/taut)
93 : ELSE
94 : factor = 0.0_dp
95 : END IF
96 114586 : rr = rng_stream%next()
97 114586 : reverse = 1.0_dp
98 : ! reverse of momentum is implemented to have the correct limit to Langevin dynamics for ndeg=1
99 : ! condition: rr < -SQRT(ndeg*kk*factor/(sigma*(1.0_dp-factor)))
100 114586 : IF ((rr*rr*sigma*(1.0_dp - factor)) > (ndeg*kk*factor) .AND. rr <= 0.0_dp) reverse = -1.0_dp
101 : ! for ndeg/=1, the reverse of momentum is not necessary. in principles, it should be there.
102 : ! in practice, it is better to skip it to avoid unnecessary slowing down of the dynamics in the small taut regime
103 : ! anyway, this should not affect the final ensemble
104 114586 : IF (ndeg /= 1) reverse = 1.0_dp
105 : resample = kk + (1.0_dp - factor)*(sigma*(sumnoises(ndeg - 1, rng_stream) + rr**2)/REAL(ndeg, KIND=dp) - kk) &
106 114586 : + 2.0_dp*rr*SQRT(kk*sigma/ndeg*(1.0_dp - factor)*factor)
107 :
108 114586 : resample = MAX(0.0_dp, resample)
109 114586 : my_res = reverse*SQRT(resample/kk)
110 : END IF
111 114684 : END FUNCTION rescaling_factor
112 :
113 : ! **************************************************************************************************
114 : !> \brief returns the sum of n independent gaussian noises squared
115 : !> (i.e. equivalent to summing the square of the return values of nn calls to gasdev)
116 : !> \param nn ...
117 : !> \param rng_stream ...
118 : !> \return ...
119 : !> \date 09.2007
120 : !> \author Teo - University of Zurich
121 : ! **************************************************************************************************
122 114586 : FUNCTION sumnoises(nn, rng_stream) RESULT(sum_gauss)
123 : INTEGER, INTENT(IN) :: nn
124 : TYPE(rng_stream_type), INTENT(INOUT) :: rng_stream
125 : REAL(KIND=dp) :: sum_gauss
126 :
127 : INTEGER :: i
128 :
129 114586 : sum_gauss = 0.0_dp
130 3923352 : DO i = 1, nn
131 3923352 : sum_gauss = sum_gauss + rng_stream%next()**2
132 : END DO
133 :
134 114586 : END FUNCTION sumnoises
135 :
136 : END MODULE csvr_system_utils
|
