Branch data Line data Source code
1 : : !!****m* ABINIT/m_predict_neb
2 : : !! NAME
3 : : !! m_predict_neb
4 : : !!
5 : : !! FUNCTION
6 : : !! This module implement the Nudged Elastic Band method (several variants)
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2012-2026 ABINIT group (MT,QD)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_predict_neb
23 : :
24 : : use defs_basis
25 : : use m_splines
26 : : use m_mep
27 : : use m_abicore
28 : : use m_errors
29 : : use m_xmpi
30 : :
31 : : use defs_abitypes, only : MPI_type
32 : : use m_results_img, only : results_img_type, gather_array_img, scatter_array_img, &
33 : : & get_geometry_img
34 : : use m_geometry, only : mkradim, metric, fcart2gred
35 : : use m_matrix, only : matr3inv
36 : :
37 : : implicit none
38 : :
39 : : private
40 : : !!***
41 : :
42 : : public :: predict_neb
43 : : !!***
44 : :
45 : : contains
46 : : !!***
47 : :
48 : : !!****f* ABINIT/predict_neb
49 : : !! NAME
50 : : !! predict_neb
51 : : !!
52 : : !! FUNCTION
53 : : !! Given the past history of images, predict the new set of images using Nudged Elastic Band method.
54 : : !! No change of acell, rprim and vel at present.
55 : : !!
56 : : !! INPUTS
57 : : !! itimimage=time index for image propagation (itimimage+1 is to be predicted here)
58 : : !! itimimage_eff=time index in the history
59 : : !! list_dynimage(nimage)=list of dynamical images. The non-dynamical ones will not change.
60 : : !! Example: in the NEB of string method, one expect the two end images to be fixed.
61 : : !! This is quite useful when ground states of the A and B states is known
62 : : !! mpi_enreg=MPI-parallelisation information
63 : : !! natom=dimension of vel_timimage and xred_timimage
64 : : !! ndynimage=number of dynamical images
65 : : !! nimage=number of images (on current proc)
66 : : !! nimage_tot=total number of images
67 : : !! ntimimage_stored=number of time steps stored in the history
68 : : !!
69 : : !! OUTPUT
70 : : !!
71 : : !! SIDE EFFECTS
72 : : !! mep_param=several parameters for Minimal Energy Path (MEP) search
73 : : !! results_img(ntimimage_stored,nimage)=datastructure that holds the history of previous computations.
74 : : !! results_img(:,:)%acell(3)
75 : : !! at input, history of the values of acell for all images
76 : : !! at output, the predicted values of acell for all images
77 : : !! results_img(:,:)%results_gs
78 : : !! at input, history of the values of energies and forces for all images
79 : : !! results_img(:,:)%rprim(3,3)
80 : : !! at input, history of the values of rprim for all images
81 : : !! at output, the predicted values of rprim for all images
82 : : !! results_img(:,:)%vel(3,natom)
83 : : !! at input, history of the values of vel for all images
84 : : !! at output, the predicted values of vel for all images
85 : : !! results_img(:,:)%vel_cell(3,3)
86 : : !! at input, history of the values of vel_cell for all images
87 : : !! at output, the predicted values of vel_cell for all images
88 : : !! results_img(:,:)%xred(3,natom)
89 : : !! at input, history of the values of xred for all images
90 : : !! at output, the predicted values of xred for all images
91 : : !!
92 : : !! SOURCE
93 : :
94 : 42 : subroutine predict_neb(itimimage,itimimage_eff,list_dynimage,mep_param,mpi_enreg,natom,&
95 : 42 : & ndynimage,nimage,nimage_tot,ntimimage_stored,results_img)
96 : :
97 : : !Arguments ------------------------------------
98 : : !scalars
99 : : integer,intent(in) :: itimimage,itimimage_eff,natom,ndynimage
100 : : integer,intent(in) :: nimage,nimage_tot,ntimimage_stored
101 : : type(mep_type),intent(inout) :: mep_param
102 : : type(MPI_type),intent(in) :: mpi_enreg
103 : : !arrays
104 : : integer,intent(in) :: list_dynimage(ndynimage)
105 : : type(results_img_type),intent(inout) :: results_img(nimage,ntimimage_stored)
106 : :
107 : : !Local variables-------------------------------
108 : : !scalars
109 : : integer :: ierr,ii,jj,iimage,iimage_min,iimage_max,natom_eff,next_itimimage
110 : : logical :: test_minus_one,test_plus_one,use_reduced_coord
111 : : real(dp) :: dvmax,dvmin,ecur,emax,emin,eref,f_para1,f_para2
112 : : character(len=100) :: msg
113 : : !arrays
114 : : integer,parameter :: voigt(3,3)=reshape([1,6,3,6,2,4,3,4,3],[3,3])
115 : : real(dp),parameter :: identity_real(3,3)=reshape([one,zero,zero,zero,one,zero,zero,zero,one],[3,3])
116 : : real(dp) :: mat3_1(3,3),mat3_2(3,3)
117 : 42 : real(dp),allocatable :: acell(:,:),buffer(:,:),buffer_all(:,:)
118 : 42 : real(dp),allocatable :: coordif(:,:,:),dimage(:),spring(:)
119 : 42 : real(dp),allocatable :: rprim(:,:,:),rprimd_start(:,:,:),rprimd_start_inv(:,:,:)
120 : 42 : real(dp),allocatable :: fcart(:,:,:),xcart(:,:,:),xred(:,:,:),strainfact_jj(:)
121 : 42 : real(dp),allocatable :: strten(:,:),strten_mat(:,:,:),rmet(:,:,:),ucvol(:),pressure(:)
122 : 42 : real(dp),allocatable :: tangent(:,:,:),vect(:,:)
123 : 42 : real(dp),allocatable,target :: etotal(:),neb_forces(:,:,:),rprimd(:,:,:)
124 : 42 : real(dp),allocatable,target :: fcart_eff(:,:,:),xcart_eff(:,:,:),xred_eff(:,:,:)
125 : 42 : real(dp),pointer :: coord_(:,:,:),etotal_all(:),fcart_eff_all(:,:,:), neb_forces_all(:,:,:)
126 : 42 : real(dp),pointer :: rprimd_all(:,:,:),xcart_eff_all(:,:,:),xred_eff_all(:,:,:)
127 : :
128 : : ! *************************************************************************
129 : :
130 : : !Check options
131 [ - + ]: 9 : if (mep_param%neb_cell_algo/=NEB_CELL_ALGO_NONE.and.&
132 : : & mep_param%mep_solver/=MEP_SOLVER_STEEPEST) then
133 : 0 : msg='Variable cell NEB only allowed with steepest descent algo!'
134 : 0 : ABI_ERROR(msg)
135 : : end if
136 : :
137 : : !In case of variable-cell, 3 additional "fake" atoms are added to unit cell vectors
138 [ + + ]: 42 : natom_eff=natom ; if (mep_param%neb_cell_algo/=NEB_CELL_ALGO_NONE) natom_eff=natom_eff+3
139 : :
140 : : !VC-NEB uses reduced coordinates
141 : 42 : use_reduced_coord=.false.
142 [ - + ]: 42 : if (mep_param%neb_cell_algo==NEB_CELL_ALGO_VCNEB) use_reduced_coord=.true.
143 : :
144 [ - + + - : 126 : ABI_MALLOC(acell,(3,nimage))
+ - + - -
+ - + ]
145 [ - + + - : 126 : ABI_MALLOC(rprim,(3,3,nimage))
+ - + - -
+ - + ]
146 [ - + + - : 168 : ABI_MALLOC(xred_eff,(3,natom_eff,nimage))
- + + - +
- + - - +
- + ]
147 : :
148 : : !Parallelism over images: only one process per image of the cell
149 [ + - ]: 42 : if (mpi_enreg%me_cell==0) then
150 : :
151 : : ! Retrieve positions and forces
152 [ - + + - : 126 : ABI_MALLOC(etotal,(nimage))
+ - + - -
+ - + ]
153 [ - + + - : 126 : ABI_MALLOC(xcart_eff,(3,natom_eff,nimage))
- + + - -
+ - + ]
154 [ - + + - : 126 : ABI_MALLOC(fcart_eff,(3,natom_eff,nimage))
- + + - -
+ - + ]
155 [ - + + - : 84 : ABI_MALLOC(rprimd,(3,3,nimage))
- + - + ]
156 [ - + + - : 84 : ABI_MALLOC(strainfact_jj,(nimage))
- + - + ]
157 [ + + ]: 316 : strainfact_jj(:)=one
158 [ - + + - : 84 : ABI_MALLOC(rprimd_start,(3,3,nimage))
- + - + ]
159 [ + + ]: 316 : do iimage=1,nimage
160 : 274 : ii=mpi_enreg%my_imgtab(iimage)
161 [ + + + + ]: 3604 : rprimd_start(:,:,iimage)=mep_param%rprimd_start(:,:,ii)
162 : : end do
163 [ - + + - : 168 : ABI_MALLOC(xred,(3,natom,nimage))
- + + - +
- + - - +
- + ]
164 [ - + + - : 126 : ABI_MALLOC(xcart,(3,natom,nimage))
- + + - -
+ - + ]
165 [ - + + - : 126 : ABI_MALLOC(fcart,(3,natom,nimage))
- + + - -
+ - + ]
166 [ - + + - : 126 : ABI_MALLOC(strten,(6,nimage))
+ - + - -
+ - + ]
167 : 42 : call get_geometry_img(results_img(:,itimimage_eff),etotal,natom,nimage,fcart,rprimd,strten,xcart,xred)
168 [ + + + + : 2508 : xred_eff(1:3,1:natom,1:nimage)=xred(1:3,1:natom,1:nimage)
+ + ]
169 [ + + + + : 2508 : xcart_eff(1:3,1:natom,1:nimage)=xcart(1:3,1:natom,1:nimage)
+ + ]
170 [ + + + + : 2508 : fcart_eff(1:3,1:natom,1:nimage)=fcart(1:3,1:natom,1:nimage)
+ + ]
171 : :
172 : : ! Retrieve unit cell vectors and derivatives (forces) in case of variable-cell NEB
173 : : ! => Retrieve pressure, volume, rmet, strten, ...
174 [ + + ]: 42 : if (mep_param%neb_cell_algo/=NEB_CELL_ALGO_NONE) then
175 [ - + + - : 18 : ABI_MALLOC(rmet,(3,3,nimage))
- + - + ]
176 [ - + + - : 18 : ABI_MALLOC(ucvol,(nimage))
- + - + ]
177 [ - + + - : 18 : ABI_MALLOC(pressure,(nimage))
- + - + ]
178 [ - + + - : 18 : ABI_MALLOC(strten_mat,(3,3,nimage))
- + - + ]
179 [ - + + - : 18 : ABI_MALLOC(rprimd_start_inv,(3,3,nimage))
- + - + ]
180 [ + + ]: 72 : pressure(1:nimage)=-(strten(1,1:nimage)+strten(2,1:nimage)+strten(3,1:nimage))*third
181 [ + + ]: 72 : do iimage=1,nimage
182 [ + + + + ]: 819 : xcart_eff(1:3,natom+1:natom+3,iimage)=zero
183 [ + + + + ]: 819 : xred_eff(1:3,natom+1:natom+3,iimage)=zero
184 : 63 : call metric(mat3_1,mat3_2,-1,rmet(:,:,iimage),rprimd(:,:,iimage),ucvol(iimage))
185 [ + + + + ]: 819 : do jj=1,3; do ii=1,3
186 : 756 : strten_mat(ii,jj,iimage)=strten(voigt(ii,jj),iimage)
187 : : end do; end do
188 : 63 : call matr3inv(rprimd_start(:,:,iimage),mat3_1(:,:))
189 [ + + + + ]: 828 : rprimd_start_inv(:,:,iimage) = transpose(mat3_1(:,:))
190 : : end do
191 : :
192 : : ! ==== Generalized Solid-State NEB (GSS-NEB)
193 : : ! See: Sheppard, Xiao, Chemelewski, Johnson, Henkelman, J. Chem. Phys. 136, 074103 (2012)
194 [ + - ]: 9 : if (mep_param%neb_cell_algo==NEB_CELL_ALGO_GSSNEB) then
195 [ + + ]: 72 : strainfact_jj(1:nimage)=(ucvol(1:nimage)**third)*(natom**sixth)
196 [ + + ]: 72 : do iimage=1,nimage
197 [ + + + + ]: 819 : mat3_1(1:3,1:3)=rprimd(1:3,1:3,iimage)-rprimd_start(1:3,1:3,iimage)
198 [ + + + + : 4032 : xcart_eff(1:3,natom+1:natom+3,iimage)=strainfact_jj(iimage)*matmul(rprimd_start_inv(:,:,iimage),mat3_1)
+ + + + +
+ + + +
+ ]
199 [ + + + + ]: 819 : strten_mat(1:3,1:3,iimage)=strten_mat(1:3,1:3,iimage)+pressure(iimage)*identity_real(1:3,1:3)
200 : : fcart_eff(1:3,natom+1:natom+3,iimage)= &
201 [ + + + + ]: 828 : & -(ucvol(iimage)/strainfact_jj(iimage))*strten_mat(1:3,1:3,iimage)
202 : : end do
203 : : end if
204 : :
205 : : ! ==== Variable-cell NEB (VC-NEB, in reduced coordinates) - At present, doesnt work
206 : : ! See: Qian, Dong, Zhou, Tian, Oganov, Wang, Comp. Phys. Comm. 184, 2111 (2013)
207 [ - + ]: 9 : if (mep_param%neb_cell_algo==NEB_CELL_ALGO_VCNEB) then
208 [ # # ]: 0 : do iimage=1,nimage
209 : : xred_eff(1:3,natom+1:natom+3,iimage)= &
210 [ # # # # : 0 : & matmul(rprimd(1:3,1:3,iimage),rprimd_start_inv(1:3,1:3,iimage))-identity_real(1:3,1:3)
# # # # #
# # # #
# ]
211 : 0 : call fcart2gred(fcart(1:3,1:natom,iimage),fcart_eff(1:3,1:natom,iimage),rprimd(1:3,1:3,iimage),natom)
212 [ # # # # : 0 : fcart_eff(1:3,1:natom,iimage)=matmul(rmet(1:3,1:3,iimage),fcart_eff(1:3,1:natom,iimage))
# # # # #
# # # # #
# # # # #
# # # ]
213 [ # # # # ]: 0 : mat3_1(1:3,1:3)=identity_real(1:3,1:3)+transpose(xred_eff(1:3,natom+1:natom+3,iimage))
214 : 0 : call matr3inv(mat3_1,mat3_2)
215 [ # # # # ]: 0 : mat3_1=transpose(mat3_2)
216 [ # # # # ]: 0 : mat3_2(1:3,1:3)=(strten_mat(1:3,1:3,iimage)+pressure(iimage))*identity_real(1:3,1:3)*ucvol(iimage)
217 [ # # # # : 0 : fcart_eff(1:3,natom+1:natom+3,iimage)=-matmul(mat3_2,mat3_1)
# # # # #
# ]
218 : : end do
219 : : end if
220 : :
221 [ - + ]: 9 : ABI_FREE(rmet)
222 [ - + ]: 9 : ABI_FREE(ucvol)
223 [ - + ]: 9 : ABI_FREE(pressure)
224 [ - + ]: 9 : ABI_FREE(strten_mat)
225 [ - + ]: 9 : ABI_FREE(rprimd_start_inv)
226 : : end if
227 : :
228 [ - + ]: 42 : ABI_FREE(xred)
229 [ - + ]: 42 : ABI_FREE(xcart)
230 [ - + ]: 42 : ABI_FREE(fcart)
231 [ - + ]: 42 : ABI_FREE(strten)
232 : :
233 : : ! Array containing effective NEB forces (F_ortho+F_spring)
234 [ - + + - : 126 : ABI_MALLOC(neb_forces,(3,natom_eff,nimage))
- + + - -
+ - + ]
235 : :
236 : : ! Parallelism: gather data of all images
237 [ - + ]: 42 : if (mpi_enreg%paral_img==1) then
238 [ # # # # : 0 : ABI_MALLOC(buffer,(9*natom_eff+10,nimage))
# # # # #
# # # # #
# # # # ]
239 [ # # # # : 0 : ABI_MALLOC(buffer_all,(9*natom_eff+10,nimage_tot))
# # # # #
# # # # #
# # ]
240 [ # # # # : 0 : ABI_MALLOC(etotal_all,(nimage_tot))
# # # # #
# # # ]
241 [ # # # # : 0 : ABI_MALLOC(xcart_eff_all,(3,natom_eff,nimage_tot))
# # # # #
# # # # #
# # ]
242 [ # # # # : 0 : ABI_MALLOC(fcart_eff_all,(3,natom_eff,nimage_tot))
# # # # #
# # # ]
243 [ # # # # : 0 : ABI_MALLOC(xred_eff_all,(3,natom_eff,nimage_tot))
# # # # #
# # # ]
244 [ # # # # : 0 : ABI_MALLOC(rprimd_all,(3,3,nimage_tot))
# # # # #
# # # ]
245 [ # # # # : 0 : ABI_MALLOC(neb_forces_all,(3,natom_eff,nimage_tot))
# # # # #
# # # ]
246 [ # # # # ]: 0 : buffer=zero;ii=0
247 [ # # ]: 0 : buffer(ii+1 ,1:nimage)=etotal(1:nimage);ii=ii+1
248 [ # # # # : 0 : buffer(ii+1:ii+9 ,1:nimage)=reshape(rprimd,(/9,nimage/));ii=ii+9
# # ]
249 [ # # # # : 0 : buffer(ii+1:ii+3*natom_eff,1:nimage)=reshape(xcart_eff,(/3*natom_eff,nimage/));ii=ii+3*natom_eff
# # ]
250 [ # # # # : 0 : buffer(ii+1:ii+3*natom_eff,1:nimage)=reshape(xred_eff,(/3*natom_eff,nimage/));ii=ii+3*natom_eff
# # ]
251 [ # # # # : 0 : buffer(ii+1:ii+3*natom_eff,1:nimage)=reshape(fcart_eff,(/3*natom_eff,nimage/));ii=ii+3*natom_eff
# # ]
252 : 0 : call gather_array_img(buffer,buffer_all,mpi_enreg,allgather=.true.)
253 : 0 : ii=0
254 [ # # ]: 0 : etotal_all(:)=buffer_all(ii+1,1:nimage_tot);ii=ii+1
255 [ # # ]: 0 : rprimd_all(:,:,:)=reshape(buffer_all(ii+1:ii+9,1:nimage_tot),(/3,3,nimage_tot/));ii=ii+9
256 [ # # ]: 0 : xcart_eff_all(:,:,:)=reshape(buffer_all(ii+1:ii+3*natom_eff,1:nimage_tot),(/3,natom_eff,nimage_tot/));ii=ii+3*natom_eff
257 [ # # ]: 0 : xred_eff_all(:,:,:)=reshape(buffer_all(ii+1:ii+3*natom_eff,1:nimage_tot),(/3,natom_eff,nimage_tot/));ii=ii+3*natom_eff
258 [ # # ]: 0 : fcart_eff_all(:,:,:)=reshape(buffer_all(ii+1:ii+3*natom_eff,1:nimage_tot),(/3,natom_eff,nimage_tot/));ii=ii+3*natom_eff
259 [ # # ]: 0 : ABI_FREE(buffer)
260 [ # # ]: 0 : ABI_FREE(buffer_all)
261 : : else
262 : 42 : etotal_all => etotal
263 : 42 : rprimd_all => rprimd
264 : 42 : xcart_eff_all => xcart_eff
265 : 42 : xred_eff_all => xred_eff
266 : 42 : fcart_eff_all => fcart_eff
267 : 42 : neb_forces_all => neb_forces
268 : : end if
269 : :
270 : : ! coordif is the vector between two images
271 : : ! dimage is the distance between two images
272 : : ! tangent is the tangent at image i
273 [ - + + - : 168 : ABI_MALLOC(coordif,(3,natom_eff,nimage_tot))
- + + - +
- + - - +
- + ]
274 [ - + + - : 126 : ABI_MALLOC(tangent,(3,natom_eff,nimage_tot))
- + + - -
+ - + ]
275 [ - + + - : 126 : ABI_MALLOC(dimage,(nimage_tot))
+ - + - -
+ - + ]
276 : :
277 : : ! Compute distances between images
278 [ + + + + ]: 486 : coordif(:,:,1)=zero;dimage(1)=zero
279 [ - + ]: 42 : coord_ => xcart_eff_all ; if (use_reduced_coord) coord_ => xred_eff_all
280 [ + + ]: 274 : do iimage=2,nimage_tot
281 [ + + + + ]: 2778 : coordif(:,:,iimage)=coord_(:,:,iimage)-coord_(:,:,iimage-1)
282 : : end do
283 : :
284 [ + + ]: 274 : do iimage=2,nimage_tot
285 : 274 : dimage(iimage)=mep_img_norm(coordif(:,:,iimage))
286 : : end do
287 : :
288 : : ! Compute tangent (not normalized)
289 [ + + + + ]: 486 : tangent(:,:,1)=zero
290 [ + + + + ]: 486 : tangent(:,:,nimage_tot)=zero
291 : : ! === Original definition of tangent
292 [ + + ]: 42 : if (mep_param%neb_algo==NEB_ALGO_STANDARD) then
293 [ + + ]: 45 : do iimage=2,nimage_tot-1
294 : : tangent(:,:,iimage)=coordif(:,:,iimage )/dimage(iimage) &
295 [ + + + + ]: 333 : & +coordif(:,:,iimage+1)/dimage(iimage+1)
296 : : end do
297 : : ! === Improved tangent (J. Chem. Phys. 113, 9978 (2000) [[cite:Henkelman2000]])
298 : : else
299 [ + + ]: 187 : do iimage=2,nimage_tot-1
300 : 154 : test_minus_one=(etotal_all(iimage-1)<etotal_all(iimage))
301 : 154 : test_plus_one =(etotal_all(iimage )<etotal_all(iimage+1))
302 [ + + ]: 187 : if (test_minus_one.and.test_plus_one)then ! V_i-1 < V_i < V_i+1
303 [ + + + + ]: 729 : tangent(:,:,iimage)=coordif(:,:,iimage+1)
304 [ + + ]: 93 : else if ((.not.test_minus_one).and.(.not.test_plus_one))then ! V_i-1 > V_i > V_i+1
305 [ + + + + ]: 666 : tangent(:,:,iimage)=coordif(:,:,iimage)
306 : : else ! V_i-1 < V_i > V_i+1 OR V_i-1 > V_i < V_i+1
307 : : dvmax=max(abs(etotal_all(iimage+1)-etotal_all(iimage)),&
308 : 39 : & abs(etotal_all(iimage-1)-etotal_all(iimage)))
309 : : dvmin=min(abs(etotal_all(iimage+1)-etotal_all(iimage)),&
310 : 39 : & abs(etotal_all(iimage-1)-etotal_all(iimage)))
311 [ + + ]: 39 : if (etotal_all(iimage+1)>etotal_all(iimage-1)) then ! V_i+1 > V_i-1
312 : : tangent(:,:,iimage)=coordif(:,:,iimage+1)*dvmax &
313 [ + + + + ]: 189 : & +coordif(:,:,iimage )*dvmin
314 : : else ! V_i+1 < V_i-1
315 : : tangent(:,:,iimage)=coordif(:,:,iimage+1)*dvmin &
316 [ + + + + ]: 342 : & +coordif(:,:,iimage )*dvmax
317 : : end if
318 : : end if
319 : : end do
320 : : end if
321 : :
322 : : ! Normalize tangent
323 [ + + ]: 232 : do iimage=2,nimage_tot-1
324 [ + + + + ]: 2292 : tangent(:,:,iimage)=tangent(:,:,iimage)/mep_img_norm(tangent(:,:,iimage))
325 : : end do
326 : :
327 : : ! Compute spring constant(s)
328 [ - + + - : 84 : ABI_MALLOC(spring,(nimage_tot))
- + - + ]
329 [ + - ]: 42 : if (abs(mep_param%neb_spring(2)-mep_param%neb_spring(1))<tol8) then ! Constant spring
330 [ + + ]: 316 : spring(:)=mep_param%neb_spring(1)
331 : : else ! Variable spring
332 [ # # # # : 0 : emax=maxval(etotal_all(:))
# # # # ]
333 : 0 : eref=max(etotal_all(1),etotal_all(nimage_tot))
334 [ # # ]: 0 : spring(:)=mep_param%neb_spring(1)
335 [ # # ]: 0 : do iimage=2,nimage_tot-1
336 : 0 : ecur=max(etotal_all(iimage-1),etotal_all(iimage))
337 [ # # ]: 0 : if (ecur<eref) then
338 : 0 : spring(iimage)=mep_param%neb_spring(1)
339 : : else
340 : : spring(iimage)=mep_param%neb_spring(2) &
341 : : & -(mep_param%neb_spring(2)-mep_param%neb_spring(1)) &
342 : 0 : & *(emax-ecur)/(emax-eref)
343 : : end if
344 : : end do
345 : : end if
346 : :
347 : : ! CI-NEB: determine image(s) with maximal energy
348 : 42 : iimage_min=-1;iimage_max=-1
349 [ + + + + ]: 42 : if (mep_param%neb_algo==NEB_ALGO_CINEB.and.itimimage>=mep_param%cineb_start) then
350 : 7 : emin=min(etotal_all(1),etotal_all(nimage_tot))
351 : 7 : emax=max(etotal_all(1),etotal_all(nimage_tot))
352 [ + + ]: 42 : do iimage=2,nimage_tot-1
353 : 35 : if (etotal_all(iimage)<emin) then
354 : : iimage_min=iimage;emin=etotal_all(iimage)
355 : : end if
356 [ + + ]: 42 : if (etotal_all(iimage)>emax) then
357 : 21 : iimage_max=iimage;emax=etotal_all(iimage)
358 : : end if
359 : : end do
360 : : end if
361 : :
362 : : ! Compute NEB forces
363 [ + - + + : 930 : neb_forces_all(:,:,1)=fcart_eff_all(:,:,1)
+ + + + +
+ ]
364 [ + + + + : 930 : neb_forces_all(:,:,nimage_tot)=fcart_eff_all(:,:,nimage_tot)
+ + + + ]
365 [ + + ]: 232 : do iimage=2,nimage_tot-1
366 : : ! === Standard NEB
367 [ + + ]: 232 : if (iimage/=iimage_max) then
368 : : ! if (iimage/=iimage_min.and.iimage/=iimage_max) then
369 : 183 : f_para1=mep_img_dotp(fcart_eff_all(:,:,iimage),tangent(:,:,iimage))
370 [ + + ]: 183 : if (mep_param%neb_algo==NEB_ALGO_STANDARD) then ! Original NEB algo
371 [ - + + - : 108 : ABI_MALLOC(vect,(3,natom_eff))
+ - + - -
+ - + ]
372 [ + + + + ]: 324 : vect(:,:)=spring(iimage+1)*coordif(:,:,iimage+1)-spring(iimage)*coordif(:,:,iimage)
373 : 36 : f_para2=mep_img_dotp(vect(:,:),tangent(:,:,iimage))
374 [ - + ]: 36 : ABI_FREE(vect)
375 : : else ! Modification from J. Chem. Phys. 113, 9978 (2000) [[cite:Henkelman2000]]
376 : 147 : f_para2=spring(iimage+1)*dimage(iimage+1)-spring(iimage)*dimage(iimage)
377 : : end if
378 : : neb_forces_all(:,:,iimage)=fcart_eff_all(:,:,iimage) & ! F_ortho
379 : : & -f_para1*tangent(:,:,iimage) & ! F_ortho
380 [ + + + + : 4191 : & +f_para2*tangent(:,:,iimage) ! F_spring
+ + + + ]
381 : : ! === CI-NEB for the image(s) with maximal energy
382 : : else
383 : 7 : f_para1=mep_img_dotp(fcart_eff_all(:,:,iimage),tangent(:,:,iimage))
384 [ + + + + : 119 : neb_forces_all(:,:,iimage)=fcart_eff_all(:,:,iimage)-two*f_para1*tangent(:,:,iimage)
+ + + + ]
385 : : end if
386 : : end do
387 : :
388 : : ! Parallelism: distribute NEB forces
389 [ - + ]: 42 : if (mpi_enreg%paral_img==1) then
390 [ # # ]: 0 : do iimage=1,nimage
391 : 0 : ii=mpi_enreg%my_imgtab(iimage)
392 [ # # # # : 0 : neb_forces(:,:,iimage)=neb_forces_all(:,:,ii)
# # # # ]
393 : : end do
394 : : end if
395 : :
396 : : ! Compute new atomic positions in each cell
397 [ + - ]: 42 : if (mep_param%mep_solver==MEP_SOLVER_STEEPEST) then ! Steepest-descent
398 : : call mep_steepest(neb_forces,list_dynimage,mep_param,natom,natom_eff,&
399 : : & ndynimage,nimage,rprimd,xcart_eff,xred_eff,&
400 : : & strainfact=strainfact_jj,rprimd_start=rprimd_start,&
401 : 42 : & use_reduced_coord=use_reduced_coord)
402 [ # # ]: 0 : else if (mep_param%mep_solver==MEP_SOLVER_QUICKMIN) then ! Quick-min
403 : : call mep_qmin(neb_forces,itimimage,list_dynimage,mep_param,natom,ndynimage, &
404 : 0 : & nimage,rprimd,xcart_eff,xred_eff)
405 [ # # ]: 0 : else if (mep_param%mep_solver==MEP_SOLVER_LBFGS) then ! Local BFGS
406 : : call mep_lbfgs(neb_forces,itimimage,list_dynimage,mep_param,natom,ndynimage,&
407 : 0 : & nimage,rprimd,xcart_eff,xred_eff)
408 [ # # ]: 0 : else if (mep_param%mep_solver==MEP_SOLVER_GBFGS) then ! Global BFGS
409 : : call mep_gbfgs(neb_forces,itimimage,list_dynimage,mep_param,mpi_enreg,natom,ndynimage,&
410 : 0 : & nimage,nimage_tot,rprimd,xcart_eff,xred_eff)
411 : : else
412 : 0 : ABI_BUG("Inconsistent solver !")
413 : : end if
414 : :
415 : : ! Compute new acell and rprim from new rprimd
416 : : !if (mep_param%neb_cell_algo/=NEB_CELL_ALGO_NONE) then
417 [ + + ]: 316 : do iimage=1,nimage
418 : 316 : call mkradim(acell(:,iimage),rprim(:,:,iimage),rprimd(:,:,iimage))
419 : : end do
420 : : !end if
421 : :
422 : : ! Free memory
423 [ - + ]: 42 : ABI_FREE(spring)
424 [ - + ]: 42 : ABI_FREE(coordif)
425 [ - + ]: 42 : ABI_FREE(tangent)
426 : 42 : ABI_FREE(dimage)
427 [ - + ]: 42 : if (mpi_enreg%paral_img==1) then
428 [ # # ]: 0 : ABI_FREE(etotal_all)
429 [ # # ]: 0 : ABI_FREE(rprimd_all)
430 [ # # ]: 0 : ABI_FREE(xred_eff_all)
431 [ # # ]: 0 : ABI_FREE(xcart_eff_all)
432 [ # # ]: 0 : ABI_FREE(fcart_eff_all)
433 [ # # ]: 0 : ABI_FREE(neb_forces_all)
434 : : end if
435 : :
436 [ - + ]: 42 : ABI_FREE(neb_forces)
437 [ - + ]: 42 : ABI_FREE(etotal)
438 [ - + ]: 42 : ABI_FREE(rprimd)
439 [ - + ]: 42 : ABI_FREE(xcart_eff)
440 [ - + ]: 42 : ABI_FREE(fcart_eff)
441 [ - + ]: 42 : ABI_FREE(strainfact_jj)
442 [ - + ]: 42 : ABI_FREE(rprimd_start)
443 : :
444 : : end if ! mpi_enreg%me_cell==0
445 : :
446 : : !Store acell, rprim, xred and vel for the new iteration
447 : 42 : call xmpi_bcast(xred_eff,0,mpi_enreg%comm_cell,ierr)
448 : 42 : call xmpi_bcast(acell,0,mpi_enreg%comm_cell,ierr)
449 : 42 : call xmpi_bcast(rprim,0,mpi_enreg%comm_cell,ierr)
450 : 42 : next_itimimage=itimimage_eff+1
451 [ - + ]: 42 : if (next_itimimage>ntimimage_stored) next_itimimage=1
452 [ + + ]: 316 : do iimage=1,nimage
453 [ + + + + ]: 2466 : results_img(iimage,next_itimimage)%xred(:,1:natom)=xred_eff(:,1:natom,iimage)
454 [ + + ]: 274 : if (mep_param%neb_cell_algo==NEB_CELL_ALGO_NONE) then
455 [ + + ]: 844 : results_img(iimage,next_itimimage)%acell(:) =results_img(iimage,itimimage_eff)%acell(:)
456 [ + + + + ]: 2743 : results_img(iimage,next_itimimage)%rprim(:,:) =results_img(iimage,itimimage_eff)%rprim(:,:)
457 : : else
458 [ + + ]: 252 : results_img(iimage,next_itimimage)%acell(:) =acell(:,iimage)
459 [ + + + + ]: 819 : results_img(iimage,next_itimimage)%rprim(:,:) =rprim(:,:,iimage)
460 : : end if
461 [ + + + + ]: 2466 : results_img(iimage,next_itimimage)%vel(:,:) =results_img(iimage,itimimage_eff)%vel(:,:)
462 [ + + + + ]: 3604 : results_img(iimage,next_itimimage)%vel_cell(:,:)=results_img(iimage,itimimage_eff)%vel_cell(:,:)
463 : : end do
464 : :
465 [ - + ]: 42 : ABI_FREE(xred_eff)
466 [ - + ]: 42 : ABI_FREE(acell)
467 [ - + ]: 42 : ABI_FREE(rprim)
468 : :
469 [ + + - + : 84 : end subroutine predict_neb
- + - + -
+ - + ]
470 : : !!***
471 : :
472 : : end module m_predict_neb
473 : : !!***
|
