!{\src2tex{textfont=tt}}
!!****f* ABINIT/pred_bfgs
!! 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 rprim) the 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 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).
!! Might be better than ionmov=2 for few degrees of freedom (less
!! than 3 or 4)
!!
!! COPYRIGHT
!! Copyright (C) 1998-2014 ABINIT group (DCA, XG, GMR, JCC, SE)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors,
!! see ~abinit/doc/developers/contributors.txt .
!!
!! 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,hessinit,hessupdt,hist2var,metric,mkrdim,var2hist
!!      xfh_recover_new,xfpack_f2vout,xfpack_vin2x,xfpack_x2vin,xred2xcart
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "abi_common.h"

subroutine pred_bfgs(ab_mover,ab_xfh,forstr,hist,ionmov,itime,zDEBUG,iexit)

 use defs_basis
 use m_profiling
 use m_abimover
 use m_abihist

 use m_bfgs, only : hessinit, hessupdt, brdene

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'pred_bfgs'
 use interfaces_41_geometry
 use interfaces_45_geomoptim, except_this_one => pred_bfgs
!End of the abilint section

implicit none

!Arguments ------------------------------------
!scalars
type(abimover),intent(in)       :: ab_mover
type(ab_xfh_type),intent(inout)    :: ab_xfh
type(abihist),intent(inout) :: hist
type(abiforstr),intent(in) :: forstr
integer,intent(in) :: itime
integer,intent(in) :: ionmov
integer,intent(in) :: iexit
logical,intent(in) :: zDEBUG

!Local variables-------------------------------
!scalars
integer  :: ndim,cycl_main
integer, parameter :: npul=0
integer  :: ierr,ii,jj,kk,nitpul
real(dp),save :: ucvol0
real(dp) :: ucvol,det
real(dp) :: etotal,etotal_prev
real(dp) :: favg,alpha0

!arrays
integer,allocatable :: ipiv(:)
real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:)
real(dp),allocatable,save :: vout(:),vout_prev(:)
real(dp),allocatable,save ::vinres(:,:),vin1(:,:)
real(dp),allocatable ::  amat(:,:),amatinv(:,:),alpha(:,:)
real(dp),allocatable :: rwork(:)
real(dp),save :: acell0(3) ! Initial acell
real(dp),save :: rprimd0(3,3) ! Initial rprimd
real(dp) :: acell(3)
real(dp) :: rprimd(3,3),rprim(3,3)
real(dp) :: gprimd(3,3)
real(dp) :: gmet(3,3)
real(dp) :: rmet(3,3)
real(dp) :: residual(3,ab_mover%natom),residual_corrected(3,ab_mover%natom)
real(dp) :: xred(3,ab_mover%natom),xcart(3,ab_mover%natom)
real(dp) :: strten(6)

!***************************************************************************
!Beginning of executable session
!***************************************************************************

 if(iexit/=0)then
   if (allocated(vin))        then
     ABI_DEALLOCATE(vin)
   end if
   if (allocated(vout))       then
     ABI_DEALLOCATE(vout)
   end if
   if (allocated(vin_prev))   then
     ABI_DEALLOCATE(vin_prev)
   end if
   if (allocated(vout_prev))  then
     ABI_DEALLOCATE(vout_prev)
   end if
   if (allocated(vinres))  then
     ABI_DEALLOCATE(vinres)
   end if
   if (allocated(vin1))  then
     ABI_DEALLOCATE(vin1)
   end if
   if (allocated(hessin))     then
     ABI_DEALLOCATE(hessin) 
   end if
   return
 end if

!write(std_out,*) 'bfgs 01'
!##########################################################
!### 01. Debugging and Verbose

 if(zDEBUG)then
   write(std_out,'(a,3a,35a,42a)') ch10,('-',kk=1,3),&
&   'Debugging and Verbose for pred_bfgs',('-',kk=1,42)
   write(std_out,*) 'ionmov: ',ionmov
   write(std_out,*) 'itime:  ',itime
 end if

!write(std_out,*) 'bfgs 02'
!##########################################################
!### 02. Compute the dimension of vectors (ndim)

 ndim=3*ab_mover%natom
 if(ab_mover%optcell==1 .or.&
& ab_mover%optcell==4 .or.&
& ab_mover%optcell==5 .or.&
& ab_mover%optcell==6) ndim=ndim+1
 if(ab_mover%optcell==2 .or.&
& ab_mover%optcell==3) ndim=ndim+6
 if(ab_mover%optcell==7 .or.&
& ab_mover%optcell==8 .or.&
& ab_mover%optcell==9) ndim=ndim+3

 if(zDEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim

!write(std_out,*) 'bfgs 03'
!##########################################################
!### 03. Allocate the vectors vin, vout and hessian matrix

!Notice that vin, vout, etc could be allocated
!From a previous dataset with a different ndim
 if(itime==1)then
   if (allocated(vin))        then
     ABI_DEALLOCATE(vin)
   end if
   if (allocated(vout))       then
     ABI_DEALLOCATE(vout)
   end if
   if (allocated(vin_prev))   then
     ABI_DEALLOCATE(vin_prev)
   end if
   if (allocated(vout_prev))  then
     ABI_DEALLOCATE(vout_prev)
   end if
    if (allocated(vinres))  then
     ABI_DEALLOCATE(vinres)
   end if
   if (allocated(vin1))  then
     ABI_DEALLOCATE(vin1)
   end if
   if (allocated(hessin))     then
     ABI_DEALLOCATE(hessin)
   end if
   if(npul>1) then
     ABI_ALLOCATE(vinres,(npul+1,ndim))
     ABI_ALLOCATE(vin1,(npul+1,ndim))
   end if
   ABI_ALLOCATE(vin,(ndim))
   ABI_ALLOCATE(vout,(ndim))
   ABI_ALLOCATE(vin_prev,(ndim))
   ABI_ALLOCATE(vout_prev,(ndim))
   ABI_ALLOCATE(hessin,(ndim,ndim))
 end if

!write(std_out,*) 'bfgs 04'
!##########################################################
!### 04. Obtain the present values from the history

 call hist2var(acell,hist,ab_mover%natom,rprim,rprimd,xcart,xred,zDEBUG)

 strten(:)=hist%histS(:,hist%ihist)
 etotal   =hist%histE(hist%ihist)

!Fill the residual with forces (No preconditioning)
!Or the preconditioned forces
 if (ab_mover%goprecon==0)then
   residual(:,:)= hist%histXF(:,:,4,hist%ihist)
 else
   residual(:,:)= forstr%fred(:,:)
 end if

 if(zDEBUG)then
   write (std_out,*) 'residual:'
   do kk=1,ab_mover%natom
     write (std_out,*) residual(:,kk)
   end do
   write (std_out,*) 'strten:'
   write (std_out,*) strten(1:3),ch10,strten(4:6)
   write (std_out,*) 'etotal:'
   write (std_out,*) etotal
 end if

 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)

!Save initial values
 if (itime==1)then
   acell0(:)=acell(:)
   rprimd0(:,:)=rprimd(:,:)
   ucvol0=ucvol
 end if

!zDEBUG (UCVOL)
 if(zDEBUG)then
   write(std_out,*) 'Volume of cell (ucvol):',ucvol
 end if

!Get rid of mean force on whole unit cell, but only if no
!generalized constraints are in effect
 if(ab_mover%nconeq==0)then
   do ii=1,3
     favg=sum(residual(ii,:))/dble(ab_mover%natom)
     residual_corrected(ii,:)=residual(ii,:)-favg
     if(ab_mover%jellslab/=0.and.ii==3)&
&     residual_corrected(ii,:)=residual(ii,:)
   end do
 else
   residual_corrected(:,:)=residual(:,:)
 end if

!write(std_out,*) 'bfgs 05'
!##########################################################
!### 05. Fill the vectors vin and vout

!Initialize input vectors : first vin, then vout
!The values of vin from the previous iteration
!should be the same
!if (itime==1)then
 call xfpack_x2vin(acell, acell0, ab_mover%natom, ndim,&
& ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
& ab_mover%symrel, ucvol, ucvol0, vin, xred)
!end if

 call xfpack_f2vout(residual_corrected, ab_mover%natom, ndim,&
& ab_mover%optcell, ab_mover%strtarget, strten, ucvol,&
& vout)

!write(std_out,*) 'bfgs 06'
!##########################################################
!### 06. Initialize or update the hessian matrix

!Initialise the Hessian matrix using gmet
 if (itime==1)then

   call hessinit(ab_mover, hessin, gmet, ndim, ucvol)

!  ! Initialize inverse hessian with identity matrix
!  ! in cartesian coordinates, which makes use of metric tensor gmet
!  ! in reduced coordinates.
!  hessin(:,:)=zero
!  do ii=1,ab_mover%natom
!  do kk=1,3
!  do jj=1,3
!  ! Warning : implemented in reduced coordinates
!  if (ab_mover%iatfix(kk,ii)==0 .and.&
!  & ab_mover%iatfix(jj,ii)==0 )then
!  hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj)
!  end if
!  end do
!  end do
!  end do
!  if(ab_mover%optcell/=0)then
!  ! These values might lead to too large changes in some cases
!  diag=ab_mover%strprecon*30.0_dp/ucvol
!  if(ab_mover%optcell==1) diag=diag/three
!  do ii=3*ab_mover%natom+1,ndim
!  hessin(ii,ii)=diag
!  end do
!  end if

   if (ab_mover%restartxf/=0) then

     call xfh_recover_new(ab_xfh,ab_mover,acell,acell0,cycl_main,residual,&
&     hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,&
&     vin_prev,vout,vout_prev,xred)

   end if

 end if

 if(itime>1)then
!  Update the hessian matrix, by taking into account the
!  current pair (x,f) and the previous one.
   call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom,ndim,vin,&
&   vin_prev,vout,vout_prev)

 end if

!zDEBUG (vin,vout and hessin before prediction)
 if(zDEBUG)then
   write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]'
   write(std_out,*) 'vin:'
   do ii=1,ndim,3
     if (ii+2<=ndim)then
       write(std_out,*) ii,vin(ii:ii+2)
     else
       write(std_out,*) ii,vin(ii:ndim)
     end if
   end do
   write(std_out,*) 'vout:'
   do ii=1,ndim,3
     if (ii+2<=ndim)then
       write(std_out,*) ii,vout(ii:ii+2)
     else
       write(std_out,*) ii,vout(ii:ndim)
     end if
   end do
   write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim
   do kk=1,ndim
     do jj=1,ndim,3
       if (jj+2<=ndim)then
         write(std_out,*) jj,hessin(jj:jj+2,kk)
       else
         write(std_out,*) jj,hessin(jj:ndim,kk)
       end if
     end do
   end do
 end if

!write(std_out,*) 'bfgs 07'
!##########################################################
!### 07. Compute the next values

 if(ionmov==2 .or. itime==1)then

!  Previous cartesian coordinates
   vin_prev(:)=vin(:)

!  New atomic cartesian coordinates are obtained from vin, hessin
!  and vout

   do ii=1,ndim
     vin(:)=vin(:)-hessin(:,ii)*vout(ii)
   end do

!Pulay mixing for vin
   nitpul=0
   if (npul>1) then
     alpha0=1.0_dp
     nitpul=min(itime, npul)
     if (itime>npul) then
       do jj=1,npul-1
         vinres(jj,:)=vinres(jj+1,:)
         vin1(jj,:)=vin1(jj+1,:)
       end do
     end if 
     vinres(nitpul,:)=vin(:)-vin_prev(:)
     vin1(nitpul,:)=vin_prev(:)
   end if

   if (nitpul>1) then
     ABI_ALLOCATE(alpha,(nitpul,ndim))
     alpha=zero
     do kk=1,ndim
       ABI_ALLOCATE(amat,(nitpul,nitpul))
       ABI_ALLOCATE(amatinv,(nitpul,nitpul))
       amat=zero;amatinv=zero
       do ii=1,nitpul
         do jj=ii,nitpul
           amat(ii,jj)=vinres(jj,kk)*vinres(ii,kk)
           amat(jj,ii)=amat(ii,jj)
         end do
       end do
       amatinv=amat
       if (abs(vin(kk)-vin_prev(kk))<tol10) then
         alpha(:,kk)=zero
       else
         ABI_ALLOCATE(ipiv,(nitpul))
         ABI_ALLOCATE(rwork,(nitpul))
!          amatinv=1.d5*amatinv
         call dgetrf(nitpul,nitpul,amatinv,nitpul,ipiv,ierr)
         call dgetri(nitpul,amatinv,nitpul,ipiv,rwork,nitpul,ierr) 
!          amatinv=1.d5*amatinv
         ABI_DEALLOCATE(ipiv)
         ABI_DEALLOCATE(rwork)
         det=zero
         do ii=1,nitpul
           do jj=1,nitpul
             alpha(ii,kk)=alpha(ii,kk)+amatinv(jj,ii)
             det=det+amatinv(jj,ii)
           end do
         end do
         alpha(:,kk)=alpha(:,kk)/det
       end if
     end do
     ABI_DEALLOCATE(amat)
     ABI_DEALLOCATE(amatinv)
     vin(:)=vin1(nitpul,:)+alpha0*(vin1(nitpul+1,:)-vin1(nitpul,:))
     vin=zero
     do ii=1,nitpul
       vin(:)=vin(:)+ alpha(ii,:)*(vin1(ii,:))
     end do
      ABI_DEALLOCATE(alpha)
   end if

   
!  Previous atomic forces
   vout_prev(:)=vout(:)

 else
   if(ionmov==3)then

     etotal_prev=hist%histE(hist%ihist-1)
!    Here the BFGS algorithm, modified to take into account the energy
     call brdene(etotal,etotal_prev,hessin,ndim,vin,vin_prev,vout,vout_prev)
   end if

!  zDEBUG (vin,vout and hessin after prediction)
   if(zDEBUG)then
     write(std_out,*) 'Vectors vin and vout [after prediction]'
     write(std_out,*) 'vin:'
     do ii=1,ndim,3
       if (ii+2<=ndim)then
         write(std_out,*) ii,vin(ii:ii+2)
       else
         write(std_out,*) ii,vin(ii:ndim)
       end if
     end do
     write(std_out,*) 'vout:'
     do ii=1,ndim,3
       if (ii+2<=ndim)then
         write(std_out,*) ii,vout(ii:ii+2)
       else
         write(std_out,*) ii,vout(ii:ndim)
       end if
     end do
   end if


!  FIXME: this should be outside the if clause on ionmov!
!  Implement fixing of atoms : put back old values for fixed
!  components
   do kk=1,ab_mover%natom
     do jj=1,3
!      Warning : implemented in reduced coordinates
       if ( ab_mover%iatfix(jj,kk)==1) then
         vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3)
       end if
     end do
   end do
 end if

!write(std_out,*) 'bfgs 08'
!##########################################################
!### 08. Update the history with the prediction

!Increase indexes
 hist%ihist=hist%ihist+1

!Transfer vin  to xred, acell and rprim
 call xfpack_vin2x(acell, acell0, ab_mover%natom, ndim,&
& ab_mover%nsym, ab_mover%optcell, rprim, rprimd0,&
& ab_mover%symrel, ucvol, ucvol0,&
& vin, xred)

 if(ab_mover%optcell/=0)then
   call mkrdim(acell,rprim,rprimd)
   call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 end if

!Compute xcart from xred, and rprimd
 call xred2xcart(ab_mover%natom,rprimd,xcart,xred)

!Fill the history with the variables
!xcart, xred, acell, rprimd
 call var2hist(acell,hist,ab_mover%natom,rprim,rprimd,xcart,xred,zDEBUG)

 if(zDEBUG)then
   write (std_out,*) 'residual:'
   do kk=1,ab_mover%natom
     write (std_out,*) residual(:,kk)
   end do
   write (std_out,*) 'strten:'
   write (std_out,*) strten(1:3),ch10,strten(4:6)
   write (std_out,*) 'etotal:'
   write (std_out,*) etotal
 end if

 hist%histV(:,:,hist%ihist)=hist%histV(:,:,hist%ihist-1)

end subroutine pred_bfgs
!!***

