Branch data Line data Source code
1 : : !!****m* ABINIT/m_pred_fire
2 : : !! NAME
3 : : !! m_pred_fire
4 : : !!
5 : : !! FUNCTION
6 : : !! Ionmov predictors (15) FIRE algorithm
7 : : !! The fast inertial relaxation engine (FIRE) method for relaxation.
8 : : !! The method is described in Erik Bitzek, Pekka Koskinen, Franz G"ahler,
9 : : !! Michael Moseler, and Peter Gumbsch, Phys. Rev. Lett. 97, 170201 [[cite:Bitzek2006]]
10 : : !!
11 : : !! COPYRIGHT
12 : : !! Copyright (C) 1998-2026 ABINIT group (hexu)
13 : : !! This file is distributed under the terms of the
14 : : !! GNU General Public License, see ~abinit/COPYING
15 : : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : : !!
17 : : !! SOURCE
18 : :
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
41 : :
42 : : !!***
43 : :
44 : : !!****f* m_pred_fire/pred_fire
45 : : !! NAME
46 : : !! pred_fire
47 : : !!
48 : : !! FUNCTION
49 : : !!
50 : : !! IONMOV 15:
51 : : !! Given a starting point xred that is a vector of length 3*natom
52 : : !! (reduced nuclei coordinates), a velocity vector (in cartesian
53 : : !! coordinates), and unit cell parameters (acell and rprimd -
54 : : !! without velocities in the present implementation), the Verlet
55 : : !! dynamics is performed, using the gradient of the energy
56 : : !! (atomic forces and stresses) as calculated by the routine scfcv.
57 : : !!
58 : : !! At each step, the dot product of velocity and force is calculated.
59 : : !! If the v.dot.f is positive for min_downhill consecutive steps, the
60 : : !! time step dtion is increased by dtinc until it reaches dtmax. If
61 : : !! v.dot.f is negative, dtion is mulitiplied by dtdec.
62 : : !!
63 : : !! Some atoms can be kept fixed, while the propagation of unit cell
64 : : !! parameters is only performed if optcell/=0.
65 : : !! No more than ab_mover%ntime steps are performed.
66 : : !! The time step is governed by dtion, contained in ab_mover
67 : : !! (coming from dtset).
68 : : !! Returned quantities are xred, and eventually acell and rprimd
69 : : !! (new ones!).
70 : : !!
71 : :
72 : : !! INPUTS
73 : : !! ab_mover <type(abimover)> : Datatype with all the information
74 : : !! needed by the preditor
75 : : !! itime : Index of the present iteration
76 : : !! ntime : Maximal number of iterations
77 : : !! ionmov : (15) FIRE. Not used in function. just to keep the same format as bfgs.
78 : : !! zDEBUG : if true print some debugging information
79 : : !!
80 : : !! OUTPUT
81 : : !!
82 : : !! SIDE EFFECTS
83 : : !! hist <type(abihist)> : History of positions,forces
84 : : !! acell, rprimd, stresses
85 : : !!
86 : : !! SOURCE
87 : 12 : 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 : 24 : 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 : 24 : real(dp) :: xred(3,ab_mover%natom),xcart(3,ab_mover%natom)
135 : : ! velocity are saved
136 : 24 : real(dp) :: vel(3,ab_mover%natom)
137 : 12 : 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 [ + + ]: 12 : if(iexit/=0)then
150 [ + - ]: 1 : ABI_SFREE(vin)
151 [ + - ]: 1 : ABI_SFREE(vout)
152 [ + - ]: 1 : ABI_SFREE(vin_prev)
153 [ + - ]: 1 : ABI_SFREE(vel_ioncell)
154 : : return
155 : : end if
156 : :
157 : 11 : write(std_out,*) 'FIRE 01'
158 : : !##########################################################
159 : : !### 01. Compute the dimension of vectors (ndim)
160 : :
161 : 11 : ndim=3*ab_mover%natom
162 [ - + ]: 11 : if(ab_mover%optcell==1) ndim=ndim+1
163 [ - + ]: 11 : if(ab_mover%optcell==2 .or.&
164 : 0 : & ab_mover%optcell==3) ndim=ndim+6
165 [ - + ]: 11 : if(ab_mover%optcell>=4) ndim=ndim+3
166 : :
167 : 11 : write(std_out,*) 'FIRE: ndim=', ndim
168 : 11 : 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 [ + + ]: 11 : if(itime==1)then
175 [ - + ]: 1 : ABI_SFREE(vin)
176 [ - + ]: 1 : ABI_SFREE(vout)
177 [ - + ]: 1 : ABI_SFREE(vin_prev)
178 [ - + ]: 1 : ABI_SFREE(vel_ioncell)
179 : :
180 [ - + + - : 3 : ABI_MALLOC(vin,(ndim))
+ - + - -
+ - + -
+ ]
181 [ - + + - : 2 : ABI_MALLOC(vout,(ndim))
- + - + -
+ ]
182 [ - + + - : 2 : ABI_MALLOC(vin_prev,(ndim))
- + - + -
+ ]
183 [ - + + - : 2 : ABI_MALLOC(vel_ioncell,(ndim))
- + - + -
+ ]
184 [ + + ]: 28 : vel_ioncell(:)=0.0
185 : : end if
186 : :
187 : 11 : write(std_out,*) 'FIRE 03'
188 : : !##########################################################
189 : : !### 03. Obtain the present values from the history
190 : :
191 : 11 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
192 : :
193 : 11 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
194 : :
195 [ + + ]: 44 : do ii=1,3
196 [ + + ]: 143 : rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
197 : : end do
198 : :
199 : 11 : ihist = abihist_findIndex(hist, 0)
200 : 11 : ihist_prev = abihist_findIndex(hist,-1)
201 : :
202 : 11 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
203 [ + + + + ]: 407 : fcart(:,:) =hist%fcart(:,:,hist%ihist)
204 [ + + ]: 77 : strten(:) =hist%strten(:,hist%ihist)
205 : 11 : etotal=hist%etot(hist%ihist)
206 : :
207 [ + + ]: 11 : if(itime==1) then
208 : : etotal_prev=0.0
209 : : else
210 : 10 : etotal_prev=hist%etot(ihist_prev)
211 : : endif
212 : :
213 : : !Fill the residual with forces (No preconditioning)
214 : : !Or the preconditioned forces
215 [ + - ]: 11 : if (ab_mover%goprecon==0)then
216 : 11 : call fcart2gred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom)
217 : : else
218 [ # # # # ]: 0 : residual(:,:)=forstr%gred(:,:)
219 : : end if
220 : :
221 : : !Save initial values
222 [ + + ]: 11 : if (itime==1)then
223 : 1 : acell0(:)=acell(:)
224 : 1 : rprimd0(:,:)=rprimd(:,:)
225 : 1 : ucvol0=ucvol
226 : : end if
227 : :
228 : :
229 [ - + ]: 11 : if(zDEBUG)then
230 : 0 : write (std_out,*) 'fcart:'
231 [ # # ]: 0 : do kk=1,ab_mover%natom
232 : 0 : write (std_out,*) fcart(:,kk)
233 : : end do
234 : 0 : write (std_out,*) 'vel:'
235 [ # # ]: 0 : do kk=1,ab_mover%natom
236 : 0 : write (std_out,*) vel(:,kk)
237 : : end do
238 : 0 : write (std_out,*) 'strten:'
239 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
240 : 0 : write (std_out,*) 'etotal:'
241 : 0 : 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 [ + + + + ]: 407 : residual_corrected(:,:)=residual(:,:)
248 [ + - ]: 11 : if(ab_mover%nconeq==0)then
249 [ + + ]: 44 : do kk=1,3
250 [ + + + - ]: 44 : if (kk/=3.or.ab_mover%jellslab==0) then
251 [ + + ]: 330 : favg=sum(residual_corrected(kk,:))/dble(ab_mover%natom)
252 [ + + ]: 330 : residual_corrected(kk,:)=residual_corrected(kk,:)-favg
253 : : end if
254 : : end do
255 : : end if
256 : :
257 : 11 : 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 [ + - - - : 11 : & 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 [ + - - - : 11 : & vout)
- - - + ]
275 : : ! Now vout -> -dE/dx
276 [ + + ]: 308 : vout(:) = -1.0*vout(:)
277 : :
278 : 11 : write(std_out,*) 'FIRE 05'
279 : : !##########################################################
280 : : !### 05. iniialize FIRE
281 [ + + ]: 11 : if ( itime==1 ) then
282 : 1 : ndownhill=0
283 : 1 : alpha=alpha0
284 [ + - ]: 1 : if (ab_mover%dtion>0)then
285 : 1 : dtratio = 1.0
286 : : end if
287 : : end if
288 : :
289 : 11 : write(std_out,*) 'FIRE 06'
290 : : !##########################################################
291 : : !### 06. update timestep
292 : : ! Note that vin & vout are in reduced coordinates.
293 [ + + ]: 308 : vf=sum(vel_ioncell*vout)
294 [ + + + - ]: 11 : if ( vf >= 0.0_dp .and. (etotal- etotal_prev <0.0_dp) ) then
295 : : !if ( vf >= 0.0_dp ) then
296 : 10 : ndownhill=ndownhill+1
297 : : ! mix v with the v projected on force vector.
298 : : vel_ioncell(:)=(1.0-alpha)*vel_ioncell(:) + alpha* vout * &
299 [ + + + + : 820 : & sqrt(sum(vel_ioncell*vel_ioncell)/sum(vout*vout))
+ + ]
300 [ + + ]: 10 : if ( ndownhill>min_downhill ) then
301 : 2 : dtratio = min(dtratio * dtinc, dtmax)
302 : 2 : alpha = alpha * alphadec
303 : : end if
304 : : else
305 : : ! reset downhill counter, velocity, alpha. decrease dtratio.
306 : 1 : ndownhill=0
307 [ + + ]: 28 : vel_ioncell(:)=0.0
308 : 1 : alpha=alpha0
309 : 1 : dtratio = dtratio*dtdec
310 : : endif
311 : :
312 : 11 : 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 [ - + + + ]: 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 [ - + + + ]: 319 : 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 : 11 : 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 [ + + ]: 110 : do kk=1,ab_mover%natom
339 [ + + ]: 407 : do jj=1,3
340 : : ! Warning : implemented in reduced coordinates
341 [ - + ]: 396 : if ( ab_mover%iatfix(jj,kk)==1) then
342 : 0 : 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 [ - + ]: 11 : if ( etotal - etotal_prev >0.0 ) then
350 [ # # # # ]: 0 : 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 [ + + ]: 308 : vin_prev(:)=vin(:)
357 : : ! endif
358 : :
359 : :
360 : :
361 : : !##########################################################
362 : : !### 08. update hist.
363 : :
364 : : !Increase indexes
365 : 11 : 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 [ + - - - : 11 : & vin, xred)
- - - - -
- - + ]
371 : :
372 [ - + ]: 11 : if(ab_mover%optcell/=0)then
373 : 0 : call mkrdim(acell,rprim,rprimd)
374 : 0 : 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 : 11 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
380 : 11 : ihist_prev = abihist_findIndex(hist,-1)
381 [ + + + + ]: 407 : hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
382 : :
383 [ - + ]: 11 : if(zDEBUG)then
384 : 0 : write (std_out,*) 'residual:'
385 [ # # ]: 0 : do kk=1,ab_mover%natom
386 : 0 : write (std_out,*) residual(:,kk)
387 : : end do
388 : 0 : write (std_out,*) 'strten:'
389 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
390 : 0 : write (std_out,*) 'etotal:'
391 : 0 : write (std_out,*) etotal
392 : : end if
393 : :
394 : :
395 : : end subroutine pred_fire
396 : : !!***
397 : :
398 : : end module m_pred_fire
399 : :
|
