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-2018 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 .

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 MODULE m_bfgs
29 
30  use defs_basis
31  use m_abicore
32  use m_errors
33  use m_abimover
34 
35  use m_io_tools,       only : open_file
36  use m_numeric_tools,  only : findmin
37 
38  implicit none
39 
40  private
41 
42 !public procedures
43  public :: hessinit   ! Initialize Hessian matrix
44  public :: hessupdt   ! Update the hessian matrix
45  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

PARENTS

      pred_bfgs,pred_delocint

CHILDREN

      findmin

SOURCE

381 subroutine brdene(etotal,etotal_prev,hessin,ndim,vin,vin_prev,vout,vout_prev)
382 
383 
384 !This section has been created automatically by the script Abilint (TD).
385 !Do not modify the following lines by hand.
386 #undef ABI_FUNC
387 #define ABI_FUNC 'brdene'
388 !End of the abilint section
389 
390  implicit none
391 
392 !Arguments ------------------------------------
393 !scalars
394  integer,intent(in) :: ndim
395  real(dp),intent(in) :: etotal
396  real(dp),intent(inout) :: etotal_prev
397 !arrays+
398  real(dp),intent(in) :: hessin(ndim,ndim),vout(ndim)
399  real(dp),intent(inout) :: vin(ndim),vin_prev(ndim),vout_prev(ndim)
400 
401 !Local variables-------------------------------
402 !scalars
403  integer :: idim,brd_status
404  real(dp) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_1,dedv_2,dedv_min
405  real(dp) :: dedv_predict,etotal_1,etotal_2,etotal_predict,lambda_1,lambda_2
406  real(dp) :: lambda_predict
407 !arrays
408  real(dp),allocatable :: dvin(:),vin_min(:),vout_min(:)
409 
410 !***************************************************************************
411 
412  ABI_ALLOCATE(dvin,(ndim))
413  ABI_ALLOCATE(vin_min,(ndim))
414  ABI_ALLOCATE(vout_min,(ndim))
415 
416  lambda_1=1.0_dp       ; lambda_2=0.0_dp
417  etotal_1=etotal      ; etotal_2=etotal_prev
418  dvin(:)=vin(:)-vin_prev(:)
419  dedv_1=dot_product(vout,dvin)
420  dedv_2=dot_product(vout_prev,dvin)
421  call findmin(dedv_1,dedv_2,dedv_predict,&
422 & d2edv2_1,d2edv2_2,d2edv2_predict,&
423 & etotal_1,etotal_2,etotal_predict,&
424 & lambda_1,lambda_2,lambda_predict,brd_status)
425 
426 !DEBUG : comes back to usual BFGS !
427 !lambda_predict=1.0_dp
428 !dedv_predict=dedv_1
429 !ENDDEBUG
430 
431 !Generates vin at the minimum, and an interpolated vout, modified
432 !to have the right value of dedv_predict, from findmin.
433  vin_min(:)=vin_prev(:)+lambda_predict*dvin(:)
434  vout_min(:)=vout_prev(:)+lambda_predict*(vout(:)-vout_prev(:))
435  dedv_min=dedv_2+lambda_predict*(dedv_1-dedv_2)
436 !Modify vout_min in order to impose dedv_predict
437  vout_min(:)=vout_min(:)+dvin(:)*(dedv_predict-dedv_min)/dot_product(dvin,dvin)
438 
439 !Previous cartesian coordinates
440  etotal_prev=etotal
441  vin_prev(:)=vin(:)
442 
443 !New atomic cartesian coordinates are obtained from vin, hessin and vout
444  vin(:)=vin_min(:)
445  do idim=1,ndim
446    vin(:)=vin(:)-hessin(:,idim)*vout_min(idim)
447  end do
448 
449 !Previous atomic forces
450  vout_prev(:)=vout(:)
451 
452  ABI_DEALLOCATE(dvin)
453  ABI_DEALLOCATE(vin_min)
454  ABI_DEALLOCATE(vout_min)
455 
456 end subroutine brdene

m_bfgs/hessinit [ Functions ]

[ Top ] [ m_bfgs ] [ Functions ]

NAME

 hessinit

FUNCTION

 Initiliase an 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.

PARENTS

      pred_bfgs,pred_diisrelax

CHILDREN

      findmin

SOURCE

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

PARENTS

      m_mep,pred_bfgs,pred_delocint,pred_diisrelax,xfh_recover_deloc
      xfh_recover_new

CHILDREN

      findmin

SOURCE

218 subroutine hessupdt(hessin,iatfix,natom,ndim,vin,vin_prev,vout,vout_prev, &
219 &                   nimage) ! optional argument
220 
221 
222 !This section has been created automatically by the script Abilint (TD).
223 !Do not modify the following lines by hand.
224 #undef ABI_FUNC
225 #define ABI_FUNC 'hessupdt'
226 !End of the abilint section
227 
228  implicit none
229 
230 !Arguments ------------------------------------
231 !scalars
232  integer,intent(in) :: natom,ndim
233  integer,intent(in),optional :: nimage
234 !arrays
235  integer,intent(in) :: iatfix(3,natom)
236  real(dp),intent(in) :: vin(ndim),vin_prev(ndim),vout(ndim),vout_prev(ndim)
237  real(dp),intent(inout) :: hessin(ndim,ndim)
238 
239 !Local variables-------------------------------
240 !scalars
241  integer :: iatom,idir,ii,jj,nimage_
242  real(dp) :: den1,den2,den3
243  character(len=100) :: msg
244 !arrays
245  real(dp) :: bfgs(ndim),din(ndim),dout(ndim),hdelta(ndim)
246 
247 !***************************************************************************
248 
249  nimage_=1;if (present(nimage)) nimage_=nimage
250 
251 !write(ab_out,*) 'VECTOR INPUT (vin)'
252 !do ii=1,ndim,3
253 !if (ii+2<=ndim)then
254 !write(ab_out,*) ii,vin(ii:ii+2)
255 !else
256 !write(ab_out,*) ii,vin(ii:ndim)
257 !end if
258 !end do
259 !write(ab_out,*) 'VECTOR OUTPUT (vout)'
260 !do ii=1,ndim,3
261 !if (ii+2<=ndim)then
262 !write(ab_out,*) ii,vout(ii:ii+2)
263 !else
264 !write(ab_out,*) ii,vout(ii:ndim)
265 !end if
266 !end do
267 
268 !write(ab_out,*) 'VECTOR INPUT (vin_prev)'
269 !do ii=1,ndim,3
270 !if (ii+2<=ndim)then
271 !write(ab_out,*) ii,vin(ii:ii+2)
272 !else
273 !write(ab_out,*) ii,vin(ii:ndim)
274 !end if
275 !end do
276 !write(ab_out,*) 'VECTOR OUTPUT (vout_prev)'
277 !do ii=1,ndim,3
278 !if (ii+2<=ndim)then
279 !write(ab_out,*) ii,vout(ii:ii+2)
280 !else
281 !write(ab_out,*) ii,vout(ii:ndim)
282 !end if
283 !end do
284 
285  if (mod(ndim,nimage_)/=0) then
286    msg='nimage must be a dividor of ndim !'
287    MSG_BUG(msg)
288  end if
289 
290 !Difference between new and previous vectors
291  din(:) =vin(:) -vin_prev(:)
292  dout(:)=vout(:)-vout_prev(:)
293 
294 !Implement fixing of atoms; must discard the change of forces on fixed atoms
295  do ii=1,nimage_
296    jj=3*natom*(ii-1)
297    do iatom=1,natom
298      do idir=1,3
299        if (iatfix(idir,iatom)==1) dout(idir+jj)=zero
300      end do
301      jj=jj+3
302    end do
303  end do
304 
305 !Compute approximate inverse Hessian times delta fcart
306 !hdelta=hessin*deltaf
307  hdelta(:)=zero
308  do ii=1,ndim
309    hdelta(:)=hdelta(:)+hessin(:,ii)*dout(ii)
310  end do
311 
312 !Calculation of dot products for the denominators
313  den1=zero ; den2=zero
314  do ii=1,ndim
315    den1=den1+dout(ii)*din(ii)
316    den2=den2+dout(ii)*hdelta(ii)
317  end do
318 
319 !DEBUG
320 !write(std_out,*)' hessupdt : den1,den2',den1,den2
321 !write(std_out,*)' din ',din
322 !write(std_out,*)' dout ',dout
323 !write(std_out,*)' hdelta ',hdelta
324 !ENDDEBUG
325 
326 !Denominators are multiplicative
327  den1=one/den1
328  den3=one/den2
329 
330 !Vectors which make a difference between the BROYDEN and
331 !the DAVIDON scheme.
332  bfgs(:)=den1*din(:)-den3*hdelta(:)
333 
334 !B.F.G.S. updating formula
335  do ii=1,ndim
336    do jj=1,ndim
337      hessin(ii,jj)=hessin(ii,jj) +den1*din(ii)*din(jj) &
338 &     -den3*hdelta(ii)*hdelta(jj) +den2*bfgs(ii)*bfgs(jj)
339    end do
340  end do
341 
342 end subroutine hessupdt