Chimie Théorique » OpenPath

   PROGRAM PES2D

   IMPLICIT NONE

     integer, parameter :: KINT = kind(1)

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

	   REAL(KREAL) :: E

!     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

! Variables

       INTEGER(KINT) :: i,j,iat,jat

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

       REAL(KREAL) :: xp,yp

       CHARACTER(132) :: Line

       LOGICAL :: Debug

      real(KREAL), parameter :: zero = 0.0_KREAL,   one = 1.0_KREAL

      REAL(KREAL), external :: Eint

      REAL(KREAL) :: dCH,dCN,aHCN

	  REAL(KREAL), PARAMETER :: dCHmin=1., dCHmax=2.5

	  REAL(KREAL), PARAMETER :: aHCNmin=0., aHCNmax=180.

	  REAL(KREAL) :: StepdCH, StepaHCN

	  INTEGER(KINT) :: NbStpdCH,NbStpaHCN

! For now dCN=1.16 is fixed. Next version should optimize it at each step !

	dCN=1.16

	NbStpdCH=100

	NbStpaHCN=100

	StepdCH=(dCHmax-dCHmin)/(NbStpdCH-1)

	StepaHCN=(aHCNmax-aHCNmin)/(NbStpaHCN-1)

	dCH=dCHmin

	DO I=1,NbStpdCH

	  aHCN=aHCNmin

	  DO J=1,NbStpaHCN

	    e=Eint(dCH,dCN,aHCN)

	    WRITE(*,'(1X,F10.4,1X,F10.4,1X,F8.2,1X,F15.6)') dCH,dCN,aHCN,e

	    aHCN=aHCN+StepaHCN

	  END DO

	  WRITE(*,*) 

	 dCH=dCH+StepdCH

	END DO

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

	end  

     function Eint(dCH,dCN,aHCN)

   IMPLICIT NONE

     integer, parameter :: KINT = kind(1)

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

      REAL(KREAL) ,INTENT(IN) :: dCH,dCN,aHCN

      REAL(KREAL) :: Eint

      INTEGER(KINT) :: i

      REAL(KREAL) :: dNH

      real(KREAL) :: autoA,autoeV,Pi

      parameter (autoA=0.52917715d0)

      parameter (autoeV=27.21183d0)

      real(KREAL) :: Q(3)

      real(KREAL) :: TT,gR,pR

      real(KREAL) :: mB,mC,mH

      real(KREAL) :: T,RN,RC,r1,r2,r3

      real(KREAL) :: Z1,Z12,V1

      real(KREAL) :: Z2,Z22,V2

      real(KREAL) :: Z3,Z32,V3

      real(KREAL) :: S1,S2,S3,S12,S22,S32,POLY

      real(KREAL) :: E1,E2,E3,HY1,HY2,HY3,SWITCH,V123

      pi=dacos(-1.0d0)

      dNH=sqrt(dCH**2+dCN**2-2*dCH*dCN*cos(aHCN*pi/180.))

      R1=dCH

      R2=dCN

      R3=dNH

!      R1=R1*autoA

!      R2=R2*autoA

!      R3=R3*autoA

!.....CH

      Z1=R1-1.0823

      Z12=Z1*Z1

      V1=-2.8521*(1.0+5.5297*Z1+8.7166*Z12+5.3082*Z1*Z12)*dexp(-5.5297*Z1)

!.....CN

      Z2=R2-1.1718

      Z22=Z2*Z2 

      V2=-7.9282*(1.0+5.2448*Z2+7.3416*Z22+4.9785*Z2*Z22)*dexp(-5.2448*Z2)

!.....NH

      Z3=R3-1.0370

      Z32=Z3*Z3

      V3=-3.9938*(1.0+3.0704*Z3)*dexp(-3.0704*Z3)

!     write(6,*) v1,v2,v3

!.....THREE BODY TERMS

      Z1=R1-1.9607

      Z2=R2-2.2794

      Z3=R3-1.8687

      S1=0.4436*Z1+0.6091*Z2+0.6575*Z3

      S2=-.8941*Z1+0.2498*Z2+0.3718*Z3

      S3=0.0622*Z1-0.7527*Z2+0.6554*Z3

      S12=S1*S1

      S22=S2*S2

      S32=S3*S3

      POLY=-3.0578*(1.0+1.9076*S1-0.5008*S2-0.0149*S3+0.6695*S12  &

          -1.3535*S22-1.0501*S32+0.2698*S1*S2-1.1120*S1*S3        &

          +1.9310*S2*S3-0.0877*S1*S12+0.0044*S2*S22+0.0700        &

          *S3*S32+0.0898*S12*S2-1.0186*S1*S22-0.0911*S12          &

          *S3+0.0017*S1*S32+0.4567*S22*S3-0.8840*S2*S32           &

          +0.3333*S1*S2*S3-0.0367*S12*S12+0.4821*S22*S22          &

          +0.2564*S32*S32-0.0017*S12*S1*S2-0.2278*S12*S22         &

          -0.1287*S1*S2*S22+0.1759*S1*S12*S3-0.0399*S12*S32       &

          -0.1447*S1*S3*S32-0.3147*S2*S22*S3+0.1233*S22*S32       &

          +0.3161*S2*S3*S32+0.0919*S12*S2*S3-0.0954*S1*S22*S3     &

          +0.1778*S1*S2*S32-0.1892*S22*S22*S2)

      E1=dexp(3.9742*Z1/2.)

      E2=dexp(4.3688*Z2/2.)

      E3=dexp(1.5176*Z3/2.)

      HY1=(E1-1.0/E1)/(E1+1.0/E1)

      HY2=(E2-1.0/E2)/(E2+1.0/E2)

      HY3=(E3-1.0/E3)/(E3+1.0/E3)

      SWITCH=(1.0-HY1)*(1.0-HY2)*(1.0-HY3)

      V123=SWITCH*POLY

!     write(6,*) v123,V1+V2+V3+V123

      Eint=(V1+V2+V3+V123)/autoeV

      return

      end function Eint