	
      PROGRAM PROTENTIAL
	
C      PARAMETER(NOFPLATES=6)
C	NO OF PLATES EXCEPT SURFACE

	
      PARAMETER(nofl=10)
C	NO OF ATOMS PER SIDE
	
	REAL COVERAGE
	REAL CRYSTCONST,NETCHARGE
	REAL AAA,AGG,AGA
	REAL BAA,BGG,BGA
	REAL LMDAA,LMDGG,LMDGA
	REAL MUAA,MUGG,MUGA
	REAL RAA,RGG,RGA
	REAL SAA,SGG,SGA
	REAL BTAA,BTGG,BTGA
	REAL NAA,NGG,NGA
	REAL CAA,CGG,CGA
	REAL DAA,DGG,DGA
	REAL HAA,HGG,HGA

	PARAMETER(COVERAGE=.4)
C	COVERAGE RATE OF THE SURFACE     
     
     
     
      PARAMETER(PI=3.14159265) 
      PARAMETER(CRYSTCONST=5.65315)
	PARAMETER(NETCHARGE=.087)

      PARAMETER(AAA=4773.163,  AGG=2839.821,  AGA=13287.605)
	PARAMETER(BAA=81.32,     BGG=114.786,   BGA=13.187   )
	PARAMETER(LMDAA=3.4841,  LMDGG=3.2834,  LMDGA=4.5989 )
	PARAMETER(MUAA=1.7154,   MUGG=1.7154,   MUGA=.2493   )
	PARAMETER(RAA=3.4,       RGG=3.4,       RGA=3.4      )
	PARAMETER(SAA=3.6,       SGG=3.6,       SGA=3.6      )
	PARAMETER(BTAA=.0074881, BTGG=.2358624, BTGA=.357192 )
	PARAMETER(NAA=.6087913,  NGG=3.4729041, NGA=6.31747  )
	PARAMETER(CAA=5.2731318, CGG=.0762977,  CGA=1.226302 )
	PARAMETER(DAA=.7510266,  DGG=19.796474, DGA=.790396  )
	PARAMETER(HAA=.1529235,  HGG=7.1459174, HGA=-.518489 )
C	 ABOVE IS THE PARA OF POTENTIAL
C	  SEE THE REF. FOR DETAILS
     
     	INTEGER NEIBOUR(20)
      COMMON/PROP/SUME,SUMV,SUMA
	COMMON/NEIBOUR/ NEIBOUR
      REAL KAPPA,DIST
	COMMON/KAPPA/ KAPPA
C	COMMON/IRR/ IRR
      COMMON/PICK/NPICK(800)
      COMMON/PROPA/RMAX,CUBE,CUBE2,HMAX,DELTA,
     A        NP,IRR,IACPT
	COMMON/IIDUM/IIDUM	
	COMMON/TOTPOT/TOTPOT
C	  KAPPA OF THE ERROR FUNCTION
	
     

	REAL*8 SUME,SUMV,SUMA
	REAL TOTPOT
      REAL AI(8*NOFL*NOFL*NOFL,5)
C	
C	ATOM INFORMATION
C		1ST: KIND OF ATOM (1 FOR GA),(2 FOR AS),(0 FOR NONE)
C		2ND: KINETIC ENERGY
C		3RD-5TH:POSITION
C	COMMON/AI/ AI()
	REAL POT,TEMPERATURE,CR
	INTEGER POSKEEPTIME,POTKEEPTIME,TOPTYPE,DISTYPE
C	 POTETIAL(VARIABLE)
     
     
C--------------------------------------------
C
C	MAIN ROUTINE START
C
C--------------------------------------------
V--                                       --V
C	WRITE(*,*) 'PLEASE GIVE LATTICE SIZE'
C	WRITE(*,*) 'AT LEAST 5'
C	READ *,LS
C	IF(LS.LE.4) THEN
C 	WRITE (*,*) 'ERROR!'
C	CALL EXIT()
C	ENDIF
C	WRITE(*,*) 'LATTICE SIZE SETTED AS ',LS,'BY',LS
C	PRINT *,''

	WRITE(*,*) 'PLEASE GIVE TOP ATOM TYPE:'
	WRITE(*,*) '1 FOR Ga,2 FOR As'
      READ *,TOPTYPE
	IF(TOPTYPE.EQ.1) THEN
	WRITE(*,*) 'OK,Ga SELECTED AS TOP'
	ELSEIF(TOPTYPE.EQ.2) THEN
	WRITE(*,*) 'OK,As SELECTED AS TOP'
	ELSE
	WRITE (*,*) 'ERROR!'
	CALL EXIT()
	ENDIF
	PRINT *,''

	WRITE(*,*) 'PLEASE GIVE DISTRIBUTION TYPE:'
	WRITE(*,*) '0 FOR RANDOM,1 FOR REGULAR'
      READ *,DISTYPE
	IF(DISTYPE.EQ.0) THEN
	WRITE(*,*) 'OK,TOP ATOM DISTRIBUTED RANDOMLY'
	ELSEIF(DISTYPE.EQ.1) THEN
	WRITE(*,*) 'OK,TOP ATOM DISTRIBUTED REGULARLY'
	ELSE
	WRITE (*,*) 'ERROR!'
	CALL EXIT()
	ENDIF
	PRINT *,''

	WRITE(*,*) 'PLEASE GIVE COVERAGE RATE'
	WRITE(*,*) 'FROM .1 TO .99'
	READ *,CR
	IF(CR.GE.1.OR.CR.LE.0) THEN
	WRITE (*,*) 'ERROR!'
	CALL EXIT()
	ENDIF
	WRITE(*,*) 'COVERAGE SELECTED AS',CR*100,'%'
	PRINT *,''

	WRITE(*,*) 'PLEASE GIVE TEMPERATURE IN K'
	WRITE(*,*) 'FROM 1 TO INFINITY'
	READ *,TEMPERATURE
	IF(TEMPERATURE.LE.0) THEN
	WRITE (*,*) 'ERROR!'
	CALL EXIT()
	ENDIF
	WRITE(*,*) 'TEMPERATURE SETED AS',TEMPERATURE,'K'
	PRINT *,''

	WRITE(*,*) 'PLEASE GIVE THE TOTAL RUN TIME'
	READ *,MAXKB
	IF(MAXKB.LE.0) THEN
	WRITE (*,*) 'ERROR!'
	CALL EXIT()
	ENDIF
	WRITE(*,*) 'THERE WILL BE TOTALLY ',MAXKB,'s RUN'
	WRITE(*,*) '----------SIM BEGIN----------'


     	KAPPA=SQRT(PI)/(CRYSTCONST*NOFL)
     	IRR=34567
	NP=8*NOFL*NOFL
	DELTA=.1    
	SUME=0
	SUMV=0
	SUMA=0
C	TR=.274
	TR=(1.381/1.6)*(1E-4)*TEMPERATURE
	BETA=1./TR
	KWRITE=10
	IIDUM=62777587
	IF(MAXKB.GT.1000) THEN
	POSKEEPTIME=MAXKB/100
	POTKEEPTIME=MAXKB/50
      ELSE
	POSKEEPTIME=50
	POTKEEPTIME=50
	ENDIF
C	PRINT *,POSKEEPTIME,POTKEEPTIME



	CALL SETUP(AI,TOPTYPE,DISTYPE,CR)
	
	OPEN(4,FILE='POT.DAT');
	PRINT *,TR
			CALL TOTALPOT(AI,NP)
			WRITE(4,*) 0,TOTPOT
			CALL SHOWPOS(AI,NP,0/10)

C	CALL SHOWPOS(AI,200,1)
      DO 97 KB=1,MAXKB
		IF(MOD(KB,POTKEEPTIME).EQ.0) THEN
			CALL TOTALPOT(AI,NP)
			WRITE(4,*) KB,TOTPOT
		ENDIF
		IF(MOD(KB,POSKEEPTIME).EQ.0) THEN
			CALL SHOWPOS(AI,NP,KB/10)
		ENDIF
          CALL PICKP(NP)
          CALL PMOVE(AI,TOTV,TOTE,BETA)
		PRINT *,KB,'th RUN'
C          CALL EVAL(TOTV,TOTE)
C
C === CALCULATE MEAN SQUARE DISPLACEMENT
C
D          TDIST=0.
D          DO 40 I=1,NP
C          TDIST=TDIST+DAX(I)**2+DAY(I)**2+DAZ(I)**2
D  40     CONTINUE
D          TDIST=TDIST/PART
C
C === PROPERTY CALCULATION & PRINT-OUT AT INTERVALS
C
D         IF(MOD(KB,KWRITE).NE.0) GOTO 99
D           FKB=FLOAT(NP)*SUMA
C           ENR=(SUME/FKB+CORE)
C           VIR=(SUMV/FKB+CORV)/(3.*TR)
C           PRES=TR*DR*(1.-VIR)
D            ACPT=FLOAT(IACPT)/FLOAT(KB*NP)
D            WRITE(6,49) KB,SUMA,ACPT
D49       FORMAT(1X,I6,4F9.3)
C ===  TEST FOR EQUILIBRATION

C  50     IF(IFLG.LT.1) CALL EQBRAT(KDIST,IFLG,NY,KB)
  97     CONTINUE





	OPEN(11,FILE='ALLCONFIG.DAT');
	IF(TOPTYPE.EQ.1) THEN
	WRITE(11,*) 'Ga SELECTED AS TOP'
	ELSEIF(TOPTYPE.EQ.2) THEN
	WRITE(11,*) 'As SELECTED AS TOP'
	ELSE
	WRITE (11,*) 'ERROR!'
	CALL EXIT()
	ENDIF

	IF(DISTYPE.EQ.0) THEN
	WRITE(11,*) 'TOP ATOM DISTRIBUTED RANDOMLY'
	ELSEIF(DISTYPE.EQ.1) THEN
	WRITE(11,*) 'TOP ATOM DISTRIBUTED REGULARLY'
	ELSE
	WRITE (11,*) 'ERROR!'
	CALL EXIT()
	ENDIF

	WRITE(11,*) 'COVERAGE SELECTED AS',CR*100,'%'

	WRITE(11,*) 'TEMPERATURE SETED AS',TEMPERATURE,'K'

	WRITE(11,*) 'THERE WILL BE TOTALLY ',MAXKB,'s RUN'

	WRITE(11,*) 'POT RECORDED EVERY ',POTKEEPTIME,'TIMES OF RUN'

	WRITE(11,*) 'POSITION RECORDED EVERY ',POSKEEPTIME,'TIMES OF RUN'

	WRITE (11,*) 'DATA BEGIN'

	DO 98 I=1,8*NOFL*NOFL*NOFL
		WRITE(11,*) AI(I,1),AI(I,3),AI(I,4),AI(I,5)
98	CONTINUE

	WRITE(11,*) 'ACCEPT RATE IS',FLOAT(IACPT)/FLOAT(KB*NP)
      CLOSE(11)
	PRINT *,'CONFIG RECORDED ON ALLCONFIG.DAT'
	PRINT *,'DONE'
c	PAUSE
C	PAUSE	
      END
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
C--------------------------------------------
C
C	MAIN ROUTINE END
C
C--------------------------------------------


	
C--------------------------------------------
C
C	SETUP SUBROUTINE START
C
C	 BUILD UP AN INITIAL LATTICE 
C
C--------------------------------------------
Cvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv     
      SUBROUTINE SETUP(AI,TOPTYPE,DISTYPE,CR)
C      PARAMETER(NOFPLATES=5)
C	NO OF PLATES EXCEPT SURFACE

	PARAMETER(nofl=10)
C	NO OF ATOMS PER SIDE

C	PARAMETER(COVERAGE=.5)

	INTEGER TOPTYPE,DISTYPE
	REAL CRYSTCONST,NETCHARGE,CR
	REAL AAA,AGG,AGA
	REAL BAA,BGG,BGA
	REAL LMDAA,LMDGG,LMDGA
	REAL MUAA,MUGG,MUGA
	REAL RAA,RGG,RGA
	REAL SAA,SGG,SGA
	REAL BTAA,BTGG,BTGA
	REAL NAA,NGG,NGA
	REAL CAA,CGG,CGA
	REAL DAA,DGG,DGA
	REAL HAA,HGG,HGA

      
      PARAMETER(PI=3.14159265) 
      PARAMETER(CRYSTCONST=5.65315)
	PARAMETER(NETCHARGE=.087)
		
      PARAMETER(AAA=4773.163,  AGG=2839.821,  AGA=13287.605)
	PARAMETER(BAA=81.32,     BGG=114.786,   BGA=13.187   )
	PARAMETER(LMDAA=3.4841,  LMDGG=3.2834,  LMDGA=4.5989 )
	PARAMETER(MUAA=1.7154,   MUGG=1.7154,   MUGA=.2493   )
	PARAMETER(RAA=3.4,       RGG=3.4,       RGA=3.4      )
	PARAMETER(SAA=3.6,       SGG=3.6,       SGA=3.6      )
	PARAMETER(BTAA=.0074881, BTGG=.2358624, BTGA=.357192 )
	PARAMETER(NAA=.6087913,  NGG=3.4729041, NGA=6.31747  )
	PARAMETER(CAA=5.2731318, CGG=.0762977,  CGA=1.226302 )
	PARAMETER(DAA=.7510266,  DGG=19.796474, DGA=.790396  )
	PARAMETER(HAA=.1529235,  HGG=7.1459174, HGA=-.518489 )



	REAL AI(8*NOFL*NOFL*NOFL,5)
 
C	SETUP MAIN LATTICE, NOFL*NOFL*NOFL*8 ATOMS
     
      AI(1,1)=1
	AI(1,2)=0
	AI(1,3)=0
	AI(1,4)=0
	AI(1,5)=0 
	

	DO 100 I=0,NOFL-1
	DO 100 J=0,NOFL-1
	DO 100 K=0,NOFL-1

	  AI((K+J*NOFL+I*NOFL*NOFL)*8+1,1)=1
        AI((K+J*NOFL+I*NOFL*NOFL)*8+1,3)=AI(1,3)+I*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+1,4)=AI(1,4)+J*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+1,5)=AI(1,5)+K*CRYSTCONST

	  AI((K+J*NOFL+I*NOFL*NOFL)*8+2,1)=1
        AI((K+J*NOFL+I*NOFL*NOFL)*8+2,3)=AI(1,3)+(I+.0)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+2,4)=AI(1,4)+(J+.5)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+2,5)=AI(1,5)+(K+.5)*CRYSTCONST

	  AI((K+J*NOFL+I*NOFL*NOFL)*8+3,1)=1
        AI((K+J*NOFL+I*NOFL*NOFL)*8+3,3)=AI(1,3)+(I+.5)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+3,4)=AI(1,4)+(J+.0)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+3,5)=AI(1,5)+(K+.5)*CRYSTCONST

	  AI((K+J*NOFL+I*NOFL*NOFL)*8+4,1)=1
        AI((K+J*NOFL+I*NOFL*NOFL)*8+4,3)=AI(1,3)+(I+.5)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+4,4)=AI(1,4)+(J+.5)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+4,5)=AI(1,5)+(K+.0)*CRYSTCONST

	  AI((K+J*NOFL+I*NOFL*NOFL)*8+5,1)=2
        AI((K+J*NOFL+I*NOFL*NOFL)*8+5,3)=AI(1,3)+(I+.25)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+5,4)=AI(1,4)+(J+.25)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+5,5)=AI(1,5)+(K+.25)*CRYSTCONST

	  AI((K+J*NOFL+I*NOFL*NOFL)*8+6,1)=2
        AI((K+J*NOFL+I*NOFL*NOFL)*8+6,3)=AI(1,3)+(I+.75)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+6,4)=AI(1,4)+(J+.75)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+6,5)=AI(1,5)+(K+.25)*CRYSTCONST

	  AI((K+J*NOFL+I*NOFL*NOFL)*8+7,1)=2
        AI((K+J*NOFL+I*NOFL*NOFL)*8+7,3)=AI(1,3)+(I+.75)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+7,4)=AI(1,4)+(J+.25)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+7,5)=AI(1,5)+(K+.75)*CRYSTCONST

	  AI((K+J*NOFL+I*NOFL*NOFL)*8+8,1)=2
        AI((K+J*NOFL+I*NOFL*NOFL)*8+8,3)=AI(1,3)+(I+.25)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+8,4)=AI(1,4)+(J+.75)*CRYSTCONST
	  AI((K+J*NOFL+I*NOFL*NOFL)*8+8,5)=AI(1,5)+(K+.75)*CRYSTCONST

100	CONTINUE
	CALL STSURFACE(AI,CR,TOPTYPE,DISTYPE)

	RETURN 
      END

C--------------------------------------------
C
C	SETUP SUBROUTINE END
C
C--------------------------------------------


C--------------------------------------------
C
C	SUBROUTINE CPOTENTIAL
C
C	CALCULATE POTENTIAL 
C
C--------------------------------------------
     
      SUBROUTINE CPOTENTIAL(AI,POT,SNOFATOM)

c      PARAMETER(NOFPLATES=10)
C	NO OF PLATES EXCEPT SURFACE

	REAL CRYSTCONST,NETCHARGE
	REAL AAA,AGG,AGA
	REAL BAA,BGG,BGA
	REAL LMDAA,LMDGG,LMDGA
	REAL MUAA,MUGG,MUGA
	REAL RAA,RGG,RGA
	REAL SAA,SGG,SGA
	REAL BTAA,BTGG,BTGA
	REAL NAA,NGG,NGA
	REAL CAA,CGG,CGA
	REAL DAA,DGG,DGA
	REAL HAA,HGG,HGA


	PARAMETER(nofl=10)
C	NO OF ATOMS PER SIDE

	PARAMETER(COVERAGE=.5)
C	COVERAGE RATE OF THE SURFACE     
     
     
     
      PARAMETER(PI=3.14159265) 
      PARAMETER(CRYSTCONST=5.65315)
	PARAMETER(NETCHARGE=.087)	

      PARAMETER(AAA=4773.163,  AGG=2839.821,  AGA=13287.605)
	PARAMETER(BAA=81.32,     BGG=114.786,   BGA=13.187   )
	PARAMETER(LMDAA=3.4841,  LMDGG=3.2834,  LMDGA=4.5989 )
	PARAMETER(MUAA=1.7154,   MUGG=1.7154,   MUGA=.2493   )
	PARAMETER(RAA=3.4,       RGG=3.4,       RGA=3.4      )
	PARAMETER(SAA=3.6,       SGG=3.6,       SGA=3.6      )
	PARAMETER(BTAA=.0074881, BTGG=.2358624, BTGA=.357192 )
	PARAMETER(NAA=.6087913,  NGG=3.4729041, NGA=6.31747  )
	PARAMETER(CAA=5.2731318, CGG=.0762977,  CGA=1.226302 )
	PARAMETER(DAA=.7510266,  DGG=19.796474, DGA=.790396  )
	PARAMETER(HAA=.1529235,  HGG=7.1459174, HGA=-.518489 )
	
C	 ABOVE IS THE PARA OF POTENTIAL
C	  SEE THE REF. FOR DETAILS
     
     
      REAL KAPPA
	COMMON/KAPPA/ KAPPA
	COMMON/IRR/ IRR
     	
	INTEGER SNOFATOM
	REAL AI(8*NOFL*NOFL*NOFL,5)


	REAL DISOFX,DISOFY,DISOFZ
	REAL DISOFXP,DISOFYP,DISOFZP
C	 DISPLACEMENT OF THE ATOM IN STUDYING AND ALL OTHER ATOM
	REAL HALFL,DISTANCE,DISTANCEP,LPOT,SPOT,KERCI,LPIN,SPIN
c	REAL FCOFR,GOFSITA,SITA
C	 FC(R),G(SITA),SITA
	REAL PA,PB,PLMD,PMU,PR,PS,PBT,PN,PC,PD,PH
C	HALF OF SIZE OF L,DISTANCE OF 2 ATOMS,LONG RANGE OF POTENTIAL,KERCI
	HALFL=NOFL*CRYSTCONST*.5
	LPOT=0.
	SPOT=0	
C	BELOW CALCULATE THE LONG&SHORT RANGE PART
	DO 200 N=1,8*NOFL*NOFL*NOFL
	  
C	  LPOT=0.
	  SPOT=0.	
	  KERCI=0.
	  
C	  DISOFX=AI(SNOFATOM,3)-AI(N,3)
C	  DISOFY=AI(SNOFATOM,4)-AI(N,4)	
C	  DISOFZ=AI(SNOFATOM,5)-AI(N,5)	
     	
C	  IF(DISOFX.GT.HALFL)  DISOFX=DISOFX-HALFL*2
C	  IF(DISOFX.LT.-HALFL) DISOFX=DISOFX+HALFL*2
C	  IF(DISOFY.GT.HALFL)  DISOFY=DISOFY-HALFL*2
C	  IF(DISOFY.LT.-HALFL) DISOFY=DISOFY+HALFL*2
C	  IF(DISOFZ.GT.HALFL)  DISOFZ=DISOFZ-HALFL*2
C	  IF(DISOFZ.LT.-HALFL) DISOFZ=DISOFZ+HALFL*2

C	  DISTANCE=SQRT(DISOFX*DISOFX+DISOFY*DISOFY+DISOFZ*DISOFZ)
C	  IF(DISTANCE.LT.SAA) THEN
C	  PRINT *,DISTANCE,N
C	  ENDIF
C	LONG RANGE PART	  
C					------HAVE CONSIDERED COVERAGE FACTOR

	  DISTANCE=GETDISTANCE(AI,SNOFATOM,N)

C	  PRINT *,'DISTANCE',DISTANCE
	  IF(SNOFATOM.NE.N.AND.DISTANCE.NE.0) THEN
	    LPIN=ERFC(KAPPA*DISTANCE)/DISTANCE
C		PRINT *,LPIN,LPOT
		
		IF(AI(SNOFATOM,1).EQ.AI(N,1)) THEN
	    LPOT=LPOT+LPIN
	    ELSE IF(AI(N,1).GE.1)	THEN
	    LPOT=LPOT-LPIN
	    ELSE IF(AI(N,1).EQ.0) THEN
	    LPOT=LPOT
	    ELSE 
	    PRINT *,'LPOT ERROR!'
	    ENDIF
C	 	PRINT *,'d,lpot',DISTANCE,LPOT
        ENDIF

200	CONTINUE 
	
	LPOT=0.108991*LPOT
C	 CONVERT LONG RANGE POT, TO THE UNIT OF eV.
	SPOT=SHORTPOT(AI,SNOFATOM)
      POT=LPOT+SPOT
c      PRINT *,'LPOT=',LPOT
c      PRINT *,'SPOT=',SPOT
c	PRINT *,'TOTAL POT=',POT
C	LONG&SHORT RANGE PART CALCULATE COMPLETE


	RETURN
	END

C--------------------------------------------
C
C	CPOTENTIAL SUBROUTINE END
C
C--------------------------------------------



C--------------------------------------------
C
C	SUBROUTINE RANDOM
C
C	GENERATE RANDOM NUMBER 
C	   
C	NOTE: SHOULD TAKE THE ABSOLUTE VALUE BEFORE USE.
C	(0 , 1)
C
C--------------------------------------------
     
C      SUBROUTINE RANDOM(RANX,RANY,RANZ)
C      COMMON/PROPA/RMAX,CUBE,CUBE2,HMAX,DELTA,
C     A        NP,IRR,IACPT
C     	
C	IRR=IRR*65539
C	RANX=FLOAT(IRR)*.4656613E-9
C	IRR=IRR*65539
C	RANY=FLOAT(IRR)*.4656613E-9
C	IRR=IRR*65539
C	RANZ=FLOAT(IRR)*.4656613E-9
C
C	RETURN 

C	END

C--------------------------------------------
C
C	RANDOM SUBROUTINE END
C
C--------------------------------------------


        REAL FUNCTION ERFC ( X )




C    *******************************************************************
C    ** APPROXIMATION TO THE COMPLEMENTARY ERROR FUNCTION             **
C    **                                                               **
C    ** REFERENCE:                                                    **
C    **                                                               **
C    ** ABRAMOWITZ AND STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS,    **
C    **    NATIONAL BUREAU OF STANDARDS, FORMULA 7.1.26               **
C    *******************************************************************


        REAL        A1, A2, A3, A4, A5, P

        PARAMETER ( A1 = 0.254829592, A2 = -0.284496736 )
        PARAMETER ( A3 = 1.421413741, A4 = -1.453152027 )
        PARAMETER ( A5 = 1.061405429, P  =  0.3275911   )


        REAL        T, X, XSQ, TP


        T  = 1.0 / ( 1.0 + P * X )
        XSQ = X * X

        TP = T * ( A1 + T * ( A2 + T * ( A3 + T * ( A4 + T * A5 ) ) ) )


        ERFC = TP * EXP ( -XSQ )

        RETURN


        END



C----------------------------------------------------------------------------	
C
C	FCOFR 
C
C----------------------------------------------------------------------------
 	   REAL FUNCTION FCOFR(DISTANCE,S,R)
	   REAL DISTANCE,S,R	
       PARAMETER(PI=3.14159265) 

	   IF(DISTANCE.LT.R) THEN
	    	FCOFR=1
	   ELSE IF(DISTANCE.GT.S) THEN
	        FCOFR=0
	   ELSE
	        FCOFR=.5+.5*COS(PI*(DISTANCE-R)/(S-R))
       ENDIF
	   IF(DISTANCE.LT.S) THEN
C	   PRINT *,'fc',FCOFR,DISTANCE
	   ENDIF
	   RETURN
	   	   END

C------------------------------------------------------------------------------
C
C	COSSITA
C
C------------------------------------------------------------------------------
	   REAL FUNCTION COSSITA(DX1,DY1,DZ1,DX2,DY2,DZ2)
	   REAL DX1,DY1,DZ1,DX2,DY2,DZ2

	   REAL DIS12,DIS13,DIS23

     	DIS12=SQRT(DX1*DX1+DY1*DY1+DZ1*DZ1)
       	DIS13=SQRT(DX2*DX2+DY2*DY2+DZ2*DZ2)
     	DIS23=SQRT((DX1-DX2)*(DX1-DX2)+(DY1-DY2)*(DY1-DY2)
     C	+(DZ1-DZ2)*(DZ1-DZ2))
	 
	   IF(DIS12*DIS13*DIS23.EQ.0) THEN
C	   PRINT *,'ERROR! COSSITA'
	   RETURN
       ENDIF
	
C       	PRINT *,'D12,D13,D23',DIS12,DIS13,DIS23
        COSSITA=(DIS12*DIS12+DIS13*DIS13-DIS23*DIS23)
     C	/(2*DIS12*DIS13)
C	   PRINT *,'COSSITA',COSSITA
       RETURN
	   END
	
C-------------------------------------------------------------------------------
C
C	GOFSITA
C
C-------------------------------------------------------------------------------
    	REAL FUNCTION GOFSITA(CS,C,D,H)
      	REAL CS,C,D,H

      	GOFSITA=1+C*C/(D*D)-C*C/(D*D+(H-CS)*(H-CS))
C		PRINT *,'C=',C
C		PRINT *,'D=',D
C		PRINT *,'H=',H
C		PRINT *,'GOFSITA=',GOFSITA
    		RETURN
      END
      

C---------------------------------------------------------------------------------
C
C     SHORTSHOT
C
C----------------------------------------------------------------------------------






    	REAL FUNCTION SHORTPOT(AI,SNOFATOM)


	INTEGER NEIBOUR(20)
	COMMON/NEIBOUR/ NEIBOUR
c      PARAMETER(NOFPLATES=10)
C	NO OF PLATES EXCEPT SURFACE

	REAL CRYSTCONST,NETCHARGE
	REAL AAA,AGG,AGA
	REAL BAA,BGG,BGA
	REAL LMDAA,LMDGG,LMDGA
	REAL MUAA,MUGG,MUGA
	REAL RAA,RGG,RGA
	REAL SAA,SGG,SGA
	REAL BTAA,BTGG,BTGA
	REAL NAA,NGG,NGA
	REAL CAA,CGG,CGA
	REAL DAA,DGG,DGA
	REAL HAA,HGG,HGA


	PARAMETER(nofl=10)
C	NO OF ATOMS PER SIDE

	PARAMETER(COVERAGE=.5)
C	COVERAGE RATE OF THE SURFACE     
     
     
     
      PARAMETER(PI=3.14159265) 
      PARAMETER(CRYSTCONST=5.65315)
	PARAMETER(NETCHARGE=.087)	

      PARAMETER(AAA=4773.163,  AGG=2839.821,  AGA=13287.605)
	PARAMETER(BAA=81.32,     BGG=114.786,   BGA=13.187   )
	PARAMETER(LMDAA=3.4841,  LMDGG=3.2834,  LMDGA=4.5989 )
	PARAMETER(MUAA=1.7154,   MUGG=1.7154,   MUGA=.2493   )
	PARAMETER(RAA=3.4,       RGG=3.4,       RGA=3.4      )
	PARAMETER(SAA=3.6,       SGG=3.6,       SGA=3.6      )
	PARAMETER(BTAA=.0074881, BTGG=.2358624, BTGA=.357192 )
	PARAMETER(NAA=.6087913,  NGG=3.4729041, NGA=6.31747  )
	PARAMETER(CAA=5.2731318, CGG=.0762977,  CGA=1.226302 )
	PARAMETER(DAA=.7510266,  DGG=19.796474, DGA=.790396  )
	PARAMETER(HAA=.1529235,  HGG=7.1459174, HGA=-.518489 )
	
C	 ABOVE IS THE PARA OF POTENTIAL
C	  SEE THE REF. FOR DETAILS
     
     
      REAL KAPPA
	COMMON/KAPPA/ KAPPA
	COMMON/IRR/ IRR
     	
	INTEGER SNOFATOM
	REAL AI(8*NOFL*NOFL*NOFL,5)


	REAL DISOFX,DISOFY,DISOFZ
	REAL DISOFXP,DISOFYP,DISOFZP
C	 DISPLACEMENT OF THE ATOM IN STUDYING AND ALL OTHER ATOM
	REAL HALFL,DISTANCE,DISTANCEP,LPOT,SPOT,KERCI,LPIN,SPIN
c	REAL FCOFR,GOFSITA,SITA
C	 FC(R),G(SITA),SITA
	REAL PA,PB,PLMD,PMU,PR,PS,PBT,PN,PC,PD,PH
	REAL PAP,PBP,PLMDP,PMUP,PRP,PSP,PBTP,PNP,PCP,PDP,PHP
C	HALF OF SIZE OF L,DISTANCE OF 2 ATOMS,LONG RANGE OF POTENTIAL,KERCI
	SPOT=0.
	KERCI=0.
	HALFL=NOFL*CRYSTCONST*.5


	CALL GETNEIBOUR(AI,SNOFATOM)

C	WRITE (*,*) (NEIBOUR(I),I=1,20)
     	NOFNEIBOUR=0
	DO 250 I=1,20
	IF(NEIBOUR(I).GT.0)	THEN
	NOFNEIBOUR=NOFNEIBOUR+1
	ENDIF
250	CONTINUE
C	PRINT *,'NOFNEIBOUR',NOFNEIBOUR
      DO 350 N=1,NOFNEIBOUR

	  DISOFX=AI(SNOFATOM,3)-AI(NEIBOUR(N),3)
	  DISOFY=AI(SNOFATOM,4)-AI(NEIBOUR(N),4)	
	  DISOFZ=AI(SNOFATOM,5)-AI(NEIBOUR(N),5)	
     	
	  IF(DISOFY.GT.HALFL)  DISOFY=DISOFY-HALFL*2
	  IF(DISOFY.LT.-HALFL) DISOFY=DISOFY+HALFL*2	  
	  IF(DISOFZ.GT.HALFL)  DISOFZ=DISOFZ-HALFL*2
	  IF(DISOFZ.LT.-HALFL) DISOFZ=DISOFZ+HALFL*2
	    
	  DISTANCE=SQRT(DISOFX*DISOFX+DISOFY*DISOFY+DISOFZ*DISOFZ)
	  KERCI=0.
	  DO 300 M=1,NOFNEIBOUR

		IF(M.NE.N) THEN
		DISOFXP=AI(SNOFATOM,3)-AI(NEIBOUR(M),3)
	    DISOFYP=AI(SNOFATOM,4)-AI(NEIBOUR(M),4)	
	    DISOFZP=AI(SNOFATOM,5)-AI(NEIBOUR(M),5)	
     	
	    IF(DISOFYP.GT.HALFL)  DISOFYP=DISOFYP-HALFL*2
	    IF(DISOFYP.LT.-HALFL) DISOFYP=DISOFYP+HALFL*2	  
	    IF(DISOFZP.GT.HALFL)  DISOFZP=DISOFZP-HALFL*2
	    IF(DISOFZP.LT.-HALFL) DISOFZP=DISOFZP+HALFL*2
	        
		DISTANCEP=SQRT(DISOFXP*DISOFXP+DISOFYP*DISOFYP
     C			+DISOFZP*DISOFZP)

  		IF(AI(SNOFATOM,1).EQ.1.AND.AI(NEIBOUR(M),1).EQ.1) THEN
	      PAP=AGG
		  PBP=BGG
		  PLMDP=LMDGG
		  PMUP  =MUGG
	      PRP=RGG
	      PSP=SGG
	      PBTP=BTGG
	      PNP=NGG
		  PCP=CGG
	      PDP=DGG
	      PHP=HGG
	    ELSE IF(AI(SNOFATOM,1).EQ.1.AND.AI(NEIBOUR(M),1).EQ.2) THEN
	      PAP=AGA
		  PBP=BGA
		  PLMDP=LMDGA
		  PMUP  =MUGA
	      PRP=RGA
	      PSP=SGA
	      PBTP=BTGA
	      PNP=NGA
		  PCP=CGA
	      PDP=DGA
	      PHP=HGA

	    ELSE IF(AI(SNOFATOM,1).EQ.2.AND.AI(NEIBOUR(M),1).EQ.1) THEN
	      PAP=AGA
		  PBP=BGA
		  PLMDP=LMDGA
		  PMUP  =MUGA
	      PRP=RGA
	      PSP=SGA
	      PBTP=BTGA
	      PNP=NGA
		  PCP=CGA
	      PDP=DGA
	      PHP=HGA

	    ELSE IF(AI(SNOFATOM,1).EQ.2.AND.AI(NEIBOUR(M),1).EQ.2) THEN
	      PAP=AAA
		  PBP=BAA
		  PLMDP=LMDAA
		  PMUP  =MUAA
	      PRP=RAA
	      PSP=SAA
	      PBTP=BTAA
	      PNP=NAA
		  PCP=CAA
	      PDP=DAA
	      PHP=HAA
		
		ENDIF
		      KERCI=KERCI+FCOFR(DISTANCEP,PSP,PRP)*GOFSITA(
     C	COSSITA(DISOFX,DISOFY,DISOFZ,DISOFXP,DISOFYP,DISOFZP),
     C    PCP,PDP,PHP)

		ENDIF
300		CONTINUE


		IF(AI(SNOFATOM,1).EQ.1.AND.AI(NEIBOUR(N),1).EQ.1) THEN
	      PA=AGG
		  PB=BGG
		  PLMD=LMDGG
		  PMU  =MUGG
	      PR=RGG
	      PS=SGG
	      PBT=BTGG
	      PN=NGG
		  PC=CGG
	      PD=DGG
	      PH=HGG
	    ELSE IF(AI(SNOFATOM,1).EQ.1.AND.AI(NEIBOUR(N),1).EQ.2) THEN
	      PA=AGA
		  PB=BGA
		  PLMD=LMDGA
		  PMU  =MUGA
	      PR=RGA
	      PS=SGA
	      PBT=BTGA
	      PN=NGA
		  PC=CGA
	      PD=DGA
	      PH=HGA

	    ELSE IF(AI(SNOFATOM,1).EQ.2.AND.AI(NEIBOUR(N),1).EQ.1) THEN
	      PA=AGA
		  PB=BGA
		  PLMD=LMDGA
		  PMU  =MUGA
	      PR=RGA
	      PS=SGA
	      PBT=BTGA
	      PN=NGA
		  PC=CGA
	      PD=DGA
	      PH=HGA

	    ELSE IF(AI(SNOFATOM,1).EQ.2.AND.AI(NEIBOUR(N),1).EQ.2) THEN
	      PA=AAA
		  PB=BAA
		  PLMD=LMDAA
		  PMU  =MUAA
	      PR=RAA
	      PS=SAA
	      PBT=BTAA
	      PN=NAA
		  PC=CAA
	      PD=DAA
	      PH=HAA


		
		ENDIF

	    B=(1+(PBT*KERCI)**PN)**(-.5/PN)
C		PRINT *,'PBT=',PBT
C		PRINT *,'KERCI=',KERCI
C		PRINT *,'PN=',PN
C		PRINT *,'B=',B
		SPIN=FCOFR(DISTANCE,PS,PR)*(PA*EXP(-PLMD*DISTANCE)
     C	  -B*PB*EXP(-PMU*DISTANCE))
C		PRINT *,'PA=',PA
C		PRINT *,'PLMD=',PLMD
C		PRINT *,'DISTANCE=',DISTANCE
C		PRINT *,'B=',B
C		PRINT *,'PB=',PB
CC		PRINT *,'PMU=',PMU
C		PRINT *,'SPIN=',SPIN
		SPOT=SPOT+SPIN



350	CONTINUE

	SHORTPOT=SPOT

	RETURN
      END
C------------------------------------------------------------------
C
C         GETNEIBOUR
C
C------------------------------------------------------------------

	SUBROUTINE GETNEIBOUR(AI,SNOFATOM)

	COMMON/NEIBOUR/ NEIBOUR

	PARAMETER(SAA=3.6,       SGG=3.6,       SGA=3.6      )

	PARAMETER(nofl=10)
C	NO OF ATOMS PER SIDE
      PARAMETER(CRYSTCONST=5.65315)


	REAL AI(8*NOFL*NOFL*NOFL,5)

	INTEGER SNOFATOM
	INTEGER NEIBOUR(20)
	INTEGER STARTNMB
		      
      DO 400 I=1,20
	  NEIBOUR(I)=0
400   CONTINUE
	
      STARTNMB=1
		
	DO 500 I=1,NOFL*NOFL*NOFL*8

c			CAN BE CHANGED TO THE DESIRED MOVABLE ATOMS.

	  IF(GETDISTANCE(AI,SNOFATOM,I).GT.0
     C  .AND.GETDISTANCE(AI,SNOFATOM,I).LT.SAA
     C  .AND.AI(I,1).GT.0) THEN	  	
	  
	  IF(STARTNMB.LE.20) THEN
	    NEIBOUR(STARTNMB)=I
		STARTNMB=STARTNMB+1
	  ELSE
	    PRINT *,'TOO MANY NEIBOUR !'
	  ENDIF
	  
	  ENDIF



500	CONTINUE


	RETURN
      END	
	

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C													   C
C	GETDISTANCE										   C
C													   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


	REAL FUNCTION GETDISTANCE(AI,ATOM1,ATOM2)

	PARAMETER(nofl=10)
C	NO OF ATOMS PER SIDE
      PARAMETER(CRYSTCONST=5.65315)


	REAL AI(8*NOFL*NOFL*NOFL,5)

	REAL DISOFXX,DISOFYY,DISOFZZ,HALFL	
	
      INTEGER ATOM1,ATOM2  
	HALFL=NOFL*CRYSTCONST*.5	  
	  
	  DISOFXX=AI(ATOM1,3)-AI(ATOM2,3)
	  DISOFYY=AI(ATOM1,4)-AI(ATOM2,4)	
	  DISOFZZ=AI(ATOM1,5)-AI(ATOM2,5)	
     	
C	  IF(DISOFXX.GT.HALFL)  DISOFXX=DISOFXX-HALFL*2
C	  IF(DISOFXX.LT.-HALFL) DISOFXX=DISOFXX+HALFL*2
	  IF(DISOFYY.GT.HALFL)  DISOFYY=DISOFYY-HALFL*2
	  IF(DISOFYY.LT.-HALFL) DISOFYY=DISOFYY+HALFL*2
	  IF(DISOFZZ.GT.HALFL)  DISOFZZ=DISOFZZ-HALFL*2
	  IF(DISOFZZ.LT.-HALFL) DISOFZZ=DISOFZZ+HALFL*2

	  GETDISTANCE=SQRT(DISOFXX*DISOFXX+
     C	  DISOFYY*DISOFYY+DISOFZZ*DISOFZZ)
     
      RETURN
      END


C     |---------------------------------------------------------------|
C     |      PREPARES A RANDOM PICKED PARTICLES LIST                  |
C     |---------------------------------------------------------------|
        SUBROUTINE PICKP(NP)
        COMMON/PICK/NPICK(800)
	COMMON/IIDUM/IIDUM	
        DIMENSION NCHEK(800)
        DO 610 I=1,NP
        NCHEK(I)=0
610      NPICK(I)=0
        NTOT=0
        DO 620 I=1,NP
C        CALL RANDOM(XX,YY,ZZ)
C        FLIP=FLOAT(NP)*ABS(XX)+0.4999
		FLIP=FLOAT(NP)*RAN2(IIDUM)+.4999
        IPICK=FLIP
        IF(IPICK.LE.0) GOTO 620
        IF(IPICK.GT.NP) GOTO 620
        IF(NCHEK(IPICK).EQ.0) THEN
           NTOT=NTOT+1
           NPICK(NTOT)=IPICK
           NCHEK(IPICK)=1
        ENDIF
620      CONTINUE
        DO 630 I=1,NP
        IF(NCHEK(I).EQ.0) THEN
           NTOT=NTOT+1
           NPICK(NTOT)=I
        ENDIF
630      CONTINUE
        IF(NTOT.NE.NP) THEN
        WRITE(6,640) NTOT
640     FORMAT(//5X,'NP.NE.NTOT',//5X,'NTOT =',I8//)
        STOP
        ENDIF
        RETURN
        END

      SUBROUTINE PMOVE(AI,TOTV,TOTE,BETA)
C     |--------------------------------------------------------------|
C     |     EVALUATE POTENTIAL  USING PAIRWISE ADDITIVE LENNARD-     |
C     |           JONES (6-12) INTERMOLECULAR POTENTIAL              |
C     |--------------------------------------------------------------|
	PARAMETER(nofl=10)

       COMMON/DISP/DAX(256),DAY(256),DAZ(256)
       COMMON/POS/X0(256),Y0(256),Z0(256)
       COMMON/PROP/SUME,SUMV,SUMA
       COMMON/PROPA/RMAX,CUBE,CUBE2,HMAX,DELTA,
     A        NP,IRR,IACPT
       COMMON/PICK/NPICK(800)
	COMMON/IIDUM/IIDUM	
      PARAMETER(CRYSTCONST=5.65315)


	REAL AI(8*NOFL*NOFL*NOFL,5)

	REAL HALFL	
	REAL POTOLD,POTNEW 
	

	REAL HMAX

	HALFL=NOFL*CRYSTCONST*.5	  
	HMAX = 10	
         TSUMA =0.
          TSUME=0.
          TSUMV=0.
C
C === OUTER LOOP OVER ATOMS
C
        DO 700 IX=1,NP
		
          I=NPICK(IX)


C	CONSIDER COVERAGE RATE
C	IF(AI(I,1).GT.0.AND.(MOD(I,8).EQ.1.OR.MOD(I,8).EQ.2)) THEN
	IF(AI(I,1).GT.0) THEN
C				   
C === STORE POSITION OF ATOM I
C
D          XI=X0(I)
D         YI=Y0(I)
D          ZI=Z0(I)
C
C === EVALUATE OLD ENERGY AND VIRIAL TERMS
C
          CALL CPOTENTIAL(AI,POTOLD,I)
	if(potold.lt.-50) print *,'potold',potold
C
C === RANDOM NUMBERS FOR MOVE
C
		XIOLD=AI(I,3)
		YIOLD=AI(I,4)
		ZIOLD=AI(I,5)
          
		
		
C		CALL RANDOM(XX,YY,ZZ)
C          XINEW = AI(I,3) + XX*DELTA
C          YINEW = AI(I,4) + YY*DELTA
C          ZINEW = AI(I,5) + ZZ*DELTA
C-----------------------------------------
C		ORIGINAL JUMP
C	XINEW = AI(I,3) + RAN2(IIDUM)*DELTA
C	YINEW = AI(I,4) + RAN2(IIDUM)*DELTA
C	ZINEW = AI(I,5) + RAN2(IIDUM)*DELTA
C-----------------------------------------
C		NEW JUMP

	XINEW = AI(I,3) + GETSIGN(RAN2(IIDUM))*DELTA
	YINEW = AI(I,4) + GETSIGN(RAN2(IIDUM))*DELTA
	ZINEW = AI(I,5) + GETSIGN(RAN2(IIDUM))*DELTA

c	XINEW = AI(I,3) + SIGN(1.,.5-RAN2(IIDUM))*DELTA
c	YINEW = AI(I,4) + SIGN(1.,.5-RAN2(IIDUM))*DELTA
c	ZINEW = AI(I,5) + SIGN(1.,.5-RAN2(IIDUM))*DELTA

C
C === CHECK FOR CONSERVATION OF PARTICLES IN PRIMARY CELL (P.B.C.)
C
          IF(YINEW.LE.0.) YINEW=YINEW+HALFL
          IF(YINEW.GE.HALFL*2) YINEW=YINEW-HALFL*2
          IF(ZINEW.LE.0.) ZINEW=ZINEW+HALFL
          IF(ZINEW.GE.HALFL*2) ZINEW=ZINEW-HALFL*2
C
C === CALCULATE THE PROPOSED ENERGY AND PROPERTY CHANGES
C
          
		AI(I,3)=XINEW
		AI(I,4)=YINEW
		AI(I,5)=ZINEW
		
				
		
		CALL CPOTENTIAL(AI,POTNEW,I)
      if(potnew .lt.-50)     print *,'potnew',potnew
		DELU = POTNEW - POTOLD
          DELV = VIRNEW - VIROLD
C
C === REJECT OR ACCEPT MOVE
C
          BETAU = DELU*BETA
          IF(BETAU.GT.HMAX) GO TO 699
          IF(DELU.LE.0.0) then 
		PROBL=1.0
          else
		PROBL=EXP(-BETAU)
		endif
C          CALL RANDOM(PX,PY,PZ)
C          PZ=ABS(PZ)
	PZ=RAN2(IIDUM)
          IF(PROBL.LT.PZ) GO TO 699
C
C ===  IF MOVE ACCEPTED, PROPOSED CONFIGURATION BECOMES THE NEW ONE
C
          IACPT = IACPT + 1

          TOTE = TOTE + DELU
D          TOTV = TOTV + DELV
          TSUMA = TSUMA + PROBL
          TSUME = TSUME + PROBL*TOTE
D          TSUMV = TSUMV + PROBL*TOTV
D          DAX(I) = DAX(I) - XX*DELTA
D          DAY(I) = DAY(I) - YY*DELTA
D          DAZ(I) = DAZ(I) - ZZ*DELTA
          GO TO 700
C
C ===  IF MOVE REJECTED, OLD CONFIGURATION BECOMES NEW CONFIGURATION
C
699       TSUME = TSUME + TOTE
D          TSUMV = TSUMV + TOTV
          TSUMA=TSUMA+1.
		AI(I,3)=XIOLD
		AI(I,4)=YIOLD
		AI(I,5)=ZIOLD
  
  
	ENDIF  
  700   CONTINUE
C
C === UPDATE SUM ACCUMULATORS
C
          SUME = SUME + TSUME
          SUMV = SUMV + TSUMV
          SUMA = SUMA + TSUMA
c		PRINT*,TSUME

	  RETURN
        END


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



C		Great Random Generator!



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      FUNCTION ran2(idum)
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL ran2,AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 711 j=NTAB+8,1,-1

          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
711      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END
	
     
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C		TOTAL POTENTIAL					 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC     
      SUBROUTINE TOTALPOT(AI,NP)
      PARAMETER(nofl=10)
C      COMMON/PROPA/RMAX,CUBE,CUBE2,HMAX,DELTA,
C     A        NP,IRR,IACPT
	COMMON/TOTPOT/TOTPOT
      REAL AI(8*NOFL*NOFL*NOFL,5)
	REAL IPOT
	TOTPOT=0
	NOFEXIST=0
      DO 800 I=1,NP
	IF(AI(I,1).NE.0) THEN
	NOFEXIST=NOFEXIST+1
      CALL CPOTENTIAL(AI,IPOT,I)
	ENDIF
      TOTPOT=TOTPOT+IPOT
800	CONTINUE
      WRITE (*,*) NOFEXIST,'EXIST','TOTPOT=',TOTPOT
	RETURN
	END


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C											   C
C											   C
C		GIVE THE POSITION OF ATOMS             C
C											   C
C											   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	SUBROUTINE SHOWPOS(AI,NOFINTEREST,NOOFFILE)
      PARAMETER(nofl=10)
      REAL AI(8*NOFL*NOFL*NOFL,5)
	INTEGER NOFINTEREST
	CHARACTER*5 OFILENAME
	INTEGER NOOFFILE
	CALL CONVERT(NOOFFILE,OFILENAME)
	OPEN(3,FILE='POS'//OFILENAME//'.DAT')
	DO 900 I=1,NOFINTEREST
C	IF((MOD(I,8).EQ.1.OR.MOD(I,8).EQ.2)) THEN
	IF(1) THEN
      WRITE (3,*) 
     C		I,AI(I,1),AI(I,3),AI(I,4),AI(I,5)
      ENDIF
900	CONTINUE

   		CLOSE(3)
	PRINT *,OFILENAME,'IS WRITTEN',NOOFFILE

     
      RETURN
	END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C												   C
C	CONVERT FILENAME FROM INT TO CHAR			   C
C												   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	subroutine convert(input,filename)
	character*5 filename
	integer input,i
	i=input
	do 911 j=1,5
	k=mod(i,10)
	i=(i-k)/10
	filename((6-j):(6-j))=char(k+48)
911	continue
	return
	end


	FUNCTION GETSIGN(RANNO)
	REAL RANNO,GETSIGN
	IF(RANNO.LE.1/3.) THEN
	GETSIGN=1
	ELSE IF((RANNO.GT.1./3.).AND.(RANNO.LE.2/3.)) THEN
	GETSIGN=0
	ELSE IF(RANNO.GT.2/3.) THEN
	GETSIGN=-1
	ENDIF
	RETURN
	END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC	
C														C
C		STSURFACE SUBROUTINE							C
C														C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C	TOPTYPE:1(Ga);2(As)
C	DISTYPE:0(RANDOM);1(REGULAR)
C
	SUBROUTINE STSURFACE(AI,COVERAGE,TOPTYPE,DISTYPE)

      PARAMETER(CRYSTCONST=5.65315)
      PARAMETER(nofl=10)
      REAL AI(8*NOFL*NOFL*NOFL,5)
	REAL COVERAGE
	INTEGER TOPTYPE,DISTYPE


	IF(TOPTYPE.EQ.1) THEN
		IF(DISTYPE.EQ.0) THEN
			DO 1150 I=1,NOFL*NOFL*8
			IF(AI(I,3).EQ.0) THEN
				IF(RAN2(IIDUM).GT.COVERAGE) AI(I,1)=0
			ENDIF
1150			CONTINUE

	
		ELSEIF(DISTYPE.EQ.1) THEN
		do 1160 i=1,NOFL*NOFL*8
		if(ai(i,3).eq.0) then
		if((ai(i,4).gt.sqrt(coverage)*nofl*crystconst)
     c	 .or.(ai(i,5).gt.sqrt(coverage)*nofl*crystconst)) ai(i,1)=0
		endif
1160		continue
			
		ELSE
		PRINT *,'ERR! IN SETUP'
		CALL EXIT(1)
		ENDIF

	ELSEIF(TOPTYPE.EQ.2) THEN
		
		do 1165 i=1,NOFL*NOFL*8
		if(ai(i,3).eq.0) then
		ai(i,1)=0
		endif
1165		continue


		IF(DISTYPE.EQ.0) THEN
			DO 1170 I=1,NOFL*NOFL*8
			if(mod(i,8).eq.0.or.mod(i,8).eq.5) then
				IF(RAN2(IIDUM).GT.COVERAGE) AI(I,1)=0
			ENDIF
1170			CONTINUE

	
		ELSEIF(DISTYPE.EQ.1) THEN
		do 1180 i=1,NOFL*NOFL*8
			if(mod(i,8).eq.0.or.mod(i,8).eq.5) then
		if((ai(i,4).gt.sqrt(coverage)*nofl*crystconst)
     c	 .or.(ai(i,5).gt.sqrt(coverage)*nofl*crystconst)) ai(i,1)=0
		endif
1180		continue
			
		ELSE
		PRINT *,'ERR! IN SETUP'
		CALL EXIT(1)
		ENDIF
	ELSE
		PRINT *,'ERR! IN SETUP'
		CALL EXIT(1)

	ENDIF

























C	 MAIN LATTICE SETUP COMPLETE
C	WRITE(*,*) AI(65,3),AI(65,4),AI(65,5)

C	 REMOVE THE SURFACE LATTICE TO THE COVERAGE RATE.
c	DO 150 I=1,800
c	IF(AI(I,3).EQ.0) THEN
c	 IF(RAN2(IIDUM).GT.COVERAGE) AI(I,1)=0
c	ENDIF
c150	CONTINUE

C	DO 150 I=1,NOFL*NOFL*8
C	IF(mod(i,8).eq.0.or.mod(i,8).eq.5) THEN
C	 IF(RAN2(IIDUM).GT.COVERAGE) AI(I,1)=0
C	ENDIF
C150	CONTINUE

c	do 150 i=1,800
c	if(ai(i,3).eq.0) then
c	 if((ai(i,4).gt.coverage*nofl*crystconst)
c     c		 .or.(ai(i,5).gt.coverage*nofl*crystconst)) ai(i,1)=0
c	endif
c150	continue

c	do 150 i=1,200
c	kkk=ai(i,1)
c	if(mod(i,8).eq.0.or.mod(i,8).eq.5) then
c	 if((ai(i,4).gt.coverage*nofl*crystconst)
c     c		 .or.(ai(i,5).gt.coverage*nofl*crystconst)) then 
c      ai(i,1)=0
c	endif
c      endif
c150	continue










	RETURN
	END