TABLE OF CONTENTS
ABINIT/m_pred_moldyn [ Modules ]
NAME
m_pred_moldyn
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, SE) 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_moldyn 22 23 use defs_basis 24 use m_abicore 25 use m_abimover 26 use m_abihist 27 28 use m_geometry, only : xcart2xred, xred2xcart 29 use m_predtk, only : fdtion 30 31 implicit none 32 33 private
ABINIT/pred_moldyn [ Functions ]
NAME
pred_moldyn
FUNCTION
Ionmov predictor (1) Molecular dynamics Molecular dynamics, with or without viscous damping This function should be called after the call to scfcv Updates positions, velocities and forces
INPUTS
ab_mover<type abimover>=Subset of dtset only related with | movement of ions and acell, contains: | dtion: Time step ! natom: Number of atoms | vis: viscosity | iatfix: Index of atoms and directions fixed | amass: Mass of ions icycle: Index of the internal cycle inside a time step (itime) itime: Index of time iteration zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist<type abihist>=Historical record of positions, forces, stresses, cell and energies, ncycle: Number of cycles of a particular time step
NOTES
* This routine is a predictor, it only produces new positions to be computed in the next iteration, this routine should produce not output at all * ncycle changes from 4 for the first iteration (itime==1) to 1 for (itime>1) * The arrays vec_tmp1 and vec_tmp2 are triky, they are use with different meanings, during the initialization they contains working positions and velocities that acumulated produce the first positions of itime=1, for itime>1 they will contain positions in 2 previous steps, those values are different from the values store in the history, thats the reason why we cannot simply use hist%xred to obtain those positions.
SOURCE
88 subroutine pred_moldyn(ab_mover,hist,icycle,itime,ncycle,ntime,zDEBUG,iexit) 89 90 !Arguments ------------------------------------ 91 !scalars 92 type(abimover),intent(in) :: ab_mover 93 type(abihist),intent(inout),target :: hist 94 integer,intent(in) :: icycle 95 integer,intent(inout) :: ncycle 96 integer,intent(in) :: itime 97 integer,intent(in) :: ntime 98 integer,intent(in) :: iexit 99 logical,intent(in) :: zDEBUG 100 101 !Local variables------------------------------- 102 !scalars 103 integer :: kk,jj,ihist,ihist_next,ihist_prev,ihist_prev2 104 integer :: ihist_prev4,ihist_prev5 105 real(dp) :: aa,alfa,bb,cc,x0,xm,em,vis,dx,dv 106 real(dp) :: fcart,fprev,fprev2 107 real(dp) :: xc 108 real(dp) :: vel,vnow,xnow,vprev 109 real(dp),save :: hh,time 110 !arrays 111 real(dp) :: acell(3),rprimd(3,3) 112 real(dp) :: xred(3,ab_mover%natom) 113 real(dp),allocatable :: xcart(:,:),xcart_prev(:,:) 114 real(dp),save,allocatable :: vec_tmp1(:,:) 115 real(dp),save,allocatable :: vec_tmp2(:,:) 116 real(dp), ABI_CONTIGUOUS pointer :: vel_cur(:,:),vel_next(:,:) 117 real(dp),pointer :: fcart_cur(:,:),fcart_prev(:,:),fcart_prev2(:,:) 118 119 !*************************************************************************** 120 !Beginning of executable session 121 !*************************************************************************** 122 123 if(iexit/=0)then 124 if(allocated(vec_tmp1)) then 125 ABI_FREE(vec_tmp1) 126 end if 127 if(allocated(vec_tmp2)) then 128 ABI_FREE(vec_tmp2) 129 end if 130 return 131 end if 132 133 vis= ab_mover%vis 134 !Just to avoid warnings of uninitialized variables 135 fprev=0.0_dp 136 fprev=0.0_dp 137 fprev2=0.0_dp 138 vnow=0.0_dp 139 vprev=0.0_dp 140 xnow=0.0_dp 141 142 !Those arrays contains intermediary results used with 143 !different meanings during the different time steps 144 !We need to preserv the allocation status, this is the 145 !reason to be 'SAVE' 146 if (itime==1.and.icycle==1)then 147 if(allocated(vec_tmp1)) then 148 ABI_FREE(vec_tmp1) 149 end if 150 if(allocated(vec_tmp2)) then 151 ABI_FREE(vec_tmp2) 152 end if 153 ABI_MALLOC(vec_tmp1,(3,ab_mover%natom)) 154 ABI_MALLOC(vec_tmp2,(3,ab_mover%natom)) 155 end if 156 157 !write(std_out,*) '00' 158 !########################################################## 159 !### 00. Copy from the history to the variables 160 161 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 162 163 ABI_MALLOC(xcart,(3,ab_mover%natom)) 164 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 165 166 if (itime==1.or.itime==2)then 167 ABI_MALLOC(xcart_prev,(3,ab_mover%natom)) 168 call xred2xcart(ab_mover%natom,rprimd,xcart_prev,hist%xred(:,:,1)) 169 end if 170 171 ihist = abihist_findIndex(hist, 0) 172 ihist_prev = abihist_findIndex(hist,-1) 173 ihist_prev2 = abihist_findIndex(hist,-2) 174 ihist_prev4 = abihist_findIndex(hist,-4) 175 ihist_prev5 = abihist_findIndex(hist,-5) 176 ihist_next = abihist_findIndex(hist,+1) 177 178 fcart_cur => hist%fcart(:,:,ihist) 179 if (itime==2) fcart_prev => hist%fcart(:,:,ihist_prev4) 180 if (itime==3) fcart_prev2 => hist%fcart(:,:,ihist_prev5) 181 if (itime >2.or. icycle>=2)fcart_prev => hist%fcart(:,:,ihist_prev) 182 if (itime >3.or. icycle>=3)fcart_prev2 => hist%fcart(:,:,ihist_prev2) 183 184 vel_cur => hist%vel(:,:,ihist) 185 vel_next => hist%vel(:,:,ihist_next) 186 187 !write(std_out,*) '01' 188 !########################################################## 189 !### 01. Get or compute the time step dtion 190 191 if (ab_mover%dtion>0)then 192 hh = ab_mover%dtion 193 else 194 hh=fdtion(ab_mover,itime,xcart,fcart_cur,vel_cur) 195 end if 196 197 !write(std_out,*) '02' 198 !########################################################## 199 !### 02. For all atoms and directions 200 do kk=1,ab_mover%natom 201 em=ab_mover%amass(kk) 202 do jj=1,3 203 204 ! write(std_out,*) '03' 205 ! ########################################################## 206 ! ### 03. Filling other values from history (forces and vel) 207 fcart=fcart_cur(jj,kk) 208 xc=xcart(jj,kk) 209 vel=hist%vel(jj,kk,1) 210 211 ! Previous values only after first iteration 212 if (itime>=2.or.icycle>=2) then 213 fprev=fcart_prev(jj,kk) 214 vprev=hist%vel(jj,kk,hist%ihist) 215 end if 216 if (itime>=3.or.icycle>=3) then 217 fprev2=fcart_prev2(jj,kk) 218 end if 219 220 if (itime==2)then 221 vec_tmp1(jj,kk)=xcart_prev(jj,kk) 222 vec_tmp2(jj,kk)=xcart(jj,kk) 223 end if 224 225 ! write(std_out,*) '04' 226 ! ########################################################## 227 ! ### 04. Take first the atoms that are not allowed to move along 228 ! ### this direction 229 ! ### Warning : implemented in cartesian coordinates 230 if (ab_mover%iatfix(jj,kk)==1) then 231 ! Their positions will be the same as xcart 232 xnow=xcart(jj,kk) 233 ! Their velocities are zero 234 vnow=0.0_dp 235 else 236 237 ! write(std_out,*) '05' 238 ! ########################################################## 239 ! ### 05. Initialization (itime==1): 240 ! ### 4 calls to obtain the forces are neeeded 241 ! ### The variables vec_tmp2 and vec_tmp1 from previous 242 ! ### calls are used in the following ones. 243 if(itime==1)then 244 x0=xcart_prev(jj,kk) 245 246 ! Prepare the second cycle 247 if(icycle==1)then 248 dx=hh*vel 249 dv=hh/em*(fcart-vis*vel) 250 xnow=x0+.5_dp*dx 251 vnow=vel+.5_dp*dv 252 vec_tmp2(jj,kk)=xc+sixth*dx 253 vec_tmp1(jj,kk)=vel+sixth*dv 254 else if(icycle==2)then 255 dx=hh*vprev 256 dv=hh/em*(fcart-vis*vprev) 257 xnow=x0+.5_dp*dx 258 vnow=vel+.5_dp*dv 259 vec_tmp2(jj,kk)=vec_tmp2(jj,kk)+third*dx 260 vec_tmp1(jj,kk)=vec_tmp1(jj,kk)+third*dv 261 else if(icycle==3)then 262 dx=hh*vprev 263 dv=hh/em*(fcart-vis*vprev) 264 xnow=x0+dx 265 vnow=vel+dv 266 vec_tmp2(jj,kk)=vec_tmp2(jj,kk)+third*dx 267 vec_tmp1(jj,kk)=vec_tmp1(jj,kk)+third*dv 268 else if(icycle==4)then 269 dx=hh*vprev 270 dv=hh/em*(fcart-vis*vprev) 271 xnow=vec_tmp2(jj,kk)+sixth*dx 272 vnow=vec_tmp1(jj,kk)+sixth*dv 273 end if 274 else !(itime/=1) 275 276 ! write(std_out,*) '06' 277 ! ########################################################## 278 ! ### 06. Change positions and velocities 279 ! ### These changes only applies for itime>2 280 if (itime>2)then 281 ! Uses a corrector to have better value of xnow, and 282 ! derive vnow. Only update atoms position and 283 ! velocity along its allowed directions 284 aa=fprev 285 bb=(fcart-fprev2)/(2._dp*hh) 286 cc=(fcart+fprev2-2._dp*fprev)/(2._dp*hh*hh) 287 x0=vec_tmp2(jj,kk) 288 xm=vec_tmp1(jj,kk) 289 if(abs(vis)<=1.d-8)then 290 ! NON-DAMPED DYNAMICS (Post-Code) 291 xnow=2._dp*x0-xm+hh**2/em/12._dp*& 292 & (fprev2+10._dp*fprev+fcart) 293 vnow=(bb*hh**2)/(3._dp*em)& 294 & +1.5_dp*aa*hh/em+& 295 & (5._dp/12._dp)*cc*hh**3/em& 296 & +x0/hh-xm/hh 297 else 298 ! DAMPED DYNAMICS (Post-Code) 299 alfa=exp(-vis*hh/em) 300 xnow=((-aa*hh*vis**2+0.5_dp*bb*hh**2*vis**2& 301 & -third*cc*hh**3*vis**2+em*bb*hh*vis& 302 & -em*cc*hh**2*vis-2._dp*em**2*cc*hh+x0*vis**3-xm*vis**3)*alfa& 303 & +aa*hh*vis**2-em*bb*hh*vis+third*cc*hh**3*vis**2& 304 & +2._dp*em**2*cc*hh+0.5D0*bb*hh**2*vis**2-em*cc*hh**2*vis+x0*vis**3)& 305 & /vis**3 306 vnow=(em*aa*vis**2*alfa-em*aa*vis**2+bb*hh*vis**2*em*alfa& 307 & -bb*hh*vis**2*em+cc*hh**2*vis**2*em*alfa-cc*hh**2*vis**2*em& 308 & -em**2*bb*vis*alfa+em**2*bb*vis-2._dp*em**2*cc*hh*vis*alfa+& 309 & 2._dp*em**2*cc*hh*vis+2._dp*em**3*cc*alfa-2._dp*em**3*cc+& 310 & vis**3*alfa**2*aa*hh-0.5_dp*vis**3*alfa**2*bb*hh**2+& 311 & third*vis**3*alfa**2*cc*hh**3-vis**2*& 312 & alfa**2*em*bb*hh+vis**2*alfa**2*em*cc*hh**2+& 313 & 2._dp*vis*alfa**2*em**2*cc*hh-vis**4*alfa**2*x0+& 314 & vis**4*alfa**2*xm)/vis**3/(alfa-1._dp)/em 315 316 end if !if(abs(vis)<=1.d-8) 317 318 xc=xnow 319 vec_tmp1(jj,kk)=vec_tmp2(jj,kk) 320 vec_tmp2(jj,kk)=xnow 321 else 322 vnow=vprev 323 end if !if(itime>2) 324 325 ! write(std_out,*) '07' 326 ! ########################################################## 327 ! ### 07. Correct positions 328 ! ### These changes only applies for itime>1 329 330 if(abs(vis)<=1.d-8)then 331 ! NON-DAMPED DYNAMICS (Pre-Code) 332 ! If the viscosity is too small, the equations become 333 ! ill conditioned due to rounding error so do regular 334 ! Verlet predictor Numerov corrector. 335 x0=vec_tmp2(jj,kk) 336 xm=vec_tmp1(jj,kk) 337 xnow=2._dp*x0-xm& 338 & + hh**2/em*fcart 339 else 340 ! DAMPED DYNAMICS (Pre-Code) 341 ! These equations come from solving 342 ! m*d2x/dt2+vis*dx/dt=a+b*t+c*t**2 343 ! analytically under the boundary conditions that 344 ! x(0)=x0 and x(-h)=xm, and the following is the 345 ! expression for x(h). a, b and c are determined 346 ! from our knowledge of the driving forces. 347 aa=fcart 348 bb=(fcart-fprev)/hh 349 x0=vec_tmp2(jj,kk) 350 xm=vec_tmp1(jj,kk) 351 alfa=exp(-vis*hh/em) 352 xnow=( (-aa*hh*vis**2 +0.5_dp*bb*hh**2*vis**2& 353 & +em*bb*hh*vis +x0*vis**3 -xm*vis**3)*alfa& 354 & +aa*hh*vis**2 -em*bb*hh*vis& 355 & +0.5_dp*bb*hh**2*vis**2 +x0*vis**3)/vis**3 356 ! End of choice between initialisation, damped 357 ! dynamics and non-damped dynamics 358 end if 359 360 end if !if(itime==1) 361 362 end if !if(ab_mover%iatfix(jj,kk)==1) 363 364 ! write(std_out,*) '08' 365 ! ########################################################## 366 ! ### 08. Update history 367 368 xcart(jj,kk)=xnow 369 vel_next(jj,kk)=vnow 370 371 ! write(std_out,*) '09' 372 ! ########################################################## 373 ! ### 09. End loops of atoms and directions 374 end do ! jj=1,3 375 end do ! kk=1,ab_mover%natom 376 377 !write(std_out,*) '10' 378 !########################################################## 379 !### 10. Filling history with the new values 380 381 hist%ihist = abihist_findIndex(hist,+1) 382 383 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 384 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 385 386 !Change ncycle for itime>1 387 if (icycle==4) ncycle=1 388 389 if (itime==1)then 390 time=0.0_dp 391 if (ab_mover%dtion<0)then 392 write(std_out,*) 'Time=',time 393 end if 394 end if 395 time=time+hh 396 hist%time(hist%ihist)=time 397 398 if (allocated(xcart)) then 399 ABI_FREE(xcart) 400 end if 401 if (allocated(xcart_prev)) then 402 ABI_FREE(xcart_prev) 403 end if 404 405 if (.false.) write(std_out,*) ntime 406 407 end subroutine pred_moldyn