Line data Source code
1 : !!****m* ABINIT/m_bfgs
2 : !! NAME
3 : !! m_bfgs
4 : !!
5 : !! FUNCTION
6 : !! This module provides several routines for the application of a
7 : !! Broyden-Fletcher-Goldfarb-Shanno (BFGS) minimization algorithm.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2012-2026 ABINIT group (XG,JCC)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
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
41 : !!***
42 :
43 : CONTAINS !===========================================================
44 : !!***
45 :
46 : !----------------------------------------------------------------------
47 :
48 : !!****f* m_bfgs/hessinit
49 : !! NAME
50 : !! hessinit
51 : !!
52 : !! FUNCTION
53 : !! Initialize the Hessian matrix, either from disk or using init_matrix.
54 : !! The size ndim must be greater or equal than 3 * ab_mover%natom.
55 : !!
56 : !! INPUTS
57 : !! fnameabi_hes=filename for Hessian matrix
58 : !! ab_mover = the input variables relevant for moving ions
59 : !! init_matrix(3,3)=matrix used for each atom (if iatfix = 0) for initialisation.
60 : !! ndim=size of the hessian and vectors
61 : !! ucvol=volume of the box (used when ab_mover%optcell is not null).
62 : !!
63 : !! OUTPUT
64 : !! hessin(ndim,ndim)=hessian matrix, initialised at output.
65 : !!
66 : !! SOURCE
67 :
68 286 : 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 286 : if (ndim < 3 * ab_mover%natom) then
95 : write(message, '(a,a,a)' )&
96 0 : & 'the size of the given hessian matrix is too small.', ch10, &
97 0 : & 'This is an internal error, contact ABINIT developers.'
98 0 : ABI_ERROR(message)
99 : end if
100 :
101 : !Special arrangement: if input hessian file exists, read data from there
102 286 : inquire (file=ab_mover%fnameabi_hes,iostat=ios,exist=ex)
103 286 : hess_ok=0
104 :
105 286 : if (ex) then
106 : ! Read inverse hessian data from file; format is
107 0 : if (open_file(ab_mover%fnameabi_hes,message,newunit=temp_unit,form='formatted',status='old') /= 0) then
108 0 : ABI_ERROR(message)
109 : end if
110 0 : read (temp_unit,*)
111 0 : read (temp_unit,*) ndim0
112 0 : 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 0 : & 'Tried to read inverse hessian from file',trim(ab_mover%fnameabi_hes),' but',ch10,&
116 0 : & 'ndim of that file =',ndim0,' , is not equal to input ndim =',ndim,ch10,&
117 0 : & ' => initialize inverse hessian with identity matrix.'
118 0 : ABI_WARNING(message)
119 0 : close(unit=temp_unit)
120 : else
121 : ! Read inverse hessian
122 0 : do jj=1,ndim
123 0 : read (temp_unit,*)
124 0 : read (temp_unit,*) (hessin(ii,jj),ii=1,ndim)
125 : end do
126 0 : close (unit=temp_unit)
127 0 : write(message,*)' Inverse hessian has been input from input hessian file',trim(ab_mover%fnameabi_hes)
128 0 : 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 29696 : hessin(:,:)=zero
137 936 : do iatom=1,ab_mover%natom
138 2886 : do idir1=1,3
139 8450 : do idir2=1,3
140 : ! Warning : implemented in reduced coordinates
141 7800 : if ( ab_mover%iatfix(idir1,iatom) ==0 .and. ab_mover%iatfix(idir2,iatom) ==0 )then
142 5412 : 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 286 : if(ab_mover%optcell/=0)then
148 : ! These values might lead to too large changes in some cases ...
149 93 : diag=ab_mover%strprecon*30.0_dp/ucvol
150 93 : if(ab_mover%optcell==1)diag=diag/three
151 296 : do idim=3*ab_mover%natom+1,ndim
152 296 : hessin(idim,idim)=diag
153 : end do
154 : end if
155 286 : call wrtout(std_out,'Inverse hessian has been initialized.','COLL')
156 : end if
157 :
158 286 : end subroutine hessinit
159 : !!***
160 :
161 : !----------------------------------------------------------------------
162 :
163 : !!****f* m_bfgs/hessupdt
164 : !! NAME
165 : !! hessupdt
166 : !!
167 : !! FUNCTION
168 : !! Update of the hessian matrix according to the Broyden formula.
169 : !! Could see Numerical Recipes (Fortran), 1986, page 307.
170 : !!
171 : !! INPUTS
172 : !! iatfix(3,natom)=1 for each atom fixed along specified direction, else 0
173 : !! natom=number of atoms in unit cell
174 : !! ndim=size of the hessian and vectors
175 : !! nimage= -- optional, default=1 --
176 : !! Number of images of the system described in
177 : !! vin, vin_prev, vout, vout_prev
178 : !! vin(ndim)=new input vector
179 : !! vin_prev(ndim)=previous input vector
180 : !! vout(ndim)=new output vector
181 : !! vout_prev(ndim)=previous output vector
182 : !!
183 : !! OUTPUT
184 : !! (see side effects)
185 : !!
186 : !! SIDE EFFECTS
187 : !! hessin(ndim,ndim)=hessian matrix, updated at output.
188 : !!
189 : !! SOURCE
190 :
191 652 : 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 1304 : real(dp) :: bfgs(ndim),din(ndim),dout(ndim),hdelta(ndim)
210 :
211 : !***************************************************************************
212 :
213 652 : 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 652 : if (mod(ndim,nimage_)/=0) then
250 0 : ABI_BUG('nimage must be a dividor of ndim !')
251 : end if
252 :
253 : !Difference between new and previous vectors
254 7957 : din(:) =vin(:) -vin_prev(:)
255 7957 : dout(:)=vout(:)-vout_prev(:)
256 :
257 : !Implement fixing of atoms; must discard the change of forces on fixed atoms
258 1304 : do ii=1,nimage_
259 652 : jj=3*natom*(ii-1)
260 3465 : do iatom=1,natom
261 8644 : do idir=1,3
262 8644 : if (iatfix(idir,iatom)==1) dout(idir+jj)=zero
263 : end do
264 2813 : jj=jj+3
265 : end do
266 : end do
267 :
268 : !Compute approximate inverse Hessian times delta fcart
269 : !hdelta=hessin*deltaf
270 7957 : hdelta(:)=zero
271 7957 : do ii=1,ndim
272 224620 : hdelta(:)=hdelta(:)+hessin(:,ii)*dout(ii)
273 : end do
274 :
275 : !Calculation of dot products for the denominators
276 7957 : den1=zero ; den2=zero
277 7957 : den1 = sum(dout(1:ndim) * din(1:ndim))
278 7957 : den2 = sum(dout(1:ndim) * hdelta(1:ndim))
279 :
280 : !DEBUG
281 : !write(std_out,*)' hessupdt : den1,den2',den1,den2
282 : !write(std_out,*)' din ',din
283 : !write(std_out,*)' dout ',dout
284 : !write(std_out,*)' hdelta ',hdelta
285 : !ENDDEBUG
286 :
287 : !Denominators are multiplicative
288 652 : den1=one/den1
289 652 : den3=one/den2
290 :
291 : !Vectors which make a difference between the BROYDEN and
292 : !the DAVIDON scheme.
293 7957 : bfgs(:)=den1*din(:)-den3*hdelta(:)
294 :
295 : !B.F.G.S. updating formula
296 7957 : do ii=1,ndim
297 224620 : do jj=1,ndim
298 : hessin(ii,jj)=hessin(ii,jj) +den1*din(ii)*din(jj) &
299 223968 : & -den3*hdelta(ii)*hdelta(jj) +den2*bfgs(ii)*bfgs(jj)
300 : end do
301 : end do
302 :
303 652 : end subroutine hessupdt
304 : !!***
305 :
306 : !----------------------------------------------------------------------
307 :
308 : !!****f* m_bfgs/brdene
309 : !! NAME
310 : !! brdene
311 : !!
312 : !! FUNCTION
313 : !! Update vin according to the Broyden formula, combined
314 : !! with a line minimisation that take into account the total energies.
315 : !! Also transfer vin to vin_prev, vout to vout_prev, and etotal to etotal_prev
316 : !! Could see Numerical Recipes (Fortran), 1986, page 307,
317 : !! as well as Schlegel, J. Comp. Chem. 3, 214 (1982) [[cite:Schlegel1982]].
318 : !!
319 : !! INPUTS
320 : !! etotal=new total energy (no meaning at output)
321 : !! hessin(ndim,ndim)=hessian matrix
322 : !! ndim=size of the hessian and vectors
323 : !! vout(ndim)=new output vector (no meaning at output)
324 : !!
325 : !! OUTPUT
326 : !! (see side effects)
327 : !!
328 : !! SIDE EFFECTS
329 : !! etotal_prev=previous total energy; contains input etotal at output
330 : !! vin(ndim)=new input vector; updated at output
331 : !! vin_prev(ndim)=previous input vector; contains input vin at output
332 : !! vout_prev(ndim)=previous output vector; contains input vout at output
333 : !!
334 : !! SOURCE
335 :
336 247 : subroutine brdene(etotal,etotal_prev,hessin,ndim,vin,vin_prev,vout,vout_prev)
337 :
338 : !Arguments ------------------------------------
339 : !scalars
340 : integer,intent(in) :: ndim
341 : real(dp),intent(in) :: etotal
342 : real(dp),intent(inout) :: etotal_prev
343 : !arrays+
344 : real(dp),intent(in) :: hessin(ndim,ndim),vout(ndim)
345 : real(dp),intent(inout) :: vin(ndim),vin_prev(ndim),vout_prev(ndim)
346 :
347 : !Local variables-------------------------------
348 : !scalars
349 : integer :: idim,brd_status
350 : real(dp) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_1,dedv_2,dedv_min
351 : real(dp) :: dedv_predict,etotal_1,etotal_2,etotal_predict,lambda_1,lambda_2
352 : real(dp) :: lambda_predict
353 : !arrays
354 247 : real(dp),allocatable :: dvin(:),vin_min(:),vout_min(:)
355 :
356 : !***************************************************************************
357 :
358 741 : ABI_MALLOC(dvin,(ndim))
359 494 : ABI_MALLOC(vin_min,(ndim))
360 494 : ABI_MALLOC(vout_min,(ndim))
361 :
362 247 : lambda_1=1.0_dp ; lambda_2=0.0_dp
363 247 : etotal_1=etotal ; etotal_2=etotal_prev
364 1750 : dvin(:)=vin(:)-vin_prev(:)
365 1750 : dedv_1=dot_product(vout,dvin)
366 1750 : dedv_2=dot_product(vout_prev,dvin)
367 : call findmin(dedv_1,dedv_2,dedv_predict,&
368 : & d2edv2_1,d2edv2_2,d2edv2_predict,&
369 : & etotal_1,etotal_2,etotal_predict,&
370 247 : & lambda_1,lambda_2,lambda_predict,brd_status)
371 :
372 : !DEBUG : comes back to usual BFGS !
373 : !lambda_predict=1.0_dp
374 : !dedv_predict=dedv_1
375 : !ENDDEBUG
376 :
377 : !Generates vin at the minimum, and an interpolated vout, modified
378 : !to have the right value of dedv_predict, from findmin.
379 1750 : vin_min(:)=vin_prev(:)+lambda_predict*dvin(:)
380 1750 : vout_min(:)=vout_prev(:)+lambda_predict*(vout(:)-vout_prev(:))
381 247 : dedv_min=dedv_2+lambda_predict*(dedv_1-dedv_2)
382 : !Modify vout_min in order to impose dedv_predict
383 3500 : vout_min(:)=vout_min(:)+dvin(:)*(dedv_predict-dedv_min)/dot_product(dvin,dvin)
384 :
385 : !Previous cartesian coordinates
386 247 : etotal_prev=etotal
387 1750 : vin_prev(:)=vin(:)
388 :
389 : !New atomic cartesian coordinates are obtained from vin, hessin and vout
390 1750 : vin(:)=vin_min(:)
391 1750 : do idim=1,ndim
392 10987 : vin(:)=vin(:)-hessin(:,idim)*vout_min(idim)
393 : end do
394 :
395 : !Previous atomic forces
396 1750 : vout_prev(:)=vout(:)
397 :
398 247 : ABI_FREE(dvin)
399 247 : ABI_FREE(vin_min)
400 247 : ABI_FREE(vout_min)
401 :
402 247 : end subroutine brdene
403 : !!***
404 :
405 : END MODULE m_bfgs
406 : !!***
|
