
     FUNCTION CFPOT(R,COEF,IDENT,DFPOT)
     IMPLICIT REAL*8(A-H,O-Z)
     DIMENSION COEF(10)
     DATA IWRT/0/,v1/1000.d0/,v2/5.d0/
     GO TO (10,20,30,50,110,210,220),IDENT
-------------------------------------------------------
     ...LENNARD-JONES POTENTIAL
-------------------------------------------------------
  10 EPS=COEF(1)
     RM=COEF(2)
     REDR=R/RM
     R6I=(1.0/REDR)**6
     CFPOT=EPS*(R6I-2.0)*R6I
     DFPOT=-12.0*EPS*(R6I-1.0)*R6I/R
     RETURN
-------------------------------------------------------
     ...12-6-8 POTENTIAL
-------------------------------------------------------
  20 EPS=COEF(1)
     RM=COEF(2)
      C8=COEF(3)*10000.0
     REDR=R/RM
     R8I=(1.0/R)**8.0
     RED6I=(1.0/REDR)**6.0
     CFPOT=EPS*RED6I*(RED6I-2.0)-C8*R8I
     DFPOT=-12.0*EPS*RED6I*(RED6I-1.0)/R+8.0*C8*R8I/R
     RETURN
-------------------------------------------------------
     ...MORSE POTENTIAL
-------------------------------------------------------
  30 EPS=COEF(1)
     RM=COEF(2)
     BETA=COEF(3)
  40 EXPX=EXP(-BETA*(R-RM)/RM)
     CFPOT=EPS*EXPX*(EXPX-2.0)
     DFPOT=-2.0*BETA*EPS*EXPX*(EXPX-1.0)/RM
     RETURN
-------------------------------------------------------
     ...MORSE-SPLINE-VAN DER WAALS
-------------------------------------------------------
  50 EPS=COEF(1)
     RM=COEF(2)
     BETA=COEF(3)
     C6=COEF(4)
     C8=COEF(5)
     RHO=R/RM
     RHO2=1.6
     RHO1=1.D0+0.693147180559945D0/BETA
     IF(RHO.GE.RHO2)GO TO 80
     IF(RHO.GE.RHO1)GO TO 70
  60 EXPX=EXP(-BETA*(R-RM)/RM)
     CFPOT=EPS*EXPX*(EXPX-2.0)
     DFPOT=-2.0*BETA*EPS*EXPX*(EXPX-1.0)/RM
     RETURN
  70 C=(C6+C8/(1.6*RM)**2)/(1.6*RM)**6.0
     CP=-(C6*6.0+8.0*C8/(1.6*RM)**2)/(1.6*RM)**7.0*RM
     DEL1=RHO-RHO1
     DEL2=RHO2-RHO
     EXP1=EXP(-BETA*(RHO1-1.0))
     CF1=EPS*EXP1*(EXP1-2.0)
     DF1=-2.0*BETA*EPS*EXP1*(EXP1-1.0)
     CF2=-C
     DF2=-CP
      DF1=DF1/CF1
      CF1=DLOG(-CF1)
      DF2=DF2/DF2
      CF2=DLOG(-CF2)
     RHO21=RHO2-RHO1
     A=(3.0*CF1+DF1*RHO21)/RHO21
     B=(3.0*CF2-DF2*RHO21)/RHO21
     S3=(2.0*A-B)/3.0
     S4=(2.0*B-A)/3.0
     S1=(CF1-S3*RHO21)/(RHO21**3.0)
     S2=(CF2-S4*RHO21)/(RHO21**3.0)
     CFPOT=DEL2*(S1*DEL2*DEL2+S3)+DEL1*(S2*DEL1*DEL1+S4)
     DFPOT=-3.0*S1*DEL2*DEL2-S3+3.0*S2*DEL1*DEL1+S4
      CFPOT=-DEXP(CFPOT)
      DFPOT=DFPOT*CFPOT
     DFPOT=DFPOT/RM
     RETURN
  80 R2I=(1.0/R)**2
     R6I=R2I**3
     CFPOT=-R6I*(C6+C8*R2I)
     DFPOT=6.0*C6/R+8.0*C8*R2I/R
     DFPOT=DFPOT*R6I
     RETURN
-------------------------------------------------------
     ...BUCKINGHAM-CORNER IWTH C8 AS A VARIABLE
5    EPS=COEF(1)
-------------------------------------------------------
  90 RM=COEF(2)
     BETA=4.0
     C6=COEF(4)/(EPS*RM**6.0)
     C8=COEF(3)
     B=(6.0*C6+8.0*C8)/(C6+C8-1.0)
     A=(8.0-2.0*C6)*EXP(B)/(B-8.0)
     REDR=R/RM
     RSQI=1.0/(REDR*REDR)
     R6I=RSQI**3
     IF(REDR.GT.1.0)GO TO 100
     RHO3I=(ABS(1.0/REDR-1.0))**3
     IF(REDR.GT.1.0)RHO3I=-RHO3I
     CFPOT=EPS*(A*EXP(-B*REDR)-EXP(-BETA*(RHO3I    ))*(C6+C8*RSQI)*R6I)
     DFPOT=EPS*(-A*B*EXP(-B*REDR)+EXP(-BETA*(RHO3I    ))*(6.0*C6+8.0*C8
    1*RSQI)*R6I/REDR -BETA*EXP(-BETA*(RHO3I+1.0))*(C6+C8*RSQI)*R6I*ABS(
    21.0/REDR-1.0)**2*3.0/(REDR*REDR))
     DFPOT=DFPOT/RM
     RETURN
 100 CFPOT=EPS*(A*EXP(-B*REDR)-(C6+C8*RSQI)*R6I)
     DFPOT=EPS*(-A*B*EXP(-B*REDR)+(6.0*C6+8.0*C8*RSQI)*R6I/REDR)/RM
     RETURN
-------------------------------------------------------------------
 H-F-D
   SEE JCP 71(6) 2637
-------------------------------------------------------------------
 110 CONTINUE
     ALFA=COEF(3)
     A=COEF(1)
     RM=COEF(2)
     C6=COEF(4)
     C8=COEF(5)
     C10=0.0
     X=R/RM
     IF(X.GE.1.280)GO TO 190
     B=1.28*RM
     C=B/R-1.0
     D=EXP(-C*C)
     RSQI=1.0/(R*R)
     FAC=2.0*B*C*RSQI
     DD=D*FAC
     DDD=D*(FAC*FAC-(4.0*B/R-4.0*C)*B*RSQI/R)
     GO TO 200
 190 D=1.0
     DD=0.0
     DDD=0.0
 200 SR=A*EXP(-ALFA*R)
     DSR=-ALFA*SR
     DDSR=-ALFA*DSR
     RSQI=1.0/(R*R)
     R6I=RSQI*RSQI*RSQI
     VDW=-(C6+(C8+C10*RSQI)*RSQI)*R6I
     CFPOT=SR+VDW*D
     DVDW=(6.0*C6+(8.0*C8+10.0*C10*RSQI)*RSQI)*R6I/R
     DDVDW=-(42.0*C6+(72.0*C8+110.0*C10*RSQI)*RSQI)*R6I*RSQI
     D1POT=DSR+DVDW*D+VDW*DD
     DFPOT=D1POT
     RETURN
-------------------------------------------------------------------
 MORSE-MORSE-SPLINE VAN DER WAALS
-------------------------------------------------------------------
 210 EPS=COEF(1)
     RM=COEF(2)
     BETA=COEF(3)
     BETAP=COEF(6)
     RHO=R/RM
     SIGMA=1.0D0+DLOG(0.5D0)/BETA
     IF(RHO.GT.SIGMA)GO TO 50
     RMP=SIGMA*RM/(1.0D0+DLOG(0.5D0)/BETAP)
     EPSP=EPS*BETA*RMP/(BETAP*RM)
     BETA=BETAP
     RM=RMP
     EPS=EPSP
     IF(RHO.NE.0.0)GO TO 40
     RETURN
-------------------------------------------------------------------
 ESMMSV potential
-------------------------------------------------------------------
 220 EPS=COEF(1)
     RM=COEF(2)
     BETA=COEF(3)
     BETAP=COEF(6)
     rho=r/rm
c      SIGMA=1.0D0+DLOG(0.5D0)/BETA
c      if(rho.gt.sigma)go to 50
c      RMP=SIGMA*RM/(1.0D0+DLOG(0.5D0)/BETAP)
c      EPSP=EPS*BETA*RMP/(BETAP*RM)
     A=COEF(7)
     ALPHA=COEF(8)
c      R1=(A-DLOG(v1))/ALPHA
c      R2=RMP-RMP/BETAP*DLOG(1.D0+DSQRT(1.D0+v2/EPSP))
-------------------------------------------------------------------
 check to see if r > r2 (if so use the MMSV potential)
-------------------------------------------------------------------
c      IF(R.GT.R2)GO TO 210
-------------------------------------------------------------------
 check to see if r < r1 (if so use repulsive exponential)
-------------------------------------------------------------------
c      IF(R.LT.R1)THEN
     CFPOT=EXP(A-ALPHA*R)
     DFPOT=-ALPHA*CFPOT
     RETURN
c      END IF
-------------------------------------------------------------------
 This can be reached only if r1 < r < r2
-------------------------------------------------------------------
99   a1=a-alpha*r1
     a2=(dlog(v2)-a1)/(r2-r1)
     a3=(a2+alpha)/(r2-r1)
     X2=DEXP(BETAP*(1.D0-R2/RMP))
     VR2=EPSP*X2*(X2-2.D0)
     VPR2=EPSP*BETAP/RMP*X2*2.d0*(1.D0-X2)
     a4=((VPR2/VR2-A2)/(R2-R1)-A3)/(R2-R1)
     cfpot=dexp(a1+(r-r1)*(a2+(r-r2)*(a3+(r-r1)*a4)))
     dfpot=cfpot*(A2+(R-R2)*(A3+(R-R1)*A4)+(R-R1)*(A3
    1 +(2.d0*R-R1-r2)*A4))
     RETURN
     END
