TABLE OF CONTENTS
ABINIT/m_pred_velverlet [ Modules ]
NAME
m_pred_velverlet
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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_pred_velverlet 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_abimover 28 use m_abihist 29 30 use m_geometry, only : xcart2xred, xred2xcart 31 32 implicit none 33 34 private
ABINIT/pred_velverlet [ Functions ]
NAME
pred_velverlet
FUNCTION
Velocity Verlet (VV) predictor of ionic positions (ionmov = 24). In constrast to Verlet algorithm, Velocity Verlet is a SYMPLECTIC integrator (for small enough step sizes "dtion", it better conserves the total energy and time-reversibility). VV is a second order integration scheme that requires a single evaluatoin of forces per time step. These properties make VV a good candidate integrator for use in Hybrid Monte Carlo simulation scheme.
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 time step
ntime = total number of time steps
zDEBUG = flag indicating whether to print debug info
iexit = flag indicating finilization of mover loop
hmcflag = optional argument indicating whether the predictor is called from the Hybrid Monte Carlo (HMC) routine
icycle = if hmcflag==1, then icycle providing information about number of HMC cycle is needed
ncycle = if hmcflag==1, then ncycle provides the total number of cycles within one HMC iteration
OUTPUT
hist = history of ionic positions, forces etc. is updated
SIDE EFFECTS
NOTES
This routine can be used either to simulate NVE molecular dynamics (ionmov = 24) or is called from pred_hmc routine (ionmov = 25) to perform updates of ionic positions in Hybrid Monte Carlo iterations.
SOURCE
81 subroutine pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,hmcflag,icycle,ncycle) 82 83 !Arguments ------------------------------------ 84 type(abimover),intent(in) :: ab_mover 85 type(abihist),intent(inout) :: hist 86 integer,intent(in) :: itime 87 integer,intent(in) :: ntime 88 integer,intent(in) :: iexit 89 logical,intent(in) :: zDEBUG 90 integer,intent(in),optional :: hmcflag 91 integer,intent(in),optional :: icycle 92 integer,intent(in),optional :: ncycle 93 94 !Local variables------------------------------- 95 96 integer :: ii,jj ! dummy integers for loop indexes 97 real(dp) :: epot,ekin !,ekin_tmp ! potential (electronic), kinetic (ionic) energies 98 real(dp) :: xcart(3,ab_mover%natom) ! Cartesian coordinates of all ions 99 real(dp) :: xred(3,ab_mover%natom) ! reduced coordinates of all ions 100 real(dp) :: vel(3,ab_mover%natom) ! ionic velocities in Cartesian coordinates 101 real(dp) :: fcart(3,ab_mover%natom),gred(3,ab_mover%natom) ! cartesian forces, and gradient in reduced coordinates 102 !real(dp) :: factor ! factor, indicating change of time step at last iteration 103 integer :: hmcflag_ 104 integer :: icycle_ 105 integer :: ncycle_ 106 107 real(dp) :: acell(3) ! lattice parameters 108 real(dp) :: rprimd(3,3) ! lattice vectors 109 real(dp),allocatable,save :: vel_prev(:,:) ! velocities at the end of each time step (half time step ahead of coordinates) 110 111 !*************************************************************************** 112 !Beginning of executable session 113 !*************************************************************************** 114 115 DBG_ENTER("COLL") 116 117 ABI_UNUSED((/ntime/)) 118 119 hmcflag_=0 120 if(present(hmcflag))then 121 hmcflag_=hmcflag 122 end if 123 124 icycle_ =0 125 if(present(icycle))then 126 icycle_=icycle 127 end if 128 129 ncycle_ =0 130 if(present(ncycle))then 131 ncycle_=ncycle 132 end if 133 134 135 if(iexit/=0)then 136 if (allocated(vel_prev)) then 137 ABI_FREE(vel_prev) 138 end if 139 return 140 end if 141 142 if((hmcflag_==0.and.itime==1).or.(hmcflag_==1.and.icycle_==1))then 143 if (allocated(vel_prev)) then 144 ABI_FREE(vel_prev) 145 end if 146 ABI_MALLOC(vel_prev,(3,ab_mover%natom)) 147 end if 148 149 150 151 ! Start preparation for velocity verlet, get information about current ionic positions, forces, velocities, etc. 152 153 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 154 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 155 fcart(:,:) = hist%fcart(:,:,hist%ihist) ! forces in Cartesian coordinates 156 vel(:,:) = hist%vel(:,:,hist%ihist) ! velocities of all ions, not needed in reality 157 epot = hist%etot(hist%ihist) ! electronic sub-system energy, not needed 158 ekin = hist%ekin(hist%ihist) ! kinetic energy, not needed 159 160 161 if(zDEBUG)then 162 write (std_out,*) 'velverlet step',itime 163 write (std_out,*) 'fcart:' 164 do ii=1,3 165 write (std_out,*) fcart(ii,:) 166 end do 167 write (std_out,*) 'gred:' 168 do ii=1,3 169 write (std_out,*) gred(ii,:) 170 end do 171 write (std_out,*) 'xcart:' 172 do ii=1,3 173 write (std_out,*) xcart(ii,:) 174 end do 175 write (std_out,*) 'vel:' 176 do ii=1,3 177 write (std_out,*) vel(ii,:) 178 end do 179 end if 180 181 182 if((hmcflag_==0.and.itime==1).or.(hmcflag_==1.and.icycle_==1))then 183 184 !the following breakdown of single time step in two halfs is needed for initialization. 185 !half step velocities "vel_prev" are saved to be used in the next iteration 186 !the velocities "vel" are only used to estimate kinetic energy at correct time instances 187 do ii=1,ab_mover%natom ! propagate velocities half time step forward 188 do jj=1,3 189 vel_prev(jj,ii) = vel(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii) 190 end do 191 end do 192 193 ! propagate velocities half time step forward 194 do ii=1,ab_mover%natom 195 do jj=1,3 196 vel(jj,ii) = vel_prev(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii) 197 end do 198 end do 199 ! use half-step behind velocity values to propagate coordinates one time step forward!!!! 200 do ii=1,ab_mover%natom 201 do jj=1,3 202 xcart(jj,ii) = xcart(jj,ii) + ab_mover%dtion*vel_prev(jj,ii) 203 end do 204 end do 205 ! now, at this 1st iteration, "vel_prev" correspond to a time instance half-step behind 206 ! that of "xcart" 207 208 else 209 210 !at this moment "vel_prev" is behind "xcart" by half of a time step 211 !(saved from the previous iteration) and these are the velocity values to be propagated 212 !using forces that are evaluated at the same time instance as xcart 213 do ii=1,ab_mover%natom ! propagate velocities one time step forward 214 do jj=1,3 215 vel_prev(jj,ii) = vel_prev(jj,ii) + ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii) 216 end do 217 end do 218 !now, the "vel_prev" velocities are half of a time step ahead and can be used to propagate xcart 219 220 !if((hmcflag_==0.and.itime==ntime-1).or.(hmcflag_==1.and.icycle_==ncycle_-1))then 221 ! factor=0.5_dp 222 !else 223 ! factor=one 224 !end if 225 226 do ii=1,ab_mover%natom ! propagate coordinates 227 do jj=1,3 228 xcart(jj,ii) = xcart(jj,ii) + ab_mover%dtion*vel_prev(jj,ii) 229 !xcart(jj,ii) = xcart(jj,ii) + factor*ab_mover%dtion*vel_prev(jj,ii) 230 end do 231 end do 232 !to estimate kinetic energy at the same time instance as the potential (electronic sub-system) energy 233 !propagate "vel" another half-time forward (these values are not to be used in the next time-step) 234 do ii=1,ab_mover%natom ! propagate velocities half time step forward 235 do jj=1,3 236 vel(jj,ii) = vel_prev(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii) 237 !vel(jj,ii) = vel_prev(jj,ii) +(factor-0.5_dp) * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii) 238 end do 239 end do 240 241 ekin=0.0 242 do ii=1,ab_mover%natom 243 do jj=1,3 244 ekin=ekin+0.5_dp*ab_mover%amass(ii)*vel(jj,ii)**2 245 end do 246 end do 247 !ekin_tmp=0.0 248 !do ii=1,ab_mover%natom 249 ! do jj=1,3 250 ! ekin_tmp=ekin_tmp+0.5_dp*ab_mover%amass(ii)*vel_prev(jj,ii)**2 251 ! end do 252 !end do 253 !write(238,*) itime,icycle,ekin_tmp,ekin,epot,factor 254 255 end if 256 257 !Convert new xcart to xred to set correct output values 258 !Update the history with the new coordinates, velocities, etc. 259 260 !Increase indexes 261 hist%ihist=abihist_findIndex(hist,+1) 262 263 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 264 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 265 266 hist%vel(:,:,hist%ihist)=vel(:,:) 267 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 268 269 DBG_EXIT("COLL") 270 271 end subroutine pred_velverlet