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 : !!***
|
