Branch data Line data Source code
1 : : !!****m* ABINIT/m_hmc
2 : : !! NAME
3 : : !! m_hmc
4 : : !!
5 : : !! FUNCTION
6 : : !! Auxiliary hmc functions
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2018-2026 ABINIT group (SPr)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! NOTES
15 : : !!
16 : : !! SOURCE
17 : :
18 : : #if defined HAVE_CONFIG_H
19 : : #include "config.h"
20 : : #endif
21 : :
22 : : #include "abi_common.h"
23 : :
24 : : module m_hmc
25 : :
26 : : use defs_basis
27 : : use m_abicore
28 : : use m_errors
29 : : use m_abimover
30 : : use m_io_tools
31 : :
32 : : !use m_geometry, only : xred2xcart
33 : : use m_numeric_tools, only : uniformrandom
34 : :
35 : : implicit none
36 : :
37 : : private
38 : :
39 : : ! *************************************************************************
40 : : public :: compute_kinetic_energy
41 : : public :: generate_random_velocities
42 : : public :: metropolis_check
43 : :
44 : : contains
45 : : !!***
46 : :
47 : : !!****f* ABINIT/m_hmc/compute_kinetic_energy
48 : : !! NAME
49 : : !! comute_kintic_energy
50 : : !!
51 : : !! FUNCTION
52 : : !! Computes kintetic energy
53 : : !!
54 : : !! INPUTS
55 : : !!
56 : : !! OUTPUT
57 : : !!
58 : : !! SIDE EFFECTS
59 : : !!
60 : : !! NOTES
61 : : !!
62 : : !! SOURCE
63 : :
64 : 0 : subroutine compute_kinetic_energy(ab_mover,vel,ekin)
65 : :
66 : : !Arguments ------------------------------------
67 : : type(abimover),intent(in) :: ab_mover
68 : : real(dp), intent(in) :: vel(3,ab_mover%natom) ! velocities
69 : : real(dp), intent(out) :: ekin ! output kinetic energy
70 : :
71 : : !Local variables-------------------------------
72 : : integer :: ii
73 : : !character(len=500) :: msg
74 : :
75 : : ! *************************************************************************
76 : :
77 : 0 : ekin=0.0
78 [ # # ]: 0 : do ii = 1, ab_mover%natom
79 [ # # ]: 0 : ekin = ekin + half * ab_mover%amass(ii) * DOT_PRODUCT(vel(:, ii), vel(:, ii))
80 : : end do
81 : :
82 : 0 : end subroutine compute_kinetic_energy
83 : : !!***
84 : :
85 : :
86 : :
87 : :
88 : :
89 : :
90 : : !!****f* ABINIT/m_hmc/generate_random_velocities
91 : : !! NAME
92 : : !! generate_random_velocities
93 : : !!
94 : : !! FUNCTION
95 : : !! Generate normally distributed random velocities
96 : : !!
97 : : !! INPUTS
98 : : !!
99 : : !! OUTPUT
100 : : !!
101 : : !! SIDE EFFECTS
102 : : !!
103 : : !! NOTES
104 : : !!
105 : : !! SOURCE
106 : :
107 : 0 : subroutine generate_random_velocities(ab_mover,kbtemp,seed,vel,ekin)
108 : :
109 : : !Arguments ------------------------------------
110 : : type(abimover),intent(in) :: ab_mover
111 : : integer, intent(inout):: seed
112 : : real(dp), intent(in) :: kbtemp
113 : : real(dp), intent(inout):: vel(3,ab_mover%natom) ! velocities
114 : : real(dp), intent(out) :: ekin ! output kinetic energy
115 : : !Local variables-------------------------------
116 : : integer :: ii,jj,natom
117 : : real(dp):: mtot,mvtot(3),mv2tot,factor
118 : : !character(len=500) :: msg
119 : :
120 : : ! *************************************************************************
121 : :
122 : :
123 : 0 : natom = ab_mover%natom
124 [ # # ]: 0 : mtot=sum(ab_mover%amass(:)) ! total mass to eventually get rid of total center of mass (CoM) momentum
125 : : !generate velocities from normal distribution with zero mean and correct standard deviation
126 [ # # ]: 0 : do ii=1,ab_mover%natom
127 [ # # ]: 0 : do jj=1,3
128 : 0 : vel(jj,ii)=sqrt(kbtemp/ab_mover%amass(ii))*cos(two_pi*uniformrandom(seed))
129 : 0 : vel(jj,ii)=vel(jj,ii)*sqrt(-2.0*log(uniformrandom(seed)))
130 : : end do
131 : : end do
132 : : !since number of atoms is most probably not big enough to obtain overall zero CoM momentum, shift the velocities
133 : : !and then renormalize
134 : : ! mvtot -> total momentum
135 [ # # # # ]: 0 : mvtot(:) = MATMUL(vel(1:3,1:natom), ab_mover%amass(1:natom))
136 [ # # ]: 0 : do ii=1,ab_mover%natom
137 [ # # ]: 0 : vel(:,ii)=vel(1:3,ii)-(mvtot(1:3)/mtot)
138 : : end do
139 : : !now the total cell momentum is zero
140 : : mv2tot=0.0
141 [ # # ]: 0 : do ii=1,ab_mover%natom
142 [ # # ]: 0 : mv2tot=mv2tot+ab_mover%amass(ii)* DOT_PRODUCT(vel(:,ii), vel(:,ii))
143 : : end do
144 : 0 : factor = mv2tot/(dble(3*ab_mover%natom))
145 : 0 : factor = sqrt(kbtemp/factor)
146 [ # # # # ]: 0 : vel(:,:)=vel(:,:)*factor
147 : :
148 : 0 : call compute_kinetic_energy(ab_mover,vel,ekin)
149 : :
150 : 0 : end subroutine generate_random_velocities
151 : : !!***
152 : :
153 : :
154 : :
155 : :
156 : : !!****f* ABINIT/m_hmc/metropolis_check
157 : : !! NAME
158 : : !! metropolis_check
159 : : !!
160 : : !! FUNCTION
161 : : !! Make an acceptance decision based on the energy differences
162 : : !!
163 : : !! INPUTS
164 : : !!
165 : : !! OUTPUT
166 : : !!
167 : : !! SIDE EFFECTS
168 : : !!
169 : : !! NOTES
170 : : !!
171 : : !! SOURCE
172 : :
173 : 202 : subroutine metropolis_check(seed,de,kbtemp,iacc)
174 : :
175 : : !Arguments ------------------------------------
176 : : integer, intent(inout):: seed
177 : : real(dp), intent(in) :: de
178 : : real(dp), intent(in) :: kbtemp
179 : : integer, intent(inout):: iacc
180 : :
181 : : !Local variables-------------------------------
182 : : real(dp) :: rnd
183 : : !character(len=500) :: msg
184 : :
185 : : ! *************************************************************************
186 : :
187 : 202 : iacc=0
188 : 202 : rnd=uniformrandom(seed)
189 [ + - ]: 202 : if(de<0)then
190 : 202 : iacc=1
191 : : else
192 [ # # ]: 0 : if(exp(-de/kbtemp)>rnd)then
193 : 0 : iacc=1
194 : : end if
195 : : end if
196 : :
197 : 202 : end subroutine metropolis_check
198 : : !!***
199 : :
200 : :
201 : :
202 : :
203 : : end module m_hmc
204 : : !!***
|
