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