TABLE OF CONTENTS


ABINIT/m_mep [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mep

FUNCTION

  This module provides several routines and datatypes for the
  Minimal Energy Path (MEP) search implementation.

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 .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

26 #if defined HAVE_CONFIG_H
27 #include "config.h"
28 #endif
29 
30 #include "abi_common.h"
31 
32 MODULE m_mep
33 
34  use defs_basis
35  use defs_abitypes
36  use m_abicore
37  use m_errors
38  use m_xmpi
39 
40  use m_geometry,    only : fred2fcart, fcart2fred, xcart2xred, xred2xcart, metric
41  use m_bfgs,        only : hessupdt
42  use m_results_img, only : results_img_type,gather_array_img
43 
44  implicit none
45 
46  private
47 
48 !public procedures
49  public :: mep_init
50  public :: mep_destroy
51  public :: mep_steepest
52  public :: mep_qmin
53  public :: mep_lbfgs
54  public :: mep_gbfgs
55  public :: mep_rk4
56  public :: mep_img_dotp
57  public :: mep_img_norm
58  public :: mep_img_dotp_red
59  public :: mep_img_norm_red

m_mep/mep_destroy [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_destroy

FUNCTION

  Destroy the content of a datastructure of type mep_type.

INPUTS

OUTPUT

SIDE EFFECTS

  mep_param=datastructure of type mep_type.
            several parameters for Minimal Energy Path (MEP) search.

PARENTS

      gstateimg

CHILDREN

      wrtout,xcart2xred,xred2xcart

SOURCE

189 subroutine mep_destroy(mep_param)
190 
191 
192 !This section has been created automatically by the script Abilint (TD).
193 !Do not modify the following lines by hand.
194 #undef ABI_FUNC
195 #define ABI_FUNC 'mep_destroy'
196 !End of the abilint section
197 
198  implicit none
199 
200 !Arguments ------------------------------------
201 !scalars
202  type(mep_type),intent(inout) :: mep_param
203 
204 !************************************************************************
205 
206  if (allocated(mep_param%bfgs_xprev)) then
207    ABI_DEALLOCATE(mep_param%bfgs_xprev)
208  end if
209  if (allocated(mep_param%gbfgs_hess))     then
210    ABI_DEALLOCATE(mep_param%gbfgs_hess)
211  end if
212  if (allocated(mep_param%bfgs_fprev)) then
213    ABI_DEALLOCATE(mep_param%bfgs_fprev)
214  end if
215  if (allocated(mep_param%lbfgs_hess))     then
216    ABI_DEALLOCATE(mep_param%lbfgs_hess)
217  end if
218  if (allocated(mep_param%qmin_vel))     then
219    ABI_DEALLOCATE(mep_param%qmin_vel)
220  end if
221  if (allocated(mep_param%rk4_xcart1)) then
222    ABI_DEALLOCATE(mep_param%rk4_xcart1)
223  end if
224  if (allocated(mep_param%rk4_fcart1)) then
225    ABI_DEALLOCATE(mep_param%rk4_fcart1)
226  end if
227  if (allocated(mep_param%rk4_fcart2)) then
228    ABI_DEALLOCATE(mep_param%rk4_fcart2)
229  end if
230  if (allocated(mep_param%rk4_fcart3)) then
231    ABI_DEALLOCATE(mep_param%rk4_fcart3)
232  end if
233 
234  nullify(mep_param%iatfix)
235 
236 end subroutine mep_destroy

m_mep/mep_gbfgs [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_gbfgs

FUNCTION

  Make a path (string of images) evolve according to a
  global Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm

INPUTS

  itime=time step
  list_dynimage(nimage)=list of dynamical images.
  mep_param=datastructure of type mep_type.
  mpi_enreg=MPI-parallelisation information
  mep_param=several parameters for Minimal Energy Path (MEP) search.
  natom=number of atoms
  ndynimage=number of dynamical images along the path
  nimage=number of images (including static ones)
  nimage_tot=total number of images
  rprimd(3,3,nimage)=dimensional primitive translations for each image along the path

OUTPUT

SIDE EFFECTS

  mep_param=datastructure of type mep_type.
            History for Runge-Kutta algorithm is filled up
  xcart(3,natom,nimage)=cartesian coordinates of atoms in each image along the path
                        before and after time evolution
  xred(3,natom,nimage)=reduced coordinates of atoms in each image along the path
                       before and after time evolution

NOTES

  Could see Numerical Recipes (Fortran), 1986, page 307.
  Has to work in cartesian coordinates

PARENTS

      predict_neb

CHILDREN

      wrtout,xcart2xred,xred2xcart

SOURCE

 712 subroutine mep_gbfgs(fcart,itime,list_dynimage,mep_param,mpi_enreg,natom,&
 713 &                    ndynimage,nimage,nimage_tot,rprimd,xcart,xred)
 714 
 715 
 716 !This section has been created automatically by the script Abilint (TD).
 717 !Do not modify the following lines by hand.
 718 #undef ABI_FUNC
 719 #define ABI_FUNC 'mep_gbfgs'
 720 !End of the abilint section
 721 
 722  implicit none
 723 
 724 !Arguments ------------------------------------
 725 !scalars
 726  integer,intent(in) :: itime,natom,ndynimage,nimage,nimage_tot
 727  type(mep_type),intent(inout) :: mep_param
 728  type(MPI_type),intent(in) :: mpi_enreg
 729 !arrays
 730  integer,intent(in) :: list_dynimage(ndynimage)
 731  real(dp),intent(in) :: fcart(3,natom,nimage),rprimd(3,3,nimage)
 732  real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage)
 733 !Local variables-------------------------------
 734 !scalars
 735  integer :: iatom,idynimage,ii,iimage,iimage_tot,indi,indj,ierr
 736  integer :: jdynimage,jatom,jj,mu,ndynimage_tot,nu
 737  logical :: reset
 738  real(dp),parameter :: initial_Hessian=1._dp ! in Bohr^2/Hartree
 739  real(dp) :: dot1,dot2,stepsize,ucvol
 740  character(len=500) :: msg
 741 !arrays
 742  integer,allocatable :: dynimage_tot(:),iatfix_fake(:,:),ind_dynimage_tot(:)
 743  integer,allocatable :: list_dynimage_tot(:)
 744  real(dp) :: favg(3),gprimd(3,3),gmet(3,3),rmet(3,3)
 745  real(dp),allocatable :: buffer(:,:),buffer_all(:,:),fred(:,:)
 746  real(dp),allocatable :: fcart_all(:,:,:),fcartp_all(:,:,:)
 747  real(dp),allocatable :: gmet_all(:,:,:),gprimd_all(:,:,:),rprimd_all(:,:,:)
 748  real(dp),allocatable :: xcart_all(:,:,:),xcartp_all(:,:,:),xred_old(:,:),xstep_all(:,:,:)
 749 
 750 !************************************************************************
 751 
 752 !Retrieve indexes of all dynamical images
 753  ABI_ALLOCATE(ind_dynimage_tot,(nimage_tot))
 754  if (mpi_enreg%paral_img==1) then
 755    ABI_ALLOCATE(dynimage_tot,(nimage_tot))
 756    dynimage_tot=0
 757    do idynimage=1,ndynimage
 758      iimage=list_dynimage(idynimage)
 759      iimage_tot=mpi_enreg%my_imgtab(iimage)
 760      dynimage_tot(iimage_tot)=1
 761    end do
 762    call xmpi_sum(dynimage_tot,mpi_enreg%comm_img,ierr)
 763    ndynimage_tot=count(dynimage_tot(:)>0)
 764    ABI_ALLOCATE(list_dynimage_tot,(ndynimage_tot))
 765    idynimage=0;ind_dynimage_tot(:)=-1
 766    do iimage_tot=1,nimage_tot
 767      if (dynimage_tot(iimage_tot)>0) then
 768        idynimage=idynimage+1
 769        ind_dynimage_tot(iimage_tot)=idynimage
 770        list_dynimage_tot(idynimage)=iimage_tot
 771      end if
 772    end do
 773    ABI_DEALLOCATE(dynimage_tot)
 774  else
 775    ndynimage_tot=ndynimage
 776    ABI_ALLOCATE(list_dynimage_tot,(ndynimage))
 777    ind_dynimage_tot(:)=-1
 778    do idynimage=1,ndynimage
 779      ind_dynimage_tot(list_dynimage(idynimage))=idynimage
 780      list_dynimage_tot(idynimage)=list_dynimage(idynimage)
 781    end do
 782  end if
 783 
 784 !Allocate history array (at first time step)
 785  if (itime==1) then
 786    if (allocated(mep_param%bfgs_xprev)) then
 787      ABI_DEALLOCATE(mep_param%bfgs_xprev)
 788    end if
 789    if (allocated(mep_param%bfgs_fprev)) then
 790      ABI_DEALLOCATE(mep_param%bfgs_fprev)
 791    end if
 792    if (allocated(mep_param%gbfgs_hess)) then
 793      ABI_DEALLOCATE(mep_param%gbfgs_hess)
 794    end if
 795    ABI_ALLOCATE(mep_param%bfgs_xprev,(3,natom,ndynimage))
 796    ABI_ALLOCATE(mep_param%bfgs_fprev,(3,natom,ndynimage))
 797    ABI_ALLOCATE(mep_param%gbfgs_hess,(3*natom*ndynimage_tot,3*natom*ndynimage_tot))
 798    mep_param%bfgs_xprev=zero
 799    mep_param%bfgs_fprev=zero
 800  end if
 801 
 802 !Retrieve positions and forces for all images
 803  ABI_ALLOCATE(xcart_all,(3,natom,ndynimage_tot))
 804  ABI_ALLOCATE(fcart_all,(3,natom,ndynimage_tot))
 805  ABI_ALLOCATE(xcartp_all,(3,natom,ndynimage_tot))
 806  ABI_ALLOCATE(fcartp_all,(3,natom,ndynimage_tot))
 807  ABI_ALLOCATE(rprimd_all,(3,3,ndynimage_tot))
 808  ABI_ALLOCATE(gprimd_all,(3,3,ndynimage_tot))
 809  ABI_ALLOCATE(gmet_all,(3,3,ndynimage_tot))
 810  if (mpi_enreg%paral_img==1) then
 811    ABI_ALLOCATE(buffer,(12*natom+27,nimage))
 812    ABI_ALLOCATE(buffer_all,(12*natom+27,nimage_tot))
 813    buffer=zero
 814    do idynimage=1,ndynimage
 815      iimage=list_dynimage(idynimage)
 816      call metric(gmet,gprimd,-1,rmet,rprimd(:,:,iimage),ucvol)
 817      buffer(         1:3 *natom  ,iimage)=reshape(xcart(:,:,iimage),(/3*natom/))
 818      buffer(3 *natom+1:6 *natom  ,iimage)=reshape(fcart(:,:,iimage),(/3*natom/))
 819      buffer(6 *natom+1:9 *natom  ,iimage)=reshape(mep_param%bfgs_xprev(:,:,idynimage),(/3*natom/))
 820      buffer(9 *natom+1:12*natom  ,iimage)=reshape(mep_param%bfgs_fprev(:,:,idynimage),(/3*natom/))
 821      buffer(12*natom+1:12*natom+9,iimage)=reshape(rprimd(:,:,iimage),(/9/))
 822      buffer(12*natom+10:12*natom+18,iimage)=reshape(gprimd(:,:),(/9/))
 823      buffer(12*natom+19:12*natom+27,iimage)=reshape(gmet(:,:),(/9/))
 824    end do
 825    call gather_array_img(buffer,buffer_all,mpi_enreg,allgather=.true.)
 826    do idynimage=1,ndynimage_tot
 827      iimage_tot=list_dynimage_tot(idynimage)
 828      do iatom=1,natom
 829        indi=12*(iatom-1)
 830        do ii=1,3
 831          xcart_all (ii,iatom,idynimage)= buffer_all(indi  +ii,iimage_tot)
 832          fcart_all (ii,iatom,idynimage)=-buffer_all(indi+3+ii,iimage_tot) ! use fcart=-cartesian_force
 833          xcartp_all(ii,iatom,idynimage)= buffer_all(indi+6+ii,iimage_tot)
 834          fcartp_all(ii,iatom,idynimage)= buffer_all(indi+9+ii,iimage_tot)
 835        end do
 836      end do
 837      indi=12*natom
 838      rprimd_all(1:3,1:3,idynimage)=reshape(buffer_all(indi+ 1:indi+ 9,iimage_tot),(/3,3/))
 839      gprimd_all(1:3,1:3,idynimage)=reshape(buffer_all(indi+10:indi+18,iimage_tot),(/3,3/))
 840      gmet_all  (1:3,1:3,idynimage)=reshape(buffer_all(indi+19:indi+27,iimage_tot),(/3,3/))
 841    end do
 842    ABI_DEALLOCATE(buffer)
 843    ABI_DEALLOCATE(buffer_all)
 844  else
 845    do idynimage=1,ndynimage
 846      iimage=list_dynimage(idynimage)
 847      xcart_all(:,:,idynimage)= xcart(:,:,iimage)
 848      fcart_all(:,:,idynimage)=-fcart(:,:,iimage) ! use fcart=-cartesian_force
 849      xcartp_all(:,:,idynimage)=mep_param%bfgs_xprev(:,:,idynimage)
 850      fcartp_all(:,:,idynimage)=mep_param%bfgs_fprev(:,:,idynimage)
 851      rprimd_all(:,:,idynimage)=rprimd(:,:,iimage)
 852      call metric(gmet_all(:,:,idynimage),gprimd_all(:,:,idynimage),-1,rmet,rprimd(:,:,iimage),ucvol)
 853    end do
 854  end if
 855 
 856 !Test if a reset is needed
 857  reset=.false.
 858  if (itime>1) then
 859    dot1=zero;dot2=zero
 860    do idynimage=1,ndynimage_tot
 861      dot1=dot2+mep_img_dotp(fcartp_all(:,:,idynimage),fcart_all (:,:,idynimage))
 862      dot2=dot1+mep_img_dotp(fcartp_all(:,:,idynimage),fcartp_all(:,:,idynimage))
 863    end do
 864    reset=((dot2<two*abs(dot1)).or.abs(dot2)<tol8)
 865    if (reset) then
 866      msg=' Resetting Hessian matrix.'
 867      call wrtout(std_out,msg,'COLL')
 868      call wrtout(ab_out ,msg,'COLL')
 869    end if
 870  end if
 871 
 872 !===> First step or reset: initialize the Hessian matrix
 873  if (itime==1.or.reset) then
 874    mep_param%gbfgs_hess(:,:)=zero
 875    do idynimage=1,ndynimage_tot
 876      indi=3*natom*(idynimage-1)
 877      do iatom=1,natom
 878        do mu=1,3
 879          do nu=1,3
 880            do ii=1,3
 881              do jj=1,3
 882                if (mep_param%iatfix(ii,iatom)==0.and. &
 883 &                  mep_param%iatfix(jj,iatom)==0) then
 884                    mep_param%gbfgs_hess(indi+mu,indi+nu)=mep_param%gbfgs_hess(indi+mu,indi+nu) &
 885 &                    +rprimd_all(mu,ii,idynimage)*rprimd_all(nu,jj,idynimage) &
 886 &                    *gmet_all(ii,jj,idynimage)*initial_Hessian
 887                end if
 888              end do
 889            end do
 890          end do
 891        end do
 892        indi=indi+3
 893      end do
 894    end do
 895 
 896 !===> Other steps: update the Hessian matrix
 897  else
 898 
 899 !  Impose here f-fprev=0 (cannot be done inside hessupdt in cartesian coordinates)
 900    ABI_ALLOCATE(fred,(3,natom))
 901    do idynimage=1,ndynimage_tot
 902      fcartp_all(:,:,idynimage)=fcartp_all(:,:,idynimage)-fcart_all(:,:,idynimage)
 903      call fcart2fred(fcartp_all(:,:,idynimage),fred,rprimd_all(:,:,idynimage),natom)
 904      where (mep_param%iatfix(:,:)==1) ! iatfix is defined in reduced coordinates
 905        fred(:,:)=zero
 906      end where
 907      call fred2fcart(favg,.TRUE.,fcartp_all(:,:,idynimage),fred,gprimd_all(:,:,idynimage),natom)
 908      do iatom=1,natom
 909        fcartp_all(:,iatom,idynimage)=fcartp_all(:,iatom,idynimage) &
 910 &                                   +fcart_all(:,iatom,idynimage)+favg(:)
 911      end do
 912    end do
 913    ABI_DEALLOCATE(fred)
 914 
 915 !  f-fprev=0 has already been imposed for fixed atoms:
 916 !  we call hessupdt with no fixed atom
 917    ABI_ALLOCATE(iatfix_fake,(3,natom))
 918    iatfix_fake(:,:)=0
 919    call hessupdt(mep_param%gbfgs_hess,&
 920 &                iatfix_fake,natom,3*natom*ndynimage_tot, &
 921                  xcart_all,xcartp_all,fcart_all,fcartp_all, &
 922 &                nimage=ndynimage_tot)
 923    ABI_DEALLOCATE(iatfix_fake)
 924  end if
 925 
 926 !Free memory
 927  ABI_DEALLOCATE(xcart_all)
 928  ABI_DEALLOCATE(xcartp_all)
 929  ABI_DEALLOCATE(fcartp_all)
 930  ABI_DEALLOCATE(rprimd_all)
 931  ABI_DEALLOCATE(gprimd_all)
 932  ABI_DEALLOCATE(gmet_all)
 933 
 934 !Update history
 935  do idynimage=1,ndynimage
 936    iimage=list_dynimage(idynimage)
 937    mep_param%bfgs_xprev(:,:,idynimage)=xcart(:,:,iimage)
 938    mep_param%bfgs_fprev(:,:,idynimage)=fcart(:,:,iimage)
 939  end do
 940 
 941 !Compute image step
 942  ABI_ALLOCATE(xstep_all,(3,natom,ndynimage_tot))
 943  xstep_all=zero
 944  do idynimage=1,ndynimage_tot
 945    indi=3*natom*(idynimage-1)
 946    do iatom=1,natom
 947      do ii=1,3
 948        do jdynimage=1,ndynimage_tot
 949          indj=3*natom*(jdynimage-1)
 950          do jatom=1,natom
 951            do jj=1,3
 952 !            Be careful: minus sign because fcart=-cartesian_force
 953              xstep_all(ii,iatom,idynimage)=xstep_all(ii,iatom,idynimage) &
 954 &                           -fcart_all(jj,jatom,jdynimage) &
 955 &                           *mep_param%gbfgs_hess(indi+ii,indj+jj)
 956            end do
 957            indj=indj+3
 958          end do
 959        end do
 960      end do
 961      indi=indi+3
 962    end do
 963  end do
 964 
 965 !Restrict image step size
 966  stepsize=zero
 967  do idynimage=1,ndynimage_tot
 968    stepsize=stepsize+mep_img_dotp(xstep_all(:,:,idynimage),xstep_all(:,:,idynimage))
 969  end do
 970  stepsize=sqrt(stepsize)
 971  if (stepsize>=mep_param%mep_mxstep*dble(ndynimage_tot)) then
 972    xstep_all=xstep_all*mep_param%mep_mxstep*dble(ndynimage_tot)/stepsize
 973    write(msg,'(a,i3,a)') " Restricting BFGS step size."
 974    call wrtout(std_out,msg,'COLL')
 975    call wrtout(ab_out ,msg,'COLL')
 976  end if
 977 
 978 !Update positions
 979  ABI_ALLOCATE(xred_old,(3,natom))
 980  do idynimage=1,ndynimage
 981    iimage=list_dynimage(idynimage)
 982    iimage_tot=mpi_enreg%my_imgtab(iimage)
 983    xred_old(:,:)=xred(:,:,iimage)
 984    xcart(:,:,iimage)=xcart(:,:,iimage)+xstep_all(:,:,ind_dynimage_tot(iimage_tot))
 985    call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage))
 986 !  In case atom is fixed, we restore its previous value
 987    do iatom=1,natom
 988      if (any(mep_param%iatfix(:,iatom)==1)) then
 989        where(mep_param%iatfix(:,iatom)==1)
 990          xred(:,iatom,iimage)=xred_old(:,iatom)
 991        end where
 992        call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage))
 993      end if
 994    end do
 995  end do
 996  ABI_DEALLOCATE(xred_old)
 997 
 998 !Free memory
 999  ABI_DEALLOCATE(fcart_all)
1000  ABI_DEALLOCATE(xstep_all)
1001  ABI_DEALLOCATE(ind_dynimage_tot)
1002  ABI_DEALLOCATE(list_dynimage_tot)
1003 
1004 end subroutine mep_gbfgs

m_mep/mep_img_dotp [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_img_dotp

FUNCTION

  Compute the dot product of two vectors in the configuration space:
    Vect1(3,natom).Vect2(3,natom)

INPUTS

  vect1(3,natom)=input vector 1
  vect2(3,natom)=input vector 2

OUTPUT

  mep_img_dotp=dot product

PARENTS

CHILDREN

SOURCE

1211 function mep_img_dotp(vect1,vect2)
1212 
1213 
1214 !This section has been created automatically by the script Abilint (TD).
1215 !Do not modify the following lines by hand.
1216 #undef ABI_FUNC
1217 #define ABI_FUNC 'mep_img_dotp'
1218 !End of the abilint section
1219 
1220  implicit none
1221 
1222 !Arguments ------------------------------------
1223 !scalars
1224  real(dp) :: mep_img_dotp
1225 !arrays
1226  real(dp),intent(in) :: vect1(:,:),vect2(:,:)
1227 !Local variables-------------------------------
1228 !scalars
1229  integer :: size1,size2
1230 !arrays
1231 
1232 !************************************************************************
1233 
1234  size1=size(vect1,1);size2=size(vect1,2)
1235  if (size1/=size(vect2,1).or.size2/=size(vect2,2)) then
1236    MSG_BUG("Error on dimensions !")
1237  end if
1238 
1239  mep_img_dotp=sum(vect1*vect2)
1240 
1241 end function mep_img_dotp

m_mep/mep_img_dotp_red [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_img_dotp_red

FUNCTION

  Compute the dot product of two vectors in the configuration space:
    Vect1(3,natom).Vect2(3,natom)
  using reduced coordinates

INPUTS

  rmet(3,3)=metric tensor
  vect1(3,natom)=input vector 1
  vect2(3,natom)=input vector 2

OUTPUT

  mep_img_dotp_red=dot product

PARENTS

CHILDREN

SOURCE

1316 function mep_img_dotp_red(rmet,vect1,vect2)
1317 
1318 
1319 !This section has been created automatically by the script Abilint (TD).
1320 !Do not modify the following lines by hand.
1321 #undef ABI_FUNC
1322 #define ABI_FUNC 'mep_img_dotp_red'
1323 !End of the abilint section
1324 
1325  implicit none
1326 
1327 !Arguments ------------------------------------
1328 !scalars
1329  real(dp) :: mep_img_dotp_red
1330 !arrays
1331  real(dp),intent(in) :: rmet(3,3)
1332  real(dp),intent(in) :: vect1(:,:),vect2(:,:)
1333 !Local variables-------------------------------
1334 !scalars
1335  integer :: iatom,ii,jj,size1,size2
1336 
1337 !************************************************************************
1338 
1339  size1=size(vect1,1);size2=size(vect1,2)
1340  if (size1/=size(vect2,1).or.size2/=size(vect2,2).or.size1/=3) then
1341    MSG_BUG("Error on dimensions !")
1342  end if
1343 
1344  mep_img_dotp_red=zero
1345  do iatom=1,size2
1346    do ii=1,3
1347      do jj=1,3
1348        mep_img_dotp_red=mep_img_dotp_red+vect1(ii,iatom)*vect2(jj,iatom)*rmet(ii,jj)
1349      end do
1350    end do
1351  end do
1352 
1353 end function mep_img_dotp_red

m_mep/mep_img_norm [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_img_norm

FUNCTION

  Compute the norm of a vector in the configuration space:
    |Vect(3,natom)|

INPUTS

  vect(3,natom)=input vector

OUTPUT

  mep_img_norm=norm

PARENTS

CHILDREN

SOURCE

1267 function mep_img_norm(vect)
1268 
1269 
1270 !This section has been created automatically by the script Abilint (TD).
1271 !Do not modify the following lines by hand.
1272 #undef ABI_FUNC
1273 #define ABI_FUNC 'mep_img_norm'
1274 !End of the abilint section
1275 
1276  implicit none
1277 
1278 !Arguments ------------------------------------
1279 !scalars
1280  real(dp) :: mep_img_norm
1281 !arrays
1282  real(dp),intent(in) :: vect(:,:)
1283 
1284 !************************************************************************
1285 
1286  mep_img_norm=sqrt(sum(vect*vect))
1287 
1288 end function mep_img_norm

m_mep/mep_img_norm_red [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_img_norm_red

FUNCTION

  Compute the norm of a vector in the configuration space:
    |Vect(3,natom)|
  using reduced coordinates

INPUTS

  rmet(3,3)=metric tensor
  vect(3,natom)=input vector

OUTPUT

  mep_img_norm_red=norm

PARENTS

CHILDREN

SOURCE

1381 function mep_img_norm_red(rmet,vect)
1382 
1383 
1384 !This section has been created automatically by the script Abilint (TD).
1385 !Do not modify the following lines by hand.
1386 #undef ABI_FUNC
1387 #define ABI_FUNC 'mep_img_norm_red'
1388 !End of the abilint section
1389 
1390  implicit none
1391 
1392 !Arguments ------------------------------------
1393 !scalars
1394  real(dp) :: mep_img_norm_red
1395 !arrays
1396  real(dp),intent(in) :: rmet(3,3)
1397  real(dp),intent(in) :: vect(:,:)
1398 
1399 !************************************************************************
1400 
1401  mep_img_norm_red=sqrt(mep_img_dotp_red(rmet,vect,vect))
1402 
1403 end function mep_img_norm_red

m_mep/mep_init [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_init

FUNCTION

  Initialize a datastructure of type mep_type.

INPUTS

  dtset <type(dataset_type)>=all input variables in current dataset

OUTPUT

SIDE EFFECTS

  mep_param=datastructure of type mep_type.
            several parameters for Minimal Energy Path (MEP) search.

PARENTS

      gstateimg

CHILDREN

      wrtout,xcart2xred,xred2xcart

SOURCE

123 subroutine mep_init(dtset,mep_param)
124 
125 
126 !This section has been created automatically by the script Abilint (TD).
127 !Do not modify the following lines by hand.
128 #undef ABI_FUNC
129 #define ABI_FUNC 'mep_init'
130 !End of the abilint section
131 
132  implicit none
133 
134 !Arguments ------------------------------------
135 !scalars
136  type(dataset_type),target,intent(in) :: dtset
137  type(mep_type),intent(inout) :: mep_param
138 
139 !************************************************************************
140 
141  if((dtset%imgmov==1).or.(dtset%imgmov==2).or.(dtset%imgmov==5))then
142    mep_param%cineb_start  = dtset%cineb_start
143    mep_param%mep_solver   = dtset%mep_solver
144    mep_param%neb_algo     = dtset%neb_algo
145    mep_param%string_algo  = dtset%string_algo
146    mep_param%fxcartfactor = dtset%fxcartfactor
147    mep_param%mep_mxstep   = dtset%mep_mxstep
148    mep_param%neb_spring   = dtset%neb_spring
149    mep_param%iatfix       =>dtset%iatfix
150  else
151    mep_param%cineb_start  = -1
152    mep_param%mep_solver   = -1
153    mep_param%neb_algo     = -1
154    mep_param%string_algo  = -1
155    mep_param%fxcartfactor = zero
156    mep_param%mep_mxstep   = 100._dp
157    mep_param%neb_spring   = zero
158    nullify(mep_param%iatfix)
159  end if
160 
161 end subroutine mep_init

m_mep/mep_lbfgs [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_lbfgs

FUNCTION

  Make a path (string of images) evolve according to a
  local Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm

INPUTS

  itime=time step
  list_dynimage(nimage)=list of dynamical images.
  mep_param=datastructure of type mep_type.
            several parameters for Minimal Energy Path (MEP) search.
  natom=number of atoms
  ndynimage=number of dynamical images along the path
  nimage=number of images (including static ones)
  rprimd(3,3,nimage)=dimensional primitive translations for each image along the path

OUTPUT

SIDE EFFECTS

  mep_param=datastructure of type mep_type.
            History for Runge-Kutta algorithm is filled up
  xcart(3,natom,nimage)=cartesian coordinates of atoms in each image along the path
                        before and after time evolution
  xred(3,natom,nimage)=reduced coordinates of atoms in each image along the path
                       before and after time evolution

NOTES

  Could see Numerical Recipes (Fortran), 1986, page 307.

PARENTS

      predict_neb

CHILDREN

      wrtout,xcart2xred,xred2xcart

SOURCE

519 subroutine mep_lbfgs(fcart,itime,list_dynimage,mep_param,natom,ndynimage,&
520 &                    nimage,rprimd,xcart,xred)
521 
522 
523 !This section has been created automatically by the script Abilint (TD).
524 !Do not modify the following lines by hand.
525 #undef ABI_FUNC
526 #define ABI_FUNC 'mep_lbfgs'
527 !End of the abilint section
528 
529  implicit none
530 
531 !Arguments ------------------------------------
532 !scalars
533  integer,intent(in) :: itime,natom,ndynimage,nimage
534  type(mep_type),intent(inout) :: mep_param
535 !arrays
536  integer,intent(in) :: list_dynimage(ndynimage)
537  real(dp),intent(in) :: rprimd(3,3,nimage)
538  real(dp),intent(in) :: fcart(3,natom,nimage)
539  real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage)
540 !Local variables-------------------------------
541 !scalars
542  integer :: iatom,idynimage,ii,iimage,indi,indj,jatom,jj
543  logical :: reset
544  real(dp),parameter :: initial_Hessian=1._dp ! in Bohr^2/Hartree
545  real(dp) :: dot1,dot2,stepsize,ucvol
546  character(len=500) :: msg
547 !arrays
548  real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3)
549  real(dp),allocatable :: fred(:,:),xstep(:,:)
550 
551 !************************************************************************
552 
553 !Allocate history array (at first time step)
554  if (itime==1) then
555    if (allocated(mep_param%bfgs_xprev)) then
556      ABI_DEALLOCATE(mep_param%bfgs_xprev)
557    end if
558    if (allocated(mep_param%bfgs_fprev)) then
559      ABI_DEALLOCATE(mep_param%bfgs_fprev)
560    end if
561    if (allocated(mep_param%lbfgs_hess)) then
562      ABI_DEALLOCATE(mep_param%lbfgs_hess)
563    end if
564    ABI_ALLOCATE(mep_param%bfgs_xprev,(3,natom,ndynimage))
565    ABI_ALLOCATE(mep_param%bfgs_fprev,(3,natom,ndynimage))
566    ABI_ALLOCATE(mep_param%lbfgs_hess,(3*natom,3*natom,ndynimage))
567    mep_param%bfgs_xprev=zero
568    mep_param%bfgs_fprev=zero
569  end if
570 
571 !Temporary storage
572  ABI_ALLOCATE(fred,(3,natom))
573  ABI_ALLOCATE(xstep,(3,natom))
574 
575 !Loop over images
576  do idynimage=1,ndynimage
577    iimage=list_dynimage(idynimage)
578    call metric(gmet,gprimd,-1,rmet,rprimd(:,:,iimage),ucvol)
579    call fcart2fred(fcart(:,:,iimage),fred,rprimd(:,:,iimage),natom)
580 
581 !  Test if a reset is needed
582    reset=.false.
583    if (itime>1) then
584      dot1=mep_img_dotp(mep_param%bfgs_fprev(:,:,idynimage),fred)
585      dot2=mep_img_dotp(mep_param%bfgs_fprev(:,:,idynimage), &
586 &                      mep_param%bfgs_fprev(:,:,idynimage))
587 !     dot1=mep_img_dotp_red(rmet,mep_param%bfgs_fprev(:,:,idynimage),fred)
588 !     dot2=mep_img_dotp_red(rmet,mep_param%bfgs_fprev(:,:,idynimage), &
589 !&                               mep_param%bfgs_fprev(:,:,idynimage))
590      reset=((dot2<two*abs(dot1)).or.abs(dot2)<tol8)
591      if (reset) then
592        write(msg,'(a,i3,a)') " Resetting Hessian matrix for image ",iimage,"."
593        call wrtout(std_out,msg,'COLL')
594        call wrtout(ab_out ,msg,'COLL')
595      end if
596    end if
597 
598 !  ===> First step or reset: initialize the Hessian matrix (in reduced coordinates)
599    if (itime==1.or.reset) then
600      mep_param%lbfgs_hess(:,:,idynimage)=zero
601      do iatom=1,natom
602        indi=3*(iatom-1)
603        do ii=1,3
604          do jj=1,3
605            if (mep_param%iatfix(ii,iatom)==0.and. &
606 &              mep_param%iatfix(jj,iatom)==0) then
607              mep_param%lbfgs_hess(indi+ii,indi+jj,idynimage)=gmet(ii,jj)*initial_Hessian
608            end if
609          end do
610        end do
611      end do
612 
613 !  ===> Other steps: update the Hessian matrix
614    else
615      call hessupdt(mep_param%lbfgs_hess(:,:,idynimage),&
616 &                  mep_param%iatfix,natom,3*natom, &
617                    xred(:,:,iimage),mep_param%bfgs_xprev(:,:,idynimage),&
618                    fred(:,:),mep_param%bfgs_fprev(:,:,idynimage))
619    end if
620 
621 !  Update history
622    mep_param%bfgs_xprev(:,:,idynimage)=xred(:,:,iimage)
623    mep_param%bfgs_fprev(:,:,idynimage)=fred(:,:)
624 
625 !  Compute image step
626    xstep=zero
627    do iatom=1,natom
628      indi=3*(iatom-1)
629      do ii=1,3
630        do jatom=1,natom
631          indj=3*(jatom-1)
632          do jj=1,3
633            xstep(ii,iatom)=xstep(ii,iatom) &
634 &             -mep_param%lbfgs_hess(indi+ii,indj+jj,idynimage)*fred(jj,jatom)
635          end do
636        end do
637      end do
638    end do
639 
640 !  Restrict image step size
641    stepsize=mep_img_norm_red(rmet,xstep)
642    if (stepsize>=mep_param%mep_mxstep) then
643      xstep=xstep*mep_param%mep_mxstep/stepsize
644      write(msg,'(a,i3,a)') " Restricting BFGS step size of image ",iimage,"."
645      call wrtout(std_out,msg,'COLL')
646      call wrtout(ab_out ,msg,'COLL')
647    end if
648 
649 !  Update positions
650    xred(:,:,iimage)=xred(:,:,iimage)+xstep(:,:)
651 
652 !  In case atom is fixed, we restore its previous value
653    where(mep_param%iatfix(:,:)==1)
654      xred(:,:,iimage)=mep_param%bfgs_xprev(:,:,idynimage)
655    end where
656 
657    call xred2xcart(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage))
658 
659 !End loop over images
660  end do
661 
662  ABI_DEALLOCATE(fred)
663  ABI_DEALLOCATE(xstep)
664 
665 end subroutine mep_lbfgs

m_mep/mep_qmin [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_qmin

FUNCTION

  Make a path (string of images) evolve according to a quick-minimizer algorithm

INPUTS

  fcart(3,natom,nimage)=cartesian forces in each image along the path
  itime=time step
  list_dynimage(nimage)=list of dynamical images.
  mep_param=datastructure of type mep_type.
            several parameters for Minimal Energy Path (MEP) search.
  natom=number of atoms
  ndynimage=number of dynamical images along the path
  nimage=number of images (including static ones)
  results_img(nimage)=datastructure that hold data for each image
                      (positions, forces, energy, ...)
  rprimd(3,3,nimage)=dimensional primitive translations for each image along the path

OUTPUT

SIDE EFFECTS

  xcart(3,natom,nimage)=cartesian coordinates of atoms in each image along the path
                        before and after time evolution
  xred(3,natom,nimage)=reduced coordinates of atoms in each image along the path
                       before and after time evolution

PARENTS

      predict_neb

CHILDREN

      wrtout,xcart2xred,xred2xcart

SOURCE

383 subroutine mep_qmin(fcart,itime,list_dynimage,mep_param,natom,ndynimage,nimage,rprimd,xcart,xred)
384 
385 
386 !This section has been created automatically by the script Abilint (TD).
387 !Do not modify the following lines by hand.
388 #undef ABI_FUNC
389 #define ABI_FUNC 'mep_qmin'
390 !End of the abilint section
391 
392  implicit none
393 
394 !Arguments ------------------------------------
395 !scalars
396  integer,intent(in) :: itime,natom,ndynimage,nimage
397  type(mep_type),intent(inout) :: mep_param
398 !arrays
399  integer,intent(in) :: list_dynimage(ndynimage)
400  real(dp),intent(in) :: fcart(3,natom,nimage),rprimd(3,3,nimage)
401  real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage)
402 !Local variables-------------------------------
403 !scalars
404  integer :: iatom,idynimage,iimage
405  real(dp) :: stepsize,vdotf
406  character(len=500) :: msg
407 !arrays
408  real(dp) :: vel_red(3)
409  real(dp),allocatable :: xred_old(:,:),xstep(:,:)
410 
411 !***********************************************************************
412 
413 !Allocate history array (at first time step)
414  if (itime==1) then
415    if (allocated(mep_param%qmin_vel)) then
416      ABI_DEALLOCATE(mep_param%qmin_vel)
417    end if
418    ABI_ALLOCATE(mep_param%qmin_vel,(3,natom,ndynimage))
419    mep_param%qmin_vel=zero
420  end if
421 
422  ABI_ALLOCATE(xred_old,(3,natom))
423  ABI_ALLOCATE(xstep,(3,natom))
424 
425  do idynimage=1,ndynimage
426    iimage=list_dynimage(idynimage)
427    xred_old(:,:)=xred(:,:,iimage)
428 
429 !  Compute velocities
430    vdotf=mep_img_dotp(mep_param%qmin_vel(:,:,idynimage),fcart(:,:,iimage))
431    if (vdotf>=zero) then
432      mep_param%qmin_vel(:,:,idynimage)=vdotf*fcart(:,:,iimage) &
433 &                              /mep_img_norm(fcart(:,:,iimage))
434    else
435      mep_param%qmin_vel(:,:,idynimage)=zero
436      write(msg,'(a,i3,a)') " Setting velocities of image ",iimage," to zero."
437      call wrtout(std_out,msg,'COLL')
438      call wrtout(ab_out ,msg,'COLL')
439    end if
440    mep_param%qmin_vel(:,:,idynimage)=mep_param%qmin_vel(:,:,idynimage) &
441 &                   +mep_param%fxcartfactor*fcart(:,:,iimage)
442 
443 !  Compute image step
444    xstep(:,:)=mep_param%fxcartfactor*mep_param%qmin_vel(:,:,idynimage)
445    stepsize=mep_img_norm(xstep)
446    if (stepsize>=mep_param%mep_mxstep) then
447      xstep=xstep*mep_param%mep_mxstep/stepsize
448      write(msg,'(a,i3,a)') " Restricting step size of image ",iimage,"."
449      call wrtout(std_out,msg,'COLL')
450      call wrtout(ab_out ,msg,'COLL')
451    end if
452 
453 !  Update positions
454    xcart(:,:,iimage)=xcart(:,:,iimage)+xstep(:,:)
455    call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage))
456 
457 !  In case atom is fixed, we restore its previous value
458    do iatom=1,natom
459      if (any(mep_param%iatfix(:,iatom)==1)) then
460        call xcart2xred(1,rprimd(:,:,iimage),mep_param%qmin_vel(:,iatom,idynimage),vel_red)
461        where(mep_param%iatfix(:,iatom)==1)
462          xred(:,iatom,iimage)=xred_old(:,iatom)
463          vel_red(:)=zero
464        end where
465        call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage))
466        call xred2xcart(1,rprimd(:,:,iimage),mep_param%qmin_vel(:,iatom,idynimage),vel_red)
467      end if
468    end do
469 
470  end do
471 
472  ABI_DEALLOCATE(xred_old)
473  ABI_DEALLOCATE(xstep)
474 
475 end subroutine mep_qmin

m_mep/mep_rk4 [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_rk4

FUNCTION

  Make a path (string of images) evolve according to a fourfth-order Runge-Kutta algorithm

INPUTS

  itime=time step
  list_dynimage(nimage)=list of dynamical images.
  mep_param=datastructure of type mep_type.
            several parameters for Minimal Energy Path (MEP) search.
  natom=number of atoms
  ndynimage=number of dynamical images along the path
  nimage=number of images (including static ones)

OUTPUT

SIDE EFFECTS

  mep_param=datastructure of type mep_type.
            History for Runge-Kutta algorithm is filled up
  xcart(3,natom,nimage)=cartesian coordinates of atoms in each image along the path
                        before and after time evolution
                        after time evolution
  xred(3,natom,nimage)=reduced coordinates of atoms in each image along the path
                       before and after time evolution

PARENTS

      predict_string

CHILDREN

      wrtout,xcart2xred,xred2xcart

SOURCE

1044 subroutine mep_rk4(fcart,itime,list_dynimage,mep_param,natom,ndynimage,nimage,rprimd,xcart,xred)
1045 
1046 
1047 !This section has been created automatically by the script Abilint (TD).
1048 !Do not modify the following lines by hand.
1049 #undef ABI_FUNC
1050 #define ABI_FUNC 'mep_rk4'
1051 !End of the abilint section
1052 
1053  implicit none
1054 
1055 !Arguments ------------------------------------
1056 !scalars
1057  integer,intent(in) :: itime,natom,ndynimage,nimage
1058  type(mep_type),intent(inout) :: mep_param
1059 !arrays
1060  integer,intent(in) :: list_dynimage(ndynimage)
1061  real(dp),intent(in) :: rprimd(3,3,nimage)
1062  real(dp),intent(in) :: fcart(3,natom,nimage)
1063  real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage)
1064 !Local variables-------------------------------
1065 !scalars
1066  integer,save :: istep_rk=0
1067  integer :: iatom,idynimage,iimage
1068  real(dp) :: stepsize
1069  character(len=500) :: msg
1070 !arrays
1071  real(dp),allocatable :: xred_old(:,:),xstep(:,:)
1072 
1073 !************************************************************************
1074 
1075 !Step for RK4 algorithm
1076  istep_rk=mod(itime,4)
1077 
1078 !Store data according to Runge-Kutta algo step
1079  if (istep_rk==1) then
1080    if (allocated(mep_param%rk4_xcart1)) then
1081      ABI_DEALLOCATE(mep_param%rk4_xcart1)
1082    end if
1083    if (allocated(mep_param%rk4_fcart1)) then
1084      ABI_DEALLOCATE(mep_param%rk4_fcart1)
1085    end if
1086    ABI_ALLOCATE(mep_param%rk4_xcart1,(3,natom,nimage))
1087    ABI_ALLOCATE(mep_param%rk4_fcart1,(3,natom,nimage))
1088    mep_param%rk4_xcart1 = xcart
1089    mep_param%rk4_fcart1 = fcart
1090  else if (istep_rk==2) then
1091    if (allocated(mep_param%rk4_fcart2)) then
1092      ABI_DEALLOCATE(mep_param%rk4_fcart2)
1093    end if
1094    ABI_ALLOCATE(mep_param%rk4_fcart2,(3,natom,nimage))
1095    mep_param%rk4_fcart2 = fcart
1096  else if (istep_rk==3) then
1097    if (allocated(mep_param%rk4_fcart3)) then
1098      ABI_DEALLOCATE(mep_param%rk4_fcart3)
1099    end if
1100    ABI_ALLOCATE(mep_param%rk4_fcart3,(3,natom,nimage))
1101    mep_param%rk4_fcart3 = fcart
1102  end if
1103 
1104  ABI_ALLOCATE(xred_old,(3,natom))
1105  if (istep_rk==0) then
1106    ABI_ALLOCATE(xstep,(3,natom))
1107  end if
1108 
1109  do idynimage=1,ndynimage
1110    iimage=list_dynimage(idynimage)
1111    xred_old(:,:)=xred(:,:,iimage)
1112 
1113 !  Note that one uses fcart, for which the sum of forces on all atoms vanish
1114 
1115 !  Intermediate Runge-Kutta step 1
1116    if      (istep_rk==1) then
1117      xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage) &
1118 &       -half*mep_param%fxcartfactor*fcart(:,:,iimage)
1119 
1120 !  Intermediate Runge-Kutta step 2
1121    else if (istep_rk==2) then
1122      xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage) &
1123 &       -half*mep_param%fxcartfactor*fcart(:,:,iimage)
1124 
1125 !  Intermediate Runge-Kutta step 3
1126    else if (istep_rk==3) then
1127      xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage) &
1128 &       -mep_param%fxcartfactor*fcart(:,:,iimage)
1129 
1130 !  Final Runge-Kutta step
1131    else if (istep_rk==0) then
1132 !    Compute image step
1133      xstep(:,:)=third*mep_param%fxcartfactor &
1134 &      *(half*fcart(:,:,iimage) &
1135 &       +half*mep_param%rk4_fcart1(:,:,iimage) &
1136 &       +mep_param%rk4_fcart2(:,:,iimage) &
1137 &       +mep_param%rk4_fcart3(:,:,iimage))
1138      stepsize=mep_img_norm(xstep)
1139      if (stepsize>=mep_param%mep_mxstep) then
1140        xstep=xstep*mep_param%mep_mxstep/stepsize
1141        write(msg,'(a,i3,a)') " Restricting step size of image ",iimage,"."
1142        call wrtout(std_out,msg,'COLL')
1143        call wrtout(ab_out ,msg,'COLL')
1144      end if
1145 !    Update positions
1146      xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage)+xstep(:,:)
1147    end if
1148 
1149    call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage))
1150 
1151 !  In case atom is fixed, we restore its previous value
1152    do iatom=1,natom
1153      if (any(mep_param%iatfix(:,iatom)==1)) then
1154        where(mep_param%iatfix(:,iatom)==1)
1155          xred(:,iatom,iimage)=xred_old(:,iatom)
1156        end where
1157        call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage))
1158      end if
1159    end do
1160 
1161  end do
1162 
1163  ABI_DEALLOCATE(xred_old)
1164  if (istep_rk==0) then
1165    ABI_DEALLOCATE(xstep)
1166  end if
1167 
1168 !Cancel storage when final RK step has been done
1169  if (istep_rk==0) then
1170    if (allocated(mep_param%rk4_xcart1)) then
1171      ABI_DEALLOCATE(mep_param%rk4_xcart1)
1172    end if
1173    if (allocated(mep_param%rk4_fcart1)) then
1174      ABI_DEALLOCATE(mep_param%rk4_fcart1)
1175    end if
1176    if (allocated(mep_param%rk4_fcart2)) then
1177      ABI_DEALLOCATE(mep_param%rk4_fcart2)
1178    end if
1179    if (allocated(mep_param%rk4_fcart3)) then
1180      ABI_DEALLOCATE(mep_param%rk4_fcart3)
1181    end if
1182  end if
1183 
1184 end subroutine mep_rk4

m_mep/mep_steepest [ Functions ]

[ Top ] [ m_mep ] [ Functions ]

NAME

  mep_steepest

FUNCTION

  Make a path (string of images) evolve according to a steepest descent algorithm

INPUTS

  fcart(3,natom,nimage)=cartesian forces in each image along the path
  list_dynimage(nimage)=list of dynamical images.
  mep_param=datastructure of type mep_type.
            several parameters for Minimal Energy Path (MEP) search.
  natom=number of atoms
  ndynimage=number of dynamical images along the path
  nimage=number of images (including static ones)
  results_img(nimage)=datastructure that hold data for each image
                      (positions, forces, energy, ...)
  rprimd(3,3,nimage)=dimensional primitive translations for each image along the path

OUTPUT

SIDE EFFECTS

  xcart(3,natom,nimage)=cartesian coordinates of atoms in each image along the path
                        before and after time evolution
  xred(3,natom,nimage)=reduced coordinates of atoms in each image along the path
                       before and after time evolution

PARENTS

      predict_neb,predict_steepest,predict_string

CHILDREN

      wrtout,xcart2xred,xred2xcart

SOURCE

276 subroutine mep_steepest(fcart,list_dynimage,mep_param,natom,ndynimage,nimage,rprimd,xcart,xred)
277 
278 
279 !This section has been created automatically by the script Abilint (TD).
280 !Do not modify the following lines by hand.
281 #undef ABI_FUNC
282 #define ABI_FUNC 'mep_steepest'
283 !End of the abilint section
284 
285  implicit none
286 
287 !Arguments ------------------------------------
288 !scalars
289  integer,intent(in) :: natom,ndynimage,nimage
290  type(mep_type),intent(in) :: mep_param
291 !arrays
292  integer,intent(in) :: list_dynimage(ndynimage)
293  real(dp),intent(in) :: fcart(3,natom,nimage),rprimd(3,3,nimage)
294  real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage)
295 !Local variables-------------------------------
296 !scalars
297  integer :: iatom,idynimage,iimage
298  real(dp) :: stepsize
299  character(len=500) :: msg
300 !arrays
301  real(dp),allocatable :: xred_old(:,:),xstep(:,:)
302 
303 !************************************************************************
304 
305  ABI_ALLOCATE(xred_old,(3,natom))
306  ABI_ALLOCATE(xstep,(3,natom))
307 
308  do idynimage=1,ndynimage
309    iimage=list_dynimage(idynimage)
310    xred_old(:,:)=xred(:,:,iimage)
311 
312 !  Compute image step
313 !  Note that one uses fcart, for which the sum of forces on all atoms vanish
314    xstep(:,:)=mep_param%fxcartfactor*fcart(:,:,iimage)
315    stepsize=mep_img_norm(xstep)
316    if (stepsize>=mep_param%mep_mxstep) then
317      xstep=xstep*mep_param%mep_mxstep/stepsize
318      write(msg,'(a,i3,a)') " Restricting step size of image ",iimage,"."
319      call wrtout(std_out,msg,'COLL')
320      call wrtout(ab_out ,msg,'COLL')
321    end if
322 
323 !  Update positions
324    xcart(:,:,iimage)=xcart(:,:,iimage)+xstep(:,:)
325    call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage))
326 
327 !  In case atom is fixed, we restore its previous value
328    do iatom=1,natom
329      if (any(mep_param%iatfix(:,iatom)==1)) then
330        where(mep_param%iatfix(:,iatom)==1)
331          xred(:,iatom,iimage)=xred_old(:,iatom)
332        end where
333        call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage))
334      end if
335    end do
336 
337  end do
338 
339  ABI_DEALLOCATE(xred_old)
340  ABI_DEALLOCATE(xstep)
341 
342 end subroutine mep_steepest

m_mep/mep_type [ Types ]

[ Top ] [ m_mep ] [ Types ]

NAME

 mep_type

FUNCTION

 Datatype with the variables required to perform MEP search

SOURCE

71  type,public :: mep_type
72 ! Scalars
73   integer  :: cineb_start  ! Starting iteration for the CI-NEB
74   integer  :: mep_solver   ! Selection of a solver for the ODE
75   integer  :: neb_algo     ! Selection of the variant of the NEB method
76   integer  :: string_algo  ! Selection of the variant of the String Method
77   real(dp) :: fxcartfactor ! Time step for steepest descent or RK4
78   real(dp) :: mep_mxstep   ! Selection of a max. step size for the ODE
79 ! Arrays
80   integer,pointer      :: iatfix(:,:)=>null() ! Atoms to fix (this pointer is associated with dtset%iatfix)
81   real(dp)             :: neb_spring(2)       ! Spring constants for the NEB method
82   real(dp),allocatable :: bfgs_xprev(:,:,:)   ! BFGS storage (prev positions)
83   real(dp),allocatable :: bfgs_fprev(:,:,:)   ! BFGS storage (prev forces)
84   real(dp),allocatable :: gbfgs_hess(:,:)     ! global-BFGS storage (Hessian matrix)
85   real(dp),allocatable :: lbfgs_hess(:,:,:)   ! local-BFGS storage (Hessian matrix)
86   real(dp),allocatable :: qmin_vel(:,:,:)     ! Quick-min algo storage (velocities)
87   real(dp),allocatable :: rk4_xcart1(:,:,:)   ! 4th-order Runge-Kutta storage
88   real(dp),allocatable :: rk4_fcart1(:,:,:)   ! 4th-order Runge-Kutta storage
89   real(dp),allocatable :: rk4_fcart2(:,:,:)   ! 4th-order Runge-Kutta storage
90   real(dp),allocatable :: rk4_fcart3(:,:,:)   ! 4th-order Runge-Kutta storage
91  end type mep_type