TABLE OF CONTENTS

ABINIT/m_bfgs [ Modules ]

[ Top ] [ Modules ]

NAME

  m_bfgs

FUNCTION

  This module provides several routines for the application of a
  Broyden-Fletcher-Goldfarb-Shanno (BFGS) minimization algorithm.

COPYRIGHT

 Copyright (C) 2012-2024 ABINIT group (XG,JCC)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_bfgs
24 
25  use defs_basis
26  use m_abicore
27  use m_errors
28  use m_abimover
29 
30  use m_io_tools,       only : open_file
31  use m_numeric_tools,  only : findmin
32 
33  implicit none
34 
35  private
36 
37 !public procedures
38  public :: hessinit   ! Initialize Hessian matrix
39  public :: hessupdt   ! Update the hessian matrix
40  public :: brdene

m_bfgs/brdene [ Functions ]

[ Top ] [ m_bfgs ] [ Functions ]

NAME

 brdene

FUNCTION

 Update vin according to the Broyden formula, combined
 with a line minimisation that take into account the total energies.
 Also transfer vin to vin_prev, vout to vout_prev, and etotal to etotal_prev
 Could see Numerical Recipes (Fortran), 1986, page 307,
 as well as Schlegel, J. Comp. Chem. 3, 214 (1982) [[cite:Schlegel1982]].

INPUTS

  etotal=new total energy (no meaning at output)
  hessin(ndim,ndim)=hessian matrix
  ndim=size of the hessian and vectors
  vout(ndim)=new output vector (no meaning at output)

OUTPUT

  (see side effects)

SIDE EFFECTS

  etotal_prev=previous total energy; contains input etotal at output
  vin(ndim)=new input vector; updated at output
  vin_prev(ndim)=previous input vector; contains input vin at output
  vout_prev(ndim)=previous output vector; contains input vout at output

SOURCE

338 subroutine brdene(etotal,etotal_prev,hessin,ndim,vin,vin_prev,vout,vout_prev)
339 
340 !Arguments ------------------------------------
341 !scalars
342  integer,intent(in) :: ndim
343  real(dp),intent(in) :: etotal
344  real(dp),intent(inout) :: etotal_prev
345 !arrays+
346  real(dp),intent(in) :: hessin(ndim,ndim),vout(ndim)
347  real(dp),intent(inout) :: vin(ndim),vin_prev(ndim),vout_prev(ndim)
348 
349 !Local variables-------------------------------
350 !scalars
351  integer :: idim,brd_status
352  real(dp) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_1,dedv_2,dedv_min
353  real(dp) :: dedv_predict,etotal_1,etotal_2,etotal_predict,lambda_1,lambda_2
354  real(dp) :: lambda_predict
355 !arrays
356  real(dp),allocatable :: dvin(:),vin_min(:),vout_min(:)
357 
358 !***************************************************************************
359 
360  ABI_MALLOC(dvin,(ndim))
361  ABI_MALLOC(vin_min,(ndim))
362  ABI_MALLOC(vout_min,(ndim))
363 
364  lambda_1=1.0_dp       ; lambda_2=0.0_dp
365  etotal_1=etotal      ; etotal_2=etotal_prev
366  dvin(:)=vin(:)-vin_prev(:)
367  dedv_1=dot_product(vout,dvin)
368  dedv_2=dot_product(vout_prev,dvin)
369  call findmin(dedv_1,dedv_2,dedv_predict,&
370 & d2edv2_1,d2edv2_2,d2edv2_predict,&
371 & etotal_1,etotal_2,etotal_predict,&
372 & lambda_1,lambda_2,lambda_predict,brd_status)
373 
374 !DEBUG : comes back to usual BFGS !
375 !lambda_predict=1.0_dp
376 !dedv_predict=dedv_1
377 !ENDDEBUG
378 
379 !Generates vin at the minimum, and an interpolated vout, modified
380 !to have the right value of dedv_predict, from findmin.
381  vin_min(:)=vin_prev(:)+lambda_predict*dvin(:)
382  vout_min(:)=vout_prev(:)+lambda_predict*(vout(:)-vout_prev(:))
383  dedv_min=dedv_2+lambda_predict*(dedv_1-dedv_2)
384 !Modify vout_min in order to impose dedv_predict
385  vout_min(:)=vout_min(:)+dvin(:)*(dedv_predict-dedv_min)/dot_product(dvin,dvin)
386 
387 !Previous cartesian coordinates
388  etotal_prev=etotal
389  vin_prev(:)=vin(:)
390 
391 !New atomic cartesian coordinates are obtained from vin, hessin and vout
392  vin(:)=vin_min(:)
393  do idim=1,ndim
394    vin(:)=vin(:)-hessin(:,idim)*vout_min(idim)
395  end do
396 
397 !Previous atomic forces
398  vout_prev(:)=vout(:)
399 
400  ABI_FREE(dvin)
401  ABI_FREE(vin_min)
402  ABI_FREE(vout_min)
403 
404 end subroutine brdene

m_bfgs/hessinit [ Functions ]

[ Top ] [ m_bfgs ] [ Functions ]

NAME

 hessinit

FUNCTION

 Initialize the Hessian matrix, either from disk or using init_matrix.
 The size ndim must be greater or equal than 3 * ab_mover%natom.

INPUTS

  fnameabi_hes=filename for Hessian matrix
  ab_mover = the input variables relevant for moving ions
  init_matrix(3,3)=matrix used for each atom (if iatfix = 0) for initialisation.
  ndim=size of the hessian and vectors
  ucvol=volume of the box (used when ab_mover%optcell is not null).

OUTPUT

  hessin(ndim,ndim)=hessian matrix, initialised at output.

SOURCE

 68 subroutine hessinit(ab_mover, hessin, init_matrix, ndim, ucvol)
 69 
 70 !Arguments ------------------------------------
 71 !scalars
 72  integer,intent(in) :: ndim
 73  real(dp),intent(in) :: ucvol
 74  type(abimover),intent(in) :: ab_mover
 75 !arrays
 76  real(dp),intent(in) :: init_matrix(3,3)
 77  real(dp),intent(out) :: hessin(ndim,ndim)
 78 
 79 !Local variables-------------------------------
 80 !scalars
 81  integer :: hess_ok,iatom,idim,idir1,idir2,ii,ios,jj,ndim0,temp_unit
 82  real(dp) :: diag
 83  logical :: ex
 84  character(len=500) :: message
 85 
 86 ! *********************************************************************
 87 
 88 !Initialization of the inverse Hessian to the unit matrix
 89 !Much better choices are possible--this simply corresponds to
 90 !taking the first minimization step as the negative of the
 91 !gradient, with the full length of the gradient vector as
 92 !the step size.  Any spring type model would probably be a better starting guess.
 93 
 94  if (ndim < 3 * ab_mover%natom) then
 95    write(message, '(a,a,a)' )&
 96 &   'the size of the given hessian matrix is too small.', ch10, &
 97 &   'This is an internal error, contact ABINIT developers.'
 98    ABI_ERROR(message)
 99  end if
100 
101 !Special arrangement: if input hessian file exists, read data from there
102  inquire (file=ab_mover%fnameabi_hes,iostat=ios,exist=ex)
103  hess_ok=0
104 
105  if (ex) then
106    ! Read inverse hessian data from file; format is
107    if (open_file(ab_mover%fnameabi_hes,message,newunit=temp_unit,form='formatted',status='old') /= 0) then
108      ABI_ERROR(message)
109    end if
110    read (temp_unit,*)
111    read (temp_unit,*) ndim0
112    if (ndim0/=ndim) then
113 !    Cannot read data because data file natom disagrees with current job
114      write(message,'(5a,i10,a,i10,2a)')&
115 &     'Tried to read inverse hessian from file',trim(ab_mover%fnameabi_hes),' but',ch10,&
116 &     'ndim of that file =',ndim0,' , is not equal to input ndim =',ndim,ch10,&
117 &     ' => initialize inverse hessian with identity matrix.'
118      ABI_WARNING(message)
119      close(unit=temp_unit)
120    else
121      ! Read inverse hessian
122      do jj=1,ndim
123        read (temp_unit,*)
124        read (temp_unit,*) (hessin(ii,jj),ii=1,ndim)
125      end do
126      close (unit=temp_unit)
127      write(message,*)' Inverse hessian has been input from input hessian file',trim(ab_mover%fnameabi_hes)
128      call wrtout(std_out,message,'COLL')
129      hess_ok=1
130    end if
131  end if
132 
133 !If hessin was not read, initialize inverse hessian with identity matrix
134 !in cartesian coordinates, which makes use of metric tensor gmet in reduced coordinates.
135  if(hess_ok==0)then
136    hessin(:,:)=zero
137    do iatom=1,ab_mover%natom
138      do idir1=1,3
139        do idir2=1,3
140 !        Warning : implemented in reduced coordinates
141          if ( ab_mover%iatfix(idir1,iatom) ==0 .and. ab_mover%iatfix(idir2,iatom) ==0 )then
142            hessin(idir1+3*(iatom-1),idir2+3*(iatom-1))=init_matrix(idir1,idir2)
143          end if
144        end do
145      end do
146    end do
147    if(ab_mover%optcell/=0)then
148 !    These values might lead to too large changes in some cases ...
149      diag=ab_mover%strprecon*30.0_dp/ucvol
150      if(ab_mover%optcell==1)diag=diag/three
151      do idim=3*ab_mover%natom+1,ndim
152        hessin(idim,idim)=diag
153      end do
154    end if
155    call wrtout(std_out,'Inverse hessian has been initialized.','COLL')
156  end if
157 
158 end subroutine hessinit

m_bfgs/hessupdt [ Functions ]

[ Top ] [ m_bfgs ] [ Functions ]

NAME

 hessupdt

FUNCTION

 Update of the hessian matrix according to the Broyden formula.
 Could see Numerical Recipes (Fortran), 1986, page 307.

INPUTS

  iatfix(3,natom)=1 for each atom fixed along specified direction, else 0
  natom=number of atoms in unit cell
  ndim=size of the hessian and vectors
  nimage= -- optional, default=1 --
         Number of images of the system described in
         vin, vin_prev, vout, vout_prev
  vin(ndim)=new input vector
  vin_prev(ndim)=previous input vector
  vout(ndim)=new output vector
  vout_prev(ndim)=previous output vector

OUTPUT

  (see side effects)

SIDE EFFECTS

  hessin(ndim,ndim)=hessian matrix, updated at output.

SOURCE

191 subroutine hessupdt(hessin,iatfix,natom,ndim,vin,vin_prev,vout,vout_prev, &
192 &                   nimage) ! optional argument
193 
194 !Arguments ------------------------------------
195 !scalars
196  integer,intent(in) :: natom,ndim
197  integer,intent(in),optional :: nimage
198 !arrays
199  integer,intent(in) :: iatfix(3,natom)
200  real(dp),intent(in) :: vin(ndim),vin_prev(ndim),vout(ndim),vout_prev(ndim)
201  real(dp),intent(inout) :: hessin(ndim,ndim)
202 
203 !Local variables-------------------------------
204 !scalars
205  integer :: iatom,idir,ii,jj,nimage_
206  real(dp) :: den1,den2,den3
207  !character(len=500) :: msg
208 !arrays
209  real(dp) :: bfgs(ndim),din(ndim),dout(ndim),hdelta(ndim)
210 
211 !***************************************************************************
212 
213  nimage_=1;if (present(nimage)) nimage_=nimage
214 
215 !write(ab_out,*) 'VECTOR INPUT (vin)'
216 !do ii=1,ndim,3
217 !if (ii+2<=ndim)then
218 !write(ab_out,*) ii,vin(ii:ii+2)
219 !else
220 !write(ab_out,*) ii,vin(ii:ndim)
221 !end if
222 !end do
223 !write(ab_out,*) 'VECTOR OUTPUT (vout)'
224 !do ii=1,ndim,3
225 !if (ii+2<=ndim)then
226 !write(ab_out,*) ii,vout(ii:ii+2)
227 !else
228 !write(ab_out,*) ii,vout(ii:ndim)
229 !end if
230 !end do
231 
232 !write(ab_out,*) 'VECTOR INPUT (vin_prev)'
233 !do ii=1,ndim,3
234 !if (ii+2<=ndim)then
235 !write(ab_out,*) ii,vin(ii:ii+2)
236 !else
237 !write(ab_out,*) ii,vin(ii:ndim)
238 !end if
239 !end do
240 !write(ab_out,*) 'VECTOR OUTPUT (vout_prev)'
241 !do ii=1,ndim,3
242 !if (ii+2<=ndim)then
243 !write(ab_out,*) ii,vout(ii:ii+2)
244 !else
245 !write(ab_out,*) ii,vout(ii:ndim)
246 !end if
247 !end do
248 
249  if (mod(ndim,nimage_)/=0) then
250    ABI_BUG('nimage must be a dividor of ndim !')
251  end if
252 
253 !Difference between new and previous vectors
254  din(:) =vin(:) -vin_prev(:)
255  dout(:)=vout(:)-vout_prev(:)
256 
257 !Implement fixing of atoms; must discard the change of forces on fixed atoms
258  do ii=1,nimage_
259    jj=3*natom*(ii-1)
260    do iatom=1,natom
261      do idir=1,3
262        if (iatfix(idir,iatom)==1) dout(idir+jj)=zero
263      end do
264      jj=jj+3
265    end do
266  end do
267 
268 !Compute approximate inverse Hessian times delta fcart
269 !hdelta=hessin*deltaf
270  hdelta(:)=zero
271  do ii=1,ndim
272    hdelta(:)=hdelta(:)+hessin(:,ii)*dout(ii)
273  end do
274 
275 !Calculation of dot products for the denominators
276  den1=zero ; den2=zero
277  do ii=1,ndim
278    den1=den1+dout(ii)*din(ii)
279    den2=den2+dout(ii)*hdelta(ii)
280  end do
281 
282 !DEBUG
283 !write(std_out,*)' hessupdt : den1,den2',den1,den2
284 !write(std_out,*)' din ',din
285 !write(std_out,*)' dout ',dout
286 !write(std_out,*)' hdelta ',hdelta
287 !ENDDEBUG
288 
289 !Denominators are multiplicative
290  den1=one/den1
291  den3=one/den2
292 
293 !Vectors which make a difference between the BROYDEN and
294 !the DAVIDON scheme.
295  bfgs(:)=den1*din(:)-den3*hdelta(:)
296 
297 !B.F.G.S. updating formula
298  do ii=1,ndim
299    do jj=1,ndim
300      hessin(ii,jj)=hessin(ii,jj) +den1*din(ii)*din(jj) &
301 &     -den3*hdelta(ii)*hdelta(jj) +den2*bfgs(ii)*bfgs(jj)
302    end do
303  end do
304 
305 end subroutine hessupdt