C
C	ABAQUS user material subroutine for the consitutitive behavior of bulk metallic
C	glasses, using the Spaepen's free volume modle. The code development can be found
C	in the paper:
C		Y.F. Gao, "An implicit finite element method for simulating inhomogeneous
C		deformation and shear bands of amorphous alloys based on the free-volume
C		model," Modelling Simul. Mater. Sci. Eng. 14, 1329-1345, 2006.
C
C	The Fortran code and the help document are last modified on May 10, 2007.
C
C		Yanfei Gao, Univ. Tennessee and Oak Ridge National Lab.
C
      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
     2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,
     3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
     4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 MATERL
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),
     4 DFGRD0(3,3),DFGRD1(3,3)
C
      DIMENSION EELAS(6),EPLAS(6),FLOW(6),FINCR(2)
	DIMENSION SSDEVN(NTENS),SSDEVNP(NTENS),SSDEVNT(NTENS),DSNDEV(NTENS)
      PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
      DATA NEWTON,TOLER/25,1.D-6/
C
C -----------------------------------------------------------
C	Please refer to ABAQUS manual for the meaning of these arrays
C -----------------------------------------------------------
C     PROPS(1) - E
C	PROPS(2) - NU
C     PROPS(3) - Vi
C     PROPS(4) - nD
C	PROPS(5) - ALPHA
C	PROPS(6) - BETA
C	PROPS(7) - SIGMA0
C
C	STATEV(1) - Vf
C	STATEV(2) - FLAG (0 FOR INITIAL VALUE)
C	STATEV(3) - dVf (INIITAL GUESS taken from dVf in previous increment)
C	STATEV(4) - dSMISES (INITIAL GUESS taken from dVf in previous increment)
C -----------------------------------------------------------
C
      IF (NDI.NE.3) THEN
         WRITE(6,1)
 1       FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',
     1          'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')
C	which means that this code is applicable for CPE4, CAX4, C3D8, etc.
      ENDIF
C
C     INITIAL VALUES
C
	IF (STATEV(2).EQ.0) THEN
C		Shear bands may be initiated from the perturbation sites of the initial
C		free volume field, or from the stress nonuniformity
	   STATEV(1)=PROPS(3)
C	   IF ((NOEL.EQ.4866).OR.(NOEL.EQ.5396)) THEN
C	      STATEV(1)=PROPS(3)*1.002D0
C	   ENDIF
	   STATEV(2)=1
	   STATEV(3)=0.D0
	   STATEV(4)=0.D0
	ENDIF
C
	EMOD=PROPS(1)
	ENU=PROPS(2)
      IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499
	EBULK3=EMOD/(ONE-TWO*ENU)
      EG2=EMOD/(ONE+ENU)
      EG=EG2/TWO
      ELAM=(EBULK3-EG2)/THREE
	ND=PROPS(4)
	ALPHA=PROPS(5)
	BETA=PROPS(6)
	SIGMA0=PROPS(7)
C
	VF=STATEV(1)
	DVF=STATEV(3)
	DSMISES=STATEV(4)
C
C	STEP (i) -- compute deviatoric stress
C
	SMISESN=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2))+
	1	(STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3))+
	1	(STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1))
	DO 10 K1=NDI+1,NTENS
		SMISESN=SMISESN+SIX*STRESS(K1)*STRESS(K1)
 10	CONTINUE
	SMISESN=SQRT(SMISESN/TWO)
	SSTRACEN=STRESS(1)+STRESS(2)+STRESS(3)
	DO 20 K1=1,NDI
		SSDEVN(K1)=STRESS(K1)-SSTRACEN/THREE
 20	CONTINUE
	DO 30 K1=NDI+1,NTENS
		SSDEVN(K1)=STRESS(K1)
 30	CONTINUE
C	
C	STEP (ii) -- compute deviatoric strain
C
	DO 40 K1=1,NDI
		DSNDEV(K1)=DSTRAN(K1)-(DSTRAN(1)+DSTRAN(2)+DSTRAN(3))/THREE
 40	CONTINUE
      DO 50 K1=NDI+1,NTENS
		DSNDEV(K1)=DSTRAN(K1)
 50	CONTINUE
C	
C	J2 INVARIANT OF DSNDEV
C
	DSNMIS=(DSNDEV(1)-DSNDEV(2))*(DSNDEV(1)-DSNDEV(2))+
	1	(DSNDEV(2)-DSNDEV(3))*(DSNDEV(2)-DSNDEV(3))+
	1	(DSNDEV(3)-DSNDEV(1))*(DSNDEV(3)-DSNDEV(1))
	DO 55 K1=NDI+1,NTENS
		DSNMIS=DSNMIS+SIX*DSNDEV(K1)*DSNDEV(K1)
 55	CONTINUE
C
C	IF DSNMIS=0, WE CANNOT PROCEED, E.G., BETASTAR=0/0 IN STEP (vi).
C	JUMP TO THE	ELASTIC DEFORMAITON CASE.
C
	IF (DSNMIS.EQ.0.D0) GOTO 200
C
C	STEP (iii) -- compute trial stress, Eq. (19)
C
	DO 60 K1=1,NDI
		SSDEVNT(K1)=SSDEVN(K1)+EG2*DSNDEV(K1)
 60	CONTINUE
	DO 70 K1=NDI+1,NTENS
		SSDEVNT(K1)=SSDEVN(K1)+EG*DSNDEV(K1)
 70	CONTINUE
	SMISEST=(SSDEVNT(1)-SSDEVNT(2))*(SSDEVNT(1)-SSDEVNT(2))+
	1	(SSDEVNT(2)-SSDEVNT(3))*(SSDEVNT(2)-SSDEVNT(3))+
	1	(SSDEVNT(3)-SSDEVNT(1))*(SSDEVNT(3)-SSDEVNT(1))
	DO 80 K1=NDI+1,NTENS
		SMISEST=SMISEST+SIX*SSDEVNT(K1)*SSDEVNT(K1)
 80	CONTINUE
	SMISEST=SQRT(SMISEST/TWO)
C
C	STEP (iv) -- compute free-volume and Mises-stress increments, Appendix A
C
	CALL DERIVATIVES(VF,SMISESN,DVF,DSMISES,DTIME,ALPHA,BETA,EMOD,ENU,
     1 ND,SIGMA0,SMISEST,FINCR)
C
	DVFNEW=DVF+FINCR(1)
	DSMNEW=DSMISES+FINCR(2)
C
	DO 90 KEWTON=1,NEWTON
C		NEED TO MODIFY THE FOLLOWING IF CONDITION, AT LEAST NORMALIZE 
	 IF ((ABS(DVFNEW-DVF).LT.TOLER).AND.(ABS(DSMNEW-DSMISES).LT.TOLER))
	1    GOTO 100
	   DVF=DVFNEW
	   DSMISES=DSMNEW
	   CALL DERIVATIVES(VF,SMISESN,DVF,DSMISES,DTIME,ALPHA,BETA,EMOD,
     1    ENU,ND,SIGMA0,SMISEST,FINCR)
	   DVFNEW=DVF+FINCR(1)
	   DSMNEW=DSMISES+FINCR(2)
 90	CONTINUE
	WRITE(6,2) NEWTON
 2    FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
     1        'CONVERGE AFTER ',I3,' ITERATIONS')
	PNEWDT=0.4D0
 100	CONTINUE
C
C	UPDATE STATE VARIABLES
C
	STATEV(1)=VF+DVF
	STATEV(3)=DVF
	STATEV(4)=DSMISES
C
C	STEP (v) -- compute Mises-strain increment, Eq. (23)
C
	SMISESNP=SMISESN+DSMISES
	DSNEFF=(SMISEST-SMISESNP)*2.D0*(1.D0+ENU)/3.D0/EMOD
C
C	STEP (vi) -- compute Eq. (21) and (22)
C	IF (SMISEST.EQ.0) BETASTAR=1.0D0
C	
	BETASTAR=1.D0-(SMISEST-SMISESNP)/SMISEST
	DO 110 K1=1,NDI
		SSDEVNP(K1)=BETASTAR*SSDEVNT(K1)
 110	CONTINUE
	DO 120 K1=NDI+1,NTENS
		SSDEVNP(K1)=BETASTAR*SSDEVNT(K1)
 120	CONTINUE
C	WRITE(*,*) BETASTAR
C
C	STEP (vii) -- compute the updated stress tensor, Eqs. (14) and (15)
C
	SSTRACENP=SSTRACEN+EBULK3*(DSTRAN(1)+DSTRAN(2)+DSTRAN(3))
	DO 130 K1=1,NDI
		STRESS(K1)=SSDEVNP(K1)+SSTRACENP/THREE
 130	CONTINUE
	DO 140 K1=NDI+1,NTENS
		STRESS(K1)=SSDEVNP(K1)
 140	CONTINUE
C
C	STEP (viii) -- compute psi in Eq. (B7)
C
	TEMP1=DEXP(SMISESNP/SIGMA0)
	TEMP2=DEXP(-SMISESNP/SIGMA0)
	TEMPCOSH=0.5*(TEMP1+TEMP2)
	TEMPSINH=0.5*(TEMP1-TEMP2)
	C1=DTIME/ALPHA*DEXP(-1.D0/STATEV(1))*THREE*(ONE-ENU)/EMOD/BETA
	C2=C1*TEMPSINH/STATEV(1)
	C3=C1*SIGMA0*(TEMPCOSH-1.D0)/STATEV(1)/STATEV(1)
	FAI=C2/(1-STATEV(3)/STATEV(1)/STATEV(1)+C3)
	C1=EMOD*DTIME/(1.D0+ENU)*DEXP(-1.D0/STATEV(1))
	C2=TEMPCOSH/SIGMA0+FAI*TEMPSINH/STATEV(1)/STATEV(1)
	PSI=C1*C2/(1+C1*C2)
C
C	STEP (ix) -- compute the elastic-plastic consistent tangent
C
	C1=THREE*EG*(1-BETASTAR-PSI)/SMISESNP/SMISESNP
	DO 150 K1=1,NDI
		DDSDDE(K1,K1)=TWO*EMOD/THREE/(1.D0+ENU)*BETASTAR+
	1		C1*SSDEVNP(K1)*SSDEVNP(K1)+EBULK3/THREE
 150	CONTINUE
	DDSDDE(1,2)=-EG2/THREE*BETASTAR+EBULK3/THREE+
	1	C1*SSDEVNP(1)*SSDEVNP(2)
	DDSDDE(2,1)=DDSDDE(1,2)
	DDSDDE(1,3)=-EG2/THREE*BETASTAR+EBULK3/THREE+
	1	C1*SSDEVNP(1)*SSDEVNP(3)
	DDSDDE(3,1)=DDSDDE(1,3)
	DDSDDE(2,3)=-EG2/THREE*BETASTAR+EBULK3/THREE+
	1	C1*SSDEVNP(2)*SSDEVNP(3)
	DDSDDE(3,2)=DDSDDE(2,3)
	DO 160 K1=1,NDI
		DO 170 K2=NDI+1,NTENS
			DDSDDE(K1,K2)=C1*SSDEVNP(K1)*SSDEVNP(K2)
			DDSDDE(K2,K1)=DDSDDE(K1,K2)
 170		CONTINUE
 160	CONTINUE
	DO 180 K1=NDI+1,NTENS
		DDSDDE(K1,K1)=EG*BETASTAR+C1*SSDEVNP(K1)*SSDEVNP(K1)
 180	CONTINUE
	IF (NTENS.EQ.6) THEN
		DDSDDE(4,5)=C1*SSDEVNP(4)*SSDEVNP(5)
		DDSDDE(5,4)=DDSDDE(4,5)
		DDSDDE(4,6)=C1*SSDEVNP(4)*SSDEVNP(6)
		DDSDDE(6,4)=DDSDDE(4,6)
		DDSDDE(5,6)=C1*SSDEVNP(5)*SSDEVNP(6)
		DDSDDE(6,5)=DDSDDE(5,6)
	ELSEIF (NTENS.EQ.4) THEN
C		DDSDDE(4,4) HAS BEEN GIVEN
	ELSE
		WRITE(6,3) 
 3		FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',
     1          'ELEMENTS WITH ONE OR THREE SHEAR COMPONENTS')
	ENDIF
C
	GOTO 300
C
 200	CONTINUE
C	
C	ELASTIC STIFFNESS MATRIX
C
      DO 210 K1=1,NDI
        DO 220 K2=1,NDI
           DDSDDE(K2,K1)=ELAM
 220    CONTINUE
        DDSDDE(K1,K1)=EG2+ELAM
 210  CONTINUE
      DO 230 K1=NDI+1,NTENS
        DDSDDE(K1,K1)=EG
 230  CONTINUE

 300	CONTINUE
C	WRITE(*,*) STRESS(6),DDSDDE(6,6)
C
      RETURN
      END
C
C	Refer to Appendix A for the following subroutine
C
	SUBROUTINE DERIVATIVES(VF,TAO,DVF,DTAO,DTIME,ALPHA,BETA,EMOD,ENU,
     1 ND,SIGMA0,SMISEST,FINCR)
C	
	INCLUDE 'ABA_PARAM.INC'
C
	DIMENSION FINCR(2)
C
	TEMP1=DEXP((TAO+DTAO)/SIGMA0)
	TEMP2=DEXP((-TAO-DTAO)/SIGMA0)
	TEMPCOSH=0.5*(TEMP1+TEMP2)
	TEMPSINH=0.5*(TEMP1-TEMP2)
C
	C1=(TEMPCOSH-1.D0)*3.D0*(1-ENU)*SIGMA0/EMOD/BETA
	C1=C1/(VF+DVF)-1.D0/ND
	F1=DVF-DTIME/ALPHA*DEXP(-1.D0/(VF+DVF))*C1
	C1=DEXP(-1.D0/(VF+DVF))*TEMPSINH
	F2=SMISEST-TAO-DTAO-EMOD/(1.D0+ENU)*DTIME*C1
C
	C1=(TEMPCOSH-1.D0)/EMOD/BETA*SIGMA0*3.D0*(1.D0-ENU)
	C1=C1*(1.D0/(VF+DVF)-1.D0)-1.D0/ND
	C2=1.D0/(VF+DVF)/(VF+DVF)*C1
	DJACOB11=1.D0-DTIME/ALPHA*DEXP(-1.D0/(VF+DVF))*C2
	C1=3.D0*(1.D0-ENU)/(VF+DVF)/EMOD/BETA*TEMPSINH
	DJACOB12=-DTIME/ALPHA*DEXP(-1.D0/(VF+DVF))*C1
	C1=TEMPSINH/(VF+DVF)/(VF+DVF)
	DJACOB21=-EMOD/(1.D0+ENU)*DTIME*DEXP(-1.D0/(VF+DVF))*C1
	C1=DEXP(-1.D0/(VF+DVF))*TEMPCOSH/SIGMA0
	DJACOB22=-1.D0-EMOD/(1.D0+ENU)*DTIME*C1
C
	DETJACO=DJACOB11*DJACOB22-DJACOB12*DJACOB21
	FINCR(1)=-1/DETJACO*(DJACOB22*F1-DJACOB12*F2)
	FINCR(2)=-1/DETJACO*(-DJACOB21*F1+DJACOB11*F2)
C
      RETURN
      END