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