TABLE OF CONTENTS
ABINIT/m_pred_fire [ Modules ]
NAME
m_pred_fire
FUNCTION
Ionmov predictors (15) FIRE algorithm The fast inertial relaxation engine (FIRE) method for relaxation. The method is described in Erik Bitzek, Pekka Koskinen, Franz G"ahler, Michael Moseler, and Peter Gumbsch, Phys. Rev. Lett. 97, 170201 [[cite:Bitzek2006]]
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (hexu) 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
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 module m_pred_fire 26 use defs_basis 27 use m_abicore 28 use m_abimover 29 use m_abihist 30 use m_xfpack 31 use m_geometry, only : mkrdim, fcart2gred, metric, xred2xcart 32 use m_errors, only: unused_var 33 34 implicit none 35 36 private 37 38 public :: pred_fire 39 40 contains
m_pred_fire/pred_fire [ Functions ]
[ Top ] [ m_pred_fire ] [ Functions ]
NAME
pred_fire
FUNCTION
IONMOV 15: 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. At each step, the dot product of velocity and force is calculated. If the v.dot.f is positive for min_downhill consecutive steps, the time step dtion is increased by dtinc until it reaches dtmax. If v.dot.f is negative, dtion is mulitiplied by dtdec. 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!).
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 : (15) FIRE. Not used in function. just to keep the same format as bfgs. zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
SOURCE
87 subroutine pred_fire(ab_mover, ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit) 88 89 !Arguments ------------------------------------ 90 !scalars 91 type(abimover),intent(in) :: ab_mover 92 type(ab_xfh_type),intent(inout) :: ab_xfh 93 type(abiforstr),intent(in) :: forstr 94 type(abihist),intent(inout) :: hist 95 integer, intent(in) :: ionmov 96 integer,intent(in) :: itime 97 integer,intent(in) :: iexit 98 logical,intent(in) :: zDEBUG 99 100 101 !Local variables------------------------------- 102 !scalars 103 integer :: ihist,ihist_prev,ndim 104 integer, parameter :: min_downhill=4 105 integer :: ii,jj,kk 106 real(dp),save :: ucvol0 107 real(dp) :: ucvol 108 real(dp) :: etotal,etotal_prev 109 real(dp) :: favg 110 ! time step, damping factor initially dtion 111 real(dp),save :: dtratio, alpha 112 ! dtinc: increment of dtratio 113 ! dtdec: decrement of dtratio 114 ! dtmax: maximum allowd value of dtratio 115 ! alphadec: decrement of alpha 116 ! alpha0: initial value of alpha 117 ! mixold: if energy goes up, linear mix old and new coordinates. mixold 118 real(dp), parameter :: dtinc=1.1, dtdec=0.5, dtmax=10.0 119 real(dp), parameter :: alphadec=0.99, alpha0=0.2, mixold=0.3 120 ! v.dot.f 121 real(dp) :: vf 122 ! number of v.dot.f >0 123 integer, save :: ndownhill 124 ! reset_lattice: whether to reset lattice if energy goes up. 125 logical, parameter :: reset_lattice = .true. 126 127 !arrays 128 real(dp) :: gprimd(3,3) 129 real(dp) :: gmet(3,3) 130 real(dp) :: rmet(3,3) 131 real(dp) :: fcart(3, ab_mover%natom) 132 real(dp) :: acell(3),strten(6), acell0(3) 133 real(dp) :: rprim(3,3),rprimd(3,3), rprimd0(3,3) 134 real(dp) :: xred(3,ab_mover%natom),xcart(3,ab_mover%natom) 135 ! velocity are saved 136 real(dp) :: vel(3,ab_mover%natom) 137 real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom) 138 real(dp),allocatable, save :: vin(:), vout(:) 139 real(dp), allocatable, save:: vin_prev(:) 140 ! velocity but correspoing to vin&vout, for ion&cell relaxation 141 real(dp),allocatable,save :: vel_ioncell(:) 142 143 ABI_UNUSED((/ionmov, ab_xfh%mxfh/)) 144 145 !*************************************************************************** 146 !Beginning of executable session 147 !*************************************************************************** 148 149 if(iexit/=0)then 150 ABI_SFREE(vin) 151 ABI_SFREE(vout) 152 ABI_SFREE(vin_prev) 153 ABI_SFREE(vel_ioncell) 154 return 155 end if 156 157 write(std_out,*) 'FIRE 01' 158 !########################################################## 159 !### 01. Compute the dimension of vectors (ndim) 160 161 ndim=3*ab_mover%natom 162 if(ab_mover%optcell==1) ndim=ndim+1 163 if(ab_mover%optcell==2 .or.& 164 & ab_mover%optcell==3) ndim=ndim+6 165 if(ab_mover%optcell>=4) ndim=ndim+3 166 167 write(std_out,*) 'FIRE: ndim=', ndim 168 write(std_out,*) 'FIRE 02' 169 !########################################################## 170 !### 02. Allocate the vectors vin 171 172 !Notice that vin, vout, etc could be allocated 173 !From a previous dataset with a different ndim 174 if(itime==1)then 175 ABI_SFREE(vin) 176 ABI_SFREE(vout) 177 ABI_SFREE(vin_prev) 178 ABI_SFREE(vel_ioncell) 179 180 ABI_MALLOC(vin,(ndim)) 181 ABI_MALLOC(vout,(ndim)) 182 ABI_MALLOC(vin_prev,(ndim)) 183 ABI_MALLOC(vel_ioncell,(ndim)) 184 vel_ioncell(:)=0.0 185 end if 186 187 write(std_out,*) 'FIRE 03' 188 !########################################################## 189 !### 03. Obtain the present values from the history 190 191 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 192 193 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 194 195 do ii=1,3 196 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 197 end do 198 199 ihist = abihist_findIndex(hist, 0) 200 ihist_prev = abihist_findIndex(hist,-1) 201 202 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 203 fcart(:,:) =hist%fcart(:,:,hist%ihist) 204 strten(:) =hist%strten(:,hist%ihist) 205 etotal=hist%etot(hist%ihist) 206 207 if(itime==1) then 208 etotal_prev=0.0 209 else 210 etotal_prev=hist%etot(ihist_prev) 211 endif 212 213 !Fill the residual with forces (No preconditioning) 214 !Or the preconditioned forces 215 if (ab_mover%goprecon==0)then 216 call fcart2gred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom) 217 else 218 residual(:,:)=forstr%gred(:,:) 219 end if 220 221 !Save initial values 222 if (itime==1)then 223 acell0(:)=acell(:) 224 rprimd0(:,:)=rprimd(:,:) 225 ucvol0=ucvol 226 end if 227 228 229 if(zDEBUG)then 230 write (std_out,*) 'fcart:' 231 do kk=1,ab_mover%natom 232 write (std_out,*) fcart(:,kk) 233 end do 234 write (std_out,*) 'vel:' 235 do kk=1,ab_mover%natom 236 write (std_out,*) vel(:,kk) 237 end do 238 write (std_out,*) 'strten:' 239 write (std_out,*) strten(1:3),ch10,strten(4:6) 240 write (std_out,*) 'etotal:' 241 write (std_out,*) etotal 242 end if 243 244 !Get rid of mean force on whole unit cell, but only if no 245 !generalized constraints are in effect 246 247 residual_corrected(:,:)=residual(:,:) 248 if(ab_mover%nconeq==0)then 249 do kk=1,3 250 if (kk/=3.or.ab_mover%jellslab==0) then 251 favg=sum(residual_corrected(kk,:))/dble(ab_mover%natom) 252 residual_corrected(kk,:)=residual_corrected(kk,:)-favg 253 end if 254 end do 255 end if 256 257 write(std_out,*) 'FIRE 04' 258 !########################################################## 259 !### 04. Fill the vectors vin and vout 260 261 !Initialize input vectors : first vin, then vout 262 ! transfer xred, acell, and rprim to vin 263 call xfpack_x2vin(acell, ab_mover%natom, ndim,& 264 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,& 265 & ab_mover%symrel, ucvol, ucvol0, vin, xred) 266 !end if 267 268 !transfer gred and strten to vout. 269 !Note: gred is not f in reduced co. 270 !but dE/dx 271 272 call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,& 273 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,& 274 & vout) 275 ! Now vout -> -dE/dx 276 vout(:) = -1.0*vout(:) 277 278 write(std_out,*) 'FIRE 05' 279 !########################################################## 280 !### 05. iniialize FIRE 281 if ( itime==1 ) then 282 ndownhill=0 283 alpha=alpha0 284 if (ab_mover%dtion>0)then 285 dtratio = 1.0 286 end if 287 end if 288 289 write(std_out,*) 'FIRE 06' 290 !########################################################## 291 !### 06. update timestep 292 ! Note that vin & vout are in reduced coordinates. 293 vf=sum(vel_ioncell*vout) 294 if ( vf >= 0.0_dp .and. (etotal- etotal_prev <0.0_dp) ) then 295 !if ( vf >= 0.0_dp ) then 296 ndownhill=ndownhill+1 297 ! mix v with the v projected on force vector. 298 vel_ioncell(:)=(1.0-alpha)*vel_ioncell(:) + alpha* vout * & 299 & sqrt(sum(vel_ioncell*vel_ioncell)/sum(vout*vout)) 300 if ( ndownhill>min_downhill ) then 301 dtratio = min(dtratio * dtinc, dtmax) 302 alpha = alpha * alphadec 303 end if 304 else 305 ! reset downhill counter, velocity, alpha. decrease dtratio. 306 ndownhill=0 307 vel_ioncell(:)=0.0 308 alpha=alpha0 309 dtratio = dtratio*dtdec 310 endif 311 312 write(std_out,*) 'FIRE 07' 313 !########################################################## 314 !### 07. MD step. update vel_ioncell 315 316 ! Here mass is not used: all masses=1 317 !write(std_out,*) 'FIRE vin: ', vin 318 ! update v 319 vel_ioncell = vel_ioncell + dtratio*ab_mover%dtion* vout 320 !write(std_out,*) 'FIRE vel: ', vel_ioncell 321 !write(std_out,*) 'FIRE delta x',dtratio*ab_mover%dtion* vel_ioncell 322 ! update x 323 vin = vin + dtratio*ab_mover%dtion* vel_ioncell 324 !write(std_out,*) 'FIRE vin: ', vin 325 326 ! write(std_out,*) 'FIRE vout: ', vout 327 ! write(std_out,*) 'FIRE vf: ', vf 328 ! write(std_out,*) 'FIRE etotal: ', etotal 329 ! write(std_out,*) 'FIRE etotal_prev: ', etotal_prev 330 ! write(std_out,*) 'FIRE deltaE: ',etotal-etotal_prev 331 write(std_out,*) 'FIRE ndownhill: ', ndownhill 332 ! write(std_out,*) 'FIRE dtratio: ', dtratio 333 ! write(std_out,*) 'FIRE dtion: ', ab_mover%dtion 334 335 336 !Implement fixing of atoms : put back old values for fixed 337 !components 338 do kk=1,ab_mover%natom 339 do jj=1,3 340 ! Warning : implemented in reduced coordinates 341 if ( ab_mover%iatfix(jj,kk)==1) then 342 vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3) 343 end if 344 end do 345 end do 346 347 ! reset_lattice to last step by a ratio if energy is increased. 348 ! disabled for debugging 349 if ( etotal - etotal_prev >0.0 ) then 350 vin= vin*(1-mixold)+vin_prev*mixold 351 end if 352 353 ! only set vin to vin_prev when energy decreased, so it's 354 ! possible to go back. 355 ! if (etotal - etotal_prev <0.0 ) then 356 vin_prev(:)=vin(:) 357 ! endif 358 359 360 361 !########################################################## 362 !### 08. update hist. 363 364 !Increase indexes 365 hist%ihist = abihist_findIndex(hist,+1) 366 !Transfer vin to xred, acell and rprim 367 call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,& 368 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,& 369 & ab_mover%symrel, ucvol, ucvol0,& 370 & vin, xred) 371 372 if(ab_mover%optcell/=0)then 373 call mkrdim(acell,rprim,rprimd) 374 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 375 end if 376 377 !Fill the history with the variables 378 !xcart, xred, acell, rprimd 379 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 380 ihist_prev = abihist_findIndex(hist,-1) 381 hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev) 382 383 if(zDEBUG)then 384 write (std_out,*) 'residual:' 385 do kk=1,ab_mover%natom 386 write (std_out,*) residual(:,kk) 387 end do 388 write (std_out,*) 'strten:' 389 write (std_out,*) strten(1:3),ch10,strten(4:6) 390 write (std_out,*) 'etotal:' 391 write (std_out,*) etotal 392 end if 393 394 395 end subroutine pred_fire