Chimie Théorique » Opt'n Path

  SUBROUTINE egrad_LEPS(nat,e,geom,grad) 

! This program computes the ernergy and gradien in cartesian coordinates

! for the cartesian geometry Geom

!----------------------------------------------------------------------

!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon, 

!  Centre National de la Recherche Scientifique,

!  Université Claude Bernard Lyon 1. All rights reserved.

!

!  This work is registered with the Agency for the Protection of Programs 

!  as IDDN.FR.001.100009.000.S.P.2014.000.30625

!

!  Authors: P. Fleurat-Lessard, P. Dayal

!  Contact: optnpath@gmail.com

!

! This file is part of "Opt'n Path".

!

!  "Opt'n Path" is free software: you can redistribute it and/or modify

!  it under the terms of the GNU Affero General Public License as

!  published by the Free Software Foundation, either version 3 of the License,

!  or (at your option) any later version.

!

!  "Opt'n Path" is distributed in the hope that it will be useful,

!  but WITHOUT ANY WARRANTY; without even the implied warranty of

!

!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

!  GNU Affero General Public License for more details.

!

!  You should have received a copy of the GNU Affero General Public License

!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.

!

! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,

! for commercial licensing opportunities.

!----------------------------------------------------------------------

   IMPLICIT NONE

     integer, parameter :: KINT = kind(1)

     integer, parameter :: KREAL = kind(1.0d0)

! Number of atoms

       INTEGER(KINT), INTENT(IN) :: Nat

! Input geometry

       REAL(KREAL), INTENT(IN) :: Geom(Nat,3)

! output energy and gradient

       REAL(KREAL), INTENT(OUT) :: E,grad(Nat*3)

! 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.001d0

! Variables

       INTEGER(KINT) :: i,iat,jat

       REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:)

       LOGICAL :: Debug

      REAL(KREAL), external :: ELEPS_xyz

  INTERFACE

     function valid(string) result (isValid)

       CHARACTER(*), intent(in) :: string

       logical                  :: isValid

     END function VALID

  END INTERFACE

  debug=valid('egrad_leps')

 if (debug) WRITE(*,*) '================ Entering Egrad_leps ==================='

 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))

      e=ELEPS_xyz(nat,Geom)

! We now calculate the gradients using numerical derivatives

      Grad=0.d0

      GradTmp=0.d0

      do iat=1,3

         do jat=1,nat

            xyztmp=geom

            do i=1,6

               xyztmp(jat,iat)=geom(jat,iat)+ISCoef(i)*hh

               gradTmp(jat,iat)=gradTmp(jat,iat)+IECoef(i)*ELEPS_xyz(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 "

    DO I=1,Nat

       WRITE(*,'(1X,I5,3(1X,F12.6))') I,Gradtmp(i,1:3)

    END DO

 END IF

      deallocate(xyztmp,gradTmp)

 if (debug) WRITE(*,*) '================ Exiting Egrad_leps ==================='

! ======================================================================

      end 

      function ELEPS_xyz(natoms,Xyz)

        use Path_module, only : order

        use Io_module, only : au2ev

   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) :: ELEPS_xyz

      INTEGER(KINT) :: i

      REAL(KREAL) :: rAB, rBC, rAC

      REAL(KREAL), PARAMETER :: a=0.05, b=0.30, c=0.05

      REAL(KREAL), PARAMETER :: dAB=4.746, dBC=4.747, dAC=3.445

      REAL(KREAL), PARAMETER :: r0=0.742, alpha=1.942

      REAL(KREAL) :: Qbc,QAc,Qab,Jab,Jbc,Jac

      rAB=0.

      rAC=0.

      rBC=0.

      do i=1,3

         rBC=rBC+(xyz(Order(3),i)-xyz(order(2),i))**2

         rAB=rAB+(xyz(order(1),i)-xyz(order(2),i))**2

         rAC=rAC+(xyz(order(1),i)-xyz(order(3),i))**2

      end do

      rAB=sqrt(rAB)

      rBC=sqrt(rBC)

      rAC=sqrt(rAC)

! For consistency with previous articles : rAC is constrained to be 3.742

      rAC=3.742

      Qab = dab*(1.5*exp(-2*alpha*(rAB-r0))-exp(-alpha*(rAB-r0)))/(2+2*a)

      Qbc = dbc*(1.5*exp(-2*alpha*(rBC-r0))-exp(-alpha*(rBC-r0)))/(2+2*b)

      Qac = dac*(1.5*exp (-2*alpha*(rAC-r0))-exp(-alpha*(rAC-r0)))/(2+2*c)

      Jab = dab*(exp(-2*alpha*(rAB-r0))-6*exp(-alpha*(rAB-r0)))/(4+4*a)

      Jbc = dbc*(exp(-2*alpha*(rBC-r0))-6*exp(-alpha*(rBC-r0)))/(4+4*b)

      Jac = dac*(exp(-2*alpha*(rAC-r0))-6*exp(-alpha*(rAC-r0)))/(4+4*c)

! V(x,y) = Qab(x)+Qbc(y)+Qac(rac)-sqrt(Jab(x)**2+Jbc(y)**2+Jac(x)**2-Jab(x)*Jbc(y)-Jbc(y)*Jac(rac)-Jab(x)*Jac(rac))

      ELEPS_xyz=(Qab+Qbc+Qac-sqrt(Jab**2+Jbc**2+Jac**2-Jab*Jbc-Jbc*Jac-Jab*Jac))/au2eV

      return

    end function ELEPS_xyz