TABLE OF CONTENTS
ABINIT/fdtion [ Functions ]
NAME
fdtion
FUNCTION
Compute the apropiated "dtion" from the present values of forces, velocity and viscosity
COPYRIGHT
Copyright (C) 1998-2018 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 . INPUTS (in) hist<type abihist>=Historical record of positions, forces | acell, stresses, and energies, 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 itime: Index of time iteration xcart(3,natom)= cartesian coordinates of atoms fcart(3,natom)= forces in cartesian coordinates vel(3,natom)= velocities OUTPUT (out) fdtion = time step computed
NOTES
PARENTS
pred_moldyn
CHILDREN
SOURCE
446 #if defined HAVE_CONFIG_H 447 #include "config.h" 448 #endif 449 450 function fdtion(ab_mover,itime,xcart,fcart,vel) 451 452 ! define dp,sixth,third,etc... 453 use defs_basis 454 use m_abimover 455 456 !This section has been created automatically by the script Abilint (TD). 457 !Do not modify the following lines by hand. 458 #undef ABI_FUNC 459 #define ABI_FUNC 'fdtion' 460 !End of the abilint section 461 462 implicit none 463 464 !Arguments --------------------------------------------- 465 !scalars 466 type(abimover),intent(in) :: ab_mover 467 integer,intent(in) :: itime 468 real(dp) :: fdtion 469 !arrays 470 real(dp) :: xcart(:,:),fcart(:,:),vel(:,:) 471 472 !Local variables ------------------------------ 473 !scalars 474 integer :: jj,kk 475 real(dp) :: max,min,val 476 real(dp) :: ff,xc,vv,em 477 478 !************************************************************************ 479 480 max=0 481 min=1e6 482 483 do kk=1,ab_mover%natom 484 em=ab_mover%amass(kk) 485 do jj=1,3 486 ff =fcart(jj,kk) 487 xc =xcart(jj,kk) 488 vv=vel(jj,kk) 489 490 if (vv>1e-8) then 491 val=abs(1.0_dp/vv) 492 write(std_out,*) 'vel',kk,jj,val 493 if (val>max) max=val 494 if (val<min) min=val 495 end if 496 497 if (ff>1e-8) then 498 val=sqrt(abs(2*em/ff)) 499 write(std_out,*) 'forces',kk,jj,val,em,ff 500 if (val>max) max=val 501 if (val<min) min=val 502 end if 503 504 end do 505 506 end do 507 508 write(std_out,*) "DTION max=",max 509 write(std_out,*) "DTION min=",min 510 511 if (itime==1)then 512 fdtion=min/10 513 else 514 fdtion=min/10 515 end if 516 517 end function fdtion
ABINIT/pred_moldyn [ Functions ]
NAME
pred_moldyn
FUNCTION
Ionmov predictor (1) Molecular dynamics Molecular dynamics, with or without viscous damping This function should be after the call to scfcv Updates positions, velocities and forces
COPYRIGHT
Copyright (C) 1998-2018 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 .
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.
PARENTS
mover
CHILDREN
hist2var,var2hist,xcart2xred,xred2xcart
SOURCE
61 #if defined HAVE_CONFIG_H 62 #include "config.h" 63 #endif 64 65 #include "abi_common.h" 66 67 subroutine pred_moldyn(ab_mover,hist,icycle,itime,ncycle,ntime,zDEBUG,iexit) 68 69 use defs_basis 70 use m_profiling_abi 71 use m_abimover 72 use m_abihist 73 74 use m_geometry, only : xcart2xred, xred2xcart 75 76 !This section has been created automatically by the script Abilint (TD). 77 !Do not modify the following lines by hand. 78 #undef ABI_FUNC 79 #define ABI_FUNC 'pred_moldyn' 80 use interfaces_45_geomoptim, except_this_one => pred_moldyn 81 !End of the abilint section 82 83 implicit none 84 85 !Arguments ------------------------------------ 86 !scalars 87 type(abimover),intent(in) :: ab_mover 88 type(abihist),intent(inout),target :: hist 89 integer,intent(in) :: icycle 90 integer,intent(inout) :: ncycle 91 integer,intent(in) :: itime 92 integer,intent(in) :: ntime 93 integer,intent(in) :: iexit 94 logical,intent(in) :: zDEBUG 95 96 !Local variables------------------------------- 97 !scalars 98 integer :: kk,jj,ihist,ihist_next,ihist_prev,ihist_prev2 99 integer :: ihist_prev4,ihist_prev5 100 real(dp) :: aa,alfa,bb,cc,x0,xm,em,vis,dx,dv 101 real(dp) :: fcart,fprev,fprev2 102 real(dp) :: xc 103 real(dp) :: vel,vnow,xnow,vprev 104 real(dp),save :: hh,time 105 !arrays 106 real(dp) :: acell(3),rprimd(3,3) 107 real(dp) :: xred(3,ab_mover%natom) 108 real(dp),allocatable :: xcart(:,:),xcart_prev(:,:) 109 real(dp),save,allocatable :: vec_tmp1(:,:) 110 real(dp),save,allocatable :: vec_tmp2(:,:) 111 real(dp), ABI_CONTIGUOUS pointer :: vel_cur(:,:),vel_next(:,:) 112 real(dp),pointer :: fcart_cur(:,:),fcart_prev(:,:),fcart_prev2(:,:) 113 114 !*************************************************************************** 115 !Beginning of executable session 116 !*************************************************************************** 117 118 if(iexit/=0)then 119 if(allocated(vec_tmp1)) then 120 ABI_DEALLOCATE(vec_tmp1) 121 end if 122 if(allocated(vec_tmp2)) then 123 ABI_DEALLOCATE(vec_tmp2) 124 end if 125 return 126 end if 127 128 vis= ab_mover%vis 129 !Just to avoid warnings of uninitialized variables 130 fprev=0.0_dp 131 fprev=0.0_dp 132 fprev2=0.0_dp 133 vnow=0.0_dp 134 vprev=0.0_dp 135 xnow=0.0_dp 136 137 !Those arrays contains intermediary results used with 138 !different meanings during the different time steps 139 !We need to preserv the allocation status, this is the 140 !reason to be 'SAVE' 141 if (itime==1.and.icycle==1)then 142 if(allocated(vec_tmp1)) then 143 ABI_DEALLOCATE(vec_tmp1) 144 end if 145 if(allocated(vec_tmp2)) then 146 ABI_DEALLOCATE(vec_tmp2) 147 end if 148 ABI_ALLOCATE(vec_tmp1,(3,ab_mover%natom)) 149 ABI_ALLOCATE(vec_tmp2,(3,ab_mover%natom)) 150 end if 151 152 !write(std_out,*) '00' 153 !########################################################## 154 !### 00. Copy from the history to the variables 155 156 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 157 158 ABI_ALLOCATE(xcart,(3,ab_mover%natom)) 159 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 160 161 if (itime==1.or.itime==2)then 162 ABI_ALLOCATE(xcart_prev,(3,ab_mover%natom)) 163 call xred2xcart(ab_mover%natom,rprimd,xcart_prev,hist%xred(:,:,1)) 164 end if 165 166 ihist = abihist_findIndex(hist, 0) 167 ihist_prev = abihist_findIndex(hist,-1) 168 ihist_prev2 = abihist_findIndex(hist,-2) 169 ihist_prev4 = abihist_findIndex(hist,-4) 170 ihist_prev5 = abihist_findIndex(hist,-5) 171 ihist_next = abihist_findIndex(hist,+1) 172 173 fcart_cur => hist%fcart(:,:,ihist) 174 if (itime==2) fcart_prev => hist%fcart(:,:,ihist_prev4) 175 if (itime==3) fcart_prev2 => hist%fcart(:,:,ihist_prev5) 176 if (itime >2.or. icycle>=2)fcart_prev => hist%fcart(:,:,ihist_prev) 177 if (itime >3.or. icycle>=3)fcart_prev2 => hist%fcart(:,:,ihist_prev2) 178 179 vel_cur => hist%vel(:,:,ihist) 180 vel_next => hist%vel(:,:,ihist_next) 181 182 !write(std_out,*) '01' 183 !########################################################## 184 !### 01. Get or compute de time step dtion 185 186 if (ab_mover%dtion>0)then 187 hh = ab_mover%dtion 188 else 189 hh=fdtion(ab_mover,itime,xcart,fcart_cur,vel_cur) 190 end if 191 192 !write(std_out,*) '02' 193 !########################################################## 194 !### 02. For all atoms and directions 195 do kk=1,ab_mover%natom 196 em=ab_mover%amass(kk) 197 do jj=1,3 198 199 ! write(std_out,*) '03' 200 ! ########################################################## 201 ! ### 03. Filling other values from history (forces and vel) 202 fcart=fcart_cur(jj,kk) 203 xc=xcart(jj,kk) 204 vel=hist%vel(jj,kk,1) 205 206 ! Previous values only after first iteration 207 if (itime>=2.or.icycle>=2) then 208 fprev=fcart_prev(jj,kk) 209 vprev=hist%vel(jj,kk,hist%ihist) 210 end if 211 if (itime>=3.or.icycle>=3) then 212 fprev2=fcart_prev2(jj,kk) 213 end if 214 215 if (itime==2)then 216 vec_tmp1(jj,kk)=xcart_prev(jj,kk) 217 vec_tmp2(jj,kk)=xcart(jj,kk) 218 end if 219 220 ! write(std_out,*) '04' 221 ! ########################################################## 222 ! ### 04. Take first the atoms that are not allowed to move along 223 ! ### this direction 224 ! ### Warning : implemented in cartesian coordinates 225 if (ab_mover%iatfix(jj,kk)==1) then 226 ! Their positions will be the same as xcart 227 xnow=xcart(jj,kk) 228 ! Their velocities are zero 229 vnow=0.0_dp 230 else 231 232 ! write(std_out,*) '05' 233 ! ########################################################## 234 ! ### 05. Initialization (itime==1): 235 ! ### 4 calls to obtain the forces are neeeded 236 ! ### The variables vec_tmp2 and vec_tmp1 from previous 237 ! ### calls are used in the following ones. 238 if(itime==1)then 239 x0=xcart_prev(jj,kk) 240 241 ! Prepare the second cycle 242 if(icycle==1)then 243 dx=hh*vel 244 dv=hh/em*(fcart-vis*vel) 245 xnow=x0+.5_dp*dx 246 vnow=vel+.5_dp*dv 247 vec_tmp2(jj,kk)=xc+sixth*dx 248 vec_tmp1(jj,kk)=vel+sixth*dv 249 else if(icycle==2)then 250 dx=hh*vprev 251 dv=hh/em*(fcart-vis*vprev) 252 xnow=x0+.5_dp*dx 253 vnow=vel+.5_dp*dv 254 vec_tmp2(jj,kk)=vec_tmp2(jj,kk)+third*dx 255 vec_tmp1(jj,kk)=vec_tmp1(jj,kk)+third*dv 256 else if(icycle==3)then 257 dx=hh*vprev 258 dv=hh/em*(fcart-vis*vprev) 259 xnow=x0+dx 260 vnow=vel+dv 261 vec_tmp2(jj,kk)=vec_tmp2(jj,kk)+third*dx 262 vec_tmp1(jj,kk)=vec_tmp1(jj,kk)+third*dv 263 else if(icycle==4)then 264 dx=hh*vprev 265 dv=hh/em*(fcart-vis*vprev) 266 xnow=vec_tmp2(jj,kk)+sixth*dx 267 vnow=vec_tmp1(jj,kk)+sixth*dv 268 end if 269 else !(itime/=1) 270 271 ! write(std_out,*) '06' 272 ! ########################################################## 273 ! ### 06. Change positions and velocities 274 ! ### These changes only applies for itime>2 275 if (itime>2)then 276 ! Uses a corrector to have better value of xnow, and 277 ! derive vnow. Only update atoms position and 278 ! velocity along its allowed directions 279 aa=fprev 280 bb=(fcart-fprev2)/(2._dp*hh) 281 cc=(fcart+fprev2-2._dp*fprev)/(2._dp*hh*hh) 282 x0=vec_tmp2(jj,kk) 283 xm=vec_tmp1(jj,kk) 284 if(abs(vis)<=1.d-8)then 285 ! NON-DAMPED DYNAMICS (Post-Code) 286 xnow=2._dp*x0-xm+hh**2/em/12._dp*& 287 & (fprev2+10._dp*fprev+fcart) 288 vnow=(bb*hh**2)/(3._dp*em)& 289 & +1.5_dp*aa*hh/em+& 290 & (5._dp/12._dp)*cc*hh**3/em& 291 & +x0/hh-xm/hh 292 else 293 ! DAMPED DYNAMICS (Post-Code) 294 alfa=exp(-vis*hh/em) 295 xnow=((-aa*hh*vis**2+0.5_dp*bb*hh**2*vis**2& 296 & -third*cc*hh**3*vis**2+em*bb*hh*vis& 297 & -em*cc*hh**2*vis-2._dp*em**2*cc*hh+x0*vis**3-xm*vis**3)*alfa& 298 & +aa*hh*vis**2-em*bb*hh*vis+third*cc*hh**3*vis**2& 299 & +2._dp*em**2*cc*hh+0.5D0*bb*hh**2*vis**2-em*cc*hh**2*vis+x0*vis**3)& 300 & /vis**3 301 vnow=(em*aa*vis**2*alfa-em*aa*vis**2+bb*hh*vis**2*em*alfa& 302 & -bb*hh*vis**2*em+cc*hh**2*vis**2*em*alfa-cc*hh**2*vis**2*em& 303 & -em**2*bb*vis*alfa+em**2*bb*vis-2._dp*em**2*cc*hh*vis*alfa+& 304 & 2._dp*em**2*cc*hh*vis+2._dp*em**3*cc*alfa-2._dp*em**3*cc+& 305 & vis**3*alfa**2*aa*hh-0.5_dp*vis**3*alfa**2*bb*hh**2+& 306 & third*vis**3*alfa**2*cc*hh**3-vis**2*& 307 & alfa**2*em*bb*hh+vis**2*alfa**2*em*cc*hh**2+& 308 & 2._dp*vis*alfa**2*em**2*cc*hh-vis**4*alfa**2*x0+& 309 & vis**4*alfa**2*xm)/vis**3/(alfa-1._dp)/em 310 311 end if !if(abs(vis)<=1.d-8) 312 313 xc=xnow 314 vec_tmp1(jj,kk)=vec_tmp2(jj,kk) 315 vec_tmp2(jj,kk)=xnow 316 else 317 vnow=vprev 318 end if !if(itime>2) 319 320 ! write(std_out,*) '07' 321 ! ########################################################## 322 ! ### 07. Correct positions 323 ! ### These changes only applies for itime>1 324 325 if(abs(vis)<=1.d-8)then 326 ! NON-DAMPED DYNAMICS (Pre-Code) 327 ! If the viscosity is too small, the equations become 328 ! ill conditioned due to rounding error so do regular 329 ! Verlet predictor Numerov corrector. 330 x0=vec_tmp2(jj,kk) 331 xm=vec_tmp1(jj,kk) 332 xnow=2._dp*x0-xm& 333 & + hh**2/em*fcart 334 else 335 ! DAMPED DYNAMICS (Pre-Code) 336 ! These equations come from solving 337 ! m*d2x/dt2+vis*dx/dt=a+b*t+c*t**2 338 ! analytically under the boundary conditions that 339 ! x(0)=x0 and x(-h)=xm, and the following is the 340 ! expression for x(h). a, b and c are determined 341 ! from our knowledge of the driving forces. 342 aa=fcart 343 bb=(fcart-fprev)/hh 344 x0=vec_tmp2(jj,kk) 345 xm=vec_tmp1(jj,kk) 346 alfa=exp(-vis*hh/em) 347 xnow=( (-aa*hh*vis**2 +0.5_dp*bb*hh**2*vis**2& 348 & +em*bb*hh*vis +x0*vis**3 -xm*vis**3)*alfa& 349 & +aa*hh*vis**2 -em*bb*hh*vis& 350 & +0.5_dp*bb*hh**2*vis**2 +x0*vis**3)/vis**3 351 ! End of choice between initialisation, damped 352 ! dynamics and non-damped dynamics 353 end if 354 355 end if !if(itime==1) 356 357 end if !if(ab_mover%iatfix(jj,kk)==1) 358 359 ! write(std_out,*) '08' 360 ! ########################################################## 361 ! ### 08. Update history 362 363 xcart(jj,kk)=xnow 364 vel_next(jj,kk)=vnow 365 366 ! write(std_out,*) '09' 367 ! ########################################################## 368 ! ### 09. End loops of atoms and directions 369 end do ! jj=1,3 370 end do ! kk=1,ab_mover%natom 371 372 !write(std_out,*) '10' 373 !########################################################## 374 !### 10. Filling history with the new values 375 376 hist%ihist = abihist_findIndex(hist,+1) 377 378 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 379 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 380 381 !Change ncycle for itime>1 382 if (icycle==4) ncycle=1 383 384 if (itime==1)then 385 time=0.0_dp 386 if (ab_mover%dtion<0)then 387 write(std_out,*) 'Time=',time 388 end if 389 end if 390 time=time+hh 391 hist%time(hist%ihist)=time 392 393 if (allocated(xcart)) then 394 ABI_DEALLOCATE(xcart) 395 end if 396 if (allocated(xcart_prev)) then 397 ABI_DEALLOCATE(xcart_prev) 398 end if 399 400 if (.false.) write(std_out,*) ntime 401 402 end subroutine pred_moldyn