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 3975 : 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 3975 : 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 3946 : mep_param%cineb_start = -1
157 3946 : mep_param%mep_solver = -1
158 3946 : mep_param%neb_algo = -1
159 3946 : mep_param%neb_cell_algo = -1
160 3946 : mep_param%string_algo = -1
161 3946 : mep_param%fxcartfactor = zero
162 3946 : mep_param%mep_mxstep = 100._dp
163 11838 : mep_param%neb_spring = zero
164 3946 : nullify(mep_param%iatfix)
165 3946 : nullify(mep_param%rprimd_start)
166 : end if
167 :
168 3975 : 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 3975 : subroutine mep_destroy(mep_param)
191 :
192 : !Arguments ------------------------------------
193 : !scalars
194 : type(mep_type),intent(inout) :: mep_param
195 :
196 : !************************************************************************
197 :
198 3975 : ABI_SFREE(mep_param%bfgs_xprev)
199 3975 : ABI_SFREE(mep_param%gbfgs_hess)
200 3975 : ABI_SFREE(mep_param%bfgs_fprev)
201 3975 : ABI_SFREE(mep_param%lbfgs_hess)
202 3975 : ABI_SFREE(mep_param%qmin_vel)
203 3975 : ABI_SFREE(mep_param%rk4_xcart1)
204 3975 : ABI_SFREE(mep_param%rk4_fcart1)
205 3975 : ABI_SFREE(mep_param%rk4_fcart2)
206 3975 : ABI_SFREE(mep_param%rk4_fcart3)
207 :
208 3975 : nullify(mep_param%iatfix)
209 3975 : nullify(mep_param%rprimd_start)
210 :
211 3975 : 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 : !!***
|
