Branch data Line data Source code
1 : : !!****m* ABINIT/m_mep
2 : : !! NAME
3 : : !! m_mep
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides several routines and datatypes for the
7 : : !! Minimal Energy Path (MEP) search implementation.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2012-2026 ABINIT group (MT)
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_mep
24 : :
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_errors
28 : : use m_dtset
29 : : use m_xmpi
30 : :
31 : : use defs_abitypes, only : MPI_type
32 : : use m_geometry, only : gred2fcart, fcart2gred, xcart2xred, xred2xcart, metric
33 : : use m_bfgs, only : hessupdt
34 : : use m_results_img, only : results_img_type, gather_array_img
35 : :
36 : : implicit none
37 : :
38 : : private
39 : :
40 : : !public procedures
41 : : public :: mep_init
42 : : public :: mep_destroy
43 : : public :: mep_steepest
44 : : public :: mep_qmin
45 : : public :: mep_lbfgs
46 : : public :: mep_gbfgs
47 : : public :: mep_rk4
48 : : public :: mep_img_dotp
49 : : public :: mep_img_norm
50 : : public :: mep_img_dotp_red
51 : : public :: mep_img_norm_red
52 : : !!***
53 : :
54 : : !!****t* m_mep/mep_type
55 : : !! NAME
56 : : !! mep_type
57 : : !!
58 : : !! FUNCTION
59 : : !! Datatype with the variables required to perform MEP search
60 : : !!
61 : : !! SOURCE
62 : :
63 : : type,public :: mep_type
64 : : ! Scalars
65 : : integer :: cineb_start ! Starting iteration for the CI-NEB
66 : : integer :: mep_solver ! Selection of a solver for the ODE
67 : : integer :: neb_algo ! Selection of the variant of the NEB method
68 : : integer :: neb_cell_algo ! Selection of the cell modification algorithm
69 : : integer :: string_algo ! Selection of the variant of the String Method
70 : : real(dp) :: fxcartfactor ! Time step for steepest descent or RK4
71 : : real(dp) :: mep_mxstep ! Selection of a max. step size for the ODE
72 : : ! Arrays
73 : : integer,pointer :: iatfix(:,:)=>null() ! Atoms to fix (this pointer is associated with dtset%iatfix)
74 : : real(dp) :: neb_spring(2) ! Spring constants for the NEB method
75 : : real(dp),allocatable :: bfgs_xprev(:,:,:) ! BFGS storage (prev positions)
76 : : real(dp),allocatable :: bfgs_fprev(:,:,:) ! BFGS storage (prev forces)
77 : : real(dp),allocatable :: gbfgs_hess(:,:) ! global-BFGS storage (Hessian matrix)
78 : : real(dp),allocatable :: lbfgs_hess(:,:,:) ! local-BFGS storage (Hessian matrix)
79 : : real(dp),allocatable :: qmin_vel(:,:,:) ! Quick-min algo storage (velocities)
80 : : real(dp),allocatable :: rk4_xcart1(:,:,:) ! 4th-order Runge-Kutta storage
81 : : real(dp),allocatable :: rk4_fcart1(:,:,:) ! 4th-order Runge-Kutta storage
82 : : real(dp),allocatable :: rk4_fcart2(:,:,:) ! 4th-order Runge-Kutta storage
83 : : real(dp),allocatable :: rk4_fcart3(:,:,:) ! 4th-order Runge-Kutta storage
84 : : real(dp),pointer :: rprimd_start(:,:,:) ! real space primitive translations at start of the MEP search
85 : : end type mep_type
86 : :
87 : : !Public constants
88 : : !NEB algorithms
89 : : integer, public :: NEB_ALGO_STANDARD = 0
90 : : integer, public :: NEB_ALGO_IMPROVED_TAN = 1
91 : : integer, public :: NEB_ALGO_CINEB = 2
92 : : !Variable-cell NEB algorithms
93 : : integer, public :: NEB_CELL_ALGO_NONE = 0
94 : : integer, public :: NEB_CELL_ALGO_GSSNEB = 1
95 : : integer, public :: NEB_CELL_ALGO_VCNEB = 2
96 : : !String method algorithms
97 : : integer, public :: STRING_ALGO_ORIGINAL = 0
98 : : integer, public :: STRING_ALGO_SIMPLIFIED_EQUAL = 1
99 : : integer, public :: STRING_ALGO_SIMPLIFIED_ENERGY = 2
100 : : !MEP solvers
101 : : integer, public :: MEP_SOLVER_STEEPEST = 0
102 : : integer, public :: MEP_SOLVER_QUICKMIN = 1
103 : : integer, public :: MEP_SOLVER_LBFGS = 2
104 : : integer, public :: MEP_SOLVER_GBFGS = 3
105 : : integer, public :: MEP_SOLVER_RK4 = 4
106 : : !!***
107 : :
108 : : CONTAINS
109 : :
110 : : !===========================================================
111 : : !!***
112 : :
113 : : !!****f* m_mep/mep_init
114 : : !! NAME
115 : : !! mep_init
116 : : !!
117 : : !! FUNCTION
118 : : !! Initialize a datastructure of type mep_type.
119 : : !!
120 : : !! INPUTS
121 : : !! dtset <type(dataset_type)>=all input variables in current dataset
122 : : !!
123 : : !! OUTPUT
124 : : !!
125 : : !! SIDE EFFECTS
126 : : !! mep_param=datastructure of type mep_type.
127 : : !! several parameters for Minimal Energy Path (MEP) search.
128 : : !!
129 : : !! SOURCE
130 : :
131 : 4008 : subroutine mep_init(dtset,mep_param)
132 : :
133 : : !Arguments ------------------------------------
134 : : !scalars
135 : : type(dataset_type),target,intent(in) :: dtset
136 : : type(mep_type),intent(inout) :: mep_param
137 : :
138 : : !************************************************************************
139 : :
140 [ + + ]: 4008 : if((dtset%imgmov==1).or.(dtset%imgmov==2).or.(dtset%imgmov==5))then
141 : 29 : mep_param%cineb_start = dtset%cineb_start
142 : 29 : mep_param%mep_solver = dtset%mep_solver
143 : 29 : mep_param%neb_algo = dtset%neb_algo
144 : 29 : mep_param%neb_cell_algo = dtset%neb_cell_algo
145 : 29 : mep_param%string_algo = dtset%string_algo
146 : 29 : mep_param%fxcartfactor = dtset%fxcartfactor
147 : 29 : mep_param%mep_mxstep = dtset%mep_mxstep
148 [ + + ]: 87 : mep_param%neb_spring = dtset%neb_spring
149 : 29 : mep_param%iatfix =>dtset%iatfix
150 : 29 : mep_param%rprimd_start =>dtset%rprimd_orig
151 : : !TODO: Q. DELACROIX DID PUT THIS - TO BE CHECKED
152 : : ! do ii=1,dtset%nimage
153 : : ! mep_param%rprimd_start(:,:,iimage)=dtset%rprimd_orig(:,:,1)
154 : : ! end do
155 : : else
156 : 3979 : mep_param%cineb_start = -1
157 : 3979 : mep_param%mep_solver = -1
158 : 3979 : mep_param%neb_algo = -1
159 : 3979 : mep_param%neb_cell_algo = -1
160 : 3979 : mep_param%string_algo = -1
161 : 3979 : mep_param%fxcartfactor = zero
162 : 3979 : mep_param%mep_mxstep = 100._dp
163 [ + + ]: 11937 : mep_param%neb_spring = zero
164 : 3979 : nullify(mep_param%iatfix)
165 : 3979 : nullify(mep_param%rprimd_start)
166 : : end if
167 : :
168 : 4008 : end subroutine mep_init
169 : : !!***
170 : :
171 : : !----------------------------------------------------------------------
172 : :
173 : : !!****f* m_mep/mep_destroy
174 : : !! NAME
175 : : !! mep_destroy
176 : : !!
177 : : !! FUNCTION
178 : : !! Destroy the content of a datastructure of type mep_type.
179 : : !!
180 : : !! INPUTS
181 : : !!
182 : : !! OUTPUT
183 : : !!
184 : : !! SIDE EFFECTS
185 : : !! mep_param=datastructure of type mep_type.
186 : : !! several parameters for Minimal Energy Path (MEP) search.
187 : : !!
188 : : !! SOURCE
189 : :
190 : 4008 : subroutine mep_destroy(mep_param)
191 : :
192 : : !Arguments ------------------------------------
193 : : !scalars
194 : : type(mep_type),intent(inout) :: mep_param
195 : :
196 : : !************************************************************************
197 : :
198 [ - + ]: 4008 : ABI_SFREE(mep_param%bfgs_xprev)
199 [ - + ]: 4008 : ABI_SFREE(mep_param%gbfgs_hess)
200 [ - + ]: 4008 : ABI_SFREE(mep_param%bfgs_fprev)
201 [ - + ]: 4008 : ABI_SFREE(mep_param%lbfgs_hess)
202 [ - + ]: 4008 : ABI_SFREE(mep_param%qmin_vel)
203 [ - + ]: 4008 : ABI_SFREE(mep_param%rk4_xcart1)
204 [ - + ]: 4008 : ABI_SFREE(mep_param%rk4_fcart1)
205 [ - + ]: 4008 : ABI_SFREE(mep_param%rk4_fcart2)
206 [ - + ]: 4008 : ABI_SFREE(mep_param%rk4_fcart3)
207 : :
208 : 4008 : nullify(mep_param%iatfix)
209 : 4008 : nullify(mep_param%rprimd_start)
210 : :
211 : 4008 : end subroutine mep_destroy
212 : : !!***
213 : :
214 : : !----------------------------------------------------------------------
215 : :
216 : : !!****f* m_mep/mep_steepest
217 : : !! NAME
218 : : !! mep_steepest
219 : : !!
220 : : !! FUNCTION
221 : : !! Make a path (string of images) evolve according to a steepest descent algorithm
222 : : !!
223 : : !! INPUTS
224 : : !! fcart(3,natom_eff,nimage)=cartesian forces in each image along the path
225 : : !! list_dynimage(nimage)=list of dynamical images.
226 : : !! mep_param=datastructure of type mep_type.
227 : : !! several parameters for Minimal Energy Path (MEP) search.
228 : : !! natom=number of atoms
229 : : !! natom_eff="effective" number of atoms, including possibly the unit cell vectors
230 : : !! ndynimage=number of dynamical images along the path
231 : : !! nimage=number of images (including static ones)
232 : : !! results_img(nimage)=datastructure that hold data for each image
233 : : !! (positions, forces, energy, ...)
234 : : !! rprimd(3,3,nimage)=dimensional primitive translations for each image along the path
235 : : !! [use_reduced_coord]=force the use of reduced coordinates instead of cartesian ones
236 : : !! [rprimd_start(3,3,nimage)]=real space primitive translations at start of the MEP search
237 : : !! Mandatory if natom_eff=natom+3
238 : : !! [strainfact(nimage)]=only for variable cell algorithms. Factor applied to strains
239 : : !! to align their dimension to atomic positions.
240 : : !! Only valid when use_reduced_coordinates=.false.
241 : : !! (essentially used for the GSS-NEB method)
242 : : !!!
243 : : !! OUTPUT
244 : : !!
245 : : !! SIDE EFFECTS
246 : : !! xcart(3,natom_eff,nimage)=cartesian coordinates of atoms in each image along the path
247 : : !! before and after time evolution. If natom_eff=natom+3, then the 3 last "atoms"
248 : : !! are the primitive vectors of the cell
249 : : !! xred(3,natom_eff,nimage)=reduced coordinates of atoms in each image along the path
250 : : !! before and after time evolution. If natom_eff=natom+3, then the 3 last "atoms"
251 : : !! are the primitive vectors of the cell
252 : : !!
253 : : !! SOURCE
254 : :
255 : 138 : subroutine mep_steepest(fcart,list_dynimage,mep_param,natom,natom_eff,ndynimage,nimage,rprimd,xcart,xred, &
256 : : & use_reduced_coord,rprimd_start,strainfact) ! optional arguments
257 : :
258 : : !Arguments ------------------------------------
259 : : !scalars
260 : : integer,intent(in) :: natom,natom_eff,ndynimage,nimage
261 : : logical,intent(in),optional :: use_reduced_coord
262 : : type(mep_type),intent(in) :: mep_param
263 : : !arrays
264 : : integer,intent(in) :: list_dynimage(ndynimage)
265 : : real(dp),intent(in) :: fcart(3,natom_eff,nimage)
266 : : real(dp),intent(inout) :: rprimd(3,3,nimage),xcart(3,natom_eff,nimage),xred(3,natom_eff,nimage)
267 : : real(dp),intent(in),optional :: rprimd_start(3,3,nimage),strainfact(nimage)
268 : : !Local variables-------------------------------
269 : : !scalars
270 : : integer :: iatom,idynimage,iimage
271 : : logical :: use_reduced_coord_
272 : : real(dp) :: stepsize
273 : : character(len=500) :: msg
274 : : !arrays
275 : : real(dp),parameter :: identity_real(3,3)=reshape([one,zero,zero,zero,one,zero,zero,zero,one],[3,3])
276 : : real(dp) :: mat3(3,3)
277 : 138 : real(dp),allocatable :: xred_old(:,:),xstep(:,:)
278 : :
279 : : !************************************************************************
280 : :
281 : 138 : use_reduced_coord_=.false.
282 [ + + ]: 138 : if (present(use_reduced_coord)) use_reduced_coord_=use_reduced_coord
283 : :
284 [ + + - + ]: 138 : if (natom_eff>=natom+3.and.(.not.present(rprimd_start))) then
285 : 0 : ABI_BUG("Mandatory arg missing (rprimd_start)!")
286 : : end if
287 : :
288 [ - + + - : 414 : ABI_MALLOC(xred_old,(3,natom_eff))
+ - + - -
+ - + ]
289 [ - + + - : 276 : ABI_MALLOC(xstep,(3,natom_eff))
- + - + ]
290 : :
291 [ + + ]: 804 : do idynimage=1,ndynimage
292 : 666 : iimage=list_dynimage(idynimage)
293 [ + + + + ]: 6534 : xred_old(:,:)=xred(:,:,iimage)
294 : :
295 : : ! Compute image step
296 : : ! Note that one uses fcart, for which the sum of forces on all atoms vanish
297 [ + + + + ]: 6534 : xstep(:,:)=mep_param%fxcartfactor*fcart(:,:,iimage)
298 : 666 : stepsize=mep_img_norm(xstep)
299 [ - + ]: 666 : if (stepsize>=mep_param%mep_mxstep) then
300 [ # # # # ]: 0 : xstep=xstep*mep_param%mep_mxstep/stepsize
301 : 0 : write(msg,'(a,i3,a)') " Restricting step size of image ",iimage,"."
302 : 0 : call wrtout(std_out,msg,'COLL')
303 : 0 : call wrtout(ab_out ,msg,'COLL')
304 : : end if
305 : :
306 : : ! Update positions
307 [ - + ]: 666 : if (use_reduced_coord_) then
308 [ # # # # ]: 0 : xred(:,:,iimage)=xred(:,:,iimage)+xstep(:,:)
309 : 0 : call xred2xcart(natom,rprimd(:,:,iimage),xcart(:,1:natom,iimage),xred(:,1:natom,iimage))
310 : : else
311 [ + + + + ]: 6534 : xcart(:,:,iimage)=xcart(:,:,iimage)+xstep(:,:)
312 : 666 : call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,1:natom,iimage),xred(:,1:natom,iimage))
313 : : end if
314 : :
315 : : ! Update unit cell vectors if they are included in the list of "atoms"
316 [ + + ]: 666 : if (natom_eff>=natom+3) then
317 [ - + ]: 45 : if (use_reduced_coord_) then
318 [ # # # # ]: 0 : mat3(1:3,1:3)=identity_real(1:3,1:3)+xred(:,natom+1:natom+3,iimage)
319 [ # # # # : 0 : rprimd(:,:,iimage)=matmul(mat3(:,:),rprimd_start(:,:,iimage))
# # # # #
# ]
320 : : else
321 [ + + + + ]: 585 : mat3(1:3,1:3)=xcart(1:3,natom+1:natom+3,iimage)
322 [ + - + + : 585 : if (present(strainfact)) mat3(1:3,1:3)=mat3(1:3,1:3)/strainfact(iimage)
+ + ]
323 : : rprimd(:,:,iimage)=matmul(rprimd_start(:,:,iimage),mat3(:,:)) &
324 [ + + + + : 2880 : & +rprimd_start(:,:,iimage)
+ + + + +
+ + + +
+ ]
325 : : end if
326 : : end if
327 : :
328 : : ! In case atom is fixed, we restore its previous value ; forbidden if variable cell
329 [ + + ]: 804 : if (mep_param%neb_cell_algo==NEB_CELL_ALGO_NONE) then
330 [ + + ]: 1863 : do iatom=1,natom
331 [ + + + + : 3834 : if (any(mep_param%iatfix(:,iatom)==1)) then
+ + ]
332 [ + + + + ]: 2340 : where(mep_param%iatfix(:,iatom)==1)
333 : : xred(:,iatom,iimage)=xred_old(:,iatom)
334 : : end where
335 : 585 : call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage))
336 : : end if
337 : : end do
338 : : end if
339 : :
340 : : end do
341 : :
342 [ - + ]: 138 : ABI_FREE(xred_old)
343 [ - + ]: 138 : ABI_FREE(xstep)
344 : :
345 : 138 : end subroutine mep_steepest
346 : : !!***
347 : :
348 : : !----------------------------------------------------------------------
349 : :
350 : : !!****f* m_mep/mep_qmin
351 : : !! NAME
352 : : !! mep_qmin
353 : : !!
354 : : !! FUNCTION
355 : : !! Make a path (string of images) evolve according to a quick-minimizer algorithm
356 : : !!
357 : : !! INPUTS
358 : : !! fcart(3,natom,nimage)=cartesian forces in each image along the path
359 : : !! itime=time step
360 : : !! list_dynimage(nimage)=list of dynamical images.
361 : : !! mep_param=datastructure of type mep_type.
362 : : !! several parameters for Minimal Energy Path (MEP) search.
363 : : !! natom=number of atoms
364 : : !! ndynimage=number of dynamical images along the path
365 : : !! nimage=number of images (including static ones)
366 : : !! results_img(nimage)=datastructure that hold data for each image
367 : : !! (positions, forces, energy, ...)
368 : : !! rprimd(3,3,nimage)=dimensional primitive translations for each image along the path
369 : : !!
370 : : !! OUTPUT
371 : : !!
372 : : !! SIDE EFFECTS
373 : : !! xcart(3,natom,nimage)=cartesian coordinates of atoms in each image along the path
374 : : !! before and after time evolution
375 : : !! xred(3,natom,nimage)=reduced coordinates of atoms in each image along the path
376 : : !! before and after time evolution
377 : : !!
378 : : !! SOURCE
379 : :
380 : 0 : subroutine mep_qmin(fcart,itime,list_dynimage,mep_param,natom,ndynimage,nimage,rprimd,xcart,xred)
381 : :
382 : : !Arguments ------------------------------------
383 : : !scalars
384 : : integer,intent(in) :: itime,natom,ndynimage,nimage
385 : : type(mep_type),intent(inout) :: mep_param
386 : : !arrays
387 : : integer,intent(in) :: list_dynimage(ndynimage)
388 : : real(dp),intent(in) :: fcart(3,natom,nimage),rprimd(3,3,nimage)
389 : : real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage)
390 : : !Local variables-------------------------------
391 : : !scalars
392 : : integer :: iatom,idynimage,iimage
393 : : real(dp) :: stepsize,vdotf
394 : : character(len=500) :: msg
395 : : !arrays
396 : : real(dp) :: vel_red(3)
397 : 0 : real(dp),allocatable :: xred_old(:,:),xstep(:,:)
398 : :
399 : : !***********************************************************************
400 : :
401 : : !Allocate history array (at first time step)
402 [ # # ]: 0 : if (itime==1) then
403 [ # # ]: 0 : ABI_SFREE(mep_param%qmin_vel)
404 [ # # # # : 0 : ABI_MALLOC(mep_param%qmin_vel,(3,natom,ndynimage))
# # # # #
# # # # #
# # # # ]
405 [ # # # # : 0 : mep_param%qmin_vel=zero
# # ]
406 : : end if
407 : :
408 [ # # # # : 0 : ABI_MALLOC(xred_old,(3,natom))
# # # # #
# # # ]
409 [ # # # # : 0 : ABI_MALLOC(xstep,(3,natom))
# # # # ]
410 : :
411 [ # # ]: 0 : do idynimage=1,ndynimage
412 : 0 : iimage=list_dynimage(idynimage)
413 [ # # # # ]: 0 : xred_old(:,:)=xred(:,:,iimage)
414 : :
415 : : ! Compute velocities
416 : 0 : vdotf=mep_img_dotp(mep_param%qmin_vel(:,:,idynimage),fcart(:,:,iimage))
417 [ # # ]: 0 : if (vdotf>=zero) then
418 : : mep_param%qmin_vel(:,:,idynimage)=vdotf*fcart(:,:,iimage) &
419 [ # # # # ]: 0 : & /mep_img_norm(fcart(:,:,iimage))
420 : : else
421 [ # # # # ]: 0 : mep_param%qmin_vel(:,:,idynimage)=zero
422 : 0 : write(msg,'(a,i3,a)') " Setting velocities of image ",iimage," to zero."
423 : 0 : call wrtout(std_out,msg,'COLL')
424 : 0 : call wrtout(ab_out ,msg,'COLL')
425 : : end if
426 : : mep_param%qmin_vel(:,:,idynimage)=mep_param%qmin_vel(:,:,idynimage) &
427 [ # # # # ]: 0 : & +mep_param%fxcartfactor*fcart(:,:,iimage)
428 : :
429 : : ! Compute image step
430 [ # # # # ]: 0 : xstep(:,:)=mep_param%fxcartfactor*mep_param%qmin_vel(:,:,idynimage)
431 : 0 : stepsize=mep_img_norm(xstep)
432 [ # # ]: 0 : if (stepsize>=mep_param%mep_mxstep) then
433 [ # # # # ]: 0 : xstep=xstep*mep_param%mep_mxstep/stepsize
434 : 0 : write(msg,'(a,i3,a)') " Restricting step size of image ",iimage,"."
435 : 0 : call wrtout(std_out,msg,'COLL')
436 : 0 : call wrtout(ab_out ,msg,'COLL')
437 : : end if
438 : :
439 : : ! Update positions
440 [ # # # # ]: 0 : xcart(:,:,iimage)=xcart(:,:,iimage)+xstep(:,:)
441 : 0 : call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage))
442 : :
443 : : ! In case atom is fixed, we restore its previous value
444 [ # # ]: 0 : do iatom=1,natom
445 [ # # # # : 0 : if (any(mep_param%iatfix(:,iatom)==1)) then
# # ]
446 : 0 : call xcart2xred(1,rprimd(:,:,iimage),mep_param%qmin_vel(:,iatom,idynimage),vel_red)
447 [ # # # # : 0 : where(mep_param%iatfix(:,iatom)==1)
# # # # #
# ]
448 : : xred(:,iatom,iimage)=xred_old(:,iatom)
449 : : vel_red(:)=zero
450 : : end where
451 : 0 : call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage))
452 : 0 : call xred2xcart(1,rprimd(:,:,iimage),mep_param%qmin_vel(:,iatom,idynimage),vel_red)
453 : : end if
454 : : end do
455 : :
456 : : end do
457 : :
458 [ # # ]: 0 : ABI_FREE(xred_old)
459 [ # # ]: 0 : ABI_FREE(xstep)
460 : :
461 : 0 : end subroutine mep_qmin
462 : : !!***
463 : :
464 : : !----------------------------------------------------------------------
465 : :
466 : : !!****f* m_mep/mep_lbfgs
467 : : !! NAME
468 : : !! mep_lbfgs
469 : : !!
470 : : !! FUNCTION
471 : : !! Make a path (string of images) evolve according to a
472 : : !! local Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm
473 : : !!
474 : : !! INPUTS
475 : : !! itime=time step
476 : : !! list_dynimage(nimage)=list of dynamical images.
477 : : !! mep_param=datastructure of type mep_type.
478 : : !! several parameters for Minimal Energy Path (MEP) search.
479 : : !! natom=number of atoms
480 : : !! ndynimage=number of dynamical images along the path
481 : : !! nimage=number of images (including static ones)
482 : : !! rprimd(3,3,nimage)=dimensional primitive translations for each image along the path
483 : : !!
484 : : !! OUTPUT
485 : : !!
486 : : !! SIDE EFFECTS
487 : : !! mep_param=datastructure of type mep_type.
488 : : !! History for Runge-Kutta algorithm is filled up
489 : : !! xcart(3,natom,nimage)=cartesian coordinates of atoms in each image along the path
490 : : !! before and after time evolution
491 : : !! xred(3,natom,nimage)=reduced coordinates of atoms in each image along the path
492 : : !! before and after time evolution
493 : : !!
494 : : !! NOTES
495 : : !! Could see Numerical Recipes (Fortran), 1986, page 307.
496 : : !!
497 : : !! SOURCE
498 : :
499 : 0 : subroutine mep_lbfgs(fcart,itime,list_dynimage,mep_param,natom,ndynimage,&
500 : 0 : & nimage,rprimd,xcart,xred)
501 : :
502 : : !Arguments ------------------------------------
503 : : !scalars
504 : : integer,intent(in) :: itime,natom,ndynimage,nimage
505 : : type(mep_type),intent(inout) :: mep_param
506 : : !arrays
507 : : integer,intent(in) :: list_dynimage(ndynimage)
508 : : real(dp),intent(in) :: rprimd(3,3,nimage)
509 : : real(dp),intent(in) :: fcart(3,natom,nimage)
510 : : real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage)
511 : : !Local variables-------------------------------
512 : : !scalars
513 : : integer :: iatom,idynimage,ii,iimage,indi,indj,jatom,jj
514 : : logical :: reset
515 : : real(dp),parameter :: initial_Hessian=1._dp ! in Bohr^2/Hartree
516 : : real(dp) :: dot1,dot2,stepsize,ucvol
517 : : character(len=500) :: msg
518 : : !arrays
519 : : real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3)
520 : 0 : real(dp),allocatable :: gred(:,:),xstep(:,:)
521 : :
522 : : !************************************************************************
523 : :
524 : : !Allocate history array (at first time step)
525 [ # # ]: 0 : if (itime==1) then
526 [ # # ]: 0 : ABI_SFREE(mep_param%bfgs_xprev)
527 [ # # ]: 0 : ABI_SFREE(mep_param%bfgs_fprev)
528 [ # # ]: 0 : ABI_SFREE(mep_param%lbfgs_hess)
529 [ # # # # : 0 : ABI_MALLOC(mep_param%bfgs_xprev,(3,natom,ndynimage))
# # # # #
# # # # #
# # # # ]
530 [ # # # # : 0 : ABI_MALLOC(mep_param%bfgs_fprev,(3,natom,ndynimage))
# # # # #
# # # #
# ]
531 [ # # # # : 0 : ABI_MALLOC(mep_param%lbfgs_hess,(3*natom,3*natom,ndynimage))
# # # # #
# # # # #
# # # # #
# # # #
# ]
532 [ # # # # : 0 : mep_param%bfgs_xprev=zero
# # ]
533 [ # # # # : 0 : mep_param%bfgs_fprev=zero
# # ]
534 : : end if
535 : :
536 : : !Temporary storage
537 [ # # # # : 0 : ABI_MALLOC(gred,(3,natom))
# # # # #
# # # ]
538 [ # # # # : 0 : ABI_MALLOC(xstep,(3,natom))
# # # # ]
539 : :
540 : : !Loop over images
541 [ # # ]: 0 : do idynimage=1,ndynimage
542 : 0 : iimage=list_dynimage(idynimage)
543 : 0 : call metric(gmet,gprimd,-1,rmet,rprimd(:,:,iimage),ucvol)
544 : 0 : call fcart2gred(fcart(:,:,iimage),gred,rprimd(:,:,iimage),natom)
545 : :
546 : : ! Test if a reset is needed
547 : 0 : reset=.false.
548 [ # # ]: 0 : if (itime>1) then
549 : 0 : dot1=mep_img_dotp(mep_param%bfgs_fprev(:,:,idynimage),gred)
550 : : dot2=mep_img_dotp(mep_param%bfgs_fprev(:,:,idynimage), &
551 : 0 : & mep_param%bfgs_fprev(:,:,idynimage))
552 : : ! dot1=mep_img_dotp_red(rmet,mep_param%bfgs_fprev(:,:,idynimage),gred)
553 : : ! dot2=mep_img_dotp_red(rmet,mep_param%bfgs_fprev(:,:,idynimage), &
554 : : !& mep_param%bfgs_fprev(:,:,idynimage))
555 [ # # # # ]: 0 : reset=((dot2<two*abs(dot1)).or.abs(dot2)<tol8)
556 : : if (reset) then
557 : 0 : write(msg,'(a,i3,a)') " Resetting Hessian matrix for image ",iimage,"."
558 : 0 : call wrtout(std_out,msg,'COLL')
559 : 0 : call wrtout(ab_out ,msg,'COLL')
560 : : end if
561 : : end if
562 : :
563 : : ! ===> First step or reset: initialize the Hessian matrix (in reduced coordinates)
564 [ # # # # ]: 0 : if (itime==1.or.reset) then
565 [ # # # # ]: 0 : mep_param%lbfgs_hess(:,:,idynimage)=zero
566 [ # # ]: 0 : do iatom=1,natom
567 : 0 : indi=3*(iatom-1)
568 [ # # ]: 0 : do ii=1,3
569 [ # # ]: 0 : do jj=1,3
570 [ # # # # ]: 0 : if (mep_param%iatfix(ii,iatom)==0.and. &
571 : 0 : & mep_param%iatfix(jj,iatom)==0) then
572 : 0 : mep_param%lbfgs_hess(indi+ii,indi+jj,idynimage)=gmet(ii,jj)*initial_Hessian
573 : : end if
574 : : end do
575 : : end do
576 : : end do
577 : :
578 : : ! ===> Other steps: update the Hessian matrix
579 : : else
580 : : call hessupdt(mep_param%lbfgs_hess(:,:,idynimage),&
581 : : & mep_param%iatfix,natom,3*natom, &
582 : : xred(:,:,iimage),mep_param%bfgs_xprev(:,:,idynimage),&
583 [ # # # # : 0 : gred(:,:),mep_param%bfgs_fprev(:,:,idynimage))
# # # # #
# ]
584 : : end if
585 : :
586 : : ! Update history
587 [ # # # # ]: 0 : mep_param%bfgs_xprev(:,:,idynimage)=xred(:,:,iimage)
588 [ # # # # ]: 0 : mep_param%bfgs_fprev(:,:,idynimage)=gred(:,:)
589 : :
590 : : ! Compute image step
591 [ # # # # ]: 0 : xstep=zero
592 [ # # ]: 0 : do iatom=1,natom
593 : 0 : indi=3*(iatom-1)
594 [ # # ]: 0 : do ii=1,3
595 [ # # ]: 0 : do jatom=1,natom
596 : 0 : indj=3*(jatom-1)
597 [ # # ]: 0 : do jj=1,3
598 : : xstep(ii,iatom)=xstep(ii,iatom) &
599 : 0 : & -mep_param%lbfgs_hess(indi+ii,indj+jj,idynimage)*gred(jj,jatom)
600 : : end do
601 : : end do
602 : : end do
603 : : end do
604 : :
605 : : ! Restrict image step size
606 : 0 : stepsize=mep_img_norm_red(rmet,xstep)
607 [ # # ]: 0 : if (stepsize>=mep_param%mep_mxstep) then
608 [ # # # # ]: 0 : xstep=xstep*mep_param%mep_mxstep/stepsize
609 : 0 : write(msg,'(a,i3,a)') " Restricting BFGS step size of image ",iimage,"."
610 : 0 : call wrtout(std_out,msg,'COLL')
611 : 0 : call wrtout(ab_out ,msg,'COLL')
612 : : end if
613 : :
614 : : ! Update positions
615 [ # # # # ]: 0 : xred(:,:,iimage)=xred(:,:,iimage)+xstep(:,:)
616 : :
617 : : ! In case atom is fixed, we restore its previous value
618 [ # # # # : 0 : where(mep_param%iatfix(:,:)==1)
# # ]
619 : : xred(:,:,iimage)=mep_param%bfgs_xprev(:,:,idynimage)
620 : : end where
621 : :
622 : 0 : call xred2xcart(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage))
623 : :
624 : : !End loop over images
625 : : end do
626 : :
627 [ # # ]: 0 : ABI_FREE(gred)
628 [ # # ]: 0 : ABI_FREE(xstep)
629 : :
630 : 0 : end subroutine mep_lbfgs
631 : : !!***
632 : :
633 : : !----------------------------------------------------------------------
634 : :
635 : : !!****f* m_mep/mep_gbfgs
636 : : !! NAME
637 : : !! mep_gbfgs
638 : : !!
639 : : !! FUNCTION
640 : : !! Make a path (string of images) evolve according to a
641 : : !! global Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm
642 : : !!
643 : : !! INPUTS
644 : : !! itime=time step
645 : : !! list_dynimage(nimage)=list of dynamical images.
646 : : !! mep_param=datastructure of type mep_type.
647 : : !! mpi_enreg=MPI-parallelisation information
648 : : !! mep_param=several parameters for Minimal Energy Path (MEP) search.
649 : : !! natom=number of atoms
650 : : !! ndynimage=number of dynamical images along the path
651 : : !! nimage=number of images (including static ones)
652 : : !! nimage_tot=total number of images
653 : : !! rprimd(3,3,nimage)=dimensional primitive translations for each image along the path
654 : : !!
655 : : !! OUTPUT
656 : : !!
657 : : !! SIDE EFFECTS
658 : : !! mep_param=datastructure of type mep_type.
659 : : !! History for Runge-Kutta algorithm is filled up
660 : : !! xcart(3,natom,nimage)=cartesian coordinates of atoms in each image along the path
661 : : !! before and after time evolution
662 : : !! xred(3,natom,nimage)=reduced coordinates of atoms in each image along the path
663 : : !! before and after time evolution
664 : : !!
665 : : !! NOTES
666 : : !! Could see Numerical Recipes (Fortran), 1986, page 307.
667 : : !! Has to work in cartesian coordinates
668 : : !!
669 : : !! SOURCE
670 : :
671 : 0 : subroutine mep_gbfgs(fcart,itime,list_dynimage,mep_param,mpi_enreg,natom,&
672 : 0 : & ndynimage,nimage,nimage_tot,rprimd,xcart,xred)
673 : :
674 : : !Arguments ------------------------------------
675 : : !scalars
676 : : integer,intent(in) :: itime,natom,ndynimage,nimage,nimage_tot
677 : : type(mep_type),intent(inout) :: mep_param
678 : : type(MPI_type),intent(in) :: mpi_enreg
679 : : !arrays
680 : : integer,intent(in) :: list_dynimage(ndynimage)
681 : : real(dp),intent(in) :: fcart(3,natom,nimage),rprimd(3,3,nimage)
682 : : real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage)
683 : : !Local variables-------------------------------
684 : : !scalars
685 : : integer :: iatom,idynimage,ii,iimage,iimage_tot,indi,indj,ierr
686 : : integer :: jdynimage,jatom,jj,mu,ndynimage_tot,nu
687 : : logical :: reset
688 : : real(dp),parameter :: initial_Hessian=1._dp ! in Bohr^2/Hartree
689 : : real(dp) :: dot1,dot2,stepsize,ucvol
690 : : character(len=500) :: msg
691 : : !arrays
692 : 0 : integer,allocatable :: dynimage_tot(:),iatfix_fake(:,:),ind_dynimage_tot(:)
693 : 0 : integer,allocatable :: list_dynimage_tot(:)
694 : : real(dp) :: favg(3),gprimd(3,3),gmet(3,3),rmet(3,3)
695 : 0 : real(dp),allocatable :: buffer(:,:),buffer_all(:,:),gred(:,:)
696 : 0 : real(dp),allocatable :: fcart_all(:,:,:),fcartp_all(:,:,:)
697 : 0 : real(dp),allocatable :: gmet_all(:,:,:),gprimd_all(:,:,:),rprimd_all(:,:,:)
698 : 0 : real(dp),allocatable :: xcart_all(:,:,:),xcartp_all(:,:,:),xred_old(:,:),xstep_all(:,:,:)
699 : :
700 : : !************************************************************************
701 : :
702 : : !Retrieve indexes of all dynamical images
703 [ # # # # : 0 : ABI_MALLOC(ind_dynimage_tot,(nimage_tot))
# # # # #
# # # ]
704 [ # # ]: 0 : if (mpi_enreg%paral_img==1) then
705 [ # # # # : 0 : ABI_MALLOC(dynimage_tot,(nimage_tot))
# # # # ]
706 [ # # ]: 0 : dynimage_tot=0
707 [ # # ]: 0 : do idynimage=1,ndynimage
708 : 0 : iimage=list_dynimage(idynimage)
709 : 0 : iimage_tot=mpi_enreg%my_imgtab(iimage)
710 : 0 : dynimage_tot(iimage_tot)=1
711 : : end do
712 : 0 : call xmpi_sum(dynimage_tot,mpi_enreg%comm_img,ierr)
713 [ # # # # ]: 0 : ndynimage_tot=count(dynimage_tot(:)>0)
714 [ # # # # : 0 : ABI_MALLOC(list_dynimage_tot,(ndynimage_tot))
# # # # #
# # # ]
715 [ # # ]: 0 : idynimage=0;ind_dynimage_tot(:)=-1
716 [ # # ]: 0 : do iimage_tot=1,nimage_tot
717 [ # # ]: 0 : if (dynimage_tot(iimage_tot)>0) then
718 : 0 : idynimage=idynimage+1
719 : 0 : ind_dynimage_tot(iimage_tot)=idynimage
720 : 0 : list_dynimage_tot(idynimage)=iimage_tot
721 : : end if
722 : : end do
723 [ # # ]: 0 : ABI_FREE(dynimage_tot)
724 : : else
725 : 0 : ndynimage_tot=ndynimage
726 [ # # # # : 0 : ABI_MALLOC(list_dynimage_tot,(ndynimage))
# # # # #
# # # ]
727 [ # # ]: 0 : ind_dynimage_tot(:)=-1
728 [ # # ]: 0 : do idynimage=1,ndynimage
729 : 0 : ind_dynimage_tot(list_dynimage(idynimage))=idynimage
730 : 0 : list_dynimage_tot(idynimage)=list_dynimage(idynimage)
731 : : end do
732 : : end if
733 : :
734 : : !Allocate history array (at first time step)
735 [ # # ]: 0 : if (itime==1) then
736 [ # # ]: 0 : ABI_SFREE(mep_param%bfgs_xprev)
737 [ # # ]: 0 : ABI_SFREE(mep_param%bfgs_fprev)
738 [ # # ]: 0 : ABI_SFREE(mep_param%gbfgs_hess)
739 [ # # # # : 0 : ABI_MALLOC(mep_param%bfgs_xprev,(3,natom,ndynimage))
# # # # #
# # # # #
# # # # ]
740 [ # # # # : 0 : ABI_MALLOC(mep_param%bfgs_fprev,(3,natom,ndynimage))
# # # # #
# # # #
# ]
741 [ # # # # : 0 : ABI_MALLOC(mep_param%gbfgs_hess,(3*natom*ndynimage_tot,3*natom*ndynimage_tot))
# # # # #
# # # # #
# # # # #
# ]
742 [ # # # # : 0 : mep_param%bfgs_xprev=zero
# # ]
743 [ # # # # : 0 : mep_param%bfgs_fprev=zero
# # ]
744 : : end if
745 : :
746 : : !Retrieve positions and forces for all images
747 [ # # # # : 0 : ABI_MALLOC(xcart_all,(3,natom,ndynimage_tot))
# # # # #
# # # # #
# # ]
748 [ # # # # : 0 : ABI_MALLOC(fcart_all,(3,natom,ndynimage_tot))
# # # # #
# # # ]
749 [ # # # # : 0 : ABI_MALLOC(xcartp_all,(3,natom,ndynimage_tot))
# # # # #
# # # ]
750 [ # # # # : 0 : ABI_MALLOC(fcartp_all,(3,natom,ndynimage_tot))
# # # # #
# # # ]
751 [ # # # # : 0 : ABI_MALLOC(rprimd_all,(3,3,ndynimage_tot))
# # # # #
# # # ]
752 [ # # # # : 0 : ABI_MALLOC(gprimd_all,(3,3,ndynimage_tot))
# # # # ]
753 [ # # # # : 0 : ABI_MALLOC(gmet_all,(3,3,ndynimage_tot))
# # # # ]
754 [ # # ]: 0 : if (mpi_enreg%paral_img==1) then
755 [ # # # # : 0 : ABI_MALLOC(buffer,(12*natom+27,nimage))
# # # # #
# # # # #
# # # # ]
756 [ # # # # : 0 : ABI_MALLOC(buffer_all,(12*natom+27,nimage_tot))
# # # # #
# # # # #
# # ]
757 [ # # # # ]: 0 : buffer=zero
758 [ # # ]: 0 : do idynimage=1,ndynimage
759 : 0 : iimage=list_dynimage(idynimage)
760 : 0 : call metric(gmet,gprimd,-1,rmet,rprimd(:,:,iimage),ucvol)
761 [ # # ]: 0 : buffer( 1:3 *natom ,iimage)=reshape(xcart(:,:,iimage),(/3*natom/))
762 [ # # ]: 0 : buffer(3 *natom+1:6 *natom ,iimage)=reshape(fcart(:,:,iimage),(/3*natom/))
763 [ # # ]: 0 : buffer(6 *natom+1:9 *natom ,iimage)=reshape(mep_param%bfgs_xprev(:,:,idynimage),(/3*natom/))
764 [ # # ]: 0 : buffer(9 *natom+1:12*natom ,iimage)=reshape(mep_param%bfgs_fprev(:,:,idynimage),(/3*natom/))
765 : 0 : buffer(12*natom+1:12*natom+9,iimage)=reshape(rprimd(:,:,iimage),(/9/))
766 : 0 : buffer(12*natom+10:12*natom+18,iimage)=reshape(gprimd(:,:),(/9/))
767 : 0 : buffer(12*natom+19:12*natom+27,iimage)=reshape(gmet(:,:),(/9/))
768 : : end do
769 : 0 : call gather_array_img(buffer,buffer_all,mpi_enreg,allgather=.true.)
770 [ # # ]: 0 : do idynimage=1,ndynimage_tot
771 : 0 : iimage_tot=list_dynimage_tot(idynimage)
772 [ # # ]: 0 : do iatom=1,natom
773 : 0 : indi=12*(iatom-1)
774 [ # # ]: 0 : do ii=1,3
775 : 0 : xcart_all (ii,iatom,idynimage)= buffer_all(indi +ii,iimage_tot)
776 : 0 : fcart_all (ii,iatom,idynimage)=-buffer_all(indi+3+ii,iimage_tot) ! use fcart=-cartesian_force
777 : 0 : xcartp_all(ii,iatom,idynimage)= buffer_all(indi+6+ii,iimage_tot)
778 : 0 : fcartp_all(ii,iatom,idynimage)= buffer_all(indi+9+ii,iimage_tot)
779 : : end do
780 : : end do
781 : 0 : indi=12*natom
782 [ # # # # ]: 0 : rprimd_all(1:3,1:3,idynimage)=reshape(buffer_all(indi+ 1:indi+ 9,iimage_tot),(/3,3/))
783 [ # # # # ]: 0 : gprimd_all(1:3,1:3,idynimage)=reshape(buffer_all(indi+10:indi+18,iimage_tot),(/3,3/))
784 [ # # # # ]: 0 : gmet_all (1:3,1:3,idynimage)=reshape(buffer_all(indi+19:indi+27,iimage_tot),(/3,3/))
785 : : end do
786 [ # # ]: 0 : ABI_FREE(buffer)
787 [ # # ]: 0 : ABI_FREE(buffer_all)
788 : : else
789 [ # # ]: 0 : do idynimage=1,ndynimage
790 : 0 : iimage=list_dynimage(idynimage)
791 [ # # # # ]: 0 : xcart_all(:,:,idynimage)= xcart(:,:,iimage)
792 [ # # # # ]: 0 : fcart_all(:,:,idynimage)=-fcart(:,:,iimage) ! use fcart=-cartesian_force
793 [ # # # # ]: 0 : xcartp_all(:,:,idynimage)=mep_param%bfgs_xprev(:,:,idynimage)
794 [ # # # # ]: 0 : fcartp_all(:,:,idynimage)=mep_param%bfgs_fprev(:,:,idynimage)
795 [ # # # # ]: 0 : rprimd_all(:,:,idynimage)=rprimd(:,:,iimage)
796 : 0 : call metric(gmet_all(:,:,idynimage),gprimd_all(:,:,idynimage),-1,rmet,rprimd(:,:,iimage),ucvol)
797 : : end do
798 : : end if
799 : :
800 : : !Test if a reset is needed
801 : 0 : reset=.false.
802 [ # # ]: 0 : if (itime>1) then
803 : : dot1=zero;dot2=zero
804 [ # # ]: 0 : do idynimage=1,ndynimage_tot
805 : 0 : dot1=dot2+mep_img_dotp(fcartp_all(:,:,idynimage),fcart_all (:,:,idynimage))
806 : 0 : dot2=dot1+mep_img_dotp(fcartp_all(:,:,idynimage),fcartp_all(:,:,idynimage))
807 : : end do
808 [ # # # # ]: 0 : reset=((dot2<two*abs(dot1)).or.abs(dot2)<tol8)
809 : : if (reset) then
810 : 0 : msg=' Resetting Hessian matrix.'
811 : 0 : call wrtout(std_out,msg,'COLL')
812 : 0 : call wrtout(ab_out ,msg,'COLL')
813 : : end if
814 : : end if
815 : :
816 : : !===> First step or reset: initialize the Hessian matrix
817 [ # # # # ]: 0 : if (itime==1.or.reset) then
818 [ # # # # ]: 0 : mep_param%gbfgs_hess(:,:)=zero
819 [ # # ]: 0 : do idynimage=1,ndynimage_tot
820 : 0 : indi=3*natom*(idynimage-1)
821 [ # # ]: 0 : do iatom=1,natom
822 [ # # ]: 0 : do mu=1,3
823 [ # # ]: 0 : do nu=1,3
824 [ # # ]: 0 : do ii=1,3
825 [ # # ]: 0 : do jj=1,3
826 [ # # # # ]: 0 : if (mep_param%iatfix(ii,iatom)==0.and. &
827 : 0 : & mep_param%iatfix(jj,iatom)==0) then
828 : : mep_param%gbfgs_hess(indi+mu,indi+nu)=mep_param%gbfgs_hess(indi+mu,indi+nu) &
829 : : & +rprimd_all(mu,ii,idynimage)*rprimd_all(nu,jj,idynimage) &
830 : 0 : & *gmet_all(ii,jj,idynimage)*initial_Hessian
831 : : end if
832 : : end do
833 : : end do
834 : : end do
835 : : end do
836 : 0 : indi=indi+3
837 : : end do
838 : : end do
839 : :
840 : : !===> Other steps: update the Hessian matrix
841 : : else
842 : :
843 : : ! Impose here f-fprev=0 (cannot be done inside hessupdt in cartesian coordinates)
844 [ # # # # : 0 : ABI_MALLOC(gred,(3,natom))
# # # # #
# # # ]
845 [ # # ]: 0 : do idynimage=1,ndynimage_tot
846 [ # # # # ]: 0 : fcartp_all(:,:,idynimage)=fcartp_all(:,:,idynimage)-fcart_all(:,:,idynimage)
847 : 0 : call fcart2gred(fcartp_all(:,:,idynimage),gred,rprimd_all(:,:,idynimage),natom)
848 [ # # # # : 0 : where (mep_param%iatfix(:,:)==1) ! iatfix is defined in reduced coordinates
# # ]
849 : : gred(:,:)=zero
850 : : end where
851 : 0 : call gred2fcart(favg,.TRUE.,fcartp_all(:,:,idynimage),gred,gprimd_all(:,:,idynimage),natom)
852 [ # # ]: 0 : do iatom=1,natom
853 : : fcartp_all(:,iatom,idynimage)=fcartp_all(:,iatom,idynimage) &
854 [ # # ]: 0 : & +fcart_all(:,iatom,idynimage)+favg(:)
855 : : end do
856 : : end do
857 [ # # ]: 0 : ABI_FREE(gred)
858 : :
859 : : ! f-fprev=0 has already been imposed for fixed atoms:
860 : : ! we call hessupdt with no fixed atom
861 [ # # # # : 0 : ABI_MALLOC(iatfix_fake,(3,natom))
# # # # #
# # # ]
862 [ # # # # ]: 0 : iatfix_fake(:,:)=0
863 : : call hessupdt(mep_param%gbfgs_hess,&
864 : : & iatfix_fake,natom,3*natom*ndynimage_tot, &
865 : : xcart_all,xcartp_all,fcart_all,fcartp_all, &
866 : 0 : & nimage=ndynimage_tot)
867 [ # # ]: 0 : ABI_FREE(iatfix_fake)
868 : : end if
869 : :
870 : : !Free memory
871 [ # # ]: 0 : ABI_FREE(xcart_all)
872 [ # # ]: 0 : ABI_FREE(xcartp_all)
873 [ # # ]: 0 : ABI_FREE(fcartp_all)
874 [ # # ]: 0 : ABI_FREE(rprimd_all)
875 [ # # ]: 0 : ABI_FREE(gprimd_all)
876 [ # # ]: 0 : ABI_FREE(gmet_all)
877 : :
878 : : !Update history
879 [ # # ]: 0 : do idynimage=1,ndynimage
880 : 0 : iimage=list_dynimage(idynimage)
881 [ # # # # ]: 0 : mep_param%bfgs_xprev(:,:,idynimage)=xcart(:,:,iimage)
882 [ # # # # ]: 0 : mep_param%bfgs_fprev(:,:,idynimage)=fcart(:,:,iimage)
883 : : end do
884 : :
885 : : !Compute image step
886 [ # # # # : 0 : ABI_MALLOC(xstep_all,(3,natom,ndynimage_tot))
# # # # #
# # # ]
887 [ # # # # : 0 : xstep_all=zero
# # ]
888 [ # # ]: 0 : do idynimage=1,ndynimage_tot
889 : 0 : indi=3*natom*(idynimage-1)
890 [ # # ]: 0 : do iatom=1,natom
891 [ # # ]: 0 : do ii=1,3
892 [ # # ]: 0 : do jdynimage=1,ndynimage_tot
893 : 0 : indj=3*natom*(jdynimage-1)
894 [ # # ]: 0 : do jatom=1,natom
895 [ # # ]: 0 : do jj=1,3
896 : : ! Be careful: minus sign because fcart=-cartesian_force
897 : : xstep_all(ii,iatom,idynimage)=xstep_all(ii,iatom,idynimage) &
898 : : & -fcart_all(jj,jatom,jdynimage) &
899 : 0 : & *mep_param%gbfgs_hess(indi+ii,indj+jj)
900 : : end do
901 : 0 : indj=indj+3
902 : : end do
903 : : end do
904 : : end do
905 : 0 : indi=indi+3
906 : : end do
907 : : end do
908 : :
909 : : !Restrict image step size
910 : : stepsize=zero
911 [ # # ]: 0 : do idynimage=1,ndynimage_tot
912 : 0 : stepsize=stepsize+mep_img_dotp(xstep_all(:,:,idynimage),xstep_all(:,:,idynimage))
913 : : end do
914 : 0 : stepsize=sqrt(stepsize)
915 [ # # ]: 0 : if (stepsize>=mep_param%mep_mxstep*dble(ndynimage_tot)) then
916 [ # # # # : 0 : xstep_all=xstep_all*mep_param%mep_mxstep*dble(ndynimage_tot)/stepsize
# # ]
917 : 0 : write(msg,'(a,i3,a)') " Restricting BFGS step size."
918 : 0 : call wrtout(std_out,msg,'COLL')
919 : 0 : call wrtout(ab_out ,msg,'COLL')
920 : : end if
921 : :
922 : : !Update positions
923 [ # # # # : 0 : ABI_MALLOC(xred_old,(3,natom))
# # # # #
# # # ]
924 [ # # ]: 0 : do idynimage=1,ndynimage
925 : 0 : iimage=list_dynimage(idynimage)
926 : 0 : iimage_tot=mpi_enreg%my_imgtab(iimage)
927 [ # # # # ]: 0 : xred_old(:,:)=xred(:,:,iimage)
928 [ # # # # ]: 0 : xcart(:,:,iimage)=xcart(:,:,iimage)+xstep_all(:,:,ind_dynimage_tot(iimage_tot))
929 : 0 : call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage))
930 : : ! In case atom is fixed, we restore its previous value
931 [ # # ]: 0 : do iatom=1,natom
932 [ # # # # : 0 : if (any(mep_param%iatfix(:,iatom)==1)) then
# # ]
933 [ # # # # ]: 0 : where(mep_param%iatfix(:,iatom)==1)
934 : : xred(:,iatom,iimage)=xred_old(:,iatom)
935 : : end where
936 : 0 : call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage))
937 : : end if
938 : : end do
939 : : end do
940 [ # # ]: 0 : ABI_FREE(xred_old)
941 : :
942 : : !Free memory
943 [ # # ]: 0 : ABI_FREE(fcart_all)
944 [ # # ]: 0 : ABI_FREE(xstep_all)
945 [ # # ]: 0 : ABI_FREE(ind_dynimage_tot)
946 [ # # ]: 0 : ABI_FREE(list_dynimage_tot)
947 : :
948 [ # # # # : 0 : end subroutine mep_gbfgs
# # ]
949 : : !!***
950 : :
951 : : !----------------------------------------------------------------------
952 : :
953 : : !!****f* m_mep/mep_rk4
954 : : !! NAME
955 : : !! mep_rk4
956 : : !!
957 : : !! FUNCTION
958 : : !! Make a path (string of images) evolve according to a fourfth-order Runge-Kutta algorithm
959 : : !!
960 : : !! INPUTS
961 : : !! itime=time step
962 : : !! list_dynimage(nimage)=list of dynamical images.
963 : : !! mep_param=datastructure of type mep_type.
964 : : !! several parameters for Minimal Energy Path (MEP) search.
965 : : !! natom=number of atoms
966 : : !! ndynimage=number of dynamical images along the path
967 : : !! nimage=number of images (including static ones)
968 : : !!
969 : : !! OUTPUT
970 : : !!
971 : : !! SIDE EFFECTS
972 : : !! mep_param=datastructure of type mep_type.
973 : : !! History for Runge-Kutta algorithm is filled up
974 : : !! xcart(3,natom,nimage)=cartesian coordinates of atoms in each image along the path
975 : : !! before and after time evolution
976 : : !! after time evolution
977 : : !! xred(3,natom,nimage)=reduced coordinates of atoms in each image along the path
978 : : !! before and after time evolution
979 : : !!
980 : : !! SOURCE
981 : :
982 : 0 : subroutine mep_rk4(fcart,itime,list_dynimage,mep_param,natom,ndynimage,nimage,rprimd,xcart,xred)
983 : :
984 : : !Arguments ------------------------------------
985 : : !scalars
986 : : integer,intent(in) :: itime,natom,ndynimage,nimage
987 : : type(mep_type),intent(inout) :: mep_param
988 : : !arrays
989 : : integer,intent(in) :: list_dynimage(ndynimage)
990 : : real(dp),intent(in) :: rprimd(3,3,nimage)
991 : : real(dp),intent(in) :: fcart(3,natom,nimage)
992 : : real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage)
993 : : !Local variables-------------------------------
994 : : !scalars
995 : : integer,save :: istep_rk=0
996 : : integer :: iatom,idynimage,iimage
997 : : real(dp) :: stepsize
998 : : character(len=500) :: msg
999 : : !arrays
1000 : 0 : real(dp),allocatable :: xred_old(:,:),xstep(:,:)
1001 : :
1002 : : !************************************************************************
1003 : :
1004 : : !Step for RK4 algorithm
1005 : 0 : istep_rk=mod(itime,4)
1006 : :
1007 : : !Store data according to Runge-Kutta algo step
1008 [ # # ]: 0 : if (istep_rk==1) then
1009 [ # # ]: 0 : ABI_SFREE(mep_param%rk4_xcart1)
1010 [ # # ]: 0 : ABI_SFREE(mep_param%rk4_fcart1)
1011 [ # # # # : 0 : ABI_MALLOC(mep_param%rk4_xcart1,(3,natom,nimage))
# # # # #
# # # # #
# # # # ]
1012 [ # # # # : 0 : ABI_MALLOC(mep_param%rk4_fcart1,(3,natom,nimage))
# # # # #
# # # #
# ]
1013 [ # # # # : 0 : mep_param%rk4_xcart1 = xcart
# # # # ]
1014 [ # # # # : 0 : mep_param%rk4_fcart1 = fcart
# # # # ]
1015 [ # # ]: 0 : else if (istep_rk==2) then
1016 [ # # ]: 0 : ABI_SFREE(mep_param%rk4_fcart2)
1017 [ # # # # : 0 : ABI_MALLOC(mep_param%rk4_fcart2,(3,natom,nimage))
# # # # #
# # # # #
# # # # ]
1018 [ # # # # : 0 : mep_param%rk4_fcart2 = fcart
# # # # ]
1019 [ # # ]: 0 : else if (istep_rk==3) then
1020 [ # # ]: 0 : ABI_SFREE(mep_param%rk4_fcart3)
1021 [ # # # # : 0 : ABI_MALLOC(mep_param%rk4_fcart3,(3,natom,nimage))
# # # # #
# # # # #
# # # # ]
1022 [ # # # # : 0 : mep_param%rk4_fcart3 = fcart
# # # # ]
1023 : : end if
1024 : :
1025 [ # # # # : 0 : ABI_MALLOC(xred_old,(3,natom))
# # # # #
# # # ]
1026 [ # # ]: 0 : if (istep_rk==0) then
1027 [ # # # # : 0 : ABI_MALLOC(xstep,(3,natom))
# # # # ]
1028 : : end if
1029 : :
1030 [ # # ]: 0 : do idynimage=1,ndynimage
1031 : 0 : iimage=list_dynimage(idynimage)
1032 [ # # # # ]: 0 : xred_old(:,:)=xred(:,:,iimage)
1033 : :
1034 : : ! Note that one uses fcart, for which the sum of forces on all atoms vanish
1035 : :
1036 : : ! Intermediate Runge-Kutta step 1
1037 [ # # ]: 0 : if (istep_rk==1) then
1038 : : xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage) &
1039 [ # # # # ]: 0 : & -half*mep_param%fxcartfactor*fcart(:,:,iimage)
1040 : :
1041 : : ! Intermediate Runge-Kutta step 2
1042 [ # # ]: 0 : else if (istep_rk==2) then
1043 : : xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage) &
1044 [ # # # # ]: 0 : & -half*mep_param%fxcartfactor*fcart(:,:,iimage)
1045 : :
1046 : : ! Intermediate Runge-Kutta step 3
1047 [ # # ]: 0 : else if (istep_rk==3) then
1048 : : xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage) &
1049 [ # # # # ]: 0 : & -mep_param%fxcartfactor*fcart(:,:,iimage)
1050 : :
1051 : : ! Final Runge-Kutta step
1052 [ # # ]: 0 : else if (istep_rk==0) then
1053 : : ! Compute image step
1054 : : xstep(:,:)=third*mep_param%fxcartfactor &
1055 : : & *(half*fcart(:,:,iimage) &
1056 : : & +half*mep_param%rk4_fcart1(:,:,iimage) &
1057 : : & +mep_param%rk4_fcart2(:,:,iimage) &
1058 [ # # # # ]: 0 : & +mep_param%rk4_fcart3(:,:,iimage))
1059 : 0 : stepsize=mep_img_norm(xstep)
1060 [ # # ]: 0 : if (stepsize>=mep_param%mep_mxstep) then
1061 [ # # # # ]: 0 : xstep=xstep*mep_param%mep_mxstep/stepsize
1062 : 0 : write(msg,'(a,i3,a)') " Restricting step size of image ",iimage,"."
1063 : 0 : call wrtout(std_out,msg,'COLL')
1064 : 0 : call wrtout(ab_out ,msg,'COLL')
1065 : : end if
1066 : : ! Update positions
1067 [ # # # # ]: 0 : xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage)+xstep(:,:)
1068 : : end if
1069 : :
1070 : 0 : call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage))
1071 : :
1072 : : ! In case atom is fixed, we restore its previous value
1073 [ # # ]: 0 : do iatom=1,natom
1074 [ # # # # : 0 : if (any(mep_param%iatfix(:,iatom)==1)) then
# # ]
1075 [ # # # # ]: 0 : where(mep_param%iatfix(:,iatom)==1)
1076 : : xred(:,iatom,iimage)=xred_old(:,iatom)
1077 : : end where
1078 : 0 : call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage))
1079 : : end if
1080 : : end do
1081 : :
1082 : : end do
1083 : :
1084 [ # # ]: 0 : ABI_FREE(xred_old)
1085 [ # # ]: 0 : if (istep_rk==0) then
1086 [ # # ]: 0 : ABI_FREE(xstep)
1087 : : end if
1088 : :
1089 : : !Cancel storage when final RK step has been done
1090 [ # # ]: 0 : if (istep_rk==0) then
1091 [ # # ]: 0 : ABI_SFREE(mep_param%rk4_xcart1)
1092 [ # # ]: 0 : ABI_SFREE(mep_param%rk4_fcart1)
1093 [ # # ]: 0 : ABI_SFREE(mep_param%rk4_fcart2)
1094 [ # # ]: 0 : ABI_SFREE(mep_param%rk4_fcart3)
1095 : : end if
1096 : :
1097 [ # # ]: 0 : end subroutine mep_rk4
1098 : : !!***
1099 : :
1100 : : !----------------------------------------------------------------------
1101 : :
1102 : : !!****f* m_mep/mep_img_dotp
1103 : : !! NAME
1104 : : !! mep_img_dotp
1105 : : !!
1106 : : !! FUNCTION
1107 : : !! Compute the dot product of two vectors in the configuration space:
1108 : : !! Vect1(3,natom).Vect2(3,natom)
1109 : : !!
1110 : : !! INPUTS
1111 : : !! vect1(3,natom)=input vector 1
1112 : : !! vect2(3,natom)=input vector 2
1113 : : !!
1114 : : !! OUTPUT
1115 : : !! mep_img_dotp=dot product
1116 : : !!
1117 : : !! SOURCE
1118 : :
1119 [ - + - + ]: 226 : function mep_img_dotp(vect1,vect2)
1120 : :
1121 : : !Arguments ------------------------------------
1122 : : !scalars
1123 : : real(dp) :: mep_img_dotp
1124 : : !arrays
1125 : : real(dp),intent(in) :: vect1(:,:),vect2(:,:)
1126 : : !Local variables-------------------------------
1127 : : !scalars
1128 : : integer :: size1,size2
1129 : : !arrays
1130 : :
1131 : : !************************************************************************
1132 : :
1133 : 226 : size1=size(vect1,1);size2=size(vect1,2)
1134 [ + - - + ]: 226 : if (size1/=size(vect2,1).or.size2/=size(vect2,2)) then
1135 : 0 : ABI_BUG("Error on dimensions !")
1136 : : end if
1137 : :
1138 [ + + + + ]: 2574 : mep_img_dotp=sum(vect1*vect2)
1139 : :
1140 : 226 : end function mep_img_dotp
1141 : : !!***
1142 : :
1143 : : !----------------------------------------------------------------------
1144 : :
1145 : : !!****f* m_mep/mep_img_norm
1146 : : !! NAME
1147 : : !! mep_img_norm
1148 : : !!
1149 : : !! FUNCTION
1150 : : !! Compute the norm of a vector in the configuration space:
1151 : : !! |Vect(3,natom)|
1152 : : !!
1153 : : !! INPUTS
1154 : : !! vect(3,natom)=input vector
1155 : : !!
1156 : : !! OUTPUT
1157 : : !! mep_img_norm=norm
1158 : : !!
1159 : : !! SOURCE
1160 : :
1161 [ - + ]: 1590 : function mep_img_norm(vect)
1162 : :
1163 : : !Arguments ------------------------------------
1164 : : !scalars
1165 : : real(dp) :: mep_img_norm
1166 : : !arrays
1167 : : real(dp),intent(in) :: vect(:,:)
1168 : :
1169 : : !************************************************************************
1170 : :
1171 [ + + + + ]: 16038 : mep_img_norm=sqrt(sum(vect*vect))
1172 : :
1173 : 1590 : end function mep_img_norm
1174 : : !!***
1175 : :
1176 : : !----------------------------------------------------------------------
1177 : :
1178 : : !!****f* m_mep/mep_img_dotp_red
1179 : : !! NAME
1180 : : !! mep_img_dotp_red
1181 : : !!
1182 : : !! FUNCTION
1183 : : !! Compute the dot product of two vectors in the configuration space:
1184 : : !! Vect1(3,natom).Vect2(3,natom)
1185 : : !! using reduced coordinates
1186 : : !!
1187 : : !! INPUTS
1188 : : !! rmet(3,3)=metric tensor
1189 : : !! vect1(3,natom)=input vector 1
1190 : : !! vect2(3,natom)=input vector 2
1191 : : !!
1192 : : !! OUTPUT
1193 : : !! mep_img_dotp_red=dot product
1194 : : !!
1195 : : !! SOURCE
1196 : :
1197 [ # # # # ]: 0 : function mep_img_dotp_red(rmet,vect1,vect2)
1198 : :
1199 : : !Arguments ------------------------------------
1200 : : !scalars
1201 : : real(dp) :: mep_img_dotp_red
1202 : : !arrays
1203 : : real(dp),intent(in) :: rmet(3,3)
1204 : : real(dp),intent(in) :: vect1(:,:),vect2(:,:)
1205 : : !Local variables-------------------------------
1206 : : !scalars
1207 : : integer :: iatom,size1,size2
1208 : :
1209 : : !************************************************************************
1210 : :
1211 : 0 : size1=size(vect1,1);size2=size(vect1,2)
1212 [ # # # # : 0 : if (size1/=size(vect2,1).or.size2/=size(vect2,2).or.size1/=3) then
# # ]
1213 : 0 : ABI_BUG("Error on dimensions !")
1214 : : end if
1215 : :
1216 : 0 : mep_img_dotp_red=zero
1217 [ # # ]: 0 : do iatom = 1, size2
1218 [ # # # # : 0 : mep_img_dotp_red = mep_img_dotp_red + dot_product(vect1(:, iatom), matmul(rmet, vect2(:, iatom)))
# # ]
1219 : : end do
1220 : :
1221 : 0 : end function mep_img_dotp_red
1222 : : !!***
1223 : :
1224 : : !----------------------------------------------------------------------
1225 : :
1226 : : !!****f* m_mep/mep_img_norm_red
1227 : : !! NAME
1228 : : !! mep_img_norm_red
1229 : : !!
1230 : : !! FUNCTION
1231 : : !! Compute the norm of a vector in the configuration space:
1232 : : !! |Vect(3,natom)|
1233 : : !! using reduced coordinates
1234 : : !!
1235 : : !! INPUTS
1236 : : !! rmet(3,3)=metric tensor
1237 : : !! vect(3,natom)=input vector
1238 : : !!
1239 : : !! OUTPUT
1240 : : !! mep_img_norm_red=norm
1241 : : !!
1242 : : !!
1243 : : !! SOURCE
1244 : :
1245 [ # # ]: 0 : function mep_img_norm_red(rmet,vect)
1246 : :
1247 : : !Arguments ------------------------------------
1248 : : !scalars
1249 : : real(dp) :: mep_img_norm_red
1250 : : !arrays
1251 : : real(dp),intent(in) :: rmet(3,3)
1252 : : real(dp),intent(in) :: vect(:,:)
1253 : : !************************************************************************
1254 : :
1255 : 0 : mep_img_norm_red=sqrt(mep_img_dotp_red(rmet,vect,vect))
1256 : :
1257 : 0 : end function mep_img_norm_red
1258 : : !!***
1259 : :
1260 [ # # # # : 0 : END MODULE m_mep
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
1261 : : !!***
|
