TABLE OF CONTENTS
- ABINIT/m_mep
- m_mep/mep_destroy
- m_mep/mep_gbfgs
- m_mep/mep_img_dotp
- m_mep/mep_img_dotp_red
- m_mep/mep_img_norm
- m_mep/mep_img_norm_red
- m_mep/mep_init
- m_mep/mep_lbfgs
- m_mep/mep_qmin
- m_mep/mep_rk4
- m_mep/mep_steepest
- m_mep/mep_type
ABINIT/m_mep [ Modules ]
NAME
m_mep
FUNCTION
This module provides several routines and datatypes for the Minimal Energy Path (MEP) search implementation.
COPYRIGHT
Copyright (C) 2012-2024 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
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 MODULE m_mep 28 29 use defs_basis 30 use m_abicore 31 use m_errors 32 use m_dtset 33 use m_xmpi 34 35 use defs_abitypes, only : MPI_type 36 use m_geometry, only : gred2fcart, fcart2gred, xcart2xred, xred2xcart, metric 37 use m_bfgs, only : hessupdt 38 use m_results_img, only : results_img_type, gather_array_img 39 40 implicit none 41 42 private 43 44 !public procedures 45 public :: mep_init 46 public :: mep_destroy 47 public :: mep_steepest 48 public :: mep_qmin 49 public :: mep_lbfgs 50 public :: mep_gbfgs 51 public :: mep_rk4 52 public :: mep_img_dotp 53 public :: mep_img_norm 54 public :: mep_img_dotp_red 55 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.
SOURCE
164 subroutine mep_destroy(mep_param) 165 166 !Arguments ------------------------------------ 167 !scalars 168 type(mep_type),intent(inout) :: mep_param 169 170 !************************************************************************ 171 172 ABI_SFREE(mep_param%bfgs_xprev) 173 ABI_SFREE(mep_param%gbfgs_hess) 174 ABI_SFREE(mep_param%bfgs_fprev) 175 ABI_SFREE(mep_param%lbfgs_hess) 176 ABI_SFREE(mep_param%qmin_vel) 177 ABI_SFREE(mep_param%rk4_xcart1) 178 ABI_SFREE(mep_param%rk4_fcart1) 179 ABI_SFREE(mep_param%rk4_fcart2) 180 ABI_SFREE(mep_param%rk4_fcart3) 181 182 nullify(mep_param%iatfix) 183 184 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
SOURCE
609 subroutine mep_gbfgs(fcart,itime,list_dynimage,mep_param,mpi_enreg,natom,& 610 & ndynimage,nimage,nimage_tot,rprimd,xcart,xred) 611 612 !Arguments ------------------------------------ 613 !scalars 614 integer,intent(in) :: itime,natom,ndynimage,nimage,nimage_tot 615 type(mep_type),intent(inout) :: mep_param 616 type(MPI_type),intent(in) :: mpi_enreg 617 !arrays 618 integer,intent(in) :: list_dynimage(ndynimage) 619 real(dp),intent(in) :: fcart(3,natom,nimage),rprimd(3,3,nimage) 620 real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage) 621 !Local variables------------------------------- 622 !scalars 623 integer :: iatom,idynimage,ii,iimage,iimage_tot,indi,indj,ierr 624 integer :: jdynimage,jatom,jj,mu,ndynimage_tot,nu 625 logical :: reset 626 real(dp),parameter :: initial_Hessian=1._dp ! in Bohr^2/Hartree 627 real(dp) :: dot1,dot2,stepsize,ucvol 628 character(len=500) :: msg 629 !arrays 630 integer,allocatable :: dynimage_tot(:),iatfix_fake(:,:),ind_dynimage_tot(:) 631 integer,allocatable :: list_dynimage_tot(:) 632 real(dp) :: favg(3),gprimd(3,3),gmet(3,3),rmet(3,3) 633 real(dp),allocatable :: buffer(:,:),buffer_all(:,:),gred(:,:) 634 real(dp),allocatable :: fcart_all(:,:,:),fcartp_all(:,:,:) 635 real(dp),allocatable :: gmet_all(:,:,:),gprimd_all(:,:,:),rprimd_all(:,:,:) 636 real(dp),allocatable :: xcart_all(:,:,:),xcartp_all(:,:,:),xred_old(:,:),xstep_all(:,:,:) 637 638 !************************************************************************ 639 640 !Retrieve indexes of all dynamical images 641 ABI_MALLOC(ind_dynimage_tot,(nimage_tot)) 642 if (mpi_enreg%paral_img==1) then 643 ABI_MALLOC(dynimage_tot,(nimage_tot)) 644 dynimage_tot=0 645 do idynimage=1,ndynimage 646 iimage=list_dynimage(idynimage) 647 iimage_tot=mpi_enreg%my_imgtab(iimage) 648 dynimage_tot(iimage_tot)=1 649 end do 650 call xmpi_sum(dynimage_tot,mpi_enreg%comm_img,ierr) 651 ndynimage_tot=count(dynimage_tot(:)>0) 652 ABI_MALLOC(list_dynimage_tot,(ndynimage_tot)) 653 idynimage=0;ind_dynimage_tot(:)=-1 654 do iimage_tot=1,nimage_tot 655 if (dynimage_tot(iimage_tot)>0) then 656 idynimage=idynimage+1 657 ind_dynimage_tot(iimage_tot)=idynimage 658 list_dynimage_tot(idynimage)=iimage_tot 659 end if 660 end do 661 ABI_FREE(dynimage_tot) 662 else 663 ndynimage_tot=ndynimage 664 ABI_MALLOC(list_dynimage_tot,(ndynimage)) 665 ind_dynimage_tot(:)=-1 666 do idynimage=1,ndynimage 667 ind_dynimage_tot(list_dynimage(idynimage))=idynimage 668 list_dynimage_tot(idynimage)=list_dynimage(idynimage) 669 end do 670 end if 671 672 !Allocate history array (at first time step) 673 if (itime==1) then 674 if (allocated(mep_param%bfgs_xprev)) then 675 ABI_FREE(mep_param%bfgs_xprev) 676 end if 677 if (allocated(mep_param%bfgs_fprev)) then 678 ABI_FREE(mep_param%bfgs_fprev) 679 end if 680 if (allocated(mep_param%gbfgs_hess)) then 681 ABI_FREE(mep_param%gbfgs_hess) 682 end if 683 ABI_MALLOC(mep_param%bfgs_xprev,(3,natom,ndynimage)) 684 ABI_MALLOC(mep_param%bfgs_fprev,(3,natom,ndynimage)) 685 ABI_MALLOC(mep_param%gbfgs_hess,(3*natom*ndynimage_tot,3*natom*ndynimage_tot)) 686 mep_param%bfgs_xprev=zero 687 mep_param%bfgs_fprev=zero 688 end if 689 690 !Retrieve positions and forces for all images 691 ABI_MALLOC(xcart_all,(3,natom,ndynimage_tot)) 692 ABI_MALLOC(fcart_all,(3,natom,ndynimage_tot)) 693 ABI_MALLOC(xcartp_all,(3,natom,ndynimage_tot)) 694 ABI_MALLOC(fcartp_all,(3,natom,ndynimage_tot)) 695 ABI_MALLOC(rprimd_all,(3,3,ndynimage_tot)) 696 ABI_MALLOC(gprimd_all,(3,3,ndynimage_tot)) 697 ABI_MALLOC(gmet_all,(3,3,ndynimage_tot)) 698 if (mpi_enreg%paral_img==1) then 699 ABI_MALLOC(buffer,(12*natom+27,nimage)) 700 ABI_MALLOC(buffer_all,(12*natom+27,nimage_tot)) 701 buffer=zero 702 do idynimage=1,ndynimage 703 iimage=list_dynimage(idynimage) 704 call metric(gmet,gprimd,-1,rmet,rprimd(:,:,iimage),ucvol) 705 buffer( 1:3 *natom ,iimage)=reshape(xcart(:,:,iimage),(/3*natom/)) 706 buffer(3 *natom+1:6 *natom ,iimage)=reshape(fcart(:,:,iimage),(/3*natom/)) 707 buffer(6 *natom+1:9 *natom ,iimage)=reshape(mep_param%bfgs_xprev(:,:,idynimage),(/3*natom/)) 708 buffer(9 *natom+1:12*natom ,iimage)=reshape(mep_param%bfgs_fprev(:,:,idynimage),(/3*natom/)) 709 buffer(12*natom+1:12*natom+9,iimage)=reshape(rprimd(:,:,iimage),(/9/)) 710 buffer(12*natom+10:12*natom+18,iimage)=reshape(gprimd(:,:),(/9/)) 711 buffer(12*natom+19:12*natom+27,iimage)=reshape(gmet(:,:),(/9/)) 712 end do 713 call gather_array_img(buffer,buffer_all,mpi_enreg,allgather=.true.) 714 do idynimage=1,ndynimage_tot 715 iimage_tot=list_dynimage_tot(idynimage) 716 do iatom=1,natom 717 indi=12*(iatom-1) 718 do ii=1,3 719 xcart_all (ii,iatom,idynimage)= buffer_all(indi +ii,iimage_tot) 720 fcart_all (ii,iatom,idynimage)=-buffer_all(indi+3+ii,iimage_tot) ! use fcart=-cartesian_force 721 xcartp_all(ii,iatom,idynimage)= buffer_all(indi+6+ii,iimage_tot) 722 fcartp_all(ii,iatom,idynimage)= buffer_all(indi+9+ii,iimage_tot) 723 end do 724 end do 725 indi=12*natom 726 rprimd_all(1:3,1:3,idynimage)=reshape(buffer_all(indi+ 1:indi+ 9,iimage_tot),(/3,3/)) 727 gprimd_all(1:3,1:3,idynimage)=reshape(buffer_all(indi+10:indi+18,iimage_tot),(/3,3/)) 728 gmet_all (1:3,1:3,idynimage)=reshape(buffer_all(indi+19:indi+27,iimage_tot),(/3,3/)) 729 end do 730 ABI_FREE(buffer) 731 ABI_FREE(buffer_all) 732 else 733 do idynimage=1,ndynimage 734 iimage=list_dynimage(idynimage) 735 xcart_all(:,:,idynimage)= xcart(:,:,iimage) 736 fcart_all(:,:,idynimage)=-fcart(:,:,iimage) ! use fcart=-cartesian_force 737 xcartp_all(:,:,idynimage)=mep_param%bfgs_xprev(:,:,idynimage) 738 fcartp_all(:,:,idynimage)=mep_param%bfgs_fprev(:,:,idynimage) 739 rprimd_all(:,:,idynimage)=rprimd(:,:,iimage) 740 call metric(gmet_all(:,:,idynimage),gprimd_all(:,:,idynimage),-1,rmet,rprimd(:,:,iimage),ucvol) 741 end do 742 end if 743 744 !Test if a reset is needed 745 reset=.false. 746 if (itime>1) then 747 dot1=zero;dot2=zero 748 do idynimage=1,ndynimage_tot 749 dot1=dot2+mep_img_dotp(fcartp_all(:,:,idynimage),fcart_all (:,:,idynimage)) 750 dot2=dot1+mep_img_dotp(fcartp_all(:,:,idynimage),fcartp_all(:,:,idynimage)) 751 end do 752 reset=((dot2<two*abs(dot1)).or.abs(dot2)<tol8) 753 if (reset) then 754 msg=' Resetting Hessian matrix.' 755 call wrtout(std_out,msg,'COLL') 756 call wrtout(ab_out ,msg,'COLL') 757 end if 758 end if 759 760 !===> First step or reset: initialize the Hessian matrix 761 if (itime==1.or.reset) then 762 mep_param%gbfgs_hess(:,:)=zero 763 do idynimage=1,ndynimage_tot 764 indi=3*natom*(idynimage-1) 765 do iatom=1,natom 766 do mu=1,3 767 do nu=1,3 768 do ii=1,3 769 do jj=1,3 770 if (mep_param%iatfix(ii,iatom)==0.and. & 771 & mep_param%iatfix(jj,iatom)==0) then 772 mep_param%gbfgs_hess(indi+mu,indi+nu)=mep_param%gbfgs_hess(indi+mu,indi+nu) & 773 & +rprimd_all(mu,ii,idynimage)*rprimd_all(nu,jj,idynimage) & 774 & *gmet_all(ii,jj,idynimage)*initial_Hessian 775 end if 776 end do 777 end do 778 end do 779 end do 780 indi=indi+3 781 end do 782 end do 783 784 !===> Other steps: update the Hessian matrix 785 else 786 787 ! Impose here f-fprev=0 (cannot be done inside hessupdt in cartesian coordinates) 788 ABI_MALLOC(gred,(3,natom)) 789 do idynimage=1,ndynimage_tot 790 fcartp_all(:,:,idynimage)=fcartp_all(:,:,idynimage)-fcart_all(:,:,idynimage) 791 call fcart2gred(fcartp_all(:,:,idynimage),gred,rprimd_all(:,:,idynimage),natom) 792 where (mep_param%iatfix(:,:)==1) ! iatfix is defined in reduced coordinates 793 gred(:,:)=zero 794 end where 795 call gred2fcart(favg,.TRUE.,fcartp_all(:,:,idynimage),gred,gprimd_all(:,:,idynimage),natom) 796 do iatom=1,natom 797 fcartp_all(:,iatom,idynimage)=fcartp_all(:,iatom,idynimage) & 798 & +fcart_all(:,iatom,idynimage)+favg(:) 799 end do 800 end do 801 ABI_FREE(gred) 802 803 ! f-fprev=0 has already been imposed for fixed atoms: 804 ! we call hessupdt with no fixed atom 805 ABI_MALLOC(iatfix_fake,(3,natom)) 806 iatfix_fake(:,:)=0 807 call hessupdt(mep_param%gbfgs_hess,& 808 & iatfix_fake,natom,3*natom*ndynimage_tot, & 809 xcart_all,xcartp_all,fcart_all,fcartp_all, & 810 & nimage=ndynimage_tot) 811 ABI_FREE(iatfix_fake) 812 end if 813 814 !Free memory 815 ABI_FREE(xcart_all) 816 ABI_FREE(xcartp_all) 817 ABI_FREE(fcartp_all) 818 ABI_FREE(rprimd_all) 819 ABI_FREE(gprimd_all) 820 ABI_FREE(gmet_all) 821 822 !Update history 823 do idynimage=1,ndynimage 824 iimage=list_dynimage(idynimage) 825 mep_param%bfgs_xprev(:,:,idynimage)=xcart(:,:,iimage) 826 mep_param%bfgs_fprev(:,:,idynimage)=fcart(:,:,iimage) 827 end do 828 829 !Compute image step 830 ABI_MALLOC(xstep_all,(3,natom,ndynimage_tot)) 831 xstep_all=zero 832 do idynimage=1,ndynimage_tot 833 indi=3*natom*(idynimage-1) 834 do iatom=1,natom 835 do ii=1,3 836 do jdynimage=1,ndynimage_tot 837 indj=3*natom*(jdynimage-1) 838 do jatom=1,natom 839 do jj=1,3 840 ! Be careful: minus sign because fcart=-cartesian_force 841 xstep_all(ii,iatom,idynimage)=xstep_all(ii,iatom,idynimage) & 842 & -fcart_all(jj,jatom,jdynimage) & 843 & *mep_param%gbfgs_hess(indi+ii,indj+jj) 844 end do 845 indj=indj+3 846 end do 847 end do 848 end do 849 indi=indi+3 850 end do 851 end do 852 853 !Restrict image step size 854 stepsize=zero 855 do idynimage=1,ndynimage_tot 856 stepsize=stepsize+mep_img_dotp(xstep_all(:,:,idynimage),xstep_all(:,:,idynimage)) 857 end do 858 stepsize=sqrt(stepsize) 859 if (stepsize>=mep_param%mep_mxstep*dble(ndynimage_tot)) then 860 xstep_all=xstep_all*mep_param%mep_mxstep*dble(ndynimage_tot)/stepsize 861 write(msg,'(a,i3,a)') " Restricting BFGS step size." 862 call wrtout(std_out,msg,'COLL') 863 call wrtout(ab_out ,msg,'COLL') 864 end if 865 866 !Update positions 867 ABI_MALLOC(xred_old,(3,natom)) 868 do idynimage=1,ndynimage 869 iimage=list_dynimage(idynimage) 870 iimage_tot=mpi_enreg%my_imgtab(iimage) 871 xred_old(:,:)=xred(:,:,iimage) 872 xcart(:,:,iimage)=xcart(:,:,iimage)+xstep_all(:,:,ind_dynimage_tot(iimage_tot)) 873 call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage)) 874 ! In case atom is fixed, we restore its previous value 875 do iatom=1,natom 876 if (any(mep_param%iatfix(:,iatom)==1)) then 877 where(mep_param%iatfix(:,iatom)==1) 878 xred(:,iatom,iimage)=xred_old(:,iatom) 879 end where 880 call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage)) 881 end if 882 end do 883 end do 884 ABI_FREE(xred_old) 885 886 !Free memory 887 ABI_FREE(fcart_all) 888 ABI_FREE(xstep_all) 889 ABI_FREE(ind_dynimage_tot) 890 ABI_FREE(list_dynimage_tot) 891 892 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
SOURCE
1080 function mep_img_dotp(vect1,vect2) 1081 1082 !Arguments ------------------------------------ 1083 !scalars 1084 real(dp) :: mep_img_dotp 1085 !arrays 1086 real(dp),intent(in) :: vect1(:,:),vect2(:,:) 1087 !Local variables------------------------------- 1088 !scalars 1089 integer :: size1,size2 1090 !arrays 1091 1092 !************************************************************************ 1093 1094 size1=size(vect1,1);size2=size(vect1,2) 1095 if (size1/=size(vect2,1).or.size2/=size(vect2,2)) then 1096 ABI_BUG("Error on dimensions !") 1097 end if 1098 1099 mep_img_dotp=sum(vect1*vect2) 1100 1101 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
SOURCE
1159 function mep_img_dotp_red(rmet,vect1,vect2) 1160 1161 !Arguments ------------------------------------ 1162 !scalars 1163 real(dp) :: mep_img_dotp_red 1164 !arrays 1165 real(dp),intent(in) :: rmet(3,3) 1166 real(dp),intent(in) :: vect1(:,:),vect2(:,:) 1167 !Local variables------------------------------- 1168 !scalars 1169 integer :: iatom,ii,jj,size1,size2 1170 1171 !************************************************************************ 1172 1173 size1=size(vect1,1);size2=size(vect1,2) 1174 if (size1/=size(vect2,1).or.size2/=size(vect2,2).or.size1/=3) then 1175 ABI_BUG("Error on dimensions !") 1176 end if 1177 1178 mep_img_dotp_red=zero 1179 do iatom=1,size2 1180 do ii=1,3 1181 do jj=1,3 1182 mep_img_dotp_red=mep_img_dotp_red+vect1(ii,iatom)*vect2(jj,iatom)*rmet(ii,jj) 1183 end do 1184 end do 1185 end do 1186 1187 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
SOURCE
1123 function mep_img_norm(vect) 1124 1125 !Arguments ------------------------------------ 1126 !scalars 1127 real(dp) :: mep_img_norm 1128 !arrays 1129 real(dp),intent(in) :: vect(:,:) 1130 1131 !************************************************************************ 1132 1133 mep_img_norm=sqrt(sum(vect*vect)) 1134 1135 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
SOURCE
1211 function mep_img_norm_red(rmet,vect) 1212 1213 !Arguments ------------------------------------ 1214 !scalars 1215 real(dp) :: mep_img_norm_red 1216 !arrays 1217 real(dp),intent(in) :: rmet(3,3) 1218 real(dp),intent(in) :: vect(:,:) 1219 1220 !************************************************************************ 1221 1222 mep_img_norm_red=sqrt(mep_img_dotp_red(rmet,vect,vect)) 1223 1224 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.
SOURCE
113 subroutine mep_init(dtset,mep_param) 114 115 !Arguments ------------------------------------ 116 !scalars 117 type(dataset_type),target,intent(in) :: dtset 118 type(mep_type),intent(inout) :: mep_param 119 120 !************************************************************************ 121 122 if((dtset%imgmov==1).or.(dtset%imgmov==2).or.(dtset%imgmov==5))then 123 mep_param%cineb_start = dtset%cineb_start 124 mep_param%mep_solver = dtset%mep_solver 125 mep_param%neb_algo = dtset%neb_algo 126 mep_param%string_algo = dtset%string_algo 127 mep_param%fxcartfactor = dtset%fxcartfactor 128 mep_param%mep_mxstep = dtset%mep_mxstep 129 mep_param%neb_spring = dtset%neb_spring 130 mep_param%iatfix =>dtset%iatfix 131 else 132 mep_param%cineb_start = -1 133 mep_param%mep_solver = -1 134 mep_param%neb_algo = -1 135 mep_param%string_algo = -1 136 mep_param%fxcartfactor = zero 137 mep_param%mep_mxstep = 100._dp 138 mep_param%neb_spring = zero 139 nullify(mep_param%iatfix) 140 end if 141 142 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.
SOURCE
431 subroutine mep_lbfgs(fcart,itime,list_dynimage,mep_param,natom,ndynimage,& 432 & nimage,rprimd,xcart,xred) 433 434 !Arguments ------------------------------------ 435 !scalars 436 integer,intent(in) :: itime,natom,ndynimage,nimage 437 type(mep_type),intent(inout) :: mep_param 438 !arrays 439 integer,intent(in) :: list_dynimage(ndynimage) 440 real(dp),intent(in) :: rprimd(3,3,nimage) 441 real(dp),intent(in) :: fcart(3,natom,nimage) 442 real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage) 443 !Local variables------------------------------- 444 !scalars 445 integer :: iatom,idynimage,ii,iimage,indi,indj,jatom,jj 446 logical :: reset 447 real(dp),parameter :: initial_Hessian=1._dp ! in Bohr^2/Hartree 448 real(dp) :: dot1,dot2,stepsize,ucvol 449 character(len=500) :: msg 450 !arrays 451 real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3) 452 real(dp),allocatable :: gred(:,:),xstep(:,:) 453 454 !************************************************************************ 455 456 !Allocate history array (at first time step) 457 if (itime==1) then 458 if (allocated(mep_param%bfgs_xprev)) then 459 ABI_FREE(mep_param%bfgs_xprev) 460 end if 461 if (allocated(mep_param%bfgs_fprev)) then 462 ABI_FREE(mep_param%bfgs_fprev) 463 end if 464 if (allocated(mep_param%lbfgs_hess)) then 465 ABI_FREE(mep_param%lbfgs_hess) 466 end if 467 ABI_MALLOC(mep_param%bfgs_xprev,(3,natom,ndynimage)) 468 ABI_MALLOC(mep_param%bfgs_fprev,(3,natom,ndynimage)) 469 ABI_MALLOC(mep_param%lbfgs_hess,(3*natom,3*natom,ndynimage)) 470 mep_param%bfgs_xprev=zero 471 mep_param%bfgs_fprev=zero 472 end if 473 474 !Temporary storage 475 ABI_MALLOC(gred,(3,natom)) 476 ABI_MALLOC(xstep,(3,natom)) 477 478 !Loop over images 479 do idynimage=1,ndynimage 480 iimage=list_dynimage(idynimage) 481 call metric(gmet,gprimd,-1,rmet,rprimd(:,:,iimage),ucvol) 482 call fcart2gred(fcart(:,:,iimage),gred,rprimd(:,:,iimage),natom) 483 484 ! Test if a reset is needed 485 reset=.false. 486 if (itime>1) then 487 dot1=mep_img_dotp(mep_param%bfgs_fprev(:,:,idynimage),gred) 488 dot2=mep_img_dotp(mep_param%bfgs_fprev(:,:,idynimage), & 489 & mep_param%bfgs_fprev(:,:,idynimage)) 490 ! dot1=mep_img_dotp_red(rmet,mep_param%bfgs_fprev(:,:,idynimage),gred) 491 ! dot2=mep_img_dotp_red(rmet,mep_param%bfgs_fprev(:,:,idynimage), & 492 !& mep_param%bfgs_fprev(:,:,idynimage)) 493 reset=((dot2<two*abs(dot1)).or.abs(dot2)<tol8) 494 if (reset) then 495 write(msg,'(a,i3,a)') " Resetting Hessian matrix for image ",iimage,"." 496 call wrtout(std_out,msg,'COLL') 497 call wrtout(ab_out ,msg,'COLL') 498 end if 499 end if 500 501 ! ===> First step or reset: initialize the Hessian matrix (in reduced coordinates) 502 if (itime==1.or.reset) then 503 mep_param%lbfgs_hess(:,:,idynimage)=zero 504 do iatom=1,natom 505 indi=3*(iatom-1) 506 do ii=1,3 507 do jj=1,3 508 if (mep_param%iatfix(ii,iatom)==0.and. & 509 & mep_param%iatfix(jj,iatom)==0) then 510 mep_param%lbfgs_hess(indi+ii,indi+jj,idynimage)=gmet(ii,jj)*initial_Hessian 511 end if 512 end do 513 end do 514 end do 515 516 ! ===> Other steps: update the Hessian matrix 517 else 518 call hessupdt(mep_param%lbfgs_hess(:,:,idynimage),& 519 & mep_param%iatfix,natom,3*natom, & 520 xred(:,:,iimage),mep_param%bfgs_xprev(:,:,idynimage),& 521 gred(:,:),mep_param%bfgs_fprev(:,:,idynimage)) 522 end if 523 524 ! Update history 525 mep_param%bfgs_xprev(:,:,idynimage)=xred(:,:,iimage) 526 mep_param%bfgs_fprev(:,:,idynimage)=gred(:,:) 527 528 ! Compute image step 529 xstep=zero 530 do iatom=1,natom 531 indi=3*(iatom-1) 532 do ii=1,3 533 do jatom=1,natom 534 indj=3*(jatom-1) 535 do jj=1,3 536 xstep(ii,iatom)=xstep(ii,iatom) & 537 & -mep_param%lbfgs_hess(indi+ii,indj+jj,idynimage)*gred(jj,jatom) 538 end do 539 end do 540 end do 541 end do 542 543 ! Restrict image step size 544 stepsize=mep_img_norm_red(rmet,xstep) 545 if (stepsize>=mep_param%mep_mxstep) then 546 xstep=xstep*mep_param%mep_mxstep/stepsize 547 write(msg,'(a,i3,a)') " Restricting BFGS step size of image ",iimage,"." 548 call wrtout(std_out,msg,'COLL') 549 call wrtout(ab_out ,msg,'COLL') 550 end if 551 552 ! Update positions 553 xred(:,:,iimage)=xred(:,:,iimage)+xstep(:,:) 554 555 ! In case atom is fixed, we restore its previous value 556 where(mep_param%iatfix(:,:)==1) 557 xred(:,:,iimage)=mep_param%bfgs_xprev(:,:,idynimage) 558 end where 559 560 call xred2xcart(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage)) 561 562 !End loop over images 563 end do 564 565 ABI_FREE(gred) 566 ABI_FREE(xstep) 567 568 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
SOURCE
310 subroutine mep_qmin(fcart,itime,list_dynimage,mep_param,natom,ndynimage,nimage,rprimd,xcart,xred) 311 312 !Arguments ------------------------------------ 313 !scalars 314 integer,intent(in) :: itime,natom,ndynimage,nimage 315 type(mep_type),intent(inout) :: mep_param 316 !arrays 317 integer,intent(in) :: list_dynimage(ndynimage) 318 real(dp),intent(in) :: fcart(3,natom,nimage),rprimd(3,3,nimage) 319 real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage) 320 !Local variables------------------------------- 321 !scalars 322 integer :: iatom,idynimage,iimage 323 real(dp) :: stepsize,vdotf 324 character(len=500) :: msg 325 !arrays 326 real(dp) :: vel_red(3) 327 real(dp),allocatable :: xred_old(:,:),xstep(:,:) 328 329 !*********************************************************************** 330 331 !Allocate history array (at first time step) 332 if (itime==1) then 333 if (allocated(mep_param%qmin_vel)) then 334 ABI_FREE(mep_param%qmin_vel) 335 end if 336 ABI_MALLOC(mep_param%qmin_vel,(3,natom,ndynimage)) 337 mep_param%qmin_vel=zero 338 end if 339 340 ABI_MALLOC(xred_old,(3,natom)) 341 ABI_MALLOC(xstep,(3,natom)) 342 343 do idynimage=1,ndynimage 344 iimage=list_dynimage(idynimage) 345 xred_old(:,:)=xred(:,:,iimage) 346 347 ! Compute velocities 348 vdotf=mep_img_dotp(mep_param%qmin_vel(:,:,idynimage),fcart(:,:,iimage)) 349 if (vdotf>=zero) then 350 mep_param%qmin_vel(:,:,idynimage)=vdotf*fcart(:,:,iimage) & 351 & /mep_img_norm(fcart(:,:,iimage)) 352 else 353 mep_param%qmin_vel(:,:,idynimage)=zero 354 write(msg,'(a,i3,a)') " Setting velocities of image ",iimage," to zero." 355 call wrtout(std_out,msg,'COLL') 356 call wrtout(ab_out ,msg,'COLL') 357 end if 358 mep_param%qmin_vel(:,:,idynimage)=mep_param%qmin_vel(:,:,idynimage) & 359 & +mep_param%fxcartfactor*fcart(:,:,iimage) 360 361 ! Compute image step 362 xstep(:,:)=mep_param%fxcartfactor*mep_param%qmin_vel(:,:,idynimage) 363 stepsize=mep_img_norm(xstep) 364 if (stepsize>=mep_param%mep_mxstep) then 365 xstep=xstep*mep_param%mep_mxstep/stepsize 366 write(msg,'(a,i3,a)') " Restricting step size of image ",iimage,"." 367 call wrtout(std_out,msg,'COLL') 368 call wrtout(ab_out ,msg,'COLL') 369 end if 370 371 ! Update positions 372 xcart(:,:,iimage)=xcart(:,:,iimage)+xstep(:,:) 373 call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage)) 374 375 ! In case atom is fixed, we restore its previous value 376 do iatom=1,natom 377 if (any(mep_param%iatfix(:,iatom)==1)) then 378 call xcart2xred(1,rprimd(:,:,iimage),mep_param%qmin_vel(:,iatom,idynimage),vel_red) 379 where(mep_param%iatfix(:,iatom)==1) 380 xred(:,iatom,iimage)=xred_old(:,iatom) 381 vel_red(:)=zero 382 end where 383 call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage)) 384 call xred2xcart(1,rprimd(:,:,iimage),mep_param%qmin_vel(:,iatom,idynimage),vel_red) 385 end if 386 end do 387 388 end do 389 390 ABI_FREE(xred_old) 391 ABI_FREE(xstep) 392 393 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
SOURCE
926 subroutine mep_rk4(fcart,itime,list_dynimage,mep_param,natom,ndynimage,nimage,rprimd,xcart,xred) 927 928 !Arguments ------------------------------------ 929 !scalars 930 integer,intent(in) :: itime,natom,ndynimage,nimage 931 type(mep_type),intent(inout) :: mep_param 932 !arrays 933 integer,intent(in) :: list_dynimage(ndynimage) 934 real(dp),intent(in) :: rprimd(3,3,nimage) 935 real(dp),intent(in) :: fcart(3,natom,nimage) 936 real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage) 937 !Local variables------------------------------- 938 !scalars 939 integer,save :: istep_rk=0 940 integer :: iatom,idynimage,iimage 941 real(dp) :: stepsize 942 character(len=500) :: msg 943 !arrays 944 real(dp),allocatable :: xred_old(:,:),xstep(:,:) 945 946 !************************************************************************ 947 948 !Step for RK4 algorithm 949 istep_rk=mod(itime,4) 950 951 !Store data according to Runge-Kutta algo step 952 if (istep_rk==1) then 953 if (allocated(mep_param%rk4_xcart1)) then 954 ABI_FREE(mep_param%rk4_xcart1) 955 end if 956 if (allocated(mep_param%rk4_fcart1)) then 957 ABI_FREE(mep_param%rk4_fcart1) 958 end if 959 ABI_MALLOC(mep_param%rk4_xcart1,(3,natom,nimage)) 960 ABI_MALLOC(mep_param%rk4_fcart1,(3,natom,nimage)) 961 mep_param%rk4_xcart1 = xcart 962 mep_param%rk4_fcart1 = fcart 963 else if (istep_rk==2) then 964 if (allocated(mep_param%rk4_fcart2)) then 965 ABI_FREE(mep_param%rk4_fcart2) 966 end if 967 ABI_MALLOC(mep_param%rk4_fcart2,(3,natom,nimage)) 968 mep_param%rk4_fcart2 = fcart 969 else if (istep_rk==3) then 970 if (allocated(mep_param%rk4_fcart3)) then 971 ABI_FREE(mep_param%rk4_fcart3) 972 end if 973 ABI_MALLOC(mep_param%rk4_fcart3,(3,natom,nimage)) 974 mep_param%rk4_fcart3 = fcart 975 end if 976 977 ABI_MALLOC(xred_old,(3,natom)) 978 if (istep_rk==0) then 979 ABI_MALLOC(xstep,(3,natom)) 980 end if 981 982 do idynimage=1,ndynimage 983 iimage=list_dynimage(idynimage) 984 xred_old(:,:)=xred(:,:,iimage) 985 986 ! Note that one uses fcart, for which the sum of forces on all atoms vanish 987 988 ! Intermediate Runge-Kutta step 1 989 if (istep_rk==1) then 990 xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage) & 991 & -half*mep_param%fxcartfactor*fcart(:,:,iimage) 992 993 ! Intermediate Runge-Kutta step 2 994 else if (istep_rk==2) then 995 xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage) & 996 & -half*mep_param%fxcartfactor*fcart(:,:,iimage) 997 998 ! Intermediate Runge-Kutta step 3 999 else if (istep_rk==3) then 1000 xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage) & 1001 & -mep_param%fxcartfactor*fcart(:,:,iimage) 1002 1003 ! Final Runge-Kutta step 1004 else if (istep_rk==0) then 1005 ! Compute image step 1006 xstep(:,:)=third*mep_param%fxcartfactor & 1007 & *(half*fcart(:,:,iimage) & 1008 & +half*mep_param%rk4_fcart1(:,:,iimage) & 1009 & +mep_param%rk4_fcart2(:,:,iimage) & 1010 & +mep_param%rk4_fcart3(:,:,iimage)) 1011 stepsize=mep_img_norm(xstep) 1012 if (stepsize>=mep_param%mep_mxstep) then 1013 xstep=xstep*mep_param%mep_mxstep/stepsize 1014 write(msg,'(a,i3,a)') " Restricting step size of image ",iimage,"." 1015 call wrtout(std_out,msg,'COLL') 1016 call wrtout(ab_out ,msg,'COLL') 1017 end if 1018 ! Update positions 1019 xcart(:,:,iimage)=mep_param%rk4_xcart1(:,:,iimage)+xstep(:,:) 1020 end if 1021 1022 call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage)) 1023 1024 ! In case atom is fixed, we restore its previous value 1025 do iatom=1,natom 1026 if (any(mep_param%iatfix(:,iatom)==1)) then 1027 where(mep_param%iatfix(:,iatom)==1) 1028 xred(:,iatom,iimage)=xred_old(:,iatom) 1029 end where 1030 call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage)) 1031 end if 1032 end do 1033 1034 end do 1035 1036 ABI_FREE(xred_old) 1037 if (istep_rk==0) then 1038 ABI_FREE(xstep) 1039 end if 1040 1041 !Cancel storage when final RK step has been done 1042 if (istep_rk==0) then 1043 if (allocated(mep_param%rk4_xcart1)) then 1044 ABI_FREE(mep_param%rk4_xcart1) 1045 end if 1046 if (allocated(mep_param%rk4_fcart1)) then 1047 ABI_FREE(mep_param%rk4_fcart1) 1048 end if 1049 if (allocated(mep_param%rk4_fcart2)) then 1050 ABI_FREE(mep_param%rk4_fcart2) 1051 end if 1052 if (allocated(mep_param%rk4_fcart3)) then 1053 ABI_FREE(mep_param%rk4_fcart3) 1054 end if 1055 end if 1056 1057 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
SOURCE
218 subroutine mep_steepest(fcart,list_dynimage,mep_param,natom,ndynimage,nimage,rprimd,xcart,xred) 219 220 !Arguments ------------------------------------ 221 !scalars 222 integer,intent(in) :: natom,ndynimage,nimage 223 type(mep_type),intent(in) :: mep_param 224 !arrays 225 integer,intent(in) :: list_dynimage(ndynimage) 226 real(dp),intent(in) :: fcart(3,natom,nimage),rprimd(3,3,nimage) 227 real(dp),intent(inout) :: xcart(3,natom,nimage),xred(3,natom,nimage) 228 !Local variables------------------------------- 229 !scalars 230 integer :: iatom,idynimage,iimage 231 real(dp) :: stepsize 232 character(len=500) :: msg 233 !arrays 234 real(dp),allocatable :: xred_old(:,:),xstep(:,:) 235 236 !************************************************************************ 237 238 ABI_MALLOC(xred_old,(3,natom)) 239 ABI_MALLOC(xstep,(3,natom)) 240 241 do idynimage=1,ndynimage 242 iimage=list_dynimage(idynimage) 243 xred_old(:,:)=xred(:,:,iimage) 244 245 ! Compute image step 246 ! Note that one uses fcart, for which the sum of forces on all atoms vanish 247 xstep(:,:)=mep_param%fxcartfactor*fcart(:,:,iimage) 248 stepsize=mep_img_norm(xstep) 249 if (stepsize>=mep_param%mep_mxstep) then 250 xstep=xstep*mep_param%mep_mxstep/stepsize 251 write(msg,'(a,i3,a)') " Restricting step size of image ",iimage,"." 252 call wrtout(std_out,msg,'COLL') 253 call wrtout(ab_out ,msg,'COLL') 254 end if 255 256 ! Update positions 257 xcart(:,:,iimage)=xcart(:,:,iimage)+xstep(:,:) 258 call xcart2xred(natom,rprimd(:,:,iimage),xcart(:,:,iimage),xred(:,:,iimage)) 259 260 ! In case atom is fixed, we restore its previous value 261 do iatom=1,natom 262 if (any(mep_param%iatfix(:,iatom)==1)) then 263 where(mep_param%iatfix(:,iatom)==1) 264 xred(:,iatom,iimage)=xred_old(:,iatom) 265 end where 266 call xred2xcart(1,rprimd(:,:,iimage),xcart(:,iatom,iimage),xred(:,iatom,iimage)) 267 end if 268 end do 269 270 end do 271 272 ABI_FREE(xred_old) 273 ABI_FREE(xstep) 274 275 end subroutine mep_steepest
m_mep/mep_type [ Types ]
NAME
mep_type
FUNCTION
Datatype with the variables required to perform MEP search
SOURCE
67 type,public :: mep_type 68 ! Scalars 69 integer :: cineb_start ! Starting iteration for the CI-NEB 70 integer :: mep_solver ! Selection of a solver for the ODE 71 integer :: neb_algo ! Selection of the variant of the NEB method 72 integer :: string_algo ! Selection of the variant of the String Method 73 real(dp) :: fxcartfactor ! Time step for steepest descent or RK4 74 real(dp) :: mep_mxstep ! Selection of a max. step size for the ODE 75 ! Arrays 76 integer,pointer :: iatfix(:,:)=>null() ! Atoms to fix (this pointer is associated with dtset%iatfix) 77 real(dp) :: neb_spring(2) ! Spring constants for the NEB method 78 real(dp),allocatable :: bfgs_xprev(:,:,:) ! BFGS storage (prev positions) 79 real(dp),allocatable :: bfgs_fprev(:,:,:) ! BFGS storage (prev forces) 80 real(dp),allocatable :: gbfgs_hess(:,:) ! global-BFGS storage (Hessian matrix) 81 real(dp),allocatable :: lbfgs_hess(:,:,:) ! local-BFGS storage (Hessian matrix) 82 real(dp),allocatable :: qmin_vel(:,:,:) ! Quick-min algo storage (velocities) 83 real(dp),allocatable :: rk4_xcart1(:,:,:) ! 4th-order Runge-Kutta storage 84 real(dp),allocatable :: rk4_fcart1(:,:,:) ! 4th-order Runge-Kutta storage 85 real(dp),allocatable :: rk4_fcart2(:,:,:) ! 4th-order Runge-Kutta storage 86 real(dp),allocatable :: rk4_fcart3(:,:,:) ! 4th-order Runge-Kutta storage 87 end type mep_type