Branch data 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 : : !!***
|
