TABLE OF CONTENTS


ABINIT/m_predict_neb [ Modules ]

[ Top ] [ Modules ]

NAME

  m_predict_neb

FUNCTION

COPYRIGHT

  Copyright (C) 2012-2018 ABINIT group (MT)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_predict_neb
28 
29  use defs_basis
30  use defs_abitypes
31  use m_splines
32  use m_mep
33  use m_abicore
34  use m_errors
35  use m_xmpi
36 
37  use m_results_img, only : results_img_type, gather_array_img, scatter_array_img, get_geometry_img
38 
39  implicit none
40 
41  private

ABINIT/predict_neb [ Functions ]

[ Top ] [ Functions ]

NAME

 predict_neb

FUNCTION

 Given the past history of images, predict the new set of images using Nudged Elastic Band method.
 No change of acell, rprim and vel at present.

INPUTS

 itimimage=time index for image propagation (itimimage+1 is to be predicted here)
 itimimage_eff=time index in the history
 list_dynimage(nimage)=list of dynamical images. The non-dynamical ones will not change.
       Example: in the NEB of string method, one expect the two end images to be fixed.
 This is quite useful when ground states of the A and B states is known
 mpi_enreg=MPI-parallelisation information
 natom=dimension of vel_timimage and xred_timimage
 ndynimage=number of dynamical images
 nimage=number of images (on current proc)
 nimage_tot=total number of images
 ntimimage_stored=number of time steps stored in the history

OUTPUT

SIDE EFFECTS

 mep_param=several parameters for Minimal Energy Path (MEP) search
 results_img(ntimimage_stored,nimage)=datastructure that holds the history of previous computations.
   results_img(:,:)%acell(3)
    at input, history of the values of acell for all images
    at output, the predicted values of acell for all images
   results_img(:,:)%results_gs
    at input, history of the values of energies and forces for all images
   results_img(:,:)%rprim(3,3)
    at input, history of the values of rprim for all images
    at output, the predicted values of rprim for all images
   results_img(:,:)%vel(3,natom)
    at input, history of the values of vel for all images
    at output, the predicted values of vel for all images
   results_img(:,:)%vel_cell(3,3)
    at input, history of the values of vel_cell for all images
    at output, the predicted values of vel_cell for all images
   results_img(:,:)%xred(3,natom)
    at input, history of the values of xred for all images
    at output, the predicted values of xred for all images

PARENTS

      predictimg

CHILDREN

      gather_array_img,get_geometry_img,mep_gbfgs,mep_lbfgs,mep_qmin
      mep_steepest,xmpi_bcast

SOURCE

103 subroutine predict_neb(itimimage,itimimage_eff,list_dynimage,mep_param,mpi_enreg,natom,&
104 &                      ndynimage,nimage,nimage_tot,ntimimage_stored,results_img)
105 
106 
107 !This section has been created automatically by the script Abilint (TD).
108 !Do not modify the following lines by hand.
109 #undef ABI_FUNC
110 #define ABI_FUNC 'predict_neb'
111 !End of the abilint section
112 
113  implicit none
114 
115 !Arguments ------------------------------------
116 !scalars
117  integer,intent(in) :: itimimage,itimimage_eff,natom,ndynimage
118  integer,intent(in) :: nimage,nimage_tot,ntimimage_stored
119  type(mep_type),intent(inout) :: mep_param
120  type(MPI_type),intent(in) :: mpi_enreg
121 !arrays
122  integer,intent(in)     :: list_dynimage(ndynimage)
123  type(results_img_type),intent(inout) :: results_img(nimage,ntimimage_stored)
124 
125 !Local variables-------------------------------
126 !scalars
127  integer :: ierr,ii,iimage,iimage_min,iimage_max,next_itimimage
128  real(dp) :: dvmax,dvmin,ecur,emax,emin,eref,f_para1,f_para2
129  logical :: test_minus_one,test_plus_one
130 !arrays
131  real(dp),allocatable :: buffer(:,:),buffer_all(:,:)
132  real(dp),allocatable :: coordif(:,:,:),dimage(:),rprimd(:,:,:),spring(:)
133  real(dp),allocatable :: tangent(:,:,:),vect(:,:)
134  real(dp),allocatable,target :: etotal(:),fcart(:,:,:),neb_forces(:,:,:)
135  real(dp),allocatable,target :: xcart(:,:,:),xred(:,:,:)
136  real(dp),pointer :: etotal_all(:),fcart_all(:,:,:),neb_forces_all(:,:,:)
137  real(dp),pointer :: xcart_all(:,:,:)
138 
139 ! *************************************************************************
140 
141  ABI_ALLOCATE(xred,(3,natom,nimage))
142 
143 !Parallelism over images: only one process per image of the cell
144  if (mpi_enreg%me_cell==0) then
145 
146 !  Retrieve positions and forces
147    ABI_ALLOCATE(etotal,(nimage))
148    ABI_ALLOCATE(xcart,(3,natom,nimage))
149    ABI_ALLOCATE(fcart,(3,natom,nimage))
150    ABI_ALLOCATE(rprimd,(3,3,nimage))
151    call get_geometry_img(etotal,natom,nimage,results_img(:,itimimage_eff),&
152 &   fcart,rprimd,xcart,xred)
153 
154 !  Array containing effective NEB forces (F_ortho+F_spring)
155    ABI_ALLOCATE(neb_forces,(3,natom,nimage))
156 
157 !  Parallelism: gather data of all images
158    if (mpi_enreg%paral_img==1) then
159      ABI_ALLOCATE(buffer,(6*natom+1,nimage))
160      ABI_ALLOCATE(buffer_all,(6*natom+1,nimage_tot))
161      ABI_ALLOCATE(xcart_all,(3,natom,nimage_tot))
162      ABI_ALLOCATE(fcart_all,(3,natom,nimage_tot))
163      ABI_ALLOCATE(etotal_all,(nimage_tot))
164      ABI_ALLOCATE(neb_forces_all,(3,natom,nimage_tot))
165      buffer=zero;ii=0
166      buffer(ii+1:ii+3*natom,1:nimage)=reshape(xcart,(/3*natom,nimage/));ii=3*natom
167      buffer(ii+1:ii+3*natom,1:nimage)=reshape(fcart,(/3*natom,nimage/));ii=6*natom
168      buffer(ii+1           ,1:nimage)=etotal(1:nimage);ii=0
169      call gather_array_img(buffer,buffer_all,mpi_enreg,allgather=.true.)
170      xcart_all(:,:,:)=reshape(buffer_all(ii+1:ii+3*natom,1:nimage_tot),(/3,natom,nimage_tot/));ii=3*natom
171      fcart_all(:,:,:)=reshape(buffer_all(ii+1:ii+3*natom,1:nimage_tot),(/3,natom,nimage_tot/));ii=6*natom
172      etotal_all(:)=buffer_all(ii+1,1:nimage_tot);ii=0
173      ABI_DEALLOCATE(buffer)
174      ABI_DEALLOCATE(buffer_all)
175    else
176      xcart_all      => xcart
177      fcart_all      => fcart
178      etotal_all     => etotal
179      neb_forces_all => neb_forces
180    end if
181 
182 !  coordif is the vector between two images
183 !  dimage is the distance between two images
184 !  tangent is the tangent at image i
185    ABI_ALLOCATE(coordif,(3,natom,nimage_tot))
186    ABI_ALLOCATE(tangent,(3,natom,nimage_tot))
187    ABI_ALLOCATE(dimage,(nimage_tot))
188 
189 !  Compute distances between images
190    coordif(:,:,1)=zero;dimage(1)=zero
191    do iimage=2,nimage_tot
192      coordif(:,:,iimage)=xcart_all(:,:,iimage)-xcart_all(:,:,iimage-1)
193      dimage(iimage)=mep_img_norm(coordif(:,:,iimage))
194    end do
195 
196 !  Compute tangent (not normalized)
197    tangent(:,:,1)=zero
198    tangent(:,:,nimage_tot)=zero
199 !  === Original definition of tangent
200    if (mep_param%neb_algo==0) then
201      do iimage=2,nimage_tot-1
202        tangent(:,:,iimage)=coordif(:,:,iimage  )/dimage(iimage) &
203 &       +coordif(:,:,iimage+1)/dimage(iimage+1)
204      end do
205 !    === Improved tangent (J. Chem. Phys. 113, 9978 (2000) [[cite:Henkelman2000]])
206    else
207      do iimage=2,nimage_tot-1
208        test_minus_one=(etotal_all(iimage-1)<etotal_all(iimage))
209        test_plus_one =(etotal_all(iimage  )<etotal_all(iimage+1))
210        if (test_minus_one.and.test_plus_one)then       ! V_i-1 < V_i < V_i+1
211          tangent(:,:,iimage)=coordif(:,:,iimage+1)
212        else if ((.not.test_minus_one).and.(.not.test_plus_one))then     ! V_i-1 > V_i > V_i+1
213          tangent(:,:,iimage)=coordif(:,:,iimage)
214        else                             ! V_i-1 < V_i > V_i+1  OR  V_i-1 > V_i < V_i+1
215          dvmax=max(abs(etotal_all(iimage+1)-etotal_all(iimage)),&
216 &         abs(etotal_all(iimage-1)-etotal_all(iimage)))
217          dvmin=min(abs(etotal_all(iimage+1)-etotal_all(iimage)),&
218 &         abs(etotal_all(iimage-1)-etotal_all(iimage)))
219          if (etotal_all(iimage+1)>etotal_all(iimage-1)) then   ! V_i+1 > V_i-1
220            tangent(:,:,iimage)=coordif(:,:,iimage+1)*dvmax &
221 &           +coordif(:,:,iimage  )*dvmin
222          else                                                  ! V_i+1 < V_i-1
223            tangent(:,:,iimage)=coordif(:,:,iimage+1)*dvmin &
224 &           +coordif(:,:,iimage  )*dvmax
225          end if
226        end if
227      end do
228    end if
229 
230 !  Normalize tangent
231    do iimage=2,nimage_tot-1
232      tangent(:,:,iimage)=tangent(:,:,iimage)/mep_img_norm(tangent(:,:,iimage))
233    end do
234 
235 !  Compute spring constant(s)
236    ABI_ALLOCATE(spring,(nimage_tot))
237    if (abs(mep_param%neb_spring(2)-mep_param%neb_spring(1))<tol8) then ! Constant spring
238      spring(:)=mep_param%neb_spring(1)
239    else                                                                ! Variable spring
240      emax=maxval(etotal_all(:))
241      eref=max(etotal_all(1),etotal_all(nimage_tot))
242      spring(1)=zero
243      do iimage=2,nimage_tot
244        ecur=max(etotal_all(iimage-1),etotal_all(iimage))
245        if (ecur<eref) then
246          spring(iimage)=mep_param%neb_spring(1)
247        else
248          spring(iimage)=mep_param%neb_spring(2) &
249 &         -(mep_param%neb_spring(2)-mep_param%neb_spring(1)) &
250 &         *(emax-ecur)/(emax-eref)
251        end if
252      end do
253    end if
254 
255 !  CI-NEB: determine image(s) with maximal energy
256    iimage_min=-1;iimage_max=-1
257    if (mep_param%neb_algo==2.and.itimimage>=mep_param%cineb_start) then
258      emin=min(etotal_all(1),etotal_all(nimage_tot))
259      emax=max(etotal_all(1),etotal_all(nimage_tot))
260      do iimage=2,nimage_tot-1
261        if (etotal_all(iimage)<emin) then
262          iimage_min=iimage;emin=etotal_all(iimage)
263        end if
264        if (etotal_all(iimage)>emax) then
265          iimage_max=iimage;emax=etotal_all(iimage)
266        end if
267      end do
268    end if
269 
270 !  Compute NEB forces
271    neb_forces_all(:,:,1)=fcart_all(:,:,1)
272    neb_forces_all(:,:,nimage_tot)=fcart_all(:,:,nimage_tot)
273    do iimage=2,nimage_tot-1
274 !    === Standard NEB
275      if (iimage/=iimage_min.and.iimage/=iimage_max) then
276        f_para1=mep_img_dotp(fcart_all(:,:,iimage),tangent(:,:,iimage))
277        if (mep_param%neb_algo==0) then   ! Original NEB algo
278          ABI_ALLOCATE(vect,(3,natom))
279          vect(:,:)=spring(iimage+1)*coordif(:,:,iimage+1)-spring(iimage)*coordif(:,:,iimage)
280          f_para2=mep_img_dotp(vect,tangent(:,:,iimage))
281          ABI_DEALLOCATE(vect)
282        else       ! Modification from J. Chem. Phys. 113, 9978 (2000) [[cite:Henkelman2000]]
283          f_para2=spring(iimage+1)*dimage(iimage+1)-spring(iimage)*dimage(iimage)
284        end if
285        neb_forces_all(:,:,iimage)=fcart_all(:,:,iimage) &       ! F_ortho
286 &      -f_para1*tangent(:,:,iimage) &  ! F_ortho
287 &      +f_para2*tangent(:,:,iimage)    ! F_spring
288 !      === CI-NEB for the image(s) with maximal energy
289      else
290        f_para1=mep_img_dotp(fcart_all(:,:,iimage),tangent(:,:,iimage))
291        neb_forces_all(:,:,iimage)=fcart_all(:,:,iimage)-two*f_para1*tangent(:,:,iimage)
292      end if
293    end do
294 
295 !  Parallelism: distribute NEB forces
296    if (mpi_enreg%paral_img==1) then
297      do iimage=1,nimage
298        ii=mpi_enreg%my_imgtab(iimage)
299        neb_forces(:,:,iimage)=neb_forces_all(:,:,ii)
300      end do
301    end if
302 
303 !  Compute new atomic positions in each cell
304    if      (mep_param%mep_solver==0) then ! Steepest-descent
305      call mep_steepest(neb_forces,list_dynimage,mep_param,natom,ndynimage,nimage,rprimd,xcart,xred)
306    else if (mep_param%mep_solver==1) then ! Quick-min
307      call mep_qmin(neb_forces,itimimage,list_dynimage,mep_param,natom,ndynimage, &
308 &     nimage,rprimd,xcart,xred)
309    else if (mep_param%mep_solver==2) then ! Local BFGS
310      call mep_lbfgs(neb_forces,itimimage,list_dynimage,mep_param,natom,ndynimage,&
311 &     nimage,rprimd,xcart,xred)
312    else if (mep_param%mep_solver==3) then ! Global BFGS
313      call mep_gbfgs(neb_forces,itimimage,list_dynimage,mep_param,mpi_enreg,natom,ndynimage,&
314 &     nimage,nimage_tot,rprimd,xcart,xred)
315    else
316      MSG_BUG("Inconsistent solver !")
317    end if
318 
319 !  Free memory
320    ABI_DEALLOCATE(spring)
321    ABI_DEALLOCATE(coordif)
322    ABI_DEALLOCATE(tangent)
323    ABI_DEALLOCATE(dimage)
324    if (mpi_enreg%paral_img==1)  then
325      ABI_DEALLOCATE(xcart_all)
326      ABI_DEALLOCATE(fcart_all)
327      ABI_DEALLOCATE(etotal_all)
328      ABI_DEALLOCATE(neb_forces_all)
329    end if
330 
331    ABI_DEALLOCATE(neb_forces)
332    ABI_DEALLOCATE(etotal)
333    ABI_DEALLOCATE(xcart)
334    ABI_DEALLOCATE(fcart)
335    ABI_DEALLOCATE(rprimd)
336  end if ! mpi_enreg%me_cell==0
337 
338 !Store acell, rprim, xred and vel for the new iteration
339  call xmpi_bcast(xred,0,mpi_enreg%comm_cell,ierr)
340  next_itimimage=itimimage_eff+1
341  if (next_itimimage>ntimimage_stored) next_itimimage=1
342  do iimage=1,nimage
343    results_img(iimage,next_itimimage)%xred(:,:)    =xred(:,:,iimage)
344    results_img(iimage,next_itimimage)%acell(:)     =results_img(iimage,itimimage_eff)%acell(:)
345    results_img(iimage,next_itimimage)%rprim(:,:)   =results_img(iimage,itimimage_eff)%rprim(:,:)
346    results_img(iimage,next_itimimage)%vel(:,:)     =results_img(iimage,itimimage_eff)%vel(:,:)
347    results_img(iimage,next_itimimage)%vel_cell(:,:)=results_img(iimage,itimimage_eff)%vel_cell(:,:)
348  end do
349  ABI_DEALLOCATE(xred)
350 
351 end subroutine predict_neb