TABLE OF CONTENTS


ABINIT/m_pred_fire [ Modules ]

[ Top ] [ 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-2018 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 .

PARENTS

CHILDREN

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 #include "abi_common.h"
29 
30 module m_pred_fire
31  use defs_basis
32  use m_abicore
33  use m_abimover
34  use m_abihist
35  use m_xfpack
36  use m_geometry,    only : mkrdim, fcart2fred, metric, xred2xcart
37 
38  implicit none
39  private
40  public :: pred_fire
41 
42 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

PARENTS

      mover

CHILDREN

      fcart2fred,hist2var,metric,mkrdim,var2hist
      xfh_recover_new,xfpack_f2vout,xfpack_vin2x,xfpack_x2vin

SOURCE

 96 subroutine pred_fire(ab_mover, ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit)
 97 
 98 
 99 !Arguments ------------------------------------
100 !scalars
101 
102 !This section has been created automatically by the script Abilint (TD).
103 !Do not modify the following lines by hand.
104 #undef ABI_FUNC
105 #define ABI_FUNC 'pred_fire'
106 !End of the abilint section
107 
108  type(abimover),intent(in) :: ab_mover
109  type(ab_xfh_type),intent(inout)    :: ab_xfh
110  type(abiforstr),intent(in) :: forstr
111  type(abihist),intent(inout) :: hist
112  integer, intent(in) :: ionmov
113  integer,intent(in) :: itime
114  integer,intent(in) :: iexit
115  logical,intent(in) :: zDEBUG
116 
117 
118 !Local variables-------------------------------
119 !scalars
120 integer  :: ihist,ihist_prev,ndim
121 integer, parameter :: min_downhill=4
122 integer  :: ierr,ii,jj,kk
123 real(dp),save :: ucvol0
124 real(dp) :: ucvol,det
125 real(dp) :: etotal,etotal_prev
126 real(dp) :: favg
127 ! time step, damping factor initially dtion
128 real(dp),save :: dtratio, alpha
129 ! dtinc: increment of dtratio
130 ! dtdec: decrement of dtratio
131 ! dtmax: maximum allowd value of dtratio
132 ! alphadec: decrement of alpha
133 ! alpha0: initial value of alpha
134 ! mixold: if energy goes up, linear mix old and new coordinates. mixold
135 real(dp), parameter :: dtinc=1.1, dtdec=0.5, dtmax=10.0
136 real(dp), parameter :: alphadec=0.99, alpha0=0.2, mixold=0.3
137 ! v.dot.f
138 real(dp) :: vf
139 ! number of v.dot.f >0
140 integer, save :: ndownhill
141 ! reset_lattice: whether to reset lattice if energy goes up.
142 logical, parameter :: reset_lattice = .true.
143 
144 !arrays
145 real(dp) :: gprimd(3,3)
146 real(dp) :: gmet(3,3)
147 real(dp) :: rmet(3,3)
148 real(dp) :: fcart(3, ab_mover%natom)
149 real(dp) :: acell(3),strten(6), acell_prev(3), acell0(3)
150 real(dp) :: rprim(3,3),rprimd(3,3), rprim_prev(3,3), rprimd_prev(3,3), rprimd0(3,3)
151 real(dp) :: xred(3,ab_mover%natom),xcart(3,ab_mover%natom)
152 real(dp) :: xred_prev(3,ab_mover%natom),xcart_prev(3,ab_mover%natom)
153 real(dp) :: vel_prev(3,ab_mover%natom)
154 ! velocity are saved
155 real(dp) :: vel(3,ab_mover%natom)
156 real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom)
157 real(dp),allocatable, save :: vin(:), vout(:)
158 real(dp), allocatable, save:: vin_prev(:)
159 ! velocity but correspoing to vin&vout, for ion&cell relaxation
160 real(dp),allocatable,save :: vel_ioncell(:)
161 
162 
163 
164 !***************************************************************************
165 !Beginning of executable session
166 !***************************************************************************
167 
168  if(iexit/=0)then
169    if (allocated(vin))           then
170      ABI_DEALLOCATE(vin)
171    end if
172    if (allocated(vout))          then
173      ABI_DEALLOCATE(vout)
174    end if
175    if (allocated(vin_prev))           then
176      ABI_DEALLOCATE(vin_prev)
177    end if
178    if (allocated(vel_ioncell))          then
179      ABI_DEALLOCATE(vel_ioncell)
180    end if
181    return
182  end if
183 
184  !write(std_out,*) 'FIRE 01'
185 !##########################################################
186 !### 01. Compute the dimension of vectors (ndim)
187 
188  ndim=3*ab_mover%natom
189  if(ab_mover%optcell==1 .or.&
190 & ab_mover%optcell==4 .or.&
191 & ab_mover%optcell==5 .or.&
192 & ab_mover%optcell==6) ndim=ndim+1
193  if(ab_mover%optcell==2 .or.&
194 & ab_mover%optcell==3) ndim=ndim+6
195  if(ab_mover%optcell==7 .or.&
196 & ab_mover%optcell==8 .or.&
197 & ab_mover%optcell==9) ndim=ndim+3
198 
199 !write(std_out,*) 'FIRE: ndim=', ndim
200 ! write(std_out,*) 'FIRE 02'
201 !##########################################################
202 !### 02. Allocate the vectors vin
203 
204 !Notice that vin, vout, etc could be allocated
205 !From a previous dataset with a different ndim
206  if(itime==1)then
207    if (allocated(vin))           then
208      ABI_DEALLOCATE(vin)
209    end if
210    if (allocated(vout))          then
211      ABI_DEALLOCATE(vout)
212    end if
213    if (allocated(vin_prev))           then
214      ABI_DEALLOCATE(vin_prev)
215    end if
216    if (allocated(vel_ioncell))          then
217      ABI_DEALLOCATE(vel_ioncell)
218    end if
219 
220    ABI_ALLOCATE(vin,(ndim))
221    ABI_ALLOCATE(vout,(ndim))
222    ABI_ALLOCATE(vin_prev,(ndim))
223    ABI_ALLOCATE(vel_ioncell,(ndim))
224    vel_ioncell(:)=0.0
225  end if
226 
227  !write(std_out,*) 'FIRE 03'
228 !##########################################################
229 !### 03. Obtain the present values from the history
230 
231  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
232 
233  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
234 
235  do ii=1,3
236    rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
237  end do
238 
239  ihist = abihist_findIndex(hist, 0)
240  ihist_prev  = abihist_findIndex(hist,-1)
241 
242  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
243  fcart(:,:)  =hist%fcart(:,:,hist%ihist)
244  strten(:)  =hist%strten(:,hist%ihist)
245  etotal=hist%etot(hist%ihist)
246 
247  if(itime==1) then
248      etotal_prev=0.0
249  else
250      etotal_prev=hist%etot(ihist_prev)
251  endif
252 
253 !Fill the residual with forces (No preconditioning)
254 !Or the preconditioned forces
255  if (ab_mover%goprecon==0)then
256    call fcart2fred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom)
257  else
258    residual(:,:)=forstr%fred(:,:)
259  end if
260 
261 !Save initial values
262  if (itime==1)then
263    acell0(:)=acell(:)
264    rprimd0(:,:)=rprimd(:,:)
265    ucvol0=ucvol
266  end if
267 
268 
269  if(zDEBUG)then
270    write (std_out,*) 'fcart:'
271    do kk=1,ab_mover%natom
272      write (std_out,*) fcart(:,kk)
273    end do
274    write (std_out,*) 'vel:'
275    do kk=1,ab_mover%natom
276      write (std_out,*) vel(:,kk)
277    end do
278    write (std_out,*) 'strten:'
279    write (std_out,*) strten(1:3),ch10,strten(4:6)
280    write (std_out,*) 'etotal:'
281    write (std_out,*) etotal
282  end if
283 
284 !Get rid of mean force on whole unit cell, but only if no
285 !generalized constraints are in effect
286 
287  residual_corrected(:,:)=residual(:,:)
288  if(ab_mover%nconeq==0)then
289    do kk=1,3
290      if (kk/=3.or.ab_mover%jellslab==0) then
291        favg=sum(residual_corrected(kk,:))/dble(ab_mover%natom)
292        residual_corrected(kk,:)=residual_corrected(kk,:)-favg
293      end if
294    end do
295  end if
296 
297  !write(std_out,*) 'FIRE 04'
298 !##########################################################
299 !### 04. Fill the vectors vin and vout
300 
301 !Initialize input vectors : first vin, then vout
302 ! transfer xred, acell, and rprim to vin
303 call xfpack_x2vin(acell, acell0, ab_mover%natom, ndim,&
304 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
305 & ab_mover%symrel, ucvol, ucvol0, vin, xred)
306 !end if
307 
308 !transfer fred and strten to vout. 
309 !Note: fred is not f in reduced co.
310 !but dE/dx
311 
312  call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,&
313 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,&
314 & vout)
315 ! Now vout -> -dE/dx
316 vout(:) = -1.0*vout(:)
317 
318  !write(std_out,*) 'FIRE 05'
319 !##########################################################
320 !### 05. iniialize FIRE
321 if ( itime==1 ) then
322    ndownhill=0
323    alpha=alpha0
324    if (ab_mover%dtion>0)then
325      dtratio = 1.0
326    end if
327 end if
328 
329  !write(std_out,*) 'FIRE 06'
330 !##########################################################
331 !### 06. update timestep
332 ! Note that vin & vout are in reduced coordinates.
333 vf=sum(vel_ioncell*vout)
334 if ( vf >= 0.0_dp .and. (etotal- etotal_prev <0.0_dp) ) then
335 !if ( vf >= 0.0_dp ) then
336     ndownhill=ndownhill+1
337     ! mix v with the v projected on force vector.
338     vel_ioncell(:)=(1.0-alpha)*vel_ioncell(:) + alpha* vout *  &
339 &               sqrt(sum(vel_ioncell*vel_ioncell)/sum(vout*vout))
340     if ( ndownhill>min_downhill ) then
341         dtratio = min(dtratio * dtinc, dtmax)
342         alpha = alpha * alphadec
343     end if
344 else
345     ! reset downhill counter, velocity, alpha. decrease dtratio.
346     ndownhill=0
347     vel_ioncell(:)=0.0
348     alpha=alpha0
349     dtratio = dtratio*dtdec
350 endif
351 
352  !write(std_out,*) 'FIRE 07'
353 !##########################################################
354 !### 07. MD step. update vel_ioncell
355 
356 ! Here mass is not used: all masses=1
357 !write(std_out,*) 'FIRE vin: ', vin
358 ! update v
359 vel_ioncell = vel_ioncell + dtratio*ab_mover%dtion* vout
360 !write(std_out,*) 'FIRE vel: ', vel_ioncell
361 !write(std_out,*) 'FIRE delta x',dtratio*ab_mover%dtion* vel_ioncell
362 ! update x
363 vin = vin + dtratio*ab_mover%dtion* vel_ioncell
364 !write(std_out,*) 'FIRE vin: ', vin
365    
366 !   write(std_out,*) 'FIRE vout: ', vout
367 !   write(std_out,*) 'FIRE vf: ', vf
368 !   write(std_out,*) 'FIRE etotal: ', etotal
369 !   write(std_out,*) 'FIRE etotal_prev: ', etotal_prev
370 !   write(std_out,*) 'FIRE deltaE: ',etotal-etotal_prev
371    write(std_out,*) 'FIRE ndownhill: ', ndownhill
372 !   write(std_out,*) 'FIRE dtratio: ', dtratio
373 !   write(std_out,*) 'FIRE dtion: ', ab_mover%dtion
374 
375 
376 !Implement fixing of atoms : put back old values for fixed
377 !components
378  do kk=1,ab_mover%natom
379    do jj=1,3
380 !    Warning : implemented in reduced coordinates
381      if ( ab_mover%iatfix(jj,kk)==1) then
382        vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3)
383      end if
384    end do
385  end do
386 
387 ! reset_lattice to last step by a ratio if energy is increased.
388 ! disabled for debugging
389 if ( etotal - etotal_prev >0.0 ) then
390     vin= vin*(1-mixold)+vin_prev*mixold
391 end if
392 
393 ! only set vin to vin_prev when energy decreased, so it's 
394 ! possible to go back.
395 ! if (etotal - etotal_prev <0.0 ) then
396  vin_prev(:)=vin(:)
397 ! endif
398 
399 
400 
401 !##########################################################
402 !### 08. update hist.
403 
404 !Increase indexes
405  hist%ihist = abihist_findIndex(hist,+1)
406 !Transfer vin  to xred, acell and rprim
407  call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,&
408 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
409 & ab_mover%symrel, ucvol, ucvol0,&
410 & vin, xred)
411 
412 
413  if(ab_mover%optcell/=0)then
414    call mkrdim(acell,rprim,rprimd)
415    call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
416  end if
417 
418 !Fill the history with the variables
419 !xcart, xred, acell, rprimd
420  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
421  ihist_prev = abihist_findIndex(hist,-1)
422  hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
423 
424  if(zDEBUG)then
425    write (std_out,*) 'residual:'
426    do kk=1,ab_mover%natom
427      write (std_out,*) residual(:,kk)
428    end do
429    write (std_out,*) 'strten:'
430    write (std_out,*) strten(1:3),ch10,strten(4:6)
431    write (std_out,*) 'etotal:'
432    write (std_out,*) etotal
433  end if
434 
435 
436 end subroutine pred_fire