cmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc
      subroutine FHH_SW_init()
      call prepot_FHH_SW
c F+H2 SW PES: Stark and Werner, JCP 104, 6515 (1996)
      return
      end
cmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc
      subroutine FHH_SW_pes(rs, vv)
      implicit none
      real*8 rs(3), vv
      !! Common block
      real*8 RSV(3), V, DV(3)
      COMMON /POTCM_FHH_SW/ RSV,V,DV

      rsv(2)=rs(3)
      rsv(3)=rs(1)
      rsv(1)=rs(2)
      call POTE_FHH_SW
      vv=v

      return
      end
cmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc
      SUBROUTINE PREPOT_FHH_SW
C
C     Ab initio P.E.S. for F+H2 built by Stark and Werner
C     Modifications by L. Banares and V. S ez R banos October 1992
C
      IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4(I-N)
      INTEGER*4 B(198),C(198),D(198)
      DIMENSION NPARM(198),A(198),P(12)
      COMMON /POTCM_FHH_SW/R(3),PE,DPE(3)

      !!--PESLIB--!!
      character*256 parapath
      integer*4 parapathlen
C
C..... Read the three-body-parameters
C..... (computed with the Levenberg-Marquardt-Procedure)
C..... nparm is the actual number of a parameter
C
      call ftn_getenv('PESLIB', parapath, parapathlen)
      OPEN(3,FILE=parapath(1:parapathlen)//'/para/FHH_SW/three.param')       


      I=1
   10 READ(3,*,END=11) NPARM(I),B(I),C(I),D(I),A(I)
      I=I+1
      GOTO 10
C
   11 MA=I-1
      if(ma.eq.0) stop 'no three.param card'
      CLOSE(3)
C
C.... Read the two-body-parameters p
C.... (computed with an Extended-Rydberg-Fit)
C
      !!--PESLIB--!!
      OPEN(3,FILE=parapath(1:parapathlen)//'/para/FHH_SW/two.param')


      I=1
   12 READ(3,*,END=111) P(I)
      I=I+1
      GOTO 12
  111 NP=I-1
      if(np.eq.0) stop 'no two.param card'
      CLOSE(3)
C
C.... Initialize the non-linear parameters
C
      RETURN
      ENTRY POTE_FHH_SW
C
      B1 = A(MA-5)
      B2 = A(MA-4)
      B3 = A(MA-3)
      X0 = A(MA-2)
      Y0 = A(MA-1)
      Z0 = A(MA)
      FIT = 0.0D0
      DFDX = 0.0D0
      DFDY = 0.0D0
      DFDZ = 0.0D0
C
C... The (modified) Aguado-Paniagua type functions
C... for the three-body-potential
C
      XEXPON = B1*(R(1)-X0)
      YEXPON = B2*(R(2)-Y0)
      ZEXPON = B3*(R(3)-Z0)

           EXPONX=DEXP(-XEXPON)
           EXPONY=DEXP(-YEXPON)
           EXPONZ=DEXP(-ZEXPON)

      FEX = R(1)*EXPONX
      FEY = R(2)*EXPONY
      FEZ = R(3)*EXPONZ


          DRHOX = EXPONX*(1-R(1)*B1)
          DRHOY = EXPONY*(1-R(2)*B2)
          DRHOZ = EXPONZ*(1-R(3)*B3)



      DO 1010 I=1,MA-6

        PEPEX = FEX
        PEPEY = FEY
        PEPEZ = FEZ

             IF(B(I).EQ.0) THEN
              PEPEX=1.D0
              GO TO 1018
             ENDIF
             DO 1020 J=1,B(I)-1
             PEPEX = PEPEX*FEX
 1020        CONTINUE
 1018        IF(C(I).EQ.0) THEN
              PEPEY=1.D0
              GO TO 1019
             ENDIF

             DO 1021 J=1,C(I)-1
             PEPEY = PEPEY*FEY
 1021        CONTINUE
 1019        IF(D(I).EQ.0) THEN
              PEPEZ=1.D0
              GO TO 1024
             ENDIF

            DO 1022 J=1,D(I)-1
             PEPEZ = PEPEZ*FEZ
 1022        CONTINUE



 1024    FXY = PEPEX*PEPEY*PEPEZ
         FIT = FIT + A(I)*FXY

c---- modification fex,fez are zero!!
        dftx=0.d0
        dfty=0.d0
        dftz=0.d0
        if(fex.ne.0.d0) DFTX = DRHOX*B(I)*PEPEX/FEX
        if(fey.ne.0.d0) DFTY = DRHOY*C(I)*PEPEY/FEY
        if(fez.ne.0.d0) DFTZ = DRHOZ*D(I)*PEPEZ/FEZ
c         DFTX = DRHOX*B(I)*PEPEX/FEX
c         DFTY = DRHOY*C(I)*PEPEY/FEY
c         DFTZ = DRHOZ*D(I)*PEPEZ/FEZ
c-----end of modificaiton

          DFTX = DFTX*PEPEY*PEPEZ
          DFTY = DFTY*PEPEX*PEPEZ
          DFTZ = DFTZ*PEPEX*PEPEY

          DFDX = DFDX+A(I)*DFTX
          DFDY = DFDY+A(I)*DFTY
          DFDZ = DFDZ+A(I)*DFTZ
C

 1010 CONTINUE
C
C.... Two-Body-Potential : Extended-Rydberg-Functional
C
      XR = R(1)-P(3)
      YR = R(2)-P(9)
      ZR = R(3)-P(3)
        XR2=XR*XR
        XR3=XR2*XR
        YR2=YR*YR
        YR3=YR2*YR
        ZR2=ZR*ZR
        ZR3=ZR2*ZR

      FX = DEXP(-P(2)*XR)
      FY = DEXP(-P(8)*YR)
      FZ = DEXP(-P(2)*ZR)

      UX = -P(1)*(1.0D0+P(2)*XR+P(4)*XR2+P(5)*XR3)
      UY = -P(7)*(1.0D0+P(8)*YR+P(10)*YR2+P(11)*YR3)
      UZ = -P(1)*(1.0D0+P(2)*ZR+P(4)*ZR2+P(5)*ZR3)

      XVAL=UX*FX+P(6)
      YVAL=UY*FY+P(12)
      ZVAL=UZ*FZ+P(6)

C

      DUX = -P(1)*(P(2)+P(4)*2.D0*XR+P(5)*3.D0*XR2)
      DUY = -P(7)*(P(8)+P(10)*2.D0*YR+P(11)*3.D0*YR2)
      DUZ = -P(1)*(P(2)+P(4)*2.D0*ZR+P(5)*3.D0*ZR2)
C
      DFX = -P(2)*FX
      DFY = -P(8)*FY
      DFZ = -P(2)*FZ
C
      DX = DUX*FX+UX*DFX
      DY = DUY*FY+UY*DFY
      DZ = DUZ*FZ+UZ*DFZ
C
C.... resulting Energy in au
C
      DPE(1) = DFDX+DX
      DPE(2) = DFDY+DY
      DPE(3) = DFDZ+DZ
C
      PE = FIT+XVAL+YVAL+ZVAL
C
      RETURN
      END