Branch data Line data Source code
1 : : !!****m* ABINIT/m_pred_bfgs
2 : : !! NAME
3 : : !! m_pred_bfgs
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, JCC, SE, FB)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : module m_pred_bfgs
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_abimover
26 : : use m_abihist
27 : : use m_xfpack
28 : : use m_lbfgs
29 : : use m_errors
30 : :
31 : : use m_geometry, only : mkrdim, fcart2gred, metric
32 : : use m_bfgs, only : hessinit, hessupdt, brdene
33 : :
34 : : implicit none
35 : :
36 : : private
37 : : !!***
38 : :
39 : : public :: pred_bfgs
40 : : public :: pred_lbfgs
41 : : !!***
42 : :
43 : : contains
44 : : !!***
45 : :
46 : : !!****f* ABINIT/pred_bfgs
47 : : !! NAME
48 : : !! pred_bfgs
49 : : !!
50 : : !! FUNCTION
51 : : !! Ionmov predictors (2 & 3) Broyden-Fletcher-Goldfarb-Shanno
52 : : !!
53 : : !! IONMOV 2:
54 : : !! Given a starting point xred that is a vector of length 3*natom
55 : : !! (reduced nuclei coordinates), and unit cell parameters
56 : : !! (acell and rprimd) the Broyden-Fletcher-Goldfarb-Shanno
57 : : !! minimization is performed on the total energy function, using
58 : : !! its gradient (atomic forces and stresse) as calculated
59 : : !! by the routine scfcv. Some atoms can be kept fixed,
60 : : !! while the optimization of unit cell parameters is
61 : : !! only performed if optcell/=0. The convergence requirement on
62 : : !! the atomic forces, dtset%tolmxf, allows an early exit.
63 : : !! Otherwise no more than dtset%ntime steps are performed.
64 : : !! Returned quantities are xred, and eventually acell and rprim (new ones!).
65 : : !! Could see Numerical Recipes (Fortran), 1986, page 307.
66 : : !!
67 : : !! IONMOV 3:
68 : : !! Conduct structural optimization using the Broyden-Fletcher-
69 : : !! Goldfarb-Shanno minimization (BFGS), modified to take into
70 : : !! account the total energy as well as the gradients (as in usual
71 : : !! BFGS). See the paper by Schlegel, J. Comp. Chem. 3, 214 (1982) [[cite:Schlegel1982]].
72 : : !! Might be better than ionmov=2 for few degrees of freedom (less than 3 or 4)
73 : : !!
74 : : !! INPUTS
75 : : !! ab_mover <type(abimover)> : Datatype with all the information
76 : : !! needed by the preditor
77 : : !! itime : Index of the present iteration
78 : : !! ntime : Maximal number of iterations
79 : : !! ionmov : (2 or 3) Specific kind of BFGS
80 : : !! zDEBUG : if true print some debugging information
81 : : !!
82 : : !! OUTPUT
83 : : !!
84 : : !! SIDE EFFECTS
85 : : !! hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
86 : : !!
87 : : !! SOURCE
88 : :
89 : 1244 : subroutine pred_bfgs(ab_mover,ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit)
90 : :
91 : : !Arguments ------------------------------------
92 : : !scalars
93 : : type(abimover),intent(in) :: ab_mover
94 : : type(ab_xfh_type),intent(inout) :: ab_xfh
95 : : type(abihist),intent(inout) :: hist
96 : : type(abiforstr),intent(in) :: forstr
97 : : integer,intent(in) :: itime
98 : : integer,intent(in) :: ionmov
99 : : integer,intent(in) :: iexit
100 : : logical,intent(in) :: zDEBUG
101 : :
102 : : !Local variables-------------------------------
103 : : !scalars
104 : : integer :: ihist_prev,ndim,cycl_main
105 : : integer, parameter :: npul=0
106 : : integer :: ierr,ii,jj,kk,nitpul
107 : : real(dp),save :: ucvol0
108 : : real(dp) :: ucvol,det
109 : : real(dp) :: etotal,etotal_prev
110 : : real(dp) :: favg,alpha0
111 : :
112 : : !arrays
113 : 1244 : integer,allocatable :: ipiv(:)
114 : : real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:)
115 : : real(dp),allocatable,save :: vout(:),vout_prev(:)
116 : : real(dp),allocatable,save ::vinres(:,:),vin1(:,:)
117 : 1244 : real(dp),allocatable :: amat(:,:),amatinv(:,:),alpha(:,:)
118 : 1244 : real(dp),allocatable :: rwork(:)
119 : : real(dp),save :: acell0(3) ! Initial acell
120 : : real(dp),save :: rprimd0(3,3) ! Initial rprimd
121 : : real(dp) :: acell(3)
122 : : real(dp) :: rprimd(3,3),rprim(3,3)
123 : : real(dp) :: gprimd(3,3)
124 : : real(dp) :: gmet(3,3)
125 : : real(dp) :: rmet(3,3)
126 : 2488 : real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom)
127 : 1244 : real(dp) :: xred(3,ab_mover%natom),strten(6)
128 : :
129 : : !***************************************************************************
130 : : !Beginning of executable session
131 : : !***************************************************************************
132 : :
133 [ + + ]: 1244 : if(iexit/=0)then
134 [ + + ]: 306 : ABI_SFREE(vin)
135 [ + + ]: 306 : ABI_SFREE(vout)
136 [ + + ]: 306 : ABI_SFREE(vin_prev)
137 [ + + ]: 306 : ABI_SFREE(vout_prev)
138 [ - + ]: 306 : ABI_SFREE(vinres)
139 [ - + ]: 306 : ABI_SFREE(vin1)
140 [ + + ]: 306 : ABI_SFREE(hessin)
141 : : return
142 : : end if
143 : :
144 : : !write(std_out,*) 'bfgs 01'
145 : : !##########################################################
146 : : !### 01. Debugging and Verbose
147 : :
148 [ - + ]: 938 : if(zDEBUG)then
149 [ # # # # : 0 : write(std_out,'(a,3a,35a,42a)') ch10,('-',kk=1,3),'Debugging and Verbose for pred_bfgs',('-',kk=1,42)
# # # # ]
150 : 0 : write(std_out,*) 'ionmov: ',ionmov
151 : 0 : write(std_out,*) 'itime: ',itime
152 : : end if
153 : :
154 : : !write(std_out,*) 'bfgs 02'
155 : : !##########################################################
156 : : !### 02. Compute the dimension of vectors (ndim)
157 : :
158 : 938 : ndim=3*ab_mover%natom
159 [ + + ]: 938 : if(ab_mover%optcell==1) ndim=ndim+1
160 [ + + ]: 938 : if(ab_mover%optcell==2 .or.&
161 : 129 : & ab_mover%optcell==3) ndim=ndim+6
162 [ + + ]: 938 : if(ab_mover%optcell>=4) ndim=ndim+3
163 : :
164 [ - + ]: 938 : if(zDEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim
165 : :
166 : : !write(std_out,*) 'bfgs 03'
167 : : !##########################################################
168 : : !### 03. Allocate the vectors vin, vout and hessian matrix
169 : :
170 : : !Notice that vin, vout, etc could be allocated
171 : : !From a previous dataset with a different ndim
172 [ + + ]: 938 : if(itime==1)then
173 [ - + ]: 286 : ABI_SFREE(vin)
174 [ - + ]: 286 : ABI_SFREE(vout)
175 [ - + ]: 286 : ABI_SFREE(vin_prev)
176 [ - + ]: 286 : ABI_SFREE(vout_prev)
177 [ - + ]: 286 : ABI_SFREE(vinres)
178 [ - + ]: 286 : ABI_SFREE(vin1)
179 [ - + ]: 286 : ABI_SFREE(hessin)
180 : : if(npul>1) then
181 : : ABI_MALLOC(vinres,(npul+1,ndim))
182 : : ABI_MALLOC(vin1,(npul+1,ndim))
183 : : end if
184 [ - + + - : 858 : ABI_MALLOC(vin,(ndim))
+ - + - -
+ - + -
+ ]
185 [ - + + - : 572 : ABI_MALLOC(vout,(ndim))
- + - + -
+ ]
186 [ - + + - : 572 : ABI_MALLOC(vin_prev,(ndim))
- + - + -
+ ]
187 [ - + + - : 572 : ABI_MALLOC(vout_prev,(ndim))
- + - + -
+ ]
188 [ - + + - : 1144 : ABI_MALLOC(hessin,(ndim,ndim))
- + + - +
- + - - +
- + - + ]
189 : : end if
190 : :
191 : : !write(std_out,*) 'bfgs 04'
192 : : !##########################################################
193 : : !### 04. Obtain the present values from the history
194 : :
195 : 938 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
196 [ + + ]: 3752 : do ii=1,3
197 [ + + ]: 12194 : rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
198 : : end do
199 : :
200 [ + + ]: 6566 : strten(:)=hist%strten(:,hist%ihist)
201 : 938 : etotal =hist%etot(hist%ihist)
202 : :
203 : : !Fill the residual with forces (No preconditioning)
204 : : !Or the preconditioned forces
205 [ + + ]: 938 : if (ab_mover%goprecon==0)then
206 : 933 : call fcart2gred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom)
207 : : else
208 [ + + + + ]: 45 : residual(:,:)=forstr%gred(:,:)
209 : : end if
210 : :
211 [ - + ]: 938 : if(zDEBUG)then
212 : 0 : write (std_out,*) 'residual:'
213 [ # # ]: 0 : do kk=1,ab_mover%natom
214 : 0 : write (std_out,*) residual(:,kk)
215 : : end do
216 : 0 : write (std_out,*) 'strten:'
217 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
218 : 0 : write (std_out,*) 'etotal:'
219 : 0 : write (std_out,*) etotal
220 : : end if
221 : :
222 : 938 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
223 : :
224 : : !Save initial values
225 [ + + ]: 938 : if (itime==1)then
226 : 286 : acell0(:)=acell(:)
227 : 286 : rprimd0(:,:)=rprimd(:,:)
228 : 286 : ucvol0=ucvol
229 : : end if
230 : :
231 : : !zDEBUG (UCVOL)
232 [ - + ]: 938 : if(zDEBUG)then
233 : 0 : write(std_out,*) 'Volume of cell (ucvol):',ucvol
234 : : end if
235 : :
236 : : !Get rid of mean force on whole unit cell, but only if no
237 : : !generalized constraints are in effect
238 [ + + + + ]: 12182 : residual_corrected(:,:)=residual(:,:)
239 [ + + ]: 938 : if(ab_mover%nconeq==0)then
240 [ + + ]: 3580 : do ii=1,3
241 [ + + + - ]: 3580 : if (ii/=3.or.ab_mover%jellslab==0) then
242 [ + + ]: 10860 : favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom)
243 [ + + ]: 10860 : residual_corrected(ii,:)=residual_corrected(ii,:)-favg
244 : : end if
245 : : end do
246 : : end if
247 : :
248 : : !write(std_out,*) 'bfgs 05'
249 : : !##########################################################
250 : : !### 05. Fill the vectors vin and vout
251 : :
252 : : !Initialize input vectors : first vin, then vout
253 : : !The values of vin from the previous iteration
254 : : !should be the same
255 : : !if (itime==1)then
256 : : call xfpack_x2vin(acell, ab_mover%natom, ndim,&
257 : : & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
258 [ + - - - : 938 : & ab_mover%symrel, ucvol, ucvol0, vin, xred)
- - - - -
- - + ]
259 : : !end if
260 : :
261 : : call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,&
262 : : & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,&
263 [ + - - - : 938 : & vout)
- - - + ]
264 : :
265 : : !write(std_out,*) 'bfgs 06'
266 : : !##########################################################
267 : : !### 06. Initialize or update the hessian matrix
268 : :
269 : : !Initialise the Hessian matrix using gmet
270 [ + + ]: 938 : if (itime==1)then
271 : :
272 : 286 : call hessinit(ab_mover, hessin, gmet, ndim, ucvol)
273 : :
274 : : ! ! Initialize inverse hessian with identity matrix
275 : : ! ! in cartesian coordinates, which makes use of metric tensor gmet
276 : : ! ! in reduced coordinates.
277 : : ! hessin(:,:)=zero
278 : : ! do ii=1,ab_mover%natom
279 : : ! do kk=1,3
280 : : ! do jj=1,3
281 : : ! ! Warning : implemented in reduced coordinates
282 : : ! if (ab_mover%iatfix(kk,ii)==0 .and.&
283 : : ! & ab_mover%iatfix(jj,ii)==0 )then
284 : : ! hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj)
285 : : ! end if
286 : : ! end do
287 : : ! end do
288 : : ! end do
289 : : ! if(ab_mover%optcell/=0)then
290 : : ! ! These values might lead to too large changes in some cases
291 : : ! diag=ab_mover%strprecon*30.0_dp/ucvol
292 : : ! if(ab_mover%optcell==1) diag=diag/three
293 : : ! do ii=3*ab_mover%natom+1,ndim
294 : : ! hessin(ii,ii)=diag
295 : : ! end do
296 : : ! end if
297 : :
298 [ + + ]: 286 : if (ab_mover%restartxf/=0) then
299 : :
300 : : call xfh_recover_new(ab_xfh,ab_mover,acell,cycl_main,residual,&
301 : : & hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,&
302 : 7 : & vin_prev,vout,vout_prev,xred)
303 : :
304 : : end if
305 : :
306 : : end if
307 : :
308 [ + + ]: 938 : if(itime>1)then
309 : : ! Update the hessian matrix, by taking into account the
310 : : ! current pair (x,f) and the previous one.
311 : : call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom,ndim,vin,&
312 [ + - - - : 652 : & vin_prev,vout,vout_prev)
- - - - -
+ ]
313 : :
314 : : end if
315 : :
316 : : !zDEBUG (vin,vout and hessin before prediction)
317 [ - + ]: 938 : if(zDEBUG)then
318 : 0 : write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]'
319 : 0 : write(std_out,*) 'vin:'
320 [ # # # # ]: 0 : do ii=1,ndim,3
321 [ # # ]: 0 : if (ii+2<=ndim)then
322 : 0 : write(std_out,*) ii,vin(ii:ii+2)
323 : : else
324 : 0 : write(std_out,*) ii,vin(ii:ndim)
325 : : end if
326 : : end do
327 : 0 : write(std_out,*) 'vout:'
328 [ # # # # ]: 0 : do ii=1,ndim,3
329 [ # # ]: 0 : if (ii+2<=ndim)then
330 : 0 : write(std_out,*) ii,vout(ii:ii+2)
331 : : else
332 : 0 : write(std_out,*) ii,vout(ii:ndim)
333 : : end if
334 : : end do
335 : 0 : write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim
336 [ # # ]: 0 : do kk=1,ndim
337 [ # # # # ]: 0 : do jj=1,ndim,3
338 [ # # ]: 0 : if (jj+2<=ndim)then
339 : 0 : write(std_out,*) jj,hessin(jj:jj+2,kk)
340 : : else
341 : 0 : write(std_out,*) jj,hessin(jj:ndim,kk)
342 : : end if
343 : : end do
344 : : end do
345 : : end if
346 : :
347 : : !write(std_out,*) 'bfgs 07'
348 : : !##########################################################
349 : : !### 07. Compute the next values
350 : :
351 [ + + + + ]: 938 : if(ionmov==2 .or. itime==1)then
352 : :
353 : : ! Previous cartesian coordinates
354 [ + + ]: 8646 : vin_prev(:)=vin(:)
355 : :
356 : : ! New atomic cartesian coordinates are obtained from vin, hessin
357 : : ! and vout
358 : :
359 [ + - + - : 729987 : vin(:) = vin(:) - matmul(hessin(:,:), vout(:))
+ - - + -
+ + + + +
+ + + - +
- + - + +
+ - ]
360 : :
361 : : !Pulay mixing for vin
362 : 691 : nitpul=0
363 : : if (npul>1) then
364 : : alpha0=1.0_dp
365 : : nitpul=min(itime, npul)
366 : : if (itime>npul) then
367 : : do jj=1,npul-1
368 : : vinres(jj,:)=vinres(jj+1,:)
369 : : vin1(jj,:)=vin1(jj+1,:)
370 : : end do
371 : : end if
372 : : vinres(nitpul,:)=vin(:)-vin_prev(:)
373 : : vin1(nitpul,:)=vin_prev(:)
374 : : end if
375 : :
376 : : if (nitpul>1) then
377 : : ABI_MALLOC(alpha,(nitpul,ndim))
378 : : alpha=zero
379 : : do kk=1,ndim
380 : : ABI_MALLOC(amat,(nitpul,nitpul))
381 : : ABI_MALLOC(amatinv,(nitpul,nitpul))
382 : : amat=zero;amatinv=zero
383 : : do ii=1,nitpul
384 : : do jj=ii,nitpul
385 : : amat(ii,jj)=vinres(jj,kk)*vinres(ii,kk)
386 : : amat(jj,ii)=amat(ii,jj)
387 : : end do
388 : : end do
389 : : amatinv=amat
390 : : if (abs(vin(kk)-vin_prev(kk))<tol10) then
391 : : alpha(:,kk)=zero
392 : : else
393 : : ABI_MALLOC(ipiv,(nitpul))
394 : : ABI_MALLOC(rwork,(nitpul))
395 : : ! amatinv=1.d5*amatinv
396 : : call dgetrf(nitpul,nitpul,amatinv,nitpul,ipiv,ierr)
397 : : call dgetri(nitpul,amatinv,nitpul,ipiv,rwork,nitpul,ierr)
398 : : ! amatinv=1.d5*amatinv
399 : : ABI_FREE(ipiv)
400 : : ABI_FREE(rwork)
401 : : det=zero
402 : : do ii=1,nitpul
403 : : do jj=1,nitpul
404 : : alpha(ii,kk)=alpha(ii,kk)+amatinv(jj,ii)
405 : : det=det+amatinv(jj,ii)
406 : : end do
407 : : end do
408 : : alpha(:,kk)=alpha(:,kk)/det
409 : : end if
410 : : end do
411 : : ABI_FREE(amat)
412 : : ABI_FREE(amatinv)
413 : : vin(:)=vin1(nitpul,:)+alpha0*(vin1(nitpul+1,:)-vin1(nitpul,:))
414 : : vin=zero
415 : : do ii=1,nitpul
416 : : vin(:)=vin(:)+ alpha(ii,:)*(vin1(ii,:))
417 : : end do
418 : : ABI_FREE(alpha)
419 : : end if
420 : :
421 : :
422 : : ! Previous atomic forces
423 [ + + ]: 8646 : vout_prev(:)=vout(:)
424 : :
425 : : else
426 [ + - ]: 247 : if(ionmov==3)then
427 : 247 : ihist_prev = abihist_findIndex(hist,-1)
428 : 247 : etotal_prev=hist%etot(ihist_prev)
429 : : ! Here the BFGS algorithm, modified to take into account the energy
430 : 247 : call brdene(etotal,etotal_prev,hessin,ndim,vin,vin_prev,vout,vout_prev)
431 : : end if
432 : :
433 : : ! zDEBUG (vin,vout and hessin after prediction)
434 [ - + ]: 247 : if(zDEBUG)then
435 : 0 : write(std_out,*) 'Vectors vin and vout [after prediction]'
436 : 0 : write(std_out,*) 'vin:'
437 [ # # # # ]: 0 : do ii=1,ndim,3
438 [ # # ]: 0 : if (ii+2<=ndim)then
439 : 0 : write(std_out,*) ii,vin(ii:ii+2)
440 : : else
441 : 0 : write(std_out,*) ii,vin(ii:ndim)
442 : : end if
443 : : end do
444 : 0 : write(std_out,*) 'vout:'
445 [ # # # # ]: 0 : do ii=1,ndim,3
446 [ # # ]: 0 : if (ii+2<=ndim)then
447 : 0 : write(std_out,*) ii,vout(ii:ii+2)
448 : : else
449 : 0 : write(std_out,*) ii,vout(ii:ndim)
450 : : end if
451 : : end do
452 : : end if
453 : :
454 : :
455 : : ! FIXME: this should be outside the if clause on ionmov!
456 : : ! Implement fixing of atoms : put back old values for fixed
457 : : ! components
458 [ + + ]: 731 : do kk=1,ab_mover%natom
459 [ + + ]: 2183 : do jj=1,3
460 : : ! Warning : implemented in reduced coordinates
461 [ + + ]: 1936 : if ( ab_mover%iatfix(jj,kk)==1) then
462 : 234 : vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3)
463 : : end if
464 : : end do
465 : : end do
466 : : end if
467 : :
468 : : !write(std_out,*) 'bfgs 08'
469 : : !##########################################################
470 : : !### 08. Update the history with the prediction
471 : :
472 : : !Increase indexes
473 : 938 : hist%ihist = abihist_findIndex(hist,+1)
474 : :
475 : : !Transfer vin to xred, acell and rprim
476 : : call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,&
477 : : & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
478 : : & ab_mover%symrel, ucvol, ucvol0,&
479 [ + - - - : 938 : & vin, xred)
- - - - -
- - + ]
480 : :
481 [ + + ]: 938 : if(ab_mover%optcell/=0)then
482 : 318 : call mkrdim(acell,rprim,rprimd)
483 : 318 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
484 : : end if
485 : :
486 : : !Fill the history with the variables
487 : : !xred, acell, rprimd, vel
488 : 938 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
489 : 938 : ihist_prev = abihist_findIndex(hist,-1)
490 [ + + + + ]: 12182 : hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
491 : :
492 [ - + ]: 938 : if(zDEBUG)then
493 : 0 : write (std_out,*) 'residual:'
494 [ # # ]: 0 : do kk=1,ab_mover%natom
495 : 0 : write (std_out,*) residual(:,kk)
496 : : end do
497 : 0 : write (std_out,*) 'strten:'
498 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
499 : 0 : write (std_out,*) 'etotal:'
500 : 0 : write (std_out,*) etotal
501 : : end if
502 : :
503 : : end subroutine pred_bfgs
504 : : !!***
505 : :
506 : : !!****f* ABINIT/pred_lbfgs
507 : : !! NAME
508 : : !! pred_lbfgs
509 : : !!
510 : : !! FUNCTION
511 : : !! Ionmov predictors (22) Limited-memory Broyden-Fletcher-Goldfarb-Shanno
512 : : !!
513 : : !! IONMOV 22:
514 : : !! Given a starting point xred that is a vector of length 3*natom
515 : : !! (reduced nuclei coordinates), and unit cell parameters
516 : : !! (acell and rprim) the L-Broyden-Fletcher-Goldfarb-Shanno
517 : : !! minimization is performed on the total energy function, using
518 : : !! its gradient (atomic forces and stress : gred or fcart and
519 : : !! stress) as calculated by the routine scfcv. Some atoms can be
520 : : !! kept fixed, while the optimization of unit cell parameters is
521 : : !! only performed if optcell/=0. The convergence requirement on
522 : : !! the atomic forces, dtset%tolmxf, allows an early exit.
523 : : !! Otherwise no more than dtset%ntime steps are performed.
524 : : !! Returned quantities are xred, and eventually acell and rprim (new ones!).
525 : : !! Could see MinPack on netlib.org
526 : : !!
527 : : !! INPUTS
528 : : !! ab_mover <type(abimover)> : Datatype with all the information needed by the preditor
529 : : !! itime : Index of the present iteration
530 : : !! ntime : Maximal number of iterations
531 : : !! ionmov : (22) Specific kind of BFGS
532 : : !! zDEBUG : if true print some debugging information
533 : : !!
534 : : !! OUTPUT
535 : : !!
536 : : !! SIDE EFFECTS
537 : : !! hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
538 : : !!
539 : : !! SOURCE
540 : :
541 : 12 : subroutine pred_lbfgs(ab_mover,ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit)
542 : :
543 : : !Arguments ------------------------------------
544 : : !scalars
545 : : type(abimover),intent(in) :: ab_mover
546 : : type(ab_xfh_type),intent(inout) :: ab_xfh
547 : : type(abihist),intent(inout) :: hist
548 : : type(abiforstr),intent(in) :: forstr
549 : : integer,intent(in) :: itime
550 : : integer,intent(in) :: ionmov
551 : : integer,intent(in) :: iexit
552 : : logical,intent(in) :: zDEBUG
553 : : character(len=500) :: ionmov22_errmsg
554 : :
555 : : !Local variables-------------------------------
556 : : !scalars
557 : : integer :: info,ihist_prev
558 : : integer :: ndim,cycl_main
559 : : integer, parameter :: npul=0
560 : : integer :: ii,jj,kk
561 : : real(dp),save :: ucvol0
562 : : real(dp) :: ucvol
563 : : real(dp) :: etotal
564 : : real(dp) :: favg
565 : :
566 : : !arrays
567 : 12 : real(dp),allocatable :: diag(:)
568 : : real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:)
569 : : real(dp),allocatable,save :: vout(:),vout_prev(:)
570 : : real(dp),allocatable,save ::vinres(:,:),vin1(:,:)
571 : : real(dp),save :: acell0(3) ! Initial acell
572 : : real(dp),save :: rprimd0(3,3) ! Initial rprimd
573 : : real(dp) :: acell(3)
574 : : real(dp) :: rprimd(3,3),rprim(3,3)
575 : : real(dp) :: gprimd(3,3)
576 : : real(dp) :: gmet(3,3)
577 : : real(dp) :: rmet(3,3)
578 : 24 : real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom)
579 : 12 : real(dp) :: xred(3,ab_mover%natom)
580 : : real(dp) :: strten(6)
581 : :
582 : : !***************************************************************************
583 : : !Beginning of executable session
584 : : !***************************************************************************
585 : :
586 [ + + ]: 12 : if(iexit/=0)then
587 : 1 : call lbfgs_destroy()
588 [ + - ]: 1 : ABI_SFREE(vin)
589 [ + - ]: 1 : ABI_SFREE(vout)
590 [ + - ]: 1 : ABI_SFREE(vin_prev)
591 [ + - ]: 1 : ABI_SFREE(vout_prev)
592 [ - + ]: 1 : ABI_SFREE(vinres)
593 [ - + ]: 1 : ABI_SFREE(vin1)
594 [ + - ]: 1 : ABI_SFREE(hessin)
595 : : return
596 : : end if
597 : :
598 : : !write(std_out,*) 'bfgs 01'
599 : : !##########################################################
600 : : !### 01. Debugging and Verbose
601 : :
602 [ - + ]: 11 : if(zDEBUG)then
603 [ # # # # : 0 : write(std_out,'(a,3a,35a,42a)') ch10,('-',kk=1,3),'Debugging and Verbose for pred_bfgs',('-',kk=1,42)
# # # # ]
604 : 0 : write(std_out,*) 'ionmov: ',ionmov
605 : 0 : write(std_out,*) 'itime: ',itime
606 : : end if
607 : :
608 : : !write(std_out,*) 'bfgs 02'
609 : : !##########################################################
610 : : !### 02. Compute the dimension of vectors (ndim)
611 : :
612 : 11 : ndim=3*ab_mover%natom
613 [ - + ]: 11 : if(ab_mover%optcell==1) ndim=ndim+1
614 [ - + ]: 11 : if(ab_mover%optcell==2 .or.&
615 : 0 : & ab_mover%optcell==3) ndim=ndim+6
616 [ - + ]: 11 : if(ab_mover%optcell>=4) ndim=ndim+3
617 : :
618 [ - + ]: 11 : if(zDEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim
619 : :
620 : : !write(std_out,*) 'bfgs 03'
621 : : !##########################################################
622 : : !### 03. Allocate the vectors vin, vout and hessian matrix
623 : :
624 : : !Notice that vin, vout, etc could be allocated
625 : : !From a previous dataset with a different ndim
626 [ + + ]: 11 : if(itime==1)then
627 [ - + ]: 1 : ABI_SFREE(vin)
628 [ - + ]: 1 : ABI_SFREE(vout)
629 [ - + ]: 1 : ABI_SFREE(vin_prev)
630 [ - + ]: 1 : ABI_SFREE(vout_prev)
631 [ - + ]: 1 : ABI_SFREE(vinres)
632 [ - + ]: 1 : ABI_SFREE(vin1)
633 [ - + ]: 1 : ABI_SFREE(hessin)
634 : : if(npul>1) then
635 : : ABI_MALLOC(vinres,(npul+1,ndim))
636 : : ABI_MALLOC(vin1,(npul+1,ndim))
637 : : end if
638 [ - + + - : 3 : ABI_MALLOC(vin,(ndim))
+ - + - -
+ - + -
+ ]
639 [ - + + - : 2 : ABI_MALLOC(vout,(ndim))
- + - + -
+ ]
640 [ - + + - : 2 : ABI_MALLOC(vin_prev,(ndim))
- + - + -
+ ]
641 [ - + + - : 2 : ABI_MALLOC(vout_prev,(ndim))
- + - + -
+ ]
642 [ - + + - : 4 : ABI_MALLOC(hessin,(ndim,ndim))
- + + - +
- + - - +
- + - + ]
643 : : end if
644 : :
645 : : !write(std_out,*) 'bfgs 04'
646 : : !##########################################################
647 : : !### 04. Obtain the present values from the history
648 : :
649 : 11 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
650 [ + + ]: 44 : do ii=1,3
651 [ + + ]: 143 : rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
652 : : end do
653 : :
654 [ + + ]: 77 : strten(:)=hist%strten(:,hist%ihist)
655 : 11 : etotal =hist%etot(hist%ihist)
656 : :
657 : : !Fill the residual with forces (No preconditioning)
658 : : !Or the preconditioned forces
659 [ + - ]: 11 : if (ab_mover%goprecon==0)then
660 : 11 : call fcart2gred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom)
661 : : else
662 [ # # # # ]: 0 : residual(:,:)= forstr%gred(:,:)
663 : : end if
664 : :
665 [ - + ]: 11 : if(zDEBUG)then
666 : 0 : write (std_out,*) 'residual:'
667 [ # # ]: 0 : do kk=1,ab_mover%natom
668 : 0 : write (std_out,*) residual(:,kk)
669 : : end do
670 : 0 : write (std_out,*) 'strten:'
671 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
672 : 0 : write (std_out,*) 'etotal:'
673 : 0 : write (std_out,*) etotal
674 : : end if
675 : :
676 : 11 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
677 : :
678 : : !Save initial values
679 [ + + ]: 11 : if (itime==1)then
680 : 1 : acell0(:)=acell(:)
681 : 1 : rprimd0(:,:)=rprimd(:,:)
682 : 1 : ucvol0=ucvol
683 : : end if
684 : :
685 : : !zDEBUG (UCVOL)
686 [ - + ]: 11 : if(zDEBUG)then
687 : 0 : write(std_out,*) 'Volume of cell (ucvol):',ucvol
688 : : end if
689 : :
690 : : !Get rid of mean force on whole unit cell, but only if no
691 : : !generalized constraints are in effect
692 [ + + + + ]: 407 : residual_corrected(:,:)=residual(:,:)
693 [ + - ]: 11 : if(ab_mover%nconeq==0)then
694 [ + + ]: 44 : do ii=1,3
695 [ + + + - ]: 44 : if (ii/=3.or.ab_mover%jellslab==0) then
696 [ + + ]: 330 : favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom)
697 [ + + ]: 330 : residual_corrected(ii,:)=residual_corrected(ii,:)-favg
698 : : end if
699 : : end do
700 : : end if
701 : :
702 : : !write(std_out,*) 'bfgs 05'
703 : : !##########################################################
704 : : !### 05. Fill the vectors vin and vout
705 : :
706 : : !Initialize input vectors : first vin, then vout
707 : : !The values of vin from the previous iteration
708 : : !should be the same
709 : : !if (itime==1)then
710 : : call xfpack_x2vin(acell, ab_mover%natom, ndim,&
711 : : & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
712 [ + - - - : 11 : & ab_mover%symrel, ucvol, ucvol0, vin, xred)
- - - - -
- - + ]
713 : : !end if
714 : :
715 : : call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,&
716 : : & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,&
717 [ + - - - : 11 : & vout)
- - - + ]
718 : :
719 : : !write(std_out,*) 'bfgs 06'
720 : : !##########################################################
721 : : !### 06. Initialize or update the hessian matrix
722 : :
723 : : !Initialise the Hessian matrix using gmet
724 [ + + ]: 11 : if (itime==1)then
725 : :
726 [ - + + - : 3 : ABI_MALLOC(diag,(ndim))
+ - + - -
+ - + ]
727 [ + + ]: 28 : do ii=1,3*ab_mover%natom
728 : : !diag(ii) = 1.00_dp / rprimd(MODULO(ii-1,3)+1,MODULO(ii-1,3)+1)**2
729 : 28 : diag(ii) = gmet(MODULO(ii-1,3)+1,MODULO(ii-1,3)+1)
730 : : end do
731 [ - + ]: 1 : if(ab_mover%optcell/=0)then
732 : : ! These values might lead to too large changes in some cases ...
733 [ # # ]: 0 : do ii=3*ab_mover%natom+1,ndim
734 : 0 : diag(ii) = ab_mover%strprecon*30.0_dp/ucvol
735 [ # # ]: 0 : if(ab_mover%optcell==1) diag(ii) = diag(ii) / three
736 : : end do
737 : : end if
738 : :
739 : : !call lbfgs_destroy()
740 : 1 : call lbfgs_init(ndim,5,diag)
741 [ - + ]: 1 : ABI_FREE(diag)
742 : :
743 [ - + ]: 1 : if (ab_mover%restartxf/=0) then
744 : :
745 : : call xfh_recover_new(ab_xfh,ab_mover,acell,cycl_main,residual,&
746 : : hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,&
747 : 0 : vin_prev,vout,vout_prev,xred)
748 : :
749 : : end if
750 : :
751 : : end if
752 : :
753 : : !zDEBUG (vin,vout and hessin before prediction)
754 [ - + ]: 11 : if(zDEBUG)then
755 : 0 : write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]'
756 : 0 : write(std_out,*) 'vin:'
757 [ # # # # ]: 0 : do ii=1,ndim,3
758 [ # # ]: 0 : if (ii+2<=ndim)then
759 : 0 : write(std_out,*) ii,vin(ii:ii+2)
760 : : else
761 : 0 : write(std_out,*) ii,vin(ii:ndim)
762 : : end if
763 : : end do
764 : 0 : write(std_out,*) 'vout:'
765 [ # # # # ]: 0 : do ii=1,ndim,3
766 [ # # ]: 0 : if (ii+2<=ndim)then
767 : 0 : write(std_out,*) ii,vout(ii:ii+2)
768 : : else
769 : 0 : write(std_out,*) ii,vout(ii:ndim)
770 : : end if
771 : : end do
772 : 0 : write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim
773 [ # # ]: 0 : do kk=1,ndim
774 [ # # # # ]: 0 : do jj=1,ndim,3
775 [ # # ]: 0 : if (jj+2<=ndim)then
776 : 0 : write(std_out,*) jj,hessin(jj:jj+2,kk)
777 : : else
778 : 0 : write(std_out,*) jj,hessin(jj:ndim,kk)
779 : : end if
780 : : end do
781 : : end do
782 : : end if
783 : :
784 : : !write(std_out,*) 'bfgs 07'
785 : : !##########################################################
786 : : !### 07. Compute the next values
787 : :
788 [ + + ]: 308 : vin_prev(:) = vin
789 [ + + ]: 308 : vout_prev(:) = vout
790 : 11 : info = lbfgs_execute(vin,etotal,vout)
791 : :
792 [ - + ]: 11 : if (info /= -1) then
793 : : write (ionmov22_errmsg, '(a,i0,3a)') &
794 : 0 : 'Lbfgs routine failed. Returned value: ', info,ch10, &
795 : 0 : 'Restart your calculation from last step or try a different ionmov'
796 : 0 : ABI_ERROR_CLASS(ionmov22_errmsg, "Ionmov22Error")
797 : : end if
798 : :
799 : : !zDEBUG (vin,vout after prediction)
800 [ - + ]: 11 : if(zDEBUG)then
801 : 0 : write(std_out,*) 'Vectors vin and vout [after prediction]'
802 : 0 : write(std_out,*) 'vin_prev:'
803 [ # # # # ]: 0 : do ii=1,ndim,3
804 [ # # ]: 0 : if (ii+2<=ndim)then
805 : 0 : write(std_out,*) ii,vin_prev(ii:ii+2)
806 : : else
807 : 0 : write(std_out,*) ii,vin_prev(ii:ndim)
808 : : end if
809 : : end do
810 : 0 : write(std_out,*) 'vin:'
811 [ # # # # ]: 0 : do ii=1,ndim,3
812 [ # # ]: 0 : if (ii+2<=ndim)then
813 : 0 : write(std_out,*) ii,vin(ii:ii+2)
814 : : else
815 : 0 : write(std_out,*) ii,vin(ii:ndim)
816 : : end if
817 : : end do
818 : 0 : write(std_out,*) 'vout:'
819 [ # # # # ]: 0 : do ii=1,ndim,3
820 [ # # ]: 0 : if (ii+2<=ndim)then
821 : 0 : write(std_out,*) ii,vout(ii:ii+2)
822 : : else
823 : 0 : write(std_out,*) ii,vout(ii:ndim)
824 : : end if
825 : : end do
826 : : end if
827 : :
828 : :
829 : : !Implement fixing of atoms : put back old values for fixed
830 : : !components
831 [ + + ]: 110 : do kk=1,ab_mover%natom
832 [ + + ]: 407 : do jj=1,3
833 : : ! Warning : implemented in reduced coordinates
834 [ - + ]: 396 : if ( ab_mover%iatfix(jj,kk)==1) then
835 : 0 : vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3)
836 : : end if
837 : : end do
838 : : end do
839 : :
840 : :
841 : : !write(std_out,*) 'bfgs 08'
842 : : !##########################################################
843 : : !### 08. Update the history with the prediction
844 : :
845 : : !Increase indexes
846 : 11 : hist%ihist = abihist_findIndex(hist,+1)
847 : :
848 : : !Transfer vin to xred, acell and rprim
849 : : call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,&
850 : : & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
851 : : & ab_mover%symrel, ucvol, ucvol0,&
852 [ + - - - : 11 : & vin, xred)
- - - - -
- - + ]
853 : :
854 [ - + ]: 11 : if(ab_mover%optcell/=0)then
855 : 0 : call mkrdim(acell,rprim,rprimd)
856 : 0 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
857 : : end if
858 : :
859 : : !Fill the history with the variables
860 : : !xcart, xred, acell, rprimd
861 : 11 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
862 : 11 : ihist_prev = abihist_findIndex(hist,-1)
863 [ + + + + ]: 407 : hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
864 : :
865 [ - + ]: 11 : if(zDEBUG)then
866 : 0 : write (std_out,*) 'residual:'
867 [ # # ]: 0 : do kk=1,ab_mover%natom
868 : 0 : write (std_out,*) residual(:,kk)
869 : : end do
870 : 0 : write (std_out,*) 'strten:'
871 : 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
872 : 0 : write (std_out,*) 'etotal:'
873 : 0 : write (std_out,*) etotal
874 : : end if
875 : :
876 : : end subroutine pred_lbfgs
877 : : !!***
878 : :
879 [ - + ]: 691 : end module m_pred_bfgs
880 : : !!***
|
