Statistics
| Revision:

## root / src / egrad_chamfre.f90 @ 8

 1 !*************############### Attention ###############*************  !(a) Two parameters aa and br(3,3) must be provided in main program  !(b) At first, z=0 for the outmost surface  !(c) The units of all parameters are not atomic units, but angstrom, eV etc.  !(d) Before using this routine, you must confirm that the energy is CONSISTENT with its derivative.  ! This check of consistency is VERY IMPORTANT for energy conservation in MD.  !!!  !! xa, ya,za present the positions of atoms. v is the energy of the system. dvdr present the force along the x,y,z directions.  !!  subroutine egrad_chamfre(Nat,V,Geom,Grad)  ! WAS xa,ya,za,v,dvdr)   use Path_module, only : Lat_a,Lat_b,order,renum   use Io_module   IMPLICIT NONE   INTEGER(KINT), INTENT(IN) :: Nat   REAL(KREAL), INTENT(OUT) :: V,grad(Nat*3)   REAL(KREAL), INTENT(IN) :: Geom(Nat,3)  ! Bohr --> Angstr  ! real(KREAL), parameter :: BOHR = 0.52917726D+00  !! PFL 2010 Nov. 28 ->  ! for now I keep this old F77 scheme. This has to be replaced by ALLOCATABLE  ! things  ! <- PFL 2010 Nov 28.   REAL(KREAL) :: xa(100),ya(100),za(100),dvdr(100,3)   REAL(KREAL) :: rshort,rlong,zshort,zlong   REAL(KREAL), SAVE :: rnc, rnn, rhh, rch, rnh, ann(100), anc(100), anh(100), ach(100), ahh(100)   REAL(KREAL) :: dv(100,10,3),df(100,10,3)   REAL(KREAL) :: px(9000),py(9000),pz(9000),tmp,temp,rtemp,yf   INTEGER(KINT) :: ka(9000),kx(9000),ky(9000),itx,ity,nf(100),mf(100,500)   REAL(KREAL) :: xi(3), xj(3), rij, er, eb   REAL(KREAL) :: fcut, tf   INTEGER(KINT) :: np,nq,nh,nc   INTEGER(KINT) :: i, j, nn, m, n, mid, n1, im, jm, n2,n3,n4, iat   LOGICAL, SAVE :: First=.TRUE.   LOGICAL :: debug, FExist   common/zdis/zshort,zlong   common/rdistanc/rshort,rlong   INTERFACE   function valid(string) result (isValid)   CHARACTER(*), intent(in) :: string   logical :: isValid   END function VALID     END INTERFACE  !**********************************************   debug=valid('egrad_chamfre')   if (debug) Call Header('Entering Egrad_chamfre')   v=0.0;dvdr=0.0  !******prepare the parameters of REBO potential  !rshort=3.5  !rlong=3.9  zshort=3.5  zlong=4.5  rnn=3.2 !! to adjust  rnc=3.5 !! to adjust  rnh=3.0 !! to adjust  rch=3.2 !! to adjust  rhh=3.2 !! to adjust  !rph=rlong   if (First) THEN   First=.FALSE.   INQUIRE(FILE="parameter.dat",EXIST=FExist)   if (.NOT.Fexist) THEN   WRITE(*,*) "!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!"   WRITE(*,*) "!! egrad_chamfre !!!!"   WRITE(*,*) "The file parameter.dat does not exists."   WRITE(*,*) "!! egrad_chamfre !!!!"   WRITE(*,*) "!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!"   STOP   END IF  open(1,file="parameter.dat")  do i=1,5 !Ni-Ni  read(1,*)ann(i)  enddo  do i=1,5 !Ni-C  read(1,*)anc(i)  enddo  do i=1,5 !Ni-H  read(1,*)anh(i)  enddo  do i=1,5 !H-H  read(1,*)ahh(i)  enddo  do i=1,5 !C-H  read(1,*)ach(i)  enddo  !********LONG INTERACTION **********  !do i=1,2  !read(1,*)a(i)  !enddo  close(1)  END IF  !************************************  np=45  nh=4  nc=1  nq=np+nh+nc  !************************************   DO I=1,Nat   Iat=Order(I)   xa(I)=Geom(Iat,1)   ya(I)=Geom(Iat,2)   za(I)=Geom(Iat,3)   END DO  px=0.0;py=0.0;pz=0.0  nn=0   do m=-1,1   do n=-1,1   do i=1,nq   if(m==0.and.n==0.and.i==1)mid=nn   nn=nn+1   px(nn)=xa(i)+m*lat_a(1)+n*lat_b(1)   py(nn)=ya(i)+m*lat_a(2)+n*lat_b(2)   pz(nn)=za(i)   ka(nn)=i;kx(nn)=m;ky(nn)=n   enddo   enddo   enddo  !***********************xyz file required for MS input   open(123,file="111.xyz",ACCESS='APPEND')   write(123,*)nn   write(123,*)   do i=1,nn   if(ka(i)<(np+1))then   write(123,10)'Ni',px(i),py(i),pz(i)   elseif (ka(i)np.and.ka(j+mid)np+4) then ! Ni-C     rtemp=rnc   endif   if (ka(i+mid)>np.and.ka(i+mid)np+4) then ! H-C     rtemp=rch   endif  ` if (ka(i+mid)>np.and.ka(i+mid)np.and.ka(j+mid)np.and.ka(i+mid)np+4.and.ka(j+mid)np+4.and.ka(j+mid)>np.and.ka(j+mid)0.1.and.temp