TABLE OF CONTENTS
- ABINIT/m_hmc
- ABINIT/m_hmc/compute_kinetic_energy
- ABINIT/m_hmc/generate_random_velocities
- ABINIT/m_hmc/metropolis_check
ABINIT/m_hmc [ Modules ]
NAME
m_hmc
FUNCTION
Auxiliary hmc functions
COPYRIGHT
Copyright (C) 2018-2022 ABINIT group (SPr) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
SOURCE
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
ABINIT/m_hmc/compute_kinetic_energy [ Functions ]
NAME
comute_kintic_energy
FUNCTION
Computes kintetic energy
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
64 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,jj 73 !character(len=500) :: msg 74 75 ! ************************************************************************* 76 77 ekin=0.0 78 do ii=1,ab_mover%natom 79 do jj=1,3 80 ekin=ekin+half*ab_mover%amass(ii)*vel(jj,ii)**2 81 end do 82 end do 83 84 end subroutine compute_kinetic_energy
ABINIT/m_hmc/generate_random_velocities [ Functions ]
NAME
generate_random_velocities
FUNCTION
Generate normally distributed random velocities
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
109 subroutine generate_random_velocities(ab_mover,kbtemp,seed,vel,ekin) 110 111 !Arguments ------------------------------------ 112 type(abimover),intent(in) :: ab_mover 113 integer, intent(inout):: seed 114 real(dp), intent(in) :: kbtemp 115 real(dp), intent(inout):: vel(3,ab_mover%natom) ! velocities 116 real(dp), intent(out) :: ekin ! output kinetic energy 117 !Local variables------------------------------- 118 integer :: ii,jj 119 real(dp):: mtot,mvtot(3),mv2tot,factor 120 !character(len=500) :: msg 121 122 ! ************************************************************************* 123 124 125 mtot=sum(ab_mover%amass(:)) ! total mass to eventually get rid of total center of mass (CoM) momentum 126 !generate velocities from normal distribution with zero mean and correct standard deviation 127 do ii=1,ab_mover%natom 128 do jj=1,3 129 vel(jj,ii)=sqrt(kbtemp/ab_mover%amass(ii))*cos(two_pi*uniformrandom(seed)) 130 vel(jj,ii)=vel(jj,ii)*sqrt(-2.0*log(uniformrandom(seed))) 131 end do 132 end do 133 !since number of atoms is most probably not big enough to obtain overall zero CoM momentum, shift the velocities 134 !and then renormalize 135 mvtot=zero ! total momentum 136 do ii=1,ab_mover%natom 137 do jj=1,3 138 mvtot(jj)=mvtot(jj)+ab_mover%amass(ii)*vel(jj,ii) 139 end do 140 end do 141 do ii=1,ab_mover%natom 142 do jj=1,3 143 vel(jj,ii)=vel(jj,ii)-(mvtot(jj)/mtot) 144 end do 145 end do 146 !now the total cell momentum is zero 147 mv2tot=0.0 148 do ii=1,ab_mover%natom 149 do jj=1,3 150 mv2tot=mv2tot+ab_mover%amass(ii)*vel(jj,ii)**2 151 end do 152 end do 153 factor = mv2tot/(dble(3*ab_mover%natom)) 154 factor = sqrt(kbtemp/factor) 155 vel(:,:)=vel(:,:)*factor 156 157 call compute_kinetic_energy(ab_mover,vel,ekin) 158 159 end subroutine generate_random_velocities
ABINIT/m_hmc/metropolis_check [ Functions ]
NAME
metropolis_check
FUNCTION
Make an acceptance decision based on the energy differences
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
182 subroutine metropolis_check(seed,de,kbtemp,iacc) 183 184 !Arguments ------------------------------------ 185 integer, intent(inout):: seed 186 real(dp), intent(in) :: de 187 real(dp), intent(in) :: kbtemp 188 integer, intent(inout):: iacc 189 190 !Local variables------------------------------- 191 real(dp) :: rnd 192 !character(len=500) :: msg 193 194 ! ************************************************************************* 195 196 iacc=0 197 rnd=uniformrandom(seed) 198 if(de<0)then 199 iacc=1 200 else 201 if(exp(-de/kbtemp)>rnd)then 202 iacc=1 203 end if 204 end if 205 206 end subroutine metropolis_check