Chimie Théorique » OpenPath

 subroutine egrad_ext(e,geomcart,gradcart)

  ! This routines calculates the energy and the gradient of 

  ! a molecule, using an external code

  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, dzdc, indzmat,Nom,Atome, massat, unit,ProgExe

  use Io_module

  !

  IMPLICIT NONE

  ! Energy (calculated if F300K=.F., else estimated)

  REAL(KREAL), INTENT (OUT) :: e

  ! Nb degree of freedom

  ! Geometry at which gradient is calculated (cf FActual also)

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

  ! Gradient calculated at Geom geometry

  REAL(KREAL), INTENT (OUT) :: gradcart(3*Nat)

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

  character(LCHARS) :: LINE

  logical           :: debug

  LOGICAL :: fexist, FSim

  LOGICAL, SAVE :: first=.true.

  LOGICAL, ALLOCATABLE :: Done(:)

  REAL(KREAL), SAVE :: Eold=1.e6

  REAL(KREAL) :: d, a_val, Pi

  REAL(KREAL) :: coef,x

  INTEGER(KINT) :: iat, jat, kat, i, j, n3at, absidg, absidg2

  INTEGER(KINT) :: kTmp, Istart,ITmp1,ITmp2, Idx

  INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11

  !

  CHARACTER(132) :: FileIn,FileOut

  CHARACTER(LCHARS) :: ListName, TitleTmp, CH32SVAR1

  CHARACTER(VLCHARS), SAVE :: RstrtCopy, RunCommand

  LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False.

  LOGICAL   :: FRdyn,FStopNucl, FTmp, Tchk,Tchk1

  INTEGER(kint) :: NbLists,LastIt,Nat4,NStep, Firstit

  INTEGER(KINT) :: NStepAd,NStepTmp

  INTEGER(KINT) :: IShowTmp

  REAL(KREAL)   :: FricPsiMin,FricPsiT, FricNuclTmp

  INTEGER(KINT), PARAMETER :: NbExtName=4

  INTEGER(KINT) :: ICouc

  INTEGER(KINT) :: ILine

  INTEGER(KINT) :: NPulses, PulseLen, NStepWarm,NStepAv,NStepEq,NStepTot

  REAL(KREAL) :: TempWarm,Mean,Slope,Dev,RTmp1,RTmp2

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

  LOGICAL, EXTERNAL :: valid

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

  Pi=dacos(-1.0d0)

  n3at=3*nat

  debug=valid('EGRAD')

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

  RunCommand=Trim(Adjustl(ProgExe))

  FileIn=Trim(CalcName) // Trim(ISuffix)

  FileOut=Trim(CalcName) // Trim(OSuffix)

  IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(RunCommand)

  ! we create the input file

  OPEN(IOTMP,File=FileIn)

  WRITE(IOTMP,'(1X,I10)') NAt

  WRITE(IOTMP,'(1X,A)') Coord

     DO I=1,Nat

        If (renum) THEN

           Iat=Order(I)

           WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),GeomCart(Iat,:)

        ELSE

           Iat=OrderInv(I)

           WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomCart(I,:)

        END IF

     END DO

  call system(RunCommand)

  if (debug) WRITE(*,*) 'DBG EGRAD, back from calculation'

  OPEN(IOTMP,FILE=FileOut, STATUS='old')

  ! We first search for the forces

  READ(IOTMP,*) e

  DO I=1,Nat

     Iat=I

     IF (renum) Iat=Order(I)

     READ(IOTMP,*)  GradCart(3*Iat-2:3*Iat)

  END DO

  CLOSE(IOTMP)

  if (debug) WRITE(*,*) '================  Egrad_ext Over ===================='

  RETURN

999 CONTINUE

  if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!'

  STOP

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

end subroutine egrad_ext