TABLE OF CONTENTS
ABINIT/pred_hmc [ Functions ]
NAME
pred_hmc
FUNCTION
Hybrid Monte Carlo simulation algorithm. The routine generates a markov chain of structural configurations (states characterized by ionic positions and lattice parameters) with probability of observing a certian state equal to Gibbs statistical weight (exp(-etotal/kT)/Z).
COPYRIGHT
Copyright (C) 2017-2018 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 .
INPUTS
ab_mover = Data structure containing information about
input variables related to MD, e.g dtion, masses, etc.
hist = history of ionic positions, forces,
itime = index of current iteration
icycle = index of current cycle of the iteration
ntime = total number of iterations
ncycle = total number of cycles
zDEBUG = flag indicating whether to print debug info
iexit = flag indicating finilization of mover loop
OUTPUT
hist = ionic positions, lattice parameters etc. are updated
SIDE EFFECTS
NOTES
PARENTS
mover
CHILDREN
hist2var,pred_velverlet,var2hist,xred2xcart
SOURCE
45 #if defined HAVE_CONFIG_H 46 #include "config.h" 47 #endif 48 49 #include "abi_common.h" 50 51 subroutine pred_hmc(ab_mover,hist,itime,icycle,ntime,ncycle,zDEBUG,iexit) 52 53 use defs_basis 54 use m_errors 55 use m_profiling_abi 56 use m_abimover 57 use m_abihist 58 59 use m_geometry, only : xred2xcart 60 use m_numeric_tools, only : uniformrandom 61 62 !This section has been created automatically by the script Abilint (TD). 63 !Do not modify the following lines by hand. 64 #undef ABI_FUNC 65 #define ABI_FUNC 'pred_hmc' 66 use interfaces_45_geomoptim, except_this_one => pred_hmc 67 !End of the abilint section 68 69 implicit none 70 71 !Arguments ------------------------------------ 72 type(abimover),intent(in) :: ab_mover 73 type(abihist),intent(inout) :: hist 74 integer,intent(in) :: itime 75 integer,intent(in) :: icycle 76 integer,intent(in) :: ntime 77 integer,intent(in) :: ncycle 78 integer,intent(in) :: iexit 79 logical,intent(in) :: zDEBUG 80 81 !Local variables------------------------------- 82 83 integer :: seed ! seed for rnd generator 84 integer :: ii,jj,iacc ! dummy integers for loop indexes and acceptance decision flag 85 real(dp) :: etotal,epot,ekin,de ! total, potential (electronic), kinetic (ionic) energies and energy difference between initial and proposed states 86 real(dp) :: mv2tot,factor ! dummies used for rescaling of velocities 87 real(dp) :: rnd 88 real(dp) :: xred(3,ab_mover%natom) ! reduced coordinates of all ions 89 real(dp) :: vel(3,ab_mover%natom) ! ionic velocities in Cartesian coordinates 90 real(dp) :: mvtot(3) ! total momentum of the cell used to rescale velocities 91 real(dp) :: mtot,kbtemp ! total ionic mass and target temperature in energy units 92 real(dp) :: acell(3) ! lattice parameters 93 real(dp) :: rprimd(3,3) ! lattice vectors 94 95 real(dp),save :: etotal_hmc_prev ! total energy of the initial state 96 real(dp),save :: acell_hmc_prev(3) ! 97 real(dp),save :: rprimd_hmc_prev(3,3) ! 98 real(dp),allocatable,save :: xcart_hmc_prev(:,:) ! Cart. coordinates of the ions corresponding to the initial state 99 real(dp),allocatable,save :: xred_hmc_prev(:,:) ! reduced coordinates of the ions corresponding to the initial state 100 101 102 ! ************************************************************************* 103 104 DBG_ENTER("COLL") 105 106 ! if (option/=1 .and. option/=2 ) then 107 ! write(msg,'(3a,i0)')& 108 !& 'The argument option should be 1 or 2,',ch10,& 109 !& 'however, option=',option 110 ! MSG_BUG(msg) 111 ! end if 112 ! 113 ! if (sizein<1) then 114 ! write(msg,'(3a,i0)')& 115 !& 'The argument sizein should be a positive number,',ch10,& 116 !& 'however, sizein=',sizein 117 ! MSG_ERROR(msg) 118 ! end if 119 120 DBG_EXIT("COLL") 121 122 123 if(iexit/=0)then 124 if (allocated(xcart_hmc_prev)) then 125 ABI_DEALLOCATE(xcart_hmc_prev) 126 end if 127 if (allocated(xred_hmc_prev)) then 128 ABI_DEALLOCATE(xred_hmc_prev) 129 end if 130 return 131 end if 132 133 134 !get current values of ionic positions and cell geometry and set up the target temperature 135 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 136 137 vel(:,:) = hist%vel(:,:,hist%ihist) ! velocities of all ions, not needed in reality 138 epot = hist%etot(hist%ihist) ! electronic sub-system energy, not needed 139 ekin = hist%ekin(hist%ihist) ! kinetic energy, not needed 140 141 142 kbtemp=(ab_mover%mdtemp(1)*kb_HaK) ! get target temperature in energy units 143 144 if(icycle==1)then 145 146 !write(239117,*) ' entering first cycle of HMC iteration ',itime 147 148 if(itime==1)then 149 seed=-239 150 151 if (allocated(xcart_hmc_prev)) then 152 ABI_DEALLOCATE(xcart_hmc_prev) 153 end if 154 if (allocated(xred_hmc_prev)) then 155 ABI_DEALLOCATE(xred_hmc_prev) 156 end if 157 ABI_ALLOCATE(xcart_hmc_prev,(3,ab_mover%natom)) 158 ABI_ALLOCATE(xred_hmc_prev,(3,ab_mover%natom)) 159 end if 160 161 !generate new set of velocities and get rid of the possible overall momentum 162 mtot=sum(ab_mover%amass(:)) ! total mass to eventually get rid of total center of mass (CoM) momentum 163 164 ! generate velocities from normal distribution with zero mean and correct standard deviation 165 do ii=1,ab_mover%natom 166 do jj=1,3 167 vel(jj,ii)=sqrt(kbtemp/ab_mover%amass(ii))*cos(two_pi*uniformrandom(seed)) 168 vel(jj,ii)=vel(jj,ii)*sqrt(-2.0*log(uniformrandom(seed))) 169 end do 170 end do 171 !since number of atoms is most probably not big enough to obtain overall zero CoM momentum, shift the velocities 172 !and then renormalize 173 mvtot=zero ! total momentum 174 do ii=1,ab_mover%natom 175 do jj=1,3 176 mvtot(jj)=mvtot(jj)+ab_mover%amass(ii)*vel(jj,ii) 177 end do 178 end do 179 do ii=1,ab_mover%natom 180 do jj=1,3 181 vel(jj,ii)=vel(jj,ii)-(mvtot(jj)/mtot) 182 end do 183 end do 184 !now the total cell momentum is zero 185 mv2tot=0.0 186 do ii=1,ab_mover%natom 187 do jj=1,3 188 mv2tot=mv2tot+ab_mover%amass(ii)*vel(jj,ii)**2 189 end do 190 end do 191 factor = mv2tot/(dble(3*ab_mover%natom)) 192 factor = sqrt(kbtemp/factor) 193 vel(:,:)=vel(:,:)*factor 194 hist%vel(:,:,hist%ihist)=vel(:,:) 195 196 197 198 !save the starting values of ionic positions (before an update attempt is performed) 199 call hist2var(acell_hmc_prev,hist,ab_mover%natom,rprimd_hmc_prev,xred_hmc_prev,zDEBUG) 200 call xred2xcart(ab_mover%natom,rprimd_hmc_prev,xcart_hmc_prev,xred_hmc_prev) 201 202 !also save the initial total energy 203 ekin=0.0 204 do ii=1,ab_mover%natom 205 do jj=1,3 206 ekin=ekin+half*ab_mover%amass(ii)*vel(jj,ii)**2 207 end do 208 end do 209 epot = hist%etot(hist%ihist) ! electronic sub-system energy 210 etotal_hmc_prev = epot + ekin ! total energy before an attempt of variables update is performed 211 212 213 call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) 214 215 !call hist2var(acell,hist,ab_mover%natom,rprim,rprimd,xcart,xred,zDEBUG) 216 !write(std_out,*) ' xcart_after_update:' 217 !write(std_out,*) ' ',xcart(1,:) 218 !write(std_out,*) ' ',xcart(2,:) 219 !write(std_out,*) ' ',xcart(3,:) 220 221 222 223 else if(icycle<ncycle)then 224 225 call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) 226 !call hist2var(acell,hist,ab_mover%natom,rprim,rprimd,xcart,xred,zDEBUG) 227 !vel(:,:) = hist%vel(:,:,hist%ihist) ! velocities of all ions, not needed in reality 228 229 else !icycle==ncycle 230 231 !the only thing left to do, is to compute the difference of the total energies and decide whether to accept the new state 232 vel(:,:)=hist%vel(:,:,hist%ihist) 233 ekin=0.0 234 do ii=1,ab_mover%natom 235 do jj=1,3 236 ekin=ekin+half*ab_mover%amass(ii)*vel(jj,ii)**2 237 end do 238 end do 239 epot = hist%etot(hist%ihist) ! electronic sub-system energy 240 etotal = epot + ekin 241 de = etotal - etotal_hmc_prev 242 243 iacc=0 244 rnd=uniformrandom(seed) 245 if(de<0)then 246 iacc=1 247 else 248 if(exp(-de/kbtemp)>rnd)then 249 iacc=1 250 end if 251 end if 252 !write(238,*) ' random number: ',rnd,' -de/kbtemp:',-de/kbtemp,' acceptance decision: ',iacc 253 !write(239,*) ' de: ',de,' estart: ',etotal_hmc_prev,' efin:', etotal 254 !write(239,*) 'ekin= ',ekin,' epot= ',epot,' e_start= ',etotal_hmc_prev,' e_end= ',etotal,' iacc=',iacc, ' dekT=',-de/kbtemp 255 256 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 257 !call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 258 259 if(iacc==0)then !in case the new state is not accepted, then roll back the coordinates 260 !write(std_out,*) ' the proposed state was not accepted, roll back the configuration' 261 !xcart(:,:)=xcart_hmc_prev(:,:) 262 !no need to roll back xcart any longer, everything is stored in xred now (v8.3.1) 263 xred(:,:)=xred_hmc_prev(:,:) 264 end if 265 266 hist%ihist=abihist_findIndex(hist,+1) 267 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 268 269 end if 270 271 ! note: add update of lattice vectors and parameters in case optcell/=0 272 273 end subroutine pred_hmc