Statistics
| Revision:

## root / src / egrad_test_2D.f90 @ 5

 1  SUBROUTINE egrad_test_2D(nat,e,geom,grad)  ! This program computes the energy and the gradient  ! to cheat PATH we use a 3 atoms system   IMPLICIT NONE   integer, parameter :: KINT = kind(1)   integer, parameter :: KREAL = kind(1.0d0)   INTEGER(KINT), INTENT(IN) :: Nat   REAL(KREAL), INTENT(OUT) :: E,grad(Nat*3)   REAL(KREAL), INTENT(IN) :: Geom(Nat,3)  ! Bohr --> Angstr  ! real(KREAL), parameter :: BOHR = 0.52917726D+00  !  ! Parameters to define the surface   INTEGER(KINT), DIMENSION(6), PARAMETER :: IECOEF = (/-1,9,-45,45,-9,1/)   INTEGER(KINT), DIMENSION(6), PARAMETER :: ISCOEF = (/-3,-2,-1,1,2,3/)   REAL(KREAL), PARAMETER :: hh=0.01d0   INTEGER(KINT), PARAMETER :: IOTMP=99  ! Variables   INTEGER(KINT) :: i, iat, jat,j   REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:)   LOGICAL :: Debug   LOGICAL, save :: First     REAL(KREAL) :: x,y   INTERFACE   function valid(string) result (isValid)   CHARACTER(*), intent(in) :: string   logical :: isValid   END function VALID   function E2D(natoms,Xyz)   use Path_module, only : order   use Io_module, only : IOTMP   IMPLICIT NONE   integer, parameter :: KINT = kind(1)   integer, parameter :: KREAL = kind(1.0d0)   INTEGER(KINT) ,INTENT(IN) :: natoms   REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3)   REAL(KREAL) :: E2D   END function E2D     END INTERFACE   debug=valid('egrad_test_2D')   if (debug) WRITE(*,*) '================ Entering Egrad_test ==================='   if (debug) THEN   WRITE(*,*) "Cartesian Geometry"   DO I=1,Nat   WRITE(*,'(1X,I5,3(1X,F12.6))') I,Geom(i,1:3)   END DO   END IF     ALLOCATE(XyZTmp(Nat,3),GradTmp(Nat,3))   if (first) THEN   OPEN(IOTMP,File="tt")   DO I=1, 501   DO j=1,501   x=(i-1.)/100.   y=(j-1.)/100.   xyztmp(2,1)=x   xyztmp(3,1)=y   write(IOTMP,*) x,y,E2D(nat,xyztmp)   END DO  ! we skip a blank line   WRITe(IOTMP,*)   END DO   CLOSE(IOTMP)   First=.FALSE.   END IF   e=E2D(nat,Geom)  ! We now calculate the gradients using numerical derivatives   Grad=0.d0   GradTmp=0.d0   do iat=1,1   do jat=2,nat   xyztmp=geom   do i=1,6   xyztmp(jat,iat)=geom(jat,iat)+ISCoef(i)*hh   gradTmp(jat,iat)=gradTmp(jat,iat)+IECoef(i)*E2D(nat,xyztmp)   end do   end do   end do   gradTmp=gradTmp/(60.*hh)   do iat=1,nat   grad(3*iat-2:3*iat)=gradTmp(iat,1:3)   end do   if (debug) THEN   WRITE(*,*) "Cartesian gradient GradTmp"   DO I=1,Nat   WRITE(*,'(1X,I5,3(1X,F12.6))') I,Gradtmp(i,1:3)   END DO   WRITE(*,*) "Cartesian gradient Grad"   WRITE(*,'(50(1X,F12.6))') Grad   END IF   deallocate(xyztmp,gradTmp)   if (debug) WRITE(*,*) '================ Exiting Egrad_test_2D ==================='  ! ======================================================================   end SUBROUTINE egrad_test_2D   function E2D(natoms,Xyz)   use Path_module, only : order   use Io_module, only : IOTMP   IMPLICIT NONE   integer, parameter :: KINT = kind(1)   integer, parameter :: KREAL = kind(1.0d0)   INTEGER(KINT) ,INTENT(IN) :: natoms   REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3)   REAL(KREAL) :: E,x,y, E2D   REAL(KREAL), save :: Data(501,501)   INTEGER(KINT) :: i,j   REAL(KREAL) :: dCH,dNH,dCN  ! real(KREAL) :: autoA,autoeV   real(KREAL) :: autoeV  ! parameter (autoA=0.52917715d0)   parameter (autoeV=27.21183d0)   real(KREAL) :: r1, r2   LOGICAL :: Debug   LOGICAL, SAVE :: First=.TRUE.     INTERFACE   function valid(string) result (isValid)   CHARACTER(*), intent(in) :: string   logical :: isValid   END function VALID     END INTERFACE   debug=valid('e2D')   if (First) THEN  ! we read the data   open(IOTMP,File="CCSDT.subtr.reg")   DO I=1, 501   DO j=1,501   READ(IOTMP,*) x,y,E   Data(int(x*100+0.1)+1,int(y*100+0.1)+1)=E   write(*,*) x,y,int(x*100+0.1)+1, int(y*100+0.1)+1,E   END DO  ! we skip a blank line   READ(IOTMP,*)   END DO   CLOSE(IOTMP)   First=.FALSE.   END IF   x=Xyz(2,1)   y=Xyz(3,1)   E2D=Data(int(x*100+0.1)+1,int(y*100+0.1)+1)   return   end function E2D