TABLE OF CONTENTS
ABINIT/m_pred_hmc [ Modules ]
NAME
m_pred_hmc
FUNCTION
COPYRIGHT
Copyright (C) 2017-2024 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 .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_pred_hmc 22 23 implicit none 24 25 private
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).
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
SOURCE
64 subroutine pred_hmc(ab_mover,hist,itime,icycle,ntime,ncycle,mttk_vars,zDEBUG,iexit) 65 66 use defs_basis 67 use m_errors 68 use m_abicore 69 use m_abimover 70 use m_abihist 71 use m_io_tools 72 use m_hmc 73 74 use m_geometry, only : xred2xcart 75 use m_numeric_tools, only : uniformrandom 76 use m_pred_velverlet, only : pred_velverlet 77 use m_pred_isothermal, only : pred_isothermal 78 implicit none 79 80 !Arguments ------------------------------------ 81 type(abimover),intent(in) :: ab_mover 82 type(abihist),intent(inout) :: hist 83 type(mttk_type),intent(inout) :: mttk_vars 84 integer,intent(in) :: itime 85 integer,intent(in) :: icycle 86 integer,intent(in) :: ntime 87 integer,intent(in) :: ncycle 88 integer,intent(in) :: iexit 89 logical,intent(in) :: zDEBUG 90 91 !Local variables------------------------------- 92 integer,save :: seed ! seed for rnd generator 93 integer :: iacc ! dummy integers for loop indexes and acceptance decision flag 94 real(dp) :: etotal,epot,ekin,de ! total, potential (electronic), kinetic (ionic) energies and energy difference 95 !real(dp) :: mv2tot,factor ! dummies used for rescaling of velocities 96 real(dp) :: xred(3,ab_mover%natom) ! reduced coordinates of all ions 97 real(dp) :: vel(3,ab_mover%natom) ! ionic velocities in Cartesian coordinates 98 !real(dp) :: mvtot(3) ! total momentum of the cell used to rescale velocities 99 real(dp) :: kbtemp !mtot, ! total ionic mass and target temperature in energy units 100 real(dp) :: acell(3) ! lattice parameters 101 real(dp) :: rprimd(3,3) ! lattice vectors 102 103 real(dp),save :: etotal_hmc_prev,epot_hmc_prev ! total energy of the initial state 104 real(dp),save :: strain(3,3),dstrain ! strain tensor 105 real(dp),save :: rprimd_original(3,3) ! initial lattice vectors <= itime=1,icycle=1 106 real(dp),allocatable,save :: xred_hmc_prev(:,:) ! reduced coordinates of the ions corresponding to the initial state 107 real(dp),allocatable,save :: fcart_hmc_prev(:,:) ! reduced coordinates of the ions corresponding to the initial state 108 109 logical,save :: strain_updated 110 logical,save :: xred_updated 111 integer,save :: strain_steps 112 logical :: strain_sweep 113 114 ! character(len=500) :: message 115 ! ************************************************************************* 116 117 DBG_ENTER("COLL") 118 119 ! if (option/=1 .and. option/=2 ) then 120 ! write(msg,'(3a,i0)')& 121 !& 'The argument option should be 1 or 2,',ch10,& 122 !& 'however, option=',option 123 ! ABI_BUG(msg) 124 ! end if 125 ! 126 ! if (sizein<1) then 127 ! write(msg,'(3a,i0)')& 128 !& 'The argument sizein should be a positive number,',ch10,& 129 !& 'however, sizein=',sizein 130 ! ABI_ERROR(msg) 131 ! end if 132 133 DBG_EXIT("COLL") 134 135 strain_sweep=.FALSE. 136 137 if(iexit/=0)then !icycle=ncycle and itime=ntime 138 if (allocated(xred_hmc_prev)) then 139 ABI_FREE(xred_hmc_prev) 140 end if 141 if (allocated(fcart_hmc_prev)) then 142 ABI_FREE(fcart_hmc_prev) 143 end if 144 !call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) ! this is needed to deallocate vel_prev array allocated in pred_velverlet 145 call pred_isothermal(ab_mover,hist,icycle,mttk_vars,ncycle,zDEBUG,iexit) 146 return 147 end if 148 149 150 !get current values of ionic positions and cell geometry and set up the target temperature 151 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 152 153 vel(:,:) = hist%vel(:,:,hist%ihist) ! velocities of all ions, not needed in reality 154 epot = hist%etot(hist%ihist) ! electronic sub-system energy, not needed 155 ekin = hist%ekin(hist%ihist) ! kinetic energy, not needed 156 157 kbtemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK ! correct temperature taking into account the possible heating/cooling 158 159 if(itime==1.and.icycle==1) then 160 if (allocated(xred_hmc_prev)) then 161 ABI_FREE(xred_hmc_prev) 162 end if 163 if (allocated(fcart_hmc_prev)) then 164 ABI_FREE(fcart_hmc_prev) 165 end if 166 167 ABI_MALLOC(xred_hmc_prev,(3,ab_mover%natom)) 168 ABI_MALLOC(fcart_hmc_prev,(3,ab_mover%natom)) 169 170 seed=-239 171 172 rprimd_original(:,:)=rprimd(:,:) 173 strain(:,:) = 0.0_dp 174 strain_steps=0 175 dstrain=0.001 176 177 strain_updated=.FALSE. 178 xred_updated=.FALSE. 179 end if 180 181 182 !IN CASE THE SWEEP IS FOR UPDATE OF ATOMIC COORDINATES************************************************ 183 !if(.NOT.strain_sweep) then 184 185 ! *---->* 186 ! 1 n 187 188 if (icycle==1) then 189 190 if(itime==1) then 191 iacc=1 192 etotal = epot + ekin 193 de=zero 194 else 195 etotal = epot + ekin 196 de = etotal - etotal_hmc_prev 197 call metropolis_check(seed,de,kbtemp,iacc) 198 !DEBUG 199 ! write(std_out,*)' m_pred_hmc, after metropolis_check : seed,de,kbtemp,iacc=',seed,de,kbtemp,iacc 200 !ENDDEBUG 201 end if 202 203 if(iacc==0)then !in case the new state is not accepted, then roll back the coordinates and energies 204 xred(:,:)= xred_hmc_prev(:,:) 205 epot = epot_hmc_prev 206 hist%fcart(:,:,hist%ihist) = fcart_hmc_prev(:,:) 207 else 208 xred_hmc_prev(:,:)=xred(:,:) 209 fcart_hmc_prev(:,:) = hist%fcart(:,:,hist%ihist) 210 epot_hmc_prev = epot !update reference potential energy 211 end if 212 213 ! write(message,'(2a,i7,a,i2,a,E24.16,a,E24.16,a,E24.16)') ch10,' HMC Sweep => ',itime,' iacc= ', iacc,' epot= ',& 214 !& epot,' ekin=',ekin,' de=',de 215 ! call wrtout(ab_out,message,'COLL') 216 ! call wrtout(std_out,message,'COLL') 217 218 !call generate_random_velocities(ab_mover,kbtemp,seed,vel,ekin) ! this routine also computes the new kinetic energy 219 !hist%vel(:,:,hist%ihist)=vel(:,:) 220 hist%vel(:,:,hist%ihist)=0 221 !call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 222 !etotal_hmc_prev=epot+ekin ! either old or current potential energy + new kinetic energy 223 224 !call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) ! 1 is indicating that velverlet is called from hmc routine 225 call pred_isothermal(ab_mover,hist,icycle,mttk_vars,ncycle,zDEBUG,iexit) 226 227 elseif(icycle > 1 .and. icycle <= ncycle)then !icycle/=1 228 229 !call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) ! 1 is indicating that velverlet is called from hmc routine 230 call pred_isothermal(ab_mover,hist,icycle,mttk_vars,ncycle,zDEBUG,iexit) 231 232 !end if 233 !END OF ATOMIC COORDINATES SWEEP************************************************ 234 235 ! else if(icycle>ncycle) then ! strain update 236 ! strain_updated = .TRUE. 237 ! strain_steps = strain_steps + 1 238 !! Metropolis update of lattice vectors and parameters in case optcell/=0 239 ! if(icycle==ncycle+1.and.xred_updated) then 240 ! !save rprimd_hmc_prev and total electronic energy etotal_hmc_prev 241 ! call hist2var(acell_hmc_prev,hist,ab_mover%natom,rprimd_hmc_prev,xred,zDEBUG) 242 ! etotal_hmc_prev = hist%etot(hist%ihist) 243 ! strain_hmc_prev(:,:) = strain(:,:) 244 ! 245 ! select case (ab_mover%optcell) 246 ! case (1) !volume optimization only 247 ! acell(:)=acell(:)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp)) 248 ! case (2,3,7,8,9) !full geometry optimization 249 ! !suggest new strain tensor values 250 ! do ii=1,3 251 ! do jj=ii,3 252 ! strain(ii,jj) = strain(ii,jj)+ 2.0_dp*dstrain*(uniformrandom(seed)-0.5_dp) 253 ! strain(jj,ii) = strain(ii,jj) 254 ! enddo 255 ! enddo 256 ! if(ab_mover%optcell==3) then !eliminate volume change if optcell==3 257 ! do ii=1,3 258 ! strain(ii,ii) = strain(ii,ii) -(strain(1,1)+strain(2,2)+strain(3,3)) 259 ! enddo 260 ! endif 261 ! do jj=1,3 ! sum over three lattice vectors 262 ! do ii=1,3 ! sum over Cart components 263 ! rprimd(ii,jj)=rprimd_original(ii,jj)+& 264 !& rprimd_original(1,jj)*strain(ii,1)+& 265 !& rprimd_original(2,jj)*strain(ii,2)+& 266 !& rprimd_original(3,jj)*strain(ii,3) 267 ! enddo 268 ! enddo 269 ! if(ab_mover%optcell==7) then 270 ! rprimd(:,1)=rprimd_original(:,1) 271 ! else if (ab_mover%optcell==8) then 272 ! rprimd(:,2)=rprimd_original(:,2) 273 ! else if (ab_mover%optcell==9) then 274 ! rprimd(:,3)=rprimd_original(:,3) 275 ! endif 276 ! case(4) 277 ! acell(1)=acell(1)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp)) 278 ! case(5) 279 ! acell(2)=acell(2)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp)) 280 ! case(6) 281 ! acell(3)=acell(3)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp)) 282 ! !case default 283 ! ! write(message,"(a,i0)") "Wrong value of optcell: ",ab_mover%optcell 284 ! ! ABI_ERROR(message) 285 ! end select 286 ! 287 ! !update the new suggested rprimd and or acell in the history record 288 ! hist%ihist=abihist_findIndex(hist,+1) 289 ! call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 290 ! else 291 ! 292 ! etotal = hist%etot(hist%ihist) 293 ! de = etotal - etotal_hmc_prev 294 ! 295 ! iacc=0 296 ! rnd=uniformrandom(seed) 297 ! if(de<0)then 298 ! iacc=1 299 ! else 300 ! if(exp(-de/kbtemp)>rnd)then 301 ! iacc=1 302 ! end if 303 ! end if 304 ! 305 ! if(iacc==0) then 306 ! strain(:,:)=strain_hmc_prev(:,:) 307 ! acell(:)=acell_hmc_prev(:) 308 ! else 309 ! call hist2var(acell_hmc_prev,hist,ab_mover%natom,rprimd_hmc_prev,xred,zDEBUG) 310 ! strain_hmc_prev(:,:) = strain(:,:) 311 ! etotal_hmc_prev=etotal 312 ! endif 313 ! 314 ! !suggest new acell/rprimd values depending on the optcell value 315 ! select case (ab_mover%optcell) 316 ! case (1) !volume optimization only 317 ! acell(:)=acell(:)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp)) 318 ! case (2,3,7,8,9) !full geometry optimization 319 ! !suggest new strain tensor values 320 ! do ii=1,3 321 ! do jj=ii,3 322 ! strain(ii,jj) = strain(ii,jj)+ 2.0_dp*dstrain*(uniformrandom(seed)-0.5_dp) 323 ! strain(jj,ii) = strain(ii,jj) 324 ! enddo 325 ! enddo 326 ! if(ab_mover%optcell==3) then !eliminate volume change if optcell==3 327 ! do ii=1,3 328 ! strain(ii,ii) = strain(ii,ii) -(strain(1,1)+strain(2,2)+strain(3,3)) 329 ! enddo 330 ! endif 331 ! do jj=1,3 ! sum over three lattice vectors 332 ! do ii=1,3 ! sum over Cart components 333 ! rprimd(ii,jj)=rprimd_original(ii,jj)+& 334 !& rprimd_original(1,jj)*strain(ii,1)+& 335 !& rprimd_original(2,jj)*strain(ii,2)+& 336 !& rprimd_original(3,jj)*strain(ii,3) 337 ! enddo 338 ! enddo 339 ! if(ab_mover%optcell==7) then 340 ! rprimd(:,1)=rprimd_original(:,1) 341 ! else if (ab_mover%optcell==8) then 342 ! rprimd(:,2)=rprimd_original(:,2) 343 ! else if (ab_mover%optcell==9) then 344 ! rprimd(:,3)=rprimd_original(:,3) 345 ! endif 346 ! case(4) 347 ! acell(1)=acell(1)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp)) 348 ! case(5) 349 ! acell(2)=acell(2)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp)) 350 ! case(6) 351 ! acell(3)=acell(3)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp)) 352 ! case default 353 ! ! write(message,"(a,i0)") "Wrong value of optcell: ",ab_mover%optcell 354 ! ! ABI_ERROR(message) 355 ! end select 356 ! 357 ! !update the new suggested rprimd/acell in the history record 358 ! hist%ihist=abihist_findIndex(hist,+1) 359 ! call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 360 ! 361 ! endif 362 363 end if 364 365 end subroutine pred_hmc