TABLE OF CONTENTS
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.
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 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.
PARENTS
mover,pred_hmc
CHILDREN
hist2var,var2hist,xcart2xred,xred2xcart
SOURCE
54 #if defined HAVE_CONFIG_H 55 #include "config.h" 56 #endif 57 58 #include "abi_common.h" 59 60 61 subroutine pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,hmcflag,icycle,ncycle) 62 63 use defs_basis 64 use m_errors 65 use m_profiling_abi 66 use m_abimover 67 use m_abihist 68 69 use m_geometry, only : xcart2xred, xred2xcart 70 71 !This section has been created automatically by the script Abilint (TD). 72 !Do not modify the following lines by hand. 73 #undef ABI_FUNC 74 #define ABI_FUNC 'pred_velverlet' 75 !End of the abilint section 76 77 implicit none 78 79 !Arguments ------------------------------------ 80 type(abimover),intent(in) :: ab_mover 81 type(abihist),intent(inout) :: hist 82 integer,intent(in) :: itime 83 integer,intent(in) :: ntime 84 integer,intent(in) :: iexit 85 logical,intent(in) :: zDEBUG 86 integer,intent(in),optional :: hmcflag 87 integer,intent(in),optional :: icycle 88 integer,intent(in),optional :: ncycle 89 90 !Local variables------------------------------- 91 92 integer :: ii,jj ! dummy integers for loop indexes 93 real(dp) :: epot,ekin,ekin_tmp ! potential (electronic), kinetic (ionic) energies 94 real(dp) :: xcart(3,ab_mover%natom) ! Cartesian coordinates of all ions 95 real(dp) :: xred(3,ab_mover%natom) ! reduced coordinates of all ions 96 real(dp) :: vel(3,ab_mover%natom) ! ionic velocities in Cartesian coordinates 97 real(dp) :: fcart(3,ab_mover%natom),fred(3,ab_mover%natom) ! forces, Cartesian and reduced coordinates 98 real(dp) :: factor ! factor, indicating change of time step at last iteration 99 integer :: hmcflag_ 100 integer :: icycle_ 101 integer :: ncycle_ 102 103 real(dp) :: acell(3) ! lattice parameters 104 real(dp) :: rprimd(3,3) ! lattice vectors 105 real(dp),allocatable,save :: vel_prev(:,:) ! velocities at the end of each time step (half time step ahead of coordinates) 106 107 !*************************************************************************** 108 !Beginning of executable session 109 !*************************************************************************** 110 111 DBG_ENTER("COLL") 112 113 ! if (option/=1 .and. option/=2 ) then 114 ! write(msg,'(3a,i0)')& 115 !& 'The argument option should be 1 or 2,',ch10,& 116 !& 'however, option=',option 117 ! MSG_BUG(msg) 118 ! end if 119 ! 120 ! if (sizein<1) then 121 ! write(msg,'(3a,i0)')& 122 !& 'The argument sizein should be a positive number,',ch10,& 123 !& 'however, sizein=',sizein 124 ! MSG_ERROR(msg) 125 ! end if 126 127 DBG_EXIT("COLL") 128 129 130 hmcflag_=0 131 if(present(hmcflag))then 132 hmcflag_=hmcflag 133 end if 134 135 icycle_ =0 136 if(present(icycle))then 137 icycle_=icycle 138 end if 139 140 ncycle_ =0 141 if(present(ncycle))then 142 ncycle_=ncycle 143 end if 144 145 146 if(iexit/=0)then 147 if (allocated(vel_prev)) then 148 ABI_DEALLOCATE(vel_prev) 149 end if 150 return 151 end if 152 153 if((hmcflag_==0.and.itime==1).or.(hmcflag_==1.and.icycle_==1))then 154 if (allocated(vel_prev)) then 155 ABI_DEALLOCATE(vel_prev) 156 end if 157 ABI_ALLOCATE(vel_prev,(3,ab_mover%natom)) 158 end if 159 160 161 162 ! Start preparation for velocity verlet, get information about current ionic positions, forces, velocities, etc. 163 164 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 165 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 166 fcart(:,:) = hist%fcart(:,:,hist%ihist) ! forces in Cartesian coordinates 167 vel(:,:) = hist%vel(:,:,hist%ihist) ! velocities of all ions, not needed in reality 168 epot = hist%etot(hist%ihist) ! electronic sub-system energy, not needed 169 ekin = hist%ekin(hist%ihist) ! kinetic energy, not needed 170 171 172 if(zDEBUG)then 173 write (std_out,*) 'velverlet step',itime 174 write (std_out,*) 'fcart:' 175 do ii=1,3 176 write (std_out,*) fcart(ii,:) 177 end do 178 write (std_out,*) 'fred:' 179 do ii=1,3 180 write (std_out,*) fred(ii,:) 181 end do 182 write (std_out,*) 'xcart:' 183 do ii=1,3 184 write (std_out,*) xcart(ii,:) 185 end do 186 write (std_out,*) 'vel:' 187 do ii=1,3 188 write (std_out,*) vel(ii,:) 189 end do 190 end if 191 192 193 if((hmcflag_==0.and.itime==1).or.(hmcflag_==1.and.icycle_==1))then 194 195 !the following breakdown of single time step in two halfs is needed for initialization. 196 !half step velocities "vel_prev" are saved to be used in the next iteration 197 !the velocities "vel" are only used to estimate kinetic energy at correct time instances 198 do ii=1,ab_mover%natom ! propagate velocities half time step forward 199 do jj=1,3 200 vel_prev(jj,ii) = vel(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii) 201 end do 202 end do 203 204 ! propagate velocities half time step forward 205 do ii=1,ab_mover%natom 206 do jj=1,3 207 vel(jj,ii) = vel_prev(jj,ii) + 0.5_dp * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii) 208 end do 209 end do 210 ! use half-step behind velocity values to propagate coordinates one time step forward!!!! 211 do ii=1,ab_mover%natom 212 do jj=1,3 213 xcart(jj,ii) = xcart(jj,ii) + ab_mover%dtion*vel_prev(jj,ii) 214 end do 215 end do 216 ! now, at this 1st iteration, "vel_prev" correspond to a time instance half-step behind 217 ! that of "xcart" 218 219 else 220 221 !at this moment "vel_prev" is behind "xcart" by half of a time step 222 !(saved from the previous iteration) and these are the velocity values to be propagated 223 !using forces that are evaluated at the same time instance as xcart 224 do ii=1,ab_mover%natom ! propagate velocities one time step forward 225 do jj=1,3 226 vel_prev(jj,ii) = vel_prev(jj,ii) + ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii) 227 end do 228 end do 229 !now, the "vel_prev" velocities are half of a time step ahead and can be used to propagate xcart 230 231 if((hmcflag_==0.and.itime==ntime-1).or.(hmcflag_==1.and.icycle_==ncycle_-1))then 232 factor=0.5_dp 233 else 234 factor=one 235 end if 236 237 do ii=1,ab_mover%natom ! propagate coordinates 238 do jj=1,3 239 xcart(jj,ii) = xcart(jj,ii) + factor*ab_mover%dtion*vel_prev(jj,ii) 240 end do 241 end do 242 !to estimate kinetic energy at the same time instance as the potential (electronic sub-system) energy 243 !propagate "vel" another half-time forward (these values are not to be used in the next time-step) 244 do ii=1,ab_mover%natom ! propagate velocities half time step forward 245 do jj=1,3 246 vel(jj,ii) = vel_prev(jj,ii) + (factor-0.5_dp) * ab_mover%dtion*fcart(jj,ii)/ab_mover%amass(ii) 247 end do 248 end do 249 250 ekin=0.0 251 do ii=1,ab_mover%natom 252 do jj=1,3 253 ekin=ekin+0.5_dp*ab_mover%amass(ii)*vel(jj,ii)**2 254 end do 255 end do 256 ekin_tmp=0.0 257 do ii=1,ab_mover%natom 258 do jj=1,3 259 ekin_tmp=ekin_tmp+0.5_dp*ab_mover%amass(ii)*vel_prev(jj,ii)**2 260 end do 261 end do 262 !write(238,*) itime,icycle,ekin_tmp,ekin,epot,factor 263 264 end if 265 266 !Convert new xcart to xred to set correct output values 267 !Update the history with the new coordinates, velocities, etc. 268 269 !Increase indexes 270 hist%ihist=abihist_findIndex(hist,+1) 271 272 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 273 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 274 275 hist%vel(:,:,hist%ihist)=vel(:,:) 276 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 277 278 end subroutine pred_velverlet