SUBROUTINE ACT_OPX(T, P, XEN, XFS, ENSTATITE, FERROSILITE, 
     1   LSTOP)
C 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IMPLICIT INTEGER(I-N) 
C 
      LOGICAL LSTOP, FIRST
      DATA FIRST / .TRUE. /
C
      IF(FIRST) THEN
         FIRST = .FALSE.
         EXPMIN = X02AEF(EXPMIN)/3.0
         EXPMAX = X02AFF(EXPMAX)/3.0
         EPSMCH = X02AAF(EPSMCH)
         SREPS = SQRT(EPSMCH)
      END IF
C
      R_OPX = 1.0 - 2.0*XEN
      IF((R_OPX.LT.-1.0).OR.(R_OPX.GT.1.0)) THEN
           LSTOP = .TRUE. 
           RETURN 
      END IF
      LSTOP = .FALSE. 
C
      HEX = -1870.0 ! cals
      VEX = -0.029 ! cals/bar
      HX = -450.0 ! cals
      W = 2000.0 ! cals
      WV = 0.0135 ! cals/bar
C     
      OP_OLD = 2.0
      OP = 0.0
      ITER = 0
      DO WHILE(ABS(OP-OP_OLD).GT.EPSMCH)
         ITER = ITER + 1
         IF(ITER.GT.250) THEN
            LSTOP = .TRUE.
            RETURN
         END IF  
         XMGM1 = 0.5*(1.0-R_OPX+OP)
         XFEM1 = 0.5*(1.0+R_OPX-OP)
         XMGM2 = 0.5*(1.0-R_OPX-OP)
         XFEM2 = 0.5*(1.0+R_OPX+OP)
         DGDOP = - 1.98726*T*LOG(XFEM1*XMGM2/(XMGM1*XFEM2))
     1           + (HEX+VEX*(P-1.0)) - (2.0*W-HX+WV*(P-1.0))*OP
         D2GDOP = 0.5*1.98726*T*(1.0/XFEM1+1.0/XMGM2
     1            + 1.0/XMGM1 + 1.0/XFEM2) - (2.0*W-HX+WV*(P-1.0))
         OP_OLD = OP
         OP=OP-DGDOP/D2GDOP
         OP = MAX(OP,R_OPX-1.0+EPSMCH,-R_OPX-1.0+EPSMCH)
         OP = MIN(OP,R_OPX+1.0-EPSMCH,-R_OPX+1.0-EPSMCH)
      END DO
C
      GAM1 = W*(XFEM1**2+XFEM2**2) + HX*XFEM1*XFEM2 
     1     + 0.25*WV*(P-1.0)*(1.0+R_OPX)**2
      GAM1 = GAM1/(1.98726*T)
      IF((GAM1.GE.EXPMIN).AND.(GAM1.LE.EXPMAX)) THEN
         GAM1 = EXP(GAM1)*XMGM1*XMGM2
      ELSE
         LSTOP = .TRUE.
         RETURN
      END IF
C 
      GAM2 = W*(XMGM1**2+XMGM2**2) + HX*XMGM1*XMGM2
     1     + 0.25*WV*(P-1.0)*(1.0-R_OPX)**2
      GAM2 = GAM2/(1.98726*T)
      IF((GAM2.GE.EXPMIN).AND.(GAM2.LE.EXPMAX)) THEN
         GAM2 = EXP(GAM2)*XFEM1*XFEM2
      ELSE
         LSTOP = .TRUE.
         RETURN
      END IF
C 
      ENSTATITE = SQRT(GAM1) 
      FERROSILITE = SQRT(GAM2) 
C 
      RETURN
      END