TABLE OF CONTENTS


ABINIT/m_pred_bfgs [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_bfgs

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, JCC, SE, FB)
  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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_pred_bfgs
27 
28  use defs_basis
29  use m_abicore
30  use m_abimover
31  use m_abihist
32  use m_xfpack
33  use m_lbfgs
34 
35  use m_geometry,    only : mkrdim, fcart2fred, metric
36  use m_bfgs,        only : hessinit, hessupdt, brdene
37 
38  implicit none
39 
40  private

ABINIT/pred_bfgs [ Functions ]

[ Top ] [ Functions ]

NAME

 pred_bfgs

FUNCTION

 Ionmov predictors (2 & 3) Broyden-Fletcher-Goldfarb-Shanno

 IONMOV 2:
 Given a starting point xred that is a vector of length 3*natom
 (reduced nuclei coordinates), and unit cell parameters
 (acell and rprimd) the Broyden-Fletcher-Goldfarb-Shanno
 minimization is performed on the total energy function, using
 its gradient (atomic forces and stresse) as calculated
 by the routine scfcv. Some atoms can be kept fixed,
 while the optimization of unit cell parameters is
 only performed if optcell/=0. The convergence requirement on
 the atomic forces, dtset%tolmxf,  allows an early exit.
 Otherwise no more than dtset%ntime steps are performed.
 Returned quantities are xred, and eventually acell and rprim (new ones!).
 Could see Numerical Recipes (Fortran), 1986, page 307.

 IONMOV 3:
 Conduct structural optimization using the Broyden-Fletcher-
 Goldfarb-Shanno minimization (BFGS), modified to take into
 account the total energy as well as the gradients (as in usual
 BFGS). See the paper by Schlegel, J. Comp. Chem. 3, 214 (1982) [[cite:Schlegel1982]].
 Might be better than ionmov=2 for few degrees of freedom (less than 3 or 4)

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 : (2 or 3) Specific kind of BFGS
 zDEBUG : if true print some debugging information

OUTPUT

SIDE EFFECTS

 hist <type(abihist)> : History of positions,forces acell, rprimd, stresses

PARENTS

      mover

CHILDREN

      brdene,dgetrf,dgetri,fcart2fred,hessinit,hessupdt,hist2var,metric
      mkrdim,var2hist,xfh_recover_new,xfpack_f2vout,xfpack_vin2x,xfpack_x2vin

SOURCE

100 subroutine pred_bfgs(ab_mover,ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit)
101 
102 
103 !This section has been created automatically by the script Abilint (TD).
104 !Do not modify the following lines by hand.
105 #undef ABI_FUNC
106 #define ABI_FUNC 'pred_bfgs'
107 !End of the abilint section
108 
109 implicit none
110 
111 !Arguments ------------------------------------
112 !scalars
113 type(abimover),intent(in)       :: ab_mover
114 type(ab_xfh_type),intent(inout)    :: ab_xfh
115 type(abihist),intent(inout) :: hist
116 type(abiforstr),intent(in) :: forstr
117 integer,intent(in) :: itime
118 integer,intent(in) :: ionmov
119 integer,intent(in) :: iexit
120 logical,intent(in) :: zDEBUG
121 
122 !Local variables-------------------------------
123 !scalars
124 integer  :: ihist_prev,ndim,cycl_main
125 integer, parameter :: npul=0
126 integer  :: ierr,ii,jj,kk,nitpul
127 real(dp),save :: ucvol0
128 real(dp) :: ucvol,det
129 real(dp) :: etotal,etotal_prev
130 real(dp) :: favg,alpha0
131 
132 !arrays
133 integer,allocatable :: ipiv(:)
134 real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:)
135 real(dp),allocatable,save :: vout(:),vout_prev(:)
136 real(dp),allocatable,save ::vinres(:,:),vin1(:,:)
137 real(dp),allocatable ::  amat(:,:),amatinv(:,:),alpha(:,:)
138 real(dp),allocatable :: rwork(:)
139 real(dp),save :: acell0(3) ! Initial acell
140 real(dp),save :: rprimd0(3,3) ! Initial rprimd
141 real(dp) :: acell(3)
142 real(dp) :: rprimd(3,3),rprim(3,3)
143 real(dp) :: gprimd(3,3)
144 real(dp) :: gmet(3,3)
145 real(dp) :: rmet(3,3)
146 real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom)
147 real(dp) :: xred(3,ab_mover%natom),strten(6)
148 
149 !***************************************************************************
150 !Beginning of executable session
151 !***************************************************************************
152 
153  if(iexit/=0)then
154    if (allocated(vin))        then
155      ABI_DEALLOCATE(vin)
156    end if
157    if (allocated(vout))       then
158      ABI_DEALLOCATE(vout)
159    end if
160    if (allocated(vin_prev))   then
161      ABI_DEALLOCATE(vin_prev)
162    end if
163    if (allocated(vout_prev))  then
164      ABI_DEALLOCATE(vout_prev)
165    end if
166    if (allocated(vinres))  then
167      ABI_DEALLOCATE(vinres)
168    end if
169    if (allocated(vin1))  then
170      ABI_DEALLOCATE(vin1)
171    end if
172    if (allocated(hessin))     then
173      ABI_DEALLOCATE(hessin)
174    end if
175    return
176  end if
177 
178 !write(std_out,*) 'bfgs 01'
179 !##########################################################
180 !### 01. Debugging and Verbose
181 
182  if(zDEBUG)then
183    write(std_out,'(a,3a,35a,42a)') ch10,('-',kk=1,3),&
184 &   'Debugging and Verbose for pred_bfgs',('-',kk=1,42)
185    write(std_out,*) 'ionmov: ',ionmov
186    write(std_out,*) 'itime:  ',itime
187  end if
188 
189 !write(std_out,*) 'bfgs 02'
190 !##########################################################
191 !### 02. Compute the dimension of vectors (ndim)
192 
193  ndim=3*ab_mover%natom
194  if(ab_mover%optcell==1 .or.&
195 & ab_mover%optcell==4 .or.&
196 & ab_mover%optcell==5 .or.&
197 & ab_mover%optcell==6) ndim=ndim+1
198  if(ab_mover%optcell==2 .or.&
199 & ab_mover%optcell==3) ndim=ndim+6
200  if(ab_mover%optcell==7 .or.&
201 & ab_mover%optcell==8 .or.&
202 & ab_mover%optcell==9) ndim=ndim+3
203 
204  if(zDEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim
205 
206 !write(std_out,*) 'bfgs 03'
207 !##########################################################
208 !### 03. Allocate the vectors vin, vout and hessian matrix
209 
210 !Notice that vin, vout, etc could be allocated
211 !From a previous dataset with a different ndim
212  if(itime==1)then
213    if (allocated(vin))        then
214      ABI_DEALLOCATE(vin)
215    end if
216    if (allocated(vout))       then
217      ABI_DEALLOCATE(vout)
218    end if
219    if (allocated(vin_prev))   then
220      ABI_DEALLOCATE(vin_prev)
221    end if
222    if (allocated(vout_prev))  then
223      ABI_DEALLOCATE(vout_prev)
224    end if
225    if (allocated(vinres))  then
226      ABI_DEALLOCATE(vinres)
227    end if
228    if (allocated(vin1))  then
229      ABI_DEALLOCATE(vin1)
230    end if
231    if (allocated(hessin))     then
232      ABI_DEALLOCATE(hessin)
233    end if
234    if(npul>1) then
235      ABI_ALLOCATE(vinres,(npul+1,ndim))
236      ABI_ALLOCATE(vin1,(npul+1,ndim))
237    end if
238    ABI_ALLOCATE(vin,(ndim))
239    ABI_ALLOCATE(vout,(ndim))
240    ABI_ALLOCATE(vin_prev,(ndim))
241    ABI_ALLOCATE(vout_prev,(ndim))
242    ABI_ALLOCATE(hessin,(ndim,ndim))
243  end if
244 
245 !write(std_out,*) 'bfgs 04'
246 !##########################################################
247 !### 04. Obtain the present values from the history
248 
249  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
250  do ii=1,3
251    rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
252  end do
253 
254  strten(:)=hist%strten(:,hist%ihist)
255  etotal   =hist%etot(hist%ihist)
256 
257 !Fill the residual with forces (No preconditioning)
258 !Or the preconditioned forces
259  if (ab_mover%goprecon==0)then
260    call fcart2fred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom)
261  else
262    residual(:,:)=forstr%fred(:,:)
263  end if
264 
265  if(zDEBUG)then
266    write (std_out,*) 'residual:'
267    do kk=1,ab_mover%natom
268      write (std_out,*) residual(:,kk)
269    end do
270    write (std_out,*) 'strten:'
271    write (std_out,*) strten(1:3),ch10,strten(4:6)
272    write (std_out,*) 'etotal:'
273    write (std_out,*) etotal
274  end if
275 
276  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
277 
278 !Save initial values
279  if (itime==1)then
280    acell0(:)=acell(:)
281    rprimd0(:,:)=rprimd(:,:)
282    ucvol0=ucvol
283  end if
284 
285 !zDEBUG (UCVOL)
286  if(zDEBUG)then
287    write(std_out,*) 'Volume of cell (ucvol):',ucvol
288  end if
289 
290 !Get rid of mean force on whole unit cell, but only if no
291 !generalized constraints are in effect
292  residual_corrected(:,:)=residual(:,:)
293  if(ab_mover%nconeq==0)then
294    do ii=1,3
295      if (ii/=3.or.ab_mover%jellslab==0) then
296        favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom)
297        residual_corrected(ii,:)=residual_corrected(ii,:)-favg
298      end if
299    end do
300  end if
301 
302 !write(std_out,*) 'bfgs 05'
303 !##########################################################
304 !### 05. Fill the vectors vin and vout
305 
306 !Initialize input vectors : first vin, then vout
307 !The values of vin from the previous iteration
308 !should be the same
309 !if (itime==1)then
310  call xfpack_x2vin(acell, acell0, ab_mover%natom, ndim,&
311 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
312 & ab_mover%symrel, ucvol, ucvol0, vin, xred)
313 !end if
314 
315  call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,&
316 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,&
317 & vout)
318 
319 !write(std_out,*) 'bfgs 06'
320 !##########################################################
321 !### 06. Initialize or update the hessian matrix
322 
323 !Initialise the Hessian matrix using gmet
324  if (itime==1)then
325 
326    call hessinit(ab_mover, hessin, gmet, ndim, ucvol)
327 
328 !  ! Initialize inverse hessian with identity matrix
329 !  ! in cartesian coordinates, which makes use of metric tensor gmet
330 !  ! in reduced coordinates.
331 !  hessin(:,:)=zero
332 !  do ii=1,ab_mover%natom
333 !  do kk=1,3
334 !  do jj=1,3
335 !  ! Warning : implemented in reduced coordinates
336 !  if (ab_mover%iatfix(kk,ii)==0 .and.&
337 !  & ab_mover%iatfix(jj,ii)==0 )then
338 !  hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj)
339 !  end if
340 !  end do
341 !  end do
342 !  end do
343 !  if(ab_mover%optcell/=0)then
344 !  ! These values might lead to too large changes in some cases
345 !  diag=ab_mover%strprecon*30.0_dp/ucvol
346 !  if(ab_mover%optcell==1) diag=diag/three
347 !  do ii=3*ab_mover%natom+1,ndim
348 !  hessin(ii,ii)=diag
349 !  end do
350 !  end if
351 
352    if (ab_mover%restartxf/=0) then
353 
354      call xfh_recover_new(ab_xfh,ab_mover,acell,acell0,cycl_main,residual,&
355 &     hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,&
356 &     vin_prev,vout,vout_prev,xred)
357 
358    end if
359 
360  end if
361 
362  if(itime>1)then
363 !  Update the hessian matrix, by taking into account the
364 !  current pair (x,f) and the previous one.
365    call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom,ndim,vin,&
366 &   vin_prev,vout,vout_prev)
367 
368  end if
369 
370 !zDEBUG (vin,vout and hessin before prediction)
371  if(zDEBUG)then
372    write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]'
373    write(std_out,*) 'vin:'
374    do ii=1,ndim,3
375      if (ii+2<=ndim)then
376        write(std_out,*) ii,vin(ii:ii+2)
377      else
378        write(std_out,*) ii,vin(ii:ndim)
379      end if
380    end do
381    write(std_out,*) 'vout:'
382    do ii=1,ndim,3
383      if (ii+2<=ndim)then
384        write(std_out,*) ii,vout(ii:ii+2)
385      else
386        write(std_out,*) ii,vout(ii:ndim)
387      end if
388    end do
389    write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim
390    do kk=1,ndim
391      do jj=1,ndim,3
392        if (jj+2<=ndim)then
393          write(std_out,*) jj,hessin(jj:jj+2,kk)
394        else
395          write(std_out,*) jj,hessin(jj:ndim,kk)
396        end if
397      end do
398    end do
399  end if
400 
401 !write(std_out,*) 'bfgs 07'
402 !##########################################################
403 !### 07. Compute the next values
404 
405  if(ionmov==2 .or. itime==1)then
406 
407 !  Previous cartesian coordinates
408    vin_prev(:)=vin(:)
409 
410 !  New atomic cartesian coordinates are obtained from vin, hessin
411 !  and vout
412 
413    do ii=1,ndim
414      vin(:)=vin(:)-hessin(:,ii)*vout(ii)
415    end do
416 
417 !Pulay mixing for vin
418    nitpul=0
419    if (npul>1) then
420      alpha0=1.0_dp
421      nitpul=min(itime, npul)
422      if (itime>npul) then
423        do jj=1,npul-1
424          vinres(jj,:)=vinres(jj+1,:)
425          vin1(jj,:)=vin1(jj+1,:)
426        end do
427      end if
428      vinres(nitpul,:)=vin(:)-vin_prev(:)
429      vin1(nitpul,:)=vin_prev(:)
430    end if
431 
432    if (nitpul>1) then
433      ABI_ALLOCATE(alpha,(nitpul,ndim))
434      alpha=zero
435      do kk=1,ndim
436        ABI_ALLOCATE(amat,(nitpul,nitpul))
437        ABI_ALLOCATE(amatinv,(nitpul,nitpul))
438        amat=zero;amatinv=zero
439        do ii=1,nitpul
440          do jj=ii,nitpul
441            amat(ii,jj)=vinres(jj,kk)*vinres(ii,kk)
442            amat(jj,ii)=amat(ii,jj)
443          end do
444        end do
445        amatinv=amat
446        if (abs(vin(kk)-vin_prev(kk))<tol10) then
447          alpha(:,kk)=zero
448        else
449          ABI_ALLOCATE(ipiv,(nitpul))
450          ABI_ALLOCATE(rwork,(nitpul))
451 !          amatinv=1.d5*amatinv
452          call dgetrf(nitpul,nitpul,amatinv,nitpul,ipiv,ierr)
453          call dgetri(nitpul,amatinv,nitpul,ipiv,rwork,nitpul,ierr)
454 !          amatinv=1.d5*amatinv
455          ABI_DEALLOCATE(ipiv)
456          ABI_DEALLOCATE(rwork)
457          det=zero
458          do ii=1,nitpul
459            do jj=1,nitpul
460              alpha(ii,kk)=alpha(ii,kk)+amatinv(jj,ii)
461              det=det+amatinv(jj,ii)
462            end do
463          end do
464          alpha(:,kk)=alpha(:,kk)/det
465        end if
466      end do
467      ABI_DEALLOCATE(amat)
468      ABI_DEALLOCATE(amatinv)
469      vin(:)=vin1(nitpul,:)+alpha0*(vin1(nitpul+1,:)-vin1(nitpul,:))
470      vin=zero
471      do ii=1,nitpul
472        vin(:)=vin(:)+ alpha(ii,:)*(vin1(ii,:))
473      end do
474      ABI_DEALLOCATE(alpha)
475    end if
476 
477 
478 !  Previous atomic forces
479    vout_prev(:)=vout(:)
480 
481  else
482    if(ionmov==3)then
483      ihist_prev = abihist_findIndex(hist,-1)
484      etotal_prev=hist%etot(ihist_prev)
485 !    Here the BFGS algorithm, modified to take into account the energy
486      call brdene(etotal,etotal_prev,hessin,ndim,vin,vin_prev,vout,vout_prev)
487    end if
488 
489 !  zDEBUG (vin,vout and hessin after prediction)
490    if(zDEBUG)then
491      write(std_out,*) 'Vectors vin and vout [after prediction]'
492      write(std_out,*) 'vin:'
493      do ii=1,ndim,3
494        if (ii+2<=ndim)then
495          write(std_out,*) ii,vin(ii:ii+2)
496        else
497          write(std_out,*) ii,vin(ii:ndim)
498        end if
499      end do
500      write(std_out,*) 'vout:'
501      do ii=1,ndim,3
502        if (ii+2<=ndim)then
503          write(std_out,*) ii,vout(ii:ii+2)
504        else
505          write(std_out,*) ii,vout(ii:ndim)
506        end if
507      end do
508    end if
509 
510 
511 !  FIXME: this should be outside the if clause on ionmov!
512 !  Implement fixing of atoms : put back old values for fixed
513 !  components
514    do kk=1,ab_mover%natom
515      do jj=1,3
516 !      Warning : implemented in reduced coordinates
517        if ( ab_mover%iatfix(jj,kk)==1) then
518          vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3)
519        end if
520      end do
521    end do
522  end if
523 
524 !write(std_out,*) 'bfgs 08'
525 !##########################################################
526 !### 08. Update the history with the prediction
527 
528 !Increase indexes
529  hist%ihist = abihist_findIndex(hist,+1)
530 
531 !Transfer vin  to xred, acell and rprim
532  call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,&
533 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
534 & ab_mover%symrel, ucvol, ucvol0,&
535 & vin, xred)
536 
537  if(ab_mover%optcell/=0)then
538    call mkrdim(acell,rprim,rprimd)
539    call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
540  end if
541 
542 !Fill the history with the variables
543 !xred, acell, rprimd, vel
544  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
545  ihist_prev = abihist_findIndex(hist,-1)
546  hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
547 
548  if(zDEBUG)then
549    write (std_out,*) 'residual:'
550    do kk=1,ab_mover%natom
551      write (std_out,*) residual(:,kk)
552    end do
553    write (std_out,*) 'strten:'
554    write (std_out,*) strten(1:3),ch10,strten(4:6)
555    write (std_out,*) 'etotal:'
556    write (std_out,*) etotal
557  end if
558 
559 end subroutine pred_bfgs

ABINIT/pred_lbfgs [ Functions ]

[ Top ] [ Functions ]

NAME

 pred_lbfgs

FUNCTION

 Ionmov predictors (22) Limited-memory Broyden-Fletcher-Goldfarb-Shanno

 IONMOV 22:
 Given a starting point xred that is a vector of length 3*natom
 (reduced nuclei coordinates), and unit cell parameters
 (acell and rprim) the L-Broyden-Fletcher-Goldfarb-Shanno
 minimization is performed on the total energy function, using
 its gradient (atomic forces and stress : fred or fcart and
 stress) as calculated by the routine scfcv. Some atoms can be
 kept fixed, while the optimization of unit cell parameters is
 only performed if optcell/=0. The convergence requirement on
 the atomic forces, dtset%tolmxf,  allows an early exit.
 Otherwise no more than dtset%ntime steps are performed.
 Returned quantities are xred, and eventually acell and rprim (new ones!).
 Could see MinPack on netlib.org

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 : (22) Specific kind of 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,lbfgs_destroy,lbfgs_init,metric,mkrdim,var2hist
      xfh_recover_new,xfpack_f2vout,xfpack_vin2x,xfpack_x2vin

SOURCE

604 subroutine pred_lbfgs(ab_mover,ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit)
605 
606 
607 !This section has been created automatically by the script Abilint (TD).
608 !Do not modify the following lines by hand.
609 #undef ABI_FUNC
610 #define ABI_FUNC 'pred_lbfgs'
611 !End of the abilint section
612 
613 implicit none
614 
615 !Arguments ------------------------------------
616 !scalars
617 type(abimover),intent(in)       :: ab_mover
618 type(ab_xfh_type),intent(inout)    :: ab_xfh
619 type(abihist),intent(inout) :: hist
620 type(abiforstr),intent(in) :: forstr
621 integer,intent(in) :: itime
622 integer,intent(in) :: ionmov
623 integer,intent(in) :: iexit
624 logical,intent(in) :: zDEBUG
625 
626 !Local variables-------------------------------
627 !scalars
628 integer :: info,ihist_prev
629 integer  :: ndim,cycl_main
630 integer, parameter :: npul=0
631 integer  :: ii,jj,kk
632 real(dp),save :: ucvol0
633 real(dp) :: ucvol
634 real(dp) :: etotal
635 real(dp) :: favg
636 
637 !arrays
638 real(dp),allocatable :: diag(:)
639 real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:)
640 real(dp),allocatable,save :: vout(:),vout_prev(:)
641 real(dp),allocatable,save ::vinres(:,:),vin1(:,:)
642 real(dp),save :: acell0(3) ! Initial acell
643 real(dp),save :: rprimd0(3,3) ! Initial rprimd
644 real(dp) :: acell(3)
645 real(dp) :: rprimd(3,3),rprim(3,3)
646 real(dp) :: gprimd(3,3)
647 real(dp) :: gmet(3,3)
648 real(dp) :: rmet(3,3)
649 real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom)
650 real(dp) :: xred(3,ab_mover%natom)
651 real(dp) :: strten(6)
652 
653 !***************************************************************************
654 !Beginning of executable session
655 !***************************************************************************
656 
657  if(iexit/=0)then
658    call lbfgs_destroy()
659 
660    if (allocated(vin))        then
661      ABI_DEALLOCATE(vin)
662    end if
663    if (allocated(vout))       then
664      ABI_DEALLOCATE(vout)
665    end if
666    if (allocated(vin_prev))   then
667      ABI_DEALLOCATE(vin_prev)
668    end if
669    if (allocated(vout_prev))  then
670      ABI_DEALLOCATE(vout_prev)
671    end if
672    if (allocated(vinres))  then
673      ABI_DEALLOCATE(vinres)
674    end if
675    if (allocated(vin1))  then
676      ABI_DEALLOCATE(vin1)
677    end if
678    if (allocated(hessin))     then
679      ABI_DEALLOCATE(hessin)
680    end if
681    return
682  end if
683 
684 !write(std_out,*) 'bfgs 01'
685 !##########################################################
686 !### 01. Debugging and Verbose
687 
688  if(zDEBUG)then
689    write(std_out,'(a,3a,35a,42a)') ch10,('-',kk=1,3),&
690 &   'Debugging and Verbose for pred_bfgs',('-',kk=1,42)
691    write(std_out,*) 'ionmov: ',ionmov
692    write(std_out,*) 'itime:  ',itime
693  end if
694 
695 !write(std_out,*) 'bfgs 02'
696 !##########################################################
697 !### 02. Compute the dimension of vectors (ndim)
698 
699  ndim=3*ab_mover%natom
700  if(ab_mover%optcell==1 .or.&
701 & ab_mover%optcell==4 .or.&
702 & ab_mover%optcell==5 .or.&
703 & ab_mover%optcell==6) ndim=ndim+1
704  if(ab_mover%optcell==2 .or.&
705 & ab_mover%optcell==3) ndim=ndim+6
706  if(ab_mover%optcell==7 .or.&
707 & ab_mover%optcell==8 .or.&
708 & ab_mover%optcell==9) ndim=ndim+3
709 
710  if(zDEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim
711 
712 !write(std_out,*) 'bfgs 03'
713 !##########################################################
714 !### 03. Allocate the vectors vin, vout and hessian matrix
715 
716 !Notice that vin, vout, etc could be allocated
717 !From a previous dataset with a different ndim
718  if(itime==1)then
719    if (allocated(vin))        then
720      ABI_DEALLOCATE(vin)
721    end if
722    if (allocated(vout))       then
723      ABI_DEALLOCATE(vout)
724    end if
725    if (allocated(vin_prev))   then
726      ABI_DEALLOCATE(vin_prev)
727    end if
728    if (allocated(vout_prev))  then
729      ABI_DEALLOCATE(vout_prev)
730    end if
731    if (allocated(vinres))  then
732      ABI_DEALLOCATE(vinres)
733    end if
734    if (allocated(vin1))  then
735      ABI_DEALLOCATE(vin1)
736    end if
737    if (allocated(hessin))     then
738      ABI_DEALLOCATE(hessin)
739    end if
740    if(npul>1) then
741      ABI_ALLOCATE(vinres,(npul+1,ndim))
742      ABI_ALLOCATE(vin1,(npul+1,ndim))
743    end if
744    ABI_ALLOCATE(vin,(ndim))
745    ABI_ALLOCATE(vout,(ndim))
746    ABI_ALLOCATE(vin_prev,(ndim))
747    ABI_ALLOCATE(vout_prev,(ndim))
748    ABI_ALLOCATE(hessin,(ndim,ndim))
749  end if
750 
751 !write(std_out,*) 'bfgs 04'
752 !##########################################################
753 !### 04. Obtain the present values from the history
754 
755  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
756  do ii=1,3
757    rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
758  end do
759 
760  strten(:)=hist%strten(:,hist%ihist)
761  etotal   =hist%etot(hist%ihist)
762 
763 !Fill the residual with forces (No preconditioning)
764 !Or the preconditioned forces
765  if (ab_mover%goprecon==0)then
766    call fcart2fred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom)
767  else
768    residual(:,:)= forstr%fred(:,:)
769  end if
770 
771  if(zDEBUG)then
772    write (std_out,*) 'residual:'
773    do kk=1,ab_mover%natom
774      write (std_out,*) residual(:,kk)
775    end do
776    write (std_out,*) 'strten:'
777    write (std_out,*) strten(1:3),ch10,strten(4:6)
778    write (std_out,*) 'etotal:'
779    write (std_out,*) etotal
780  end if
781 
782  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
783 
784 !Save initial values
785  if (itime==1)then
786    acell0(:)=acell(:)
787    rprimd0(:,:)=rprimd(:,:)
788    ucvol0=ucvol
789  end if
790 
791 !zDEBUG (UCVOL)
792  if(zDEBUG)then
793    write(std_out,*) 'Volume of cell (ucvol):',ucvol
794  end if
795 
796 !Get rid of mean force on whole unit cell, but only if no
797 !generalized constraints are in effect
798  residual_corrected(:,:)=residual(:,:)
799  if(ab_mover%nconeq==0)then
800    do ii=1,3
801      if (ii/=3.or.ab_mover%jellslab==0) then
802        favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom)
803        residual_corrected(ii,:)=residual_corrected(ii,:)-favg
804      end if
805    end do
806  end if
807 
808 !write(std_out,*) 'bfgs 05'
809 !##########################################################
810 !### 05. Fill the vectors vin and vout
811 
812 !Initialize input vectors : first vin, then vout
813 !The values of vin from the previous iteration
814 !should be the same
815 !if (itime==1)then
816  call xfpack_x2vin(acell, acell0, ab_mover%natom, ndim,&
817 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
818 & ab_mover%symrel, ucvol, ucvol0, vin, xred)
819 !end if
820 
821  call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,&
822 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol,&
823 & vout)
824 
825 !write(std_out,*) 'bfgs 06'
826 !##########################################################
827 !### 06. Initialize or update the hessian matrix
828 
829 !Initialise the Hessian matrix using gmet
830  if (itime==1)then
831 
832    ABI_ALLOCATE(diag,(ndim))
833    do ii=1,3*ab_mover%natom
834 !      diag(ii) = 1.00_dp / rprimd(MODULO(ii-1,3)+1,MODULO(ii-1,3)+1)**2
835      diag(ii) = gmet(MODULO(ii-1,3)+1,MODULO(ii-1,3)+1)
836    end do
837    if(ab_mover%optcell/=0)then
838 !     These values might lead to too large changes in some cases ...
839      do ii=3*ab_mover%natom+1,ndim
840        diag(ii) = ab_mover%strprecon*30.0_dp/ucvol
841        if(ab_mover%optcell==1) diag(ii) = diag(ii) / three
842      end do
843    end if
844 
845    call lbfgs_init(ndim,5,diag)
846    ABI_DEALLOCATE(diag)
847 
848    if (ab_mover%restartxf/=0) then
849 
850      call xfh_recover_new(ab_xfh,ab_mover,acell,acell0,cycl_main,residual,&
851 &     hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,&
852 &     vin_prev,vout,vout_prev,xred)
853 
854    end if
855 
856  end if
857 
858 !zDEBUG (vin,vout and hessin before prediction)
859  if(zDEBUG)then
860    write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]'
861    write(std_out,*) 'vin:'
862    do ii=1,ndim,3
863      if (ii+2<=ndim)then
864        write(std_out,*) ii,vin(ii:ii+2)
865      else
866        write(std_out,*) ii,vin(ii:ndim)
867      end if
868    end do
869    write(std_out,*) 'vout:'
870    do ii=1,ndim,3
871      if (ii+2<=ndim)then
872        write(std_out,*) ii,vout(ii:ii+2)
873      else
874        write(std_out,*) ii,vout(ii:ndim)
875      end if
876    end do
877    write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim
878    do kk=1,ndim
879      do jj=1,ndim,3
880        if (jj+2<=ndim)then
881          write(std_out,*) jj,hessin(jj:jj+2,kk)
882        else
883          write(std_out,*) jj,hessin(jj:ndim,kk)
884        end if
885      end do
886    end do
887  end if
888 
889 !write(std_out,*) 'bfgs 07'
890 !##########################################################
891 !### 07. Compute the next values
892 
893  vin_prev(:) = vin
894  vout_prev(:) = vout
895  info = lbfgs_execute(vin,etotal,vout)
896 
897 !zDEBUG (vin,vout after prediction)
898  if(zDEBUG)then
899    write(std_out,*) 'Vectors vin and vout [after prediction]'
900    write(std_out,*) 'vin_prev:'
901    do ii=1,ndim,3
902      if (ii+2<=ndim)then
903        write(std_out,*) ii,vin_prev(ii:ii+2)
904      else
905        write(std_out,*) ii,vin_prev(ii:ndim)
906      end if
907    end do
908    write(std_out,*) 'vin:'
909    do ii=1,ndim,3
910      if (ii+2<=ndim)then
911        write(std_out,*) ii,vin(ii:ii+2)
912      else
913        write(std_out,*) ii,vin(ii:ndim)
914      end if
915    end do
916    write(std_out,*) 'vout:'
917    do ii=1,ndim,3
918      if (ii+2<=ndim)then
919        write(std_out,*) ii,vout(ii:ii+2)
920      else
921        write(std_out,*) ii,vout(ii:ndim)
922      end if
923    end do
924  end if
925 
926 
927 !Implement fixing of atoms : put back old values for fixed
928 !components
929  do kk=1,ab_mover%natom
930    do jj=1,3
931 !    Warning : implemented in reduced coordinates
932      if ( ab_mover%iatfix(jj,kk)==1) then
933        vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3)
934      end if
935    end do
936  end do
937 
938 
939 !write(std_out,*) 'bfgs 08'
940 !##########################################################
941 !### 08. Update the history with the prediction
942 
943 !Increase indexes
944  hist%ihist = abihist_findIndex(hist,+1)
945 
946 !Transfer vin  to xred, acell and rprim
947  call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,&
948 & ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
949 & ab_mover%symrel, ucvol, ucvol0,&
950 & vin, xred)
951 
952  if(ab_mover%optcell/=0)then
953    call mkrdim(acell,rprim,rprimd)
954    call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
955  end if
956 
957 !Fill the history with the variables
958 !xcart, xred, acell, rprimd
959  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
960  ihist_prev = abihist_findIndex(hist,-1)
961  hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
962 
963  if(zDEBUG)then
964    write (std_out,*) 'residual:'
965    do kk=1,ab_mover%natom
966      write (std_out,*) residual(:,kk)
967    end do
968    write (std_out,*) 'strten:'
969    write (std_out,*) strten(1:3),ch10,strten(4:6)
970    write (std_out,*) 'etotal:'
971    write (std_out,*) etotal
972  end if
973 
974 end subroutine pred_lbfgs