TABLE OF CONTENTS
ABINIT/pred_verlet [ Functions ]
NAME
pred_verlet
FUNCTION
Ionmov predictors (6 & 7) Verlet algorithm IONMOV 6: Given a starting point xred that is a vector of length 3*natom (reduced nuclei coordinates), a velocity vector (in cartesian coordinates), and unit cell parameters (acell and rprimd - without velocities in the present implementation), the Verlet dynamics is performed, using the gradient of the energy (atomic forces and stresses) as calculated by the routine scfcv. Some atoms can be kept fixed, while the propagation of unit cell parameters is only performed if optcell/=0. No more than ab_mover%ntime steps are performed. The time step is governed by dtion, contained in ab_mover (coming from dtset). Returned quantities are xred, and eventually acell and rprimd (new ones!). IONMOV 7: Block every atom for which the scalar product of velocity and forces is negative, in order to reach the minimum. The convergence requirement on the atomic forces, ab_mover%tolmxf, allows an early exit.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, JCC, 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
ab_mover <type(abimover)> : Datatype with all the information needed by the preditor itime : Index of the present iteration ntime : Maximal number of iterations ionmov : (6 or 7) Specific kind of VERLET zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
NOTES
PARENTS
mover
CHILDREN
fcart2fred,hist2var,metric,mkrdim,var2hist,wrtout,xcart2xred xfpack_f2vout,xfpack_vin2x,xfpack_x2vin,xred2xcart
SOURCE
64 #if defined HAVE_CONFIG_H 65 #include "config.h" 66 #endif 67 68 #include "abi_common.h" 69 70 subroutine pred_verlet(ab_mover,hist,ionmov,itime,ntime,zDEBUG,iexit) 71 72 use defs_basis 73 use m_profiling_abi 74 use m_abimover 75 use m_abihist 76 77 use m_geometry, only : mkrdim, xcart2xred, xred2xcart, fcart2fred, metric 78 79 !This section has been created automatically by the script Abilint (TD). 80 !Do not modify the following lines by hand. 81 #undef ABI_FUNC 82 #define ABI_FUNC 'pred_verlet' 83 use interfaces_14_hidewrite 84 use interfaces_45_geomoptim, except_this_one => pred_verlet 85 !End of the abilint section 86 87 implicit none 88 89 !Arguments ------------------------------------ 90 !scalars 91 type(abimover),intent(in) :: ab_mover 92 type(abihist),intent(inout) :: hist 93 integer,intent(in) :: itime 94 integer,intent(in) :: ntime 95 integer,intent(in) :: ionmov 96 integer,intent(in) :: iexit 97 logical,intent(in) :: zDEBUG 98 99 !Local variables------------------------------- 100 !scalars 101 integer :: ii,istopped,jj,kk,ndim,nstopped 102 real(dp) :: amass_tot,etotal,diag,favg,ekin_corr,scprod,taylor,ucvol0 103 real(dp),save :: ucvol,ucvol_next 104 character(len=500) :: message 105 !arrays 106 integer :: stopped(ab_mover%natom) 107 real(dp) :: acell0(3),fcart(3,ab_mover%natom) 108 real(dp) :: fred_corrected(3,ab_mover%natom) 109 real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3), strten(6) 110 real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom) 111 real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom) 112 real(dp) :: vel(3,ab_mover%natom),vel_nexthalf(3,ab_mover%natom) 113 real(dp),save :: acell(3),acell_next(3) 114 real(dp),save :: rprimd(3,3),rprim(3,3),rprimd_next(3,3),rprim_next(3,3) 115 real(dp),allocatable,save :: hessin(:,:) 116 real(dp),allocatable,save :: vin(:),vin_prev(:),vin_next(:) 117 real(dp),allocatable,save :: vout(:),vout_prev(:),vel_prevhalf(:,:) 118 119 !*************************************************************************** 120 !Beginning of executable session 121 !*************************************************************************** 122 123 if(iexit/=0)then 124 if (allocated(vin)) then 125 ABI_DEALLOCATE(vin) 126 end if 127 if (allocated(vin_next)) then 128 ABI_DEALLOCATE(vin_next) 129 end if 130 if (allocated(vout)) then 131 ABI_DEALLOCATE(vout) 132 end if 133 if (allocated(vin_prev)) then 134 ABI_DEALLOCATE(vin_prev) 135 end if 136 if (allocated(vout_prev)) then 137 ABI_DEALLOCATE(vout_prev) 138 end if 139 if (allocated(hessin)) then 140 ABI_DEALLOCATE(hessin) 141 end if 142 if (allocated(vel_prevhalf)) then 143 ABI_DEALLOCATE(vel_prevhalf) 144 end if 145 return 146 end if 147 148 !write(std_out,*) 'verlet 01' 149 !########################################################## 150 !### 01. Compute the dimension of vectors (ndim) 151 152 ndim=3*ab_mover%natom 153 if(ab_mover%optcell==1 .or.& 154 & ab_mover%optcell==4 .or.& 155 & ab_mover%optcell==5 .or.& 156 & ab_mover%optcell==6) ndim=ndim+1 157 if(ab_mover%optcell==2 .or.& 158 & ab_mover%optcell==3) ndim=ndim+6 159 if(ab_mover%optcell==7 .or.& 160 & ab_mover%optcell==8 .or.& 161 & ab_mover%optcell==9) ndim=ndim+3 162 163 !write(std_out,*) 'verlet 02' 164 !########################################################## 165 !### 02. Allocate the vectors vin, vout and hessian matrix 166 167 !Notice that vin, vout, etc could be allocated 168 !From a previous dataset with a different ndim 169 if(itime==1)then 170 if (allocated(vin)) then 171 ABI_DEALLOCATE(vin) 172 end if 173 if (allocated(vin_next)) then 174 ABI_DEALLOCATE(vin_next) 175 end if 176 if (allocated(vout)) then 177 ABI_DEALLOCATE(vout) 178 end if 179 if (allocated(vin_prev)) then 180 ABI_DEALLOCATE(vin_prev) 181 end if 182 if (allocated(vout_prev)) then 183 ABI_DEALLOCATE(vout_prev) 184 end if 185 if (allocated(hessin)) then 186 ABI_DEALLOCATE(hessin) 187 end if 188 if (allocated(vel_prevhalf)) then 189 ABI_DEALLOCATE(vel_prevhalf) 190 end if 191 192 ABI_ALLOCATE(vin,(ndim)) 193 ABI_ALLOCATE(vin_next,(ndim)) 194 ABI_ALLOCATE(vout,(ndim)) 195 ABI_ALLOCATE(vin_prev,(ndim)) 196 ABI_ALLOCATE(vout_prev,(ndim)) 197 ABI_ALLOCATE(hessin,(ndim,ndim)) 198 ABI_ALLOCATE(vel_prevhalf,(3,ab_mover%natom)) 199 end if 200 201 !write(std_out,*) 'verlet 03' 202 !########################################################## 203 !### 03. Obtain the present values from the history 204 205 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 206 207 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 208 do ii=1,3 209 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 210 end do 211 212 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 213 fcart(:,:) =hist%fcart(:,:,hist%ihist) 214 strten(:) =hist%strten(:,hist%ihist) 215 vel(:,:) =hist%vel(:,:,hist%ihist) 216 etotal =hist%etot(hist%ihist) 217 218 if(zDEBUG)then 219 write (std_out,*) 'fcart:' 220 do kk=1,ab_mover%natom 221 write (std_out,*) fcart(:,kk) 222 end do 223 write (std_out,*) 'vel:' 224 do kk=1,ab_mover%natom 225 write (std_out,*) vel(:,kk) 226 end do 227 write (std_out,*) 'strten:' 228 write (std_out,*) strten(1:3),ch10,strten(4:6) 229 write (std_out,*) 'etotal:' 230 write (std_out,*) etotal 231 end if 232 233 acell0=acell ; ucvol0=ucvol 234 235 !Get rid of mean force on whole unit cell, but only if no 236 !generalized constraints are in effect 237 call fcart2fred(fcart,fred_corrected,rprimd,ab_mover%natom) 238 if(ab_mover%nconeq==0)then 239 amass_tot=sum(ab_mover%amass(:)) 240 do kk=1,3 241 if (kk/=3.or.ab_mover%jellslab==0) then 242 favg=sum(fred_corrected(kk,:))/dble(ab_mover%natom) 243 fred_corrected(kk,:)=fred_corrected(kk,:)-favg*ab_mover%amass(:)/amass_tot 244 end if 245 end do 246 end if 247 248 !write(std_out,*) 'verlet 04' 249 !########################################################## 250 !### 04. Fill the vectors vin and vout 251 252 !Initialize input vectors : first vin, then vout 253 call xfpack_x2vin(acell, acell0, ab_mover%natom, ndim,& 254 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd,& 255 & ab_mover%symrel, ucvol, ucvol0, vin, xred) 256 call xfpack_f2vout(fred_corrected, ab_mover%natom, ndim,& 257 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,& 258 & vout) 259 260 !write(std_out,*) 'verlet 05' 261 !########################################################## 262 !### 05. Initialize or update the hessian matrix 263 264 !Here, set up the matrix of transformation between forces and 265 !acceleration. Masses must be included here. 266 !Beside this feature, one could define 267 !a preconditioner, in which case it should 268 !be the inverse hessian, like in Broyden. This explains the 269 !name chosen for this transformation matrix. This would allow 270 !to find easily the optimal geometry with ionmov=7. 271 !The default, now implemented, corresponds to the identity matrix 272 !in cartesian coordinates, which makes use of metric tensor gmet 273 !in reduced coordinates. 274 275 !Initialise the Hessian matrix using gmet 276 if (itime==1)then 277 ! Initialize inverse hessian with identity matrix 278 ! in cartesian coordinates, which makes use of metric tensor gmet 279 ! in reduced coordinates. 280 hessin(:,:)=zero 281 do ii=1,ab_mover%natom 282 do kk=1,3 283 do jj=1,3 284 ! Warning : implemented in reduced coordinates 285 if (ab_mover%iatfix(kk,ii)==0 .and.& 286 & ab_mover%iatfix(jj,ii)==0 )then 287 hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj)/ab_mover%amass(ii) 288 end if 289 end do 290 end do 291 end do 292 if(ab_mover%optcell/=0)then 293 ! These values might lead to too large changes in some cases ... 294 diag=ab_mover%strprecon*30.0_dp/ucvol 295 if(ab_mover%optcell==1) diag=diag/three 296 do ii=3*ab_mover%natom+1,ndim 297 hessin(ii,ii)=diag 298 end do 299 end if 300 end if 301 302 !zDEBUG (vin,vout and hessin before prediction) 303 if(zDEBUG)then 304 write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]' 305 write(std_out,*) 'vin:' 306 do ii=1,ndim,3 307 if (ii+2<=ndim)then 308 write(std_out,*) ii,vin(ii:ii+2) 309 else 310 write(std_out,*) ii,vin(ii:ndim) 311 end if 312 end do 313 write(std_out,*) 'vout:' 314 do ii=1,ndim,3 315 if (ii+2<=ndim)then 316 write(std_out,*) ii,vout(ii:ii+2) 317 else 318 write(std_out,*) ii,vout(ii:ndim) 319 end if 320 end do 321 write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim 322 do kk=1,ndim 323 do jj=1,ndim,3 324 if (jj+2<=ndim)then 325 write(std_out,*) jj,hessin(jj:jj+2,kk) 326 else 327 write(std_out,*) jj,hessin(jj:ndim,kk) 328 end if 329 end do 330 end do 331 end if 332 333 !write(std_out,*) 'verlet 06' 334 !########################################################## 335 !### 06. Compute the next values 336 337 !%%% VERLET ALGORITHM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 338 339 if(zDEBUG)then 340 write(std_out,*) 'Shifted data ENTER' 341 write(std_out,*) 'acell(:)',acell(:) 342 write(std_out,*) 'rprimd(:,:)',rprimd(:,:) 343 write(std_out,*) 'ucvol',ucvol 344 write(std_out,*) 'vel_prevhalf(:,:)',vel_prevhalf(:,:) 345 write(std_out,*) 'vin_prev(:)',vin_prev(:) 346 write(std_out,*) 'vin(:)',vin(:) 347 write(std_out,*) 'xcart(:,:)',xcart(:,:) 348 write(std_out,*) 'xred(:,:)',xred(:,:) 349 end if 350 351 !Compute next atomic coordinates and cell parameters, using 352 !Verlet algorithm 353 354 !1. First propagate the position, without acceleration 355 if(itime/=1)then 356 vin_next(:)=2*vin(:)-vin_prev(:) 357 taylor=one 358 else 359 ! Initialisation : no vin_prev is available, but the ionic velocity 360 ! is available, in cartesian coordinates 361 ! Uses the velocity 362 xcart_next(:,:)=xcart(:,:)+ab_mover%dtion*vel(:,:) 363 ! Convert back to xred_next (reduced coordinates) 364 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 365 ! Impose no change of acell, ucvol, rprim, and rprimd 366 acell_next(:)=acell(:) 367 ucvol_next=ucvol 368 rprim_next(:,:)=rprim(:,:) 369 rprimd_next(:,:)=rprimd(:,:) 370 write(std_out,*) 'ucvol',ucvol 371 ! Store all these next values in vin_next 372 call xfpack_x2vin(acell_next,acell0,ab_mover%natom,& 373 & ndim,ab_mover%nsym,ab_mover%optcell,rprim_next,& 374 & rprimd,ab_mover%symrel,ucvol_next,ucvol0,& 375 & vin_next,xred_next) 376 taylor=half 377 end if 378 379 !2. Now, take into account the acceleration 380 do ii=1,ndim 381 ! Note the minus sign: the forces are minus the gradients, 382 ! contained in vout. 383 vin_next(:)=vin_next(:)-ab_mover%dtion**2*hessin(:,ii)*& 384 & vout(ii)*taylor 385 end do 386 387 !3. Implement fixing of atoms : put back old values for fixed 388 !components 389 do kk=1,ab_mover%natom 390 do jj=1,3 391 ! Warning : implemented in reduced coordinates 392 if (ab_mover%iatfix(jj,kk) == 1) then 393 vin_next(jj+(kk-1)*3)=vin(jj+(kk-1)*3) 394 end if 395 end do 396 end do 397 398 !4. Now, compute the velocity at the next half-step 399 !Get xred_next, and eventually acell_next, ucvol_next, rprim_next and 400 !rprimd_next, from vin_next 401 call xfpack_vin2x(acell_next,acell0,ab_mover%natom,ndim,& 402 & ab_mover%nsym,ab_mover%optcell,rprim_next,rprimd,& 403 & ab_mover%symrel,ucvol_next,ucvol0,vin_next,xred_next) 404 if(ab_mover%optcell/=0)then 405 call mkrdim(acell_next,rprim_next,rprimd_next) 406 call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol_next) 407 else 408 rprimd_next(:,:)=rprimd(:,:) 409 end if 410 !Convert input xred_next (reduced coordinates) to 411 !xcart_next (cartesian) 412 call xred2xcart(ab_mover%natom,rprimd_next,xcart_next,xred_next) 413 !Compute the velocity at half of the new step 414 vel_nexthalf(:,:)=(xcart_next(:,:)-xcart(:,:))/ab_mover%dtion 415 416 !5. If needed, compute the velocity at present position 417 if(itime/=1)then 418 vel(:,:)=(vel_nexthalf(:,:)+vel_prevhalf(:,:))*0.5_dp 419 end if 420 421 !%%% VERLET ALGORITHM BLOCKING ATOMS %%%%%%%%%%%%%%%%%%%%%% 422 423 !Here, stop the atoms for which the scalar product of velocity 424 !and force is negative, and recompute the kinetic energy. 425 if(ionmov==7)then 426 stopped(:)=0 427 do ii=1,ab_mover%natom 428 scprod=fcart(1,ii)*vel(1,ii)+& 429 & fcart(2,ii)*vel(2,ii)+& 430 & fcart(3,ii)*vel(3,ii) 431 if(scprod<0.0_dp .and. itime/=1)then 432 stopped(ii)=1 433 write(std_out,*) 'Stopped atom',ii 434 ! Shift the velocities of the previous half-step and current 435 ! half-step, so that the acceleration is correct but the 436 ! present velocity vanishes. 437 vel_prevhalf(:,ii)=vel_prevhalf(:,ii)-vel(:,ii) 438 vel_nexthalf(:,ii)=vel_nexthalf(:,ii)-vel(:,ii) 439 vel(:,ii)=0.0_dp 440 xcart_next(:,ii)=xcart(:,ii)+ab_mover%dtion*vel_nexthalf(:,ii) 441 end if 442 end do 443 444 if(zDEBUG)then 445 write (std_out,*) 'fcart:' 446 do kk=1,ab_mover%natom 447 write (std_out,*) fcart(:,kk) 448 end do 449 write (std_out,*) 'vel_prevhalf:' 450 do kk=1,ab_mover%natom 451 write (std_out,*) vel_prevhalf(:,kk) 452 end do 453 write (std_out,*) 'vel_nexthalf:' 454 do kk=1,ab_mover%natom 455 write (std_out,*) vel_nexthalf(:,kk) 456 end do 457 write (std_out,*) 'vel:' 458 do kk=1,ab_mover%natom 459 write (std_out,*) vel(:,kk) 460 end do 461 write (std_out,*) 'xcart_next:' 462 do kk=1,ab_mover%natom 463 write (std_out,*) xcart_next(:,kk) 464 end do 465 end if 466 467 ! Establish a list of stopped atoms 468 nstopped=sum(stopped(:)) 469 470 if(nstopped/=0)then 471 write(message,'(a)') ' List of stopped atoms (ionmov=7) :' 472 call wrtout(ab_out,message,'COLL') 473 istopped=1 474 do ii=1,ab_mover%natom 475 if(stopped(ii)==1)then 476 stopped(istopped)=ii 477 istopped=istopped+1 478 end if 479 end do 480 do ii=1,nstopped,16 481 write(message, '(16i4)' ) stopped(ii:min(ii+15,nstopped)) 482 call wrtout(ab_out,message,'COLL') 483 end do 484 ! Now, compute the corrected kinetic energy 485 ! Generate xred_next from xcart_next 486 call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next) 487 ! Store xred_next, and eventual acell_next and rprim_next in vin 488 call xfpack_x2vin(acell_next,acell0,& 489 & ab_mover%natom,ndim,ab_mover%nsym,ab_mover%optcell,& 490 & rprim_next,rprimd,& 491 & ab_mover%symrel,ucvol_next,ucvol0,vin_next,xred_next) 492 493 ekin_corr=0.0_dp 494 do ii=1,ab_mover%natom 495 do jj=1,3 496 ! Warning : the fixing of atomis is implemented in reduced 497 ! coordinates, so that this expression is wrong 498 if (ab_mover%iatfix(jj,ii) == 0) then 499 ekin_corr=ekin_corr+0.5_dp*ab_mover%amass(ii)*vel(jj,ii)**2 500 end if 501 end do 502 end do 503 ! End of test nstopped/=0 504 end if 505 506 ! End of test ionmov==7 507 end if 508 509 !write(std_out,*) 'verlet 07' 510 !########################################################## 511 !### 07. Shift the data from next values to the present 512 513 !acell(:)=acell_next(:) 514 !rprim(:,:)=rprim_next(:,:) 515 !ucvol=ucvol_next 516 vel_prevhalf(:,:)=vel_nexthalf(:,:) 517 vin_prev(:)=vin(:) 518 vin(:)=vin_next(:) 519 xcart(:,:)=xcart_next(:,:) 520 xred(:,:)=xred_next(:,:) 521 522 write(std_out,*) 'Shifted data EXIT' 523 write(std_out,*) 'acell(:)',acell(:) 524 write(std_out,*) 'rprim(:,:)',rprim(:,:) 525 write(std_out,*) 'rprimd(:,:)',rprimd(:,:) 526 write(std_out,*) 'ucvol',ucvol 527 write(std_out,*) 'vel_prevhalf(:,:)',vel_prevhalf(:,:) 528 write(std_out,*) 'vin_prev(:)',vin_prev(:) 529 write(std_out,*) 'vin(:)',vin(:) 530 write(std_out,*) 'xred(:,:)',xred(:,:) 531 532 533 !write(std_out,*) 'verlet 08' 534 !########################################################## 535 !### 08. Update the history with the prediction 536 537 !Increase indexes 538 hist%ihist=abihist_findIndex(hist,+1) 539 540 !Fill the history with the variables 541 !xred, acell, rprimd, vel 542 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 543 hist%vel(:,:,hist%ihist)=vel(:,:) 544 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 545 546 if(zDEBUG)then 547 write (std_out,*) 'vel:' 548 do kk=1,ab_mover%natom 549 write (std_out,*) vel(:,kk) 550 end do 551 end if 552 553 if (.false.) write(std_out,*) ntime 554 555 end subroutine pred_verlet