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 :
|
