*     MDBNCH                                F.ERCOLESSI  17-DEC-1988
*                                                    REV 17-MAR-1990
*                                                    REV 17-DEC-1992
*                                                    REV  9-NOV-1993
*                                                    REV  2-NOV-1994
*                                                    REV 30-NOV-1994
*
*     (ALL REVS ARE JUST COSMETIC CLEANUPS, MOSTLY TO IMPROVE STANDARD
*      CONFORMANCE.  THE PROGRAM DOES THE EXACT SAME THING SINCE THE 
*      17-DEC-1988 RELEASE).
*
*     MDBNCH IS A MOLECULAR DYNAMICS BENCHMARK.
*     THE SYSTEM SIMULATED IS GOLD, USING A MANY-BODY 'GLUE'
*     INTERACTION POTENTIAL. THREE DIFFERENT NUMBER OF PARTICLES
*     ARE USED: 256, 2048 AND 16384.
*
*     THE BENCHMARK DOES NOT REQUIRE ANY INPUT, AND CAN BE RUN
*     SIMPLY BY COMPILE-LINK-GO.
*     IT WRITES ITS RESULTS ON THE 'STANDARD OUTPUT', I.E. USING
*     FORTRAN PRINT STATEMENTS.  NO OTHER OUTPUT IS PRODUCED.
*     THAT'S ALL I/O: NO SCRATCH FILES ARE USED.
*     
*     THE BENCHMARK IS SELF-CONTAINED: NO EXTERNAL ROUTINES ARE 
*     NEEDED, WITH THE FOLLOWING EXCEPTION. 
*     A
*                DOUBLE PRECISION FUNCTION SECOND()
*
*     SHOULD BE MADE AVAILABLE AT LINK TIME FOR TIMING PURPOSES.
*
*     THE BENCHMARK CONTAINS A BLOCK DATA: MAKE SURE IT IS 
*     APPROPRIATELY PROCESSED BY YOUR COMPILER AND/OR LINKER.
*
*     THE BENCHMARK IS INTENDED TO BE RUN USING 64-BIT PRECISION FOR
*     ALL REAL VARIABLES.  TO THIS END, ALL REAL VARIABLES AND
*     CONSTANTS ARE DECLARED 'DOUBLE PRECISION', MEANING THAT
*     WE HAVE 32-BIT MACHINES IN MIND AS A TARGET.  IF YOU ARE GOING 
*     TO RUN IT ON 64-BIT MACHINES, MAKE SURE TO SPECIFY THE COMPILER 
*     OPTION TO TREAT DOUBLE PRECISION AS SINGLE.
*
*     PLEASE SEND ALL RESULTS TO FURIO@SISSA.IT, INCLUDING OUTPUT,
*     EXACT MACHINE TYPE, OPERATING SYSTEM AND COMPILER RELEASE,
*     COMPILER OPTIONS, AND ANY OTHER USEFUL INFORMATION.   THANKS.
*
*     THERE IS AN OFFICIAL MDBNCH WWW PAGE AT THE URL
*                   HTTP://WWW.SISSA.IT/FURIO/MDBNCH.HTML
*     (EVERYTHING IN LOWER CASE, THAT CANNOT BE USED HERE :-)
*     ESTABLISHED ON NOVEMBER 3, 1994.
*     THE READER IS REFERRED THERE FOR FURTHER INFORMATIONS.
*
*     NO HUMANS ALLOWED BEYOND THIS POINT.
*
      PROGRAM MDBNCH
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/CONST/TWOPI,BOLTZ
      COMMON/CNTRL/NTABLE,LISTEM,LMETHD,
     $NDIFF,NSTAT,NGOFR,NTRAJ,IVDUMP,ILIN,
     $LOCK(3,3),LILOCK(2,14),NFREED,LOCKCM,
     $NSCALE,ITPART,ITWALL
      COMMON/COUNT/NFI,LCOUNT,LISTER,KNTSTA,KNTGOR,LEP,MANYON
      COMMON/DEN/RNRHO,RCRHO,RNRHO2,RCRHO2,DR2RHO,RHO(KR),DRHO(KR)
      COMMON/GEAR/F02,F12,F32,F42,F52
      COMMON/GLUE/DMIN,DMAX,DD,UJ(KG),DUJ(KG)
      COMMON/IDENT/ELEMEN,REF,TODAY,NOW
      CHARACTER ELEMEN*10,REF*72,TODAY*10,NOW*10
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/LSTUPD/RLIST(3,NM)
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/PARAM/DELTA,DELTA2,GAMMA,VSCALE,CTRLCE
     $,CTRLMI,CTRLMA,RSQUPD,RANSQ,VMAS,BOX(3,3)
      COMMON/POT/RN,RC,RN2,RC2,DR2,VJ(KP),FJ(KP)
      COMMON/PRESS/PEXT,PAI(3,3)
      COMMON/PRINT/LNGPRT,IPRIND
      COMMON/SCRATC/DUMMY1(NM),DUMMY2(NM),DUMMY3(NM),DUMMY4(NM)
      COMMON/STATIS/FGS(NG),GRANG,FACNG,SCABY2,RESZ,DONTR,FONTR,SIG2,
     $NGS(NG),NGMAX,NZHIGH,NZLOW,MULTIP
      COMMON/SUMS/TEMPSM,TEMWSM,EKINSM,POT2SM,PGLUSM,POSTSM
     $,TOTESM,DENSSM,ALSM,VOLUSM,AREASM,HEIGSM
      COMMON/THRU/ATMASS,ECOH,R0,SPAREF
      COMMON/TIMERS/TIMSTR,TIMFIX,TIMBLD,TIMFRC,TIMINT
      EXTERNAL SECOND
      PRINT '(/,1X,2A)',
     $'     MDBNCH: A MOLECULAR DYNAMICS BENCHMARK, VERSION ',
     $'OF DECEMBER 17, 1988'
      TIMALL=SECOND()
      NBSIZE=4
      CALL MTE(NBSIZE)
      NSTEPS=1000
      NLIST=10
      METHOD=0
      SKIN=1.0D0
      NCORR=0
      NPRINT=100
      CALL MASTER(NSTEPS,NLIST,METHOD,SKIN,NCORR,NPRINT)
      NBSIZE=8
      CALL MTE(NBSIZE)
      NSTEPS=10
      NLIST=10
      METHOD=1
      SKIN=1.0D0
      NCORR=0
      NPRINT=1
      CALL MASTER(NSTEPS,NLIST,METHOD,SKIN,NCORR,NPRINT)
      NBSIZE=8
      CALL MTE(NBSIZE)
      NSTEPS=10
      NLIST=10
      METHOD=0
      SKIN=1.0D0
      NCORR=0
      NPRINT=1
      CALL MASTER(NSTEPS,NLIST,METHOD,SKIN,NCORR,NPRINT)
      NBSIZE=8
      CALL MTE(NBSIZE)
      NSTEPS=10
      NLIST=10
      METHOD=0
      SKIN=1.0D0
      NCORR=10
      NPRINT=1
      CALL MASTER(NSTEPS,NLIST,METHOD,SKIN,NCORR,NPRINT)
      NBSIZE=8
      CALL MTE(NBSIZE)
      NSTEPS=10
      NLIST=10
      METHOD=1
      SKIN=1.5D0
      NCORR=0
      NPRINT=1
      CALL MASTER(NSTEPS,NLIST,METHOD,SKIN,NCORR,NPRINT)
      NBSIZE=8
      CALL MTE(NBSIZE)
      NSTEPS=10
      NLIST=5
      METHOD=1
      SKIN=0.3D0
      NCORR=0
      NPRINT=1
      CALL MASTER(NSTEPS,NLIST,METHOD,SKIN,NCORR,NPRINT)
      NBSIZE=16
      CALL MTE(NBSIZE)
      NSTEPS=10
      NLIST=5
      METHOD=1
      SKIN=0.5D0
      NCORR=0
      NPRINT=1
      CALL MASTER(NSTEPS,NLIST,METHOD,SKIN,NCORR,NPRINT)
      TIMALL=SECOND()-TIMALL
      PRINT '(/,1X,79(''*''),/,1X,A,F12.6,A,/)',
     $'COMPLETE BENCHMARK EXECUTION TIME :',
     $TIMALL,' CP SECONDS.'
      END
      SUBROUTINE MASTER(NSTEPS,NLIST,METHOD,SKIN,NCORR,NPRINT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/CONST/TWOPI,BOLTZ
      COMMON/CNTRL/NTABLE,LISTEM,LMETHD,
     $NDIFF,NSTAT,NGOFR,NTRAJ,IVDUMP,ILIN,
     $LOCK(3,3),LILOCK(2,14),NFREED,LOCKCM,
     $NSCALE,ITPART,ITWALL
      COMMON/COUNT/NFI,LCOUNT,LISTER,KNTSTA,KNTGOR,LEP,MANYON
      COMMON/DEN/RNRHO,RCRHO,RNRHO2,RCRHO2,DR2RHO,RHO(KR),DRHO(KR)
      COMMON/GEAR/F02,F12,F32,F42,F52
      COMMON/GLUE/DMIN,DMAX,DD,UJ(KG),DUJ(KG)
      COMMON/IDENT/ELEMEN,REF,TODAY,NOW
      CHARACTER ELEMEN*10,REF*72,TODAY*10,NOW*10
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/LSTUPD/RLIST(3,NM)
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/PARAM/DELTA,DELTA2,GAMMA,VSCALE,CTRLCE
     $,CTRLMI,CTRLMA,RSQUPD,RANSQ,VMAS,BOX(3,3)
      COMMON/PBCS/HALF,PBCX,PBCY,PBCZ
      COMMON/POT/RN,RC,RN2,RC2,DR2,VJ(KP),FJ(KP)
      COMMON/PRESS/PEXT,PAI(3,3)
      COMMON/PRINT/LNGPRT,IPRIND
      COMMON/SCRATC/DUMMY1(NM),DUMMY2(NM),DUMMY3(NM),DUMMY4(NM)
      COMMON/STATIS/FGS(NG),GRANG,FACNG,SCABY2,RESZ,DONTR,FONTR,SIG2,
     $NGS(NG),NGMAX,NZHIGH,NZLOW,MULTIP
      COMMON/SUMS/TEMPSM,TEMWSM,EKINSM,POT2SM,PGLUSM,POSTSM
     $,TOTESM,DENSSM,ALSM,VOLUSM,AREASM,HEIGSM
      COMMON/THRU/ATMASS,ECOH,R0,SPAREF
      COMMON/TIMERS/TIMSTR,TIMFIX,TIMBLD,TIMFRC,TIMINT
      EXTERNAL SECOND
      PRINT '(/,1X,79(''*''),/)'
      PRINT '(1X,A,I6,A,I5,A)',
     $'MD BENCHMARK FOR',
     $MOLSA,' PARTICLES,',
     $NSTEPS,' STEPS.'
      IF(METHOD.EQ.0)THEN
      PRINT '(1X,A,I3,A,F5.2)',
     $'O(N**2) BRUTE FORCE LIST FORMATION EVERY',NLIST,
     $' WITH SKIN =',SKIN
      ELSE
      PRINT '(1X,A,I3,A,F5.2)',
     $'O(N) CELL-METHOD LIST FORMATION EVERY',NLIST,
     $' WITH SKIN =',SKIN
      ENDIF
      IF(NCORR.EQ.0)THEN
      PRINT '(1X,A)',
     $'PAIR CORRELATION FUNCTION NOT COMPUTED'
      ELSE
      PRINT '(1X,A,I3,A)',
     $'PAIR CORRELATION FUNCTION COMPUTED EVERY',NCORR,
     $' STEPS'
      ENDIF
      CALL MINIT(NLIST,METHOD,SKIN,NCORR)
      DO 1 ISTEP=1,NSTEPS
      NFI=NFI+1
      IF((MOD(ISTEP,NPRINT).EQ.0).OR.
     $(ISTEP.EQ.1).OR.(ISTEP.EQ.NSTEPS))THEN
      LNGPRT=1
      ELSE
      LNGPRT=0
      ENDIF
      CALL MSTEP
1     CONTINUE
      TIMSTP=SECOND()
      TIMCPU=TIMSTP-TIMSTR
      TIMSTE=TIMCPU-TIMFIX
      PRINT '(/,1X,I5,A,I5,A)',
     $NSTEPS,' TIME STEPS,',LCOUNT,' LIST UPDATES'
      IF(TIMFIX.NE.0.D0)
     $PRINT '(1X,F11.6,A)',TIMFIX,
     $' SEC. CP TIME FOR INITIALIZATION'
      IF(TIMBLD.NE.0.D0)
     $PRINT '(1X,F11.6,A,F10.6,A)',TIMBLD,
     $' SEC. CP TIME UPDATING THE LIST        (',
     $TIMBLD/LCOUNT,' SEC/UPD. )'
      IF(TIMFRC.NE.0.D0)
     $PRINT '(1X,F11.6,A,F10.6,A)',TIMFRC,
     $' SEC. CP TIME CALCULATING FORCES       (',
     $TIMFRC/NSTEPS,' SEC/STEP )'
      IF(TIMINT.NE.0.D0)
     $PRINT '(1X,F11.6,A,F10.6,A)',TIMINT,
     $' SEC. CP TIME FOR TIME INTEGRATION     (',
     $TIMINT/NSTEPS,' SEC/STEP )'
      TIMOTH=TIMCPU-(TIMFIX+TIMBLD+TIMFRC+TIMINT)
      IF(TIMOTH.NE.0.D0)
     $PRINT '(1X,F11.6,A,F10.6,A)',TIMOTH,
     $' SEC. CP TIME FOR OTHER TASKS          (',
     $TIMOTH/NSTEPS,' SEC/STEP )'
      IF(TIMSTE.NE.0.D0)
     $PRINT '(1X,F11.6,A,F10.6,A)',TIMSTE,
     $' SEC. CP TIME EXCLUDING INITIALIZATION (',
     $TIMSTE/NSTEPS,' SEC/STEP )'
      IF(TIMCPU.NE.0.D0)
     $PRINT '(1X,F11.6,A)',TIMCPU,
     $' SEC. TOTAL CP TIME'
      RETURN
      END
      SUBROUTINE MINIT(NLIST,METHOD,SKIN,NCORR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/CONST/TWOPI,BOLTZ
      COMMON/CNTRL/NTABLE,LISTEM,LMETHD,
     $NDIFF,NSTAT,NGOFR,NTRAJ,IVDUMP,ILIN,
     $LOCK(3,3),LILOCK(2,14),NFREED,LOCKCM,
     $NSCALE,ITPART,ITWALL
      COMMON/COUNT/NFI,LCOUNT,LISTER,KNTSTA,KNTGOR,LEP,MANYON
      COMMON/DEN/RNRHO,RCRHO,RNRHO2,RCRHO2,DR2RHO,RHO(KR),DRHO(KR)
      COMMON/GEAR/F02,F12,F32,F42,F52
      COMMON/GLUE/DMIN,DMAX,DD,UJ(KG),DUJ(KG)
      COMMON/IDENT/ELEMEN,REF,TODAY,NOW
      CHARACTER ELEMEN*10,REF*72,TODAY*10,NOW*10
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/PARAM/DELTA,DELTA2,GAMMA,VSCALE,CTRLCE
     $,CTRLMI,CTRLMA,RSQUPD,RANSQ,VMAS,BOX(3,3)
      COMMON/PBCS/HALF,PBCX,PBCY,PBCZ
      COMMON/POT/RN,RC,RN2,RC2,DR2,VJ(KP),FJ(KP)
      COMMON/PRESS/PEXT,PAI(3,3)
      COMMON/PRINT/LNGPRT,IPRIND
      COMMON/SCRATC/DUMMY1(NM),DUMMY2(NM),DUMMY3(NM),DUMMY4(NM)
      COMMON/STATIS/FGS(NG),GRANG,FACNG,SCABY2,RESZ,DONTR,FONTR,SIG2,
     $NGS(NG),NGMAX,NZHIGH,NZLOW,MULTIP
      COMMON/SUMS/TEMPSM,TEMWSM,EKINSM,POT2SM,PGLUSM,POSTSM
     $,TOTESM,DENSSM,ALSM,VOLUSM,AREASM,HEIGSM
      COMMON/THRU/ATMASS,ECOH,R0,SPAREF
      COMMON/TIMERS/TIMSTR,TIMFIX,TIMBLD,TIMFRC,TIMINT
      DIMENSION H(3,3),H1(3,3),H2(3,3),H3(3,3),H4(3,3),H5(3,3),HIN(3,3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      EQUIVALENCE(X(1,-2,1),H1(1,1))
      EQUIVALENCE(X(1,-2,2),H2(1,1))
      EQUIVALENCE(X(1,-2,3),H3(1,1))
      EQUIVALENCE(X(1,-2,4),H4(1,1))
      EQUIVALENCE(X(1,-2,5),H5(1,1))
      EQUIVALENCE(XIN(1,-2),HIN(1,1))
      DIMENSION HINPUT(3,3)
      CHARACTER OLDREF*72,YSTRDY*10,BEFORE*10
      EXTERNAL SECOND
      CHARACTER*10 VERS
      PARAMETER(VERS='62.2      ')
      PARAMETER(PI=3.141592653589793D0)
      DATA CF/1.03641D-28/
2     FORMAT(1X,A)
3     FORMAT(/,1X,A)
1     FORMAT('1')
      TIMSTR=SECOND()
      TIMFIX=0.D0
      TIMBLD=0.D0
      TIMFRC=0.D0
      TIMINT=0.D0
      TWOPI=PI+PI
      BOLTZ=11606.D0
      OLDREF=REF
      YSTRDY=TODAY
      BEFORE=NOW
      DLTOLD=DELTA
      CALL RESET(X,3*(NM+3)*5)
      CALL IRESET(NGS,NGMAX)
      CALL RESET(FGS,NGMAX)
      MANYON=1
      IBEG=0
      NTRAJ=0
      NSTAT=0
      NGOFR=NCORR
      NDIFF=0
      IPRIND=0
      LNGPRT=0
      NTABLE=0
      LISTEM=NLIST
      DELTA=0.05D0
      NSCALE=0
      VSCALE=1.D0
      ITPART=0
      ITWALL=0
      CENTRE=-3.0000D0
      TOLLER=0.0002D0
      IHCHGE=1
      SPACNG=4.15D0
      LOCKCM=0
      LOCK(1,1)=-2
      LOCK(2,1)=-1
      LOCK(3,1)=-1
      LOCK(1,2)=-1
      LOCK(2,2)=-2
      LOCK(3,2)=-1
      LOCK(1,3)=-1
      LOCK(2,3)=-1
      LOCK(3,3)=-2
      VMAS=0.5D0
      IWDAMP=0
      GAMMA=0.D0
      PEXT=0.D0
      DO 8086 I=1,3
      PAI(1,I)=0.D0
      PAI(2,I)=0.D0
      PAI(3,I)=0.D0
8086  CONTINUE
      ILIN=0
      IF(NTRAJ.EQ.0)THEN
      IVDUMP=0
      ELSE IF(NTRAJ.GT.0)THEN
      IVDUMP=1
      ELSE
      IVDUMP=0
      NTRAJ=-NTRAJ
      ENDIF
      SKIN=ABS(SKIN)
      LMETHD=METHOD
      IF(LISTEM.EQ.1)SKIN=0.D0
      IF(NSCALE.LE.0)VSCALE=1.D0
      IF(IWDAMP.EQ.0)GAMMA=0.D0
      IF(TOLLER.LT.0.D0)TOLLER=-TOLLER
      LISTER=LISTEM
      LCOUNT=0
      RSQUPD=(0.5D0*SKIN)**2
      F02=3.D0/16.D0
      F12=251.D0/360.D0
      F32=11.D0/18.D0
      F42=1.D0/6.D0
      F52=1.D0/60.D0
      DELTA2=DELTA**2
      CTRLCE=CENTRE
      CTRLMI=CENTRE-TOLLER
      CTRLMA=CENTRE+TOLLER
      LPBCSM=LPBC(1)+LPBC(2)+LPBC(3)
      HALF=0.5D0
      PBCX=LPBC(1)
      PBCY=LPBC(2)
      PBCZ=LPBC(3)
      IF(LPBCSM.GT.1)THEN
      SPCING=R0*BOX(3,3)/NPLA
      ELSE
      SPCING=R0
      ENDIF
      SIG2=1.D0/(2.D0*(SPCING/2.D0)**2)
      L=LPBC(1)+LPBC(2)+LPBC(3)
      CALL POTENT
      IF(MANYON.NE.0)THEN
      CALL DENSIT
      CALL ELGLUE
      RANSQ=(MAX(RC,RCRHO)+SKIN)**2
      ELSE
      RANSQ=(RC+SKIN)**2
      ENDIF
      IF((DELTA.NE.DLTOLD).AND.(DLTOLD.GT.0.D0))THEN
      RATIO=DELTA/DLTOLD
      RATIO2=RATIO**2
      RATIO4=RATIO2**2
      DO 1091 K=1,3
      DO 1092 I=-2,MOLSP
      X(K,I,1)=X(K,I,1)*RATIO
      X(K,I,2)=X(K,I,2)*RATIO2
      X(K,I,3)=X(K,I,3)*RATIO2*RATIO
      X(K,I,4)=X(K,I,4)*RATIO4
      X(K,I,5)=X(K,I,5)*RATIO4*RATIO
1092  CONTINUE
1091  CONTINUE
      ENDIF
      IF(IHCHGE.NE.0)THEN
      R0NEW=R0*SPACNG/SPAREF
      DO 1096 K=1,3
      DO 1096 J=1,3
1096  HINPUT(J,K)=R0NEW*BOX(J,K)
      DO 1093 K=1,3
      DO 1093 J=1,3
1093  IF(HINPUT(J,K).NE.H(J,K))GOTO 1094
      IHCHGE=0
1094  CONTINUE
      ENDIF
      IF(IHCHGE.NE.0)THEN
      DO 1095 K=1,3
      DO 1095 J=1,3
1095  H(J,K)=HINPUT(J,K)
      ENDIF
      MALOCK=0
      DO 1100 K=1,3
      DO 1100 J=1,3
      IF(LOCK(J,K).GT.MALOCK)THEN
      MALOCK=LOCK(J,K)
      ELSE IF(LOCK(J,K).LT.0)THEN
      H1(J,K)=0.D0
      H2(J,K)=0.D0
      H3(J,K)=0.D0
      H4(J,K)=0.D0
      H5(J,K)=0.D0
      IF(LOCK(J,K).EQ.-1)THEN
      H(J,K)=0.D0
      HIN(J,K)=0.D0
      IF(J.EQ.K)
     $PRINT 2,'MINIT: A DIAGONAL ELEMENT OF H IS ZERO.'
      ENDIF
      ENDIF
1100  CONTINUE
      L=1
      DO 1120 I=1,MALOCK
      M=0
      DO 1130 K=1,3
      DO 1130 J=1,3
      IF(LOCK(J,K).EQ.I)THEN
      M=M+1
      LILOCK(1,L+M)=J
      LILOCK(2,L+M)=K
      ENDIF
1130  CONTINUE
      IF(M.GT.1)THEN
      LILOCK(1,L)=M
      LILOCK(2,L)=0
      L=L+M+1
      ELSE IF(M.EQ.1)THEN
      LOCK(LILOCK(1,L+M),LILOCK(2,L+M))=0
      PRINT '(1X,A,I3)',
     $'MINIT: SINGLE ELEMENT IN LOCK CLASS NR.',I
      ELSE
      ENDIF
1120  CONTINUE
      LILOCK(1,L)=0
      LILOCK(2,L)=0
      IF(L.GT.14)PRINT '(/,1X,A,I3)',
     $'MINIT: LILOCK OVERFLOW, INDEX IS',L
      CALL SYMM(H1)
      CALL SYMM(H2)
      CALL SYMM(H3)
      CALL SYMM(H4)
      CALL SYMM(H5)
      NFREED=0
      DO 101 J=1,3
      DO 102 I=1,3
      IF(LOCK(I,J).EQ.0)NFREED=NFREED+1
102   CONTINUE
101   CONTINUE
      L=1
200   CONTINUE
      IF(LILOCK(1,L).NE.0)THEN
      NFREED=NFREED+1
      L=L+LILOCK(1,L)+1
      GOTO 200
      ENDIF
      IF(IBEG.EQ.0)THEN
      DO 98 I=-2,MOLSA
      XIN(1,I)=X0(1,I)
      XIN(2,I)=X0(2,I)
      XIN(3,I)=X0(3,I)
98    CONTINUE
      CALL RESET(FGS,NGMAX)
      CALL IRESET(NGS,NGMAX)
      NFI=0
      KNTGOR=0
      KNTSTA=0
      TEMPSM=0.D0
      TEMWSM=0.D0
      EKINSM=0.D0
      POT2SM=0.D0
      PGLUSM=0.D0
      POSTSM=0.D0
      TOTESM=0.D0
      DENSSM=0.D0
      ALSM=0.D0
      VOLUSM=0.D0
      AREASM=0.D0
      HEIGSM=0.D0
      ENDIF
      PRINT '(/,1X,2A)',
     $' STEP LP  KIN.E   POT.E   TOT.E ',
     $'  DIFFUS     PX       PY       PZ   '
      PRINT '(1X,2A)',
     $' ---- -- ------- ------- -------',
     $' -------- -------- -------- --------'
      TIMFIX=TIMFIX+(SECOND()-TIMSTR)
      RETURN
      END
      SUBROUTINE SYMM(HN)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION HN(3,3)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/CNTRL/NTABLE,LISTEM,LMETHD,
     $NDIFF,NSTAT,NGOFR,NTRAJ,IVDUMP,ILIN,
     $LOCK(3,3),LILOCK(2,14),NFREED,LOCKCM,
     $NSCALE,ITPART,ITWALL
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      DIMENSION H(3,3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      L=1
1     CONTINUE
      IF(LILOCK(1,L).NE.0)THEN
      SUM=0.D0
      DO 2 M=1,LILOCK(1,L)
      I=LILOCK(1,L+M)
      J=LILOCK(2,L+M)
      IF(H(I,J).NE.0.D0)THEN
      SUM=SUM+HN(I,J)/H(I,J)
      ELSE
      SUM=SUM+HN(I,J)
      ENDIF
2     CONTINUE
      SUM=SUM/LILOCK(1,L)
      DO 3 M=1,LILOCK(1,L)
      I=LILOCK(1,L+M)
      J=LILOCK(2,L+M)
      IF(H(I,J).NE.0.D0)THEN
      HN(I,J)=SUM*H(I,J)
      ELSE
      HN(I,J)=SUM
      ENDIF
3     CONTINUE
      L=L+LILOCK(1,L)+1
      GOTO 1
      ENDIF
      RETURN
      END
      SUBROUTINE POTENT
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/POT/RN,RC,RN2,RC2,DR2,VJ(KP),FJ(KP)
      COMMON/THRU/ATMASS,ECOH,R0,SPAREF
      CHARACTER*72 TITLE
      TITLE='TWO-BODY POTENTIAL'
      NINT=KP
      RN=1.69D0
      RC=3.7D0
      R0=0.2878207442D+01
      RHARD=RN
      RN2=RN**2
      RC2=RC**2
      DR2=(RC2-RN2)/(NINT-1)
      DO 3 I=1,KP
      RSQ=RN2+(I-1)*DR2
      R=SQRT(RSQ)
      CALL POTFUN(R , PHI , DPHI , D2PHI)
      VJ(I)=PHI
      FJ(I)=-DPHI/R
3     CONTINUE
      RETURN
      END
      SUBROUTINE DENSIT
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/DEN/RNRHO,RCRHO,RNRHO2,RCRHO2,DR2RHO,RHO(KR),DRHO(KR)
      CHARACTER*72 TITLE
      TITLE='DENSITY FUNCTION'
      NINT=KR
      RNRHO=1.69D0
      RCRHO=3.9D0
      RNRHO2=RNRHO**2
      RCRHO2=RCRHO**2
      DR2RHO=(RCRHO2-RNRHO2)/(NINT-1)
      DO 13 I=1,KR
      RSQ=RNRHO2+(I-1)*DR2RHO
      R=SQRT(RSQ)
      CALL DENFUN(R , RH , DRH , D2RH)
      RHO(I)=RH
      DRHO(I)=-DRH/R
13    CONTINUE
      RETURN
      END
      SUBROUTINE ELGLUE
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/GLUE/DMIN,DMAX,DD,UJ(KG),DUJ(KG)
      CHARACTER*72 TITLE
      TITLE='GLUE'
      NINT=KG
      DMIN=0.D0
      DMAX=20.D0
      DD=(DMAX-DMIN)/(NINT-1)
      DO 4 I=1,KG
      DENS=DMIN+(I-1)*DD
      CALL GLUFUN(DENS , U0 , U1 , U2)
      UJ(I)=U0
      DUJ(I)=U1
4     CONTINUE
      RETURN
      END
      BLOCK DATA AU053
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/DENDAT/RRD,RRB,RRC,RHOD,RHOA,
     $R1I,R2I,R3I,R1II,R2II,R3II,R2III,R3III
      COMMON/GLUDAT/DB,UB,DSW,
     $B0I,B1I,B2I,B3I,B4I,B2II,B3II,B4II,
     $B2III,B3III
      COMMON/POTDAT/D,A,RC,PHI1,PHI2,
     $A0I,A1I,A2I,A3I,A4I,
     $A0II,A1II,A2II,A3II,A4II,A5II,A6II,
     $A3III,A4III,A5III
      DATA RRD,RRB,RRC,RHOD,RHOA,
     $R1I,R2I,R3I,R1II,R2II,R3II,R2III,R3III/
     $0.2878207442141723D+01,
     $0.3500000000000000D+01,
     $0.3900000000000000D+01,
     $0.1000000000000000D+01,
     $0.0000000000000000D+00,
     $-0.6800000000000000D+00,
     $0.7500000000000000D+00,
     $-0.1333333333333333D+01,
     $-0.6800000000000000D+00,
     $0.7500000000000000D+00,
     $-0.1527241171296038D+01,
     $0.5578188675490974D+01,
     $0.6132971688727435D+01/
      DATA DB,UB,DSW/
     $0.1200000000000000D+02,
     $-0.3300000000000000D+01,
     $0.9358157767784574D+01/
      DATA B0I,B1I,B2I,B3I,B4I,B2II,B3II,B4II,B2III,B3III/
     $-0.2793388616771698D+01,
     $-0.3419999999999999D+00,
     $0.3902327808424106D-01,
     $0.7558829951858879D-02,
     $0.3090472511796849D-03,
     $0.8618226772941980D-01,
     $0.4341701445034724D-02,
     $-0.3044398779375916D-03,
     $0.8618226772941980D-01,
     $0.4325981467602070D-02/
      DATA D,A,RC,PHI1,PHI2/
     $0.2878207442141723D+01,
     $0.4070400000000000D+01,
     $0.3700000000000000D+01,
     $-0.8000000000000000D-01,
     $0.0000000000000000D+00/
      DATA A0I,A1I,A2I,A3I,A4I,A0II,A1II,A2II,A3II,A4II,A5II,A6II,
     $A3III,A4III,A5III/
     $-0.8000000000000000D-01,
     $0.0000000000000000D+00,
     $0.7619231375231362D+00,
     $-0.8333333333333333D+00,
     $-0.1211483464993159D+00,
     $-0.8000000000000000D-01,
     $0.0000000000000000D+00,
     $0.7619231375231362D+00,
     $-0.8333333333333333D+00,
     $-0.1096009851140349D+01,
     $0.2158417178555998D+01,
     $-0.9128915709636862D+00,
     $0.0000000000000000D+00,
     $0.0000000000000000D+00,
     $0.0000000000000000D+00/
      END
      SUBROUTINE DENFUN(R,RHO,DRHO,D2RHO)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/DENDAT/RRD,RRB,RRC,RHOD,RHOA,
     $R1I,R2I,R3I,R1II,R2II,R3II,R2III,R3III
      IF(R.GE.RRC)THEN
      RHO=0.D0
      DRHO=0.D0
      D2RHO=0.D0
      ELSE IF(R.GE.RRB)THEN
      X=R-RRC
      RHO=(X**2)*(R2III+X*R3III)
      DRHO=X*(2.D0*R2III+X*3.D0*R3III)
      D2RHO=2.D0*R2III+X*6.D0*R3III
      ELSE IF(R.GE.RRD)THEN
      X=R-RRD
      RHO=RHOD+X*(R1II+X*(R2II+X*R3II))
      DRHO=R1II+X*(2.D0*R2II+X*3.D0*R3II)
      D2RHO=2.D0*R2II+X*6.D0*R3II
      ELSE
      X=R-RRD
      RHO=RHOD+X*(R1I+X*(R2I+X*R3I))
      DRHO=R1I+X*(2.D0*R2I+X*3.D0*R3I)
      D2RHO=2.D0*R2I+X*6.D0*R3I
      ENDIF
      RETURN
      END
      SUBROUTINE GLUFUN(DENS,U,U1,U2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/GLUDAT/DB,UB,DSW,
     $B0I,B1I,B2I,B3I,B4I,B2II,B3II,B4II,
     $B2III,B3III
      IF(DENS.GT.DB)THEN
      X=DENS-DB
      U=UB+(X**2)*(B2III+X*B3III)
      U1=X*(2.D0*B2III+X*3.D0*B3III)
      U2=2.D0*B2III+X*6.D0*B3III
      ELSE IF(DENS.GT.DSW)THEN
      X=DENS-DB
      U=UB+(X**2)*(B2II+X*(B3II+X*B4II))
      U1=X*(2.D0*B2II+X*(3.D0*B3II+X*4.D0*B4II))
      U2=2.D0*B2II+X*(6.D0*B3II+X*12.D0*B4II)
      ELSE
      X=DENS-DSW
      U=B0I+X*(B1I+X*(B2I+X*(B3I+X*B4I)))
      U1=B1I+X*(2.D0*B2I+X*(3.D0*B3I+X*4.D0*B4I))
      U2=2.D0*B2I+X*(6.D0*B3I+X*12.D0*B4I)
      ENDIF
      RETURN
      END
      SUBROUTINE POTFUN(R,PHI,DPHI,D2PHI)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/POTDAT/D,A,RC,PHI1,PHI2,
     $A0I,A1I,A2I,A3I,A4I,
     $A0II,A1II,A2II,A3II,A4II,A5II,A6II,
     $A3III,A4III,A5III
      IF(R.GE.RC)THEN
      PHI=0.D0
      DPHI=0.D0
      D2PHI=0.D0
      ELSE IF(R.GE.A)THEN
      X=R-RC
      PHI=(X**3)*(A5III*X**2+A4III*X+A3III)
      DPHI=(X**2)*(5.D0*A5III*X**2+4.D0*A4III*X+3.D0*A3III)
      D2PHI=X*(20.D0*A5III*X**2+12.D0*A4III*X+6.D0*A3III)
      ELSE IF(R.GE.D)THEN
      X=R-D
      PHI=A0II+X*(A1II+X*(A2II+
     $X*(A3II+X*(A4II+X*(A5II+X*A6II)))))
      DPHI=A1II+X*(2.D0*A2II+X*(3.D0*A3II+
     $X*(4.D0*A4II+X*(5.D0*A5II+X*6.D0*A6II))))
      D2PHI=2.D0*A2II+X*(6.D0*A3II+X*(12.D0*A4II+
     $X*(20.D0*A5II+X*30.D0*A6II)))
      ELSE
      X=R-D
      PHI=A0I+X*(A1I+X*(A2I+X*(A3I+X*A4I)))
      DPHI=A1I+X*(2.D0*A2I+X*(3.D0*A3I+X*4.D0*A4I))
      D2PHI=2.D0*(A2I+X*(3.D0*A3I+X*6.D0*A4I))
      ENDIF
      RETURN
      END
      SUBROUTINE MSTEP
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/CNTRL/NTABLE,LISTEM,LMETHD,
     $NDIFF,NSTAT,NGOFR,NTRAJ,IVDUMP,ILIN,
     $LOCK(3,3),LILOCK(2,14),NFREED,LOCKCM,
     $NSCALE,ITPART,ITWALL
      COMMON/CONST/TWOPI,BOLTZ
      COMMON/COUNT/NFI,LCOUNT,LISTER,KNTSTA,KNTGOR,LEP,MANYON
      COMMON/DEN/RNRHO,RCRHO,RNRHO2,RCRHO2,DR2RHO,RHO(KR),DRHO(KR)
      COMMON/GEAR/F02,F12,F32,F42,F52
      COMMON/GLUE/DMIN,DMAX,DD,UJ(KG),DUJ(KG)
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/LISCOM/LIST(LL),MRKR1(NM),MRKR2(NM),LISLEN
      COMMON/LSTUPD/RLIST(3,NM)
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/MOTION/VIRKIN(3,3),VIRPOT(3,3),XNP(3,-2:NM)
      COMMON/PARAM/DELTA,DELTA2,GAMMA,VSCALE,CTRLCE
     $,CTRLMI,CTRLMA,RSQUPD,RANSQ,VMAS,BOX(3,3)
      COMMON/POT/RN,RC,RN2,RC2,DR2,VJ(KP),FJ(KP)
      COMMON/PRESS/PEXT,PAI(3,3)
      COMMON/PRINT/LNGPRT,IPRIND
      COMMON/SCRATC/DUMMY1(NM),DUMMY2(NM),DUMMY3(NM),DUMMY4(NM)
      COMMON/STATIS/FGS(NG),GRANG,FACNG,SCABY2,RESZ,DONTR,FONTR,SIG2,
     $NGS(NG),NGMAX,NZHIGH,NZLOW,MULTIP
      COMMON/SUMS/TEMPSM,TEMWSM,EKINSM,POT2SM,PGLUSM,POSTSM
     $,TOTESM,DENSSM,ALSM,VOLUSM,AREASM,HEIGSM
      COMMON/TIMERS/TIMSTR,TIMFIX,TIMBLD,TIMFRC,TIMINT
      COMMON/WALLS/HI(3,3),G(3,3),DH,AREA,VOLUME,SCM(3)
      DIMENSION H(3,3),HD(3,3),HIN(3,3)
      DIMENSION HIHD(3,3),HIHDT(3,3),HHH(3,3),HCOM(3,3)
      DIMENSION HDT(3,3),HSAV(3,3),HT(3,3)
      DIMENSION GI(3,3)
      DIMENSION VM(3,3)
      DIMENSION SMIN(3),SMAX(3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      EQUIVALENCE(X(1,-2,1),HD(1,1))
      EQUIVALENCE(XNP(1,-2),VM(1,1))
      EQUIVALENCE(XIN(1,-2),HIN(1,1))
      EXTERNAL SECOND
      SAVE TEMP,TEMW,TEMN,EKINP,EKINW,TOTE
      LOGICAL FIRSTP,UPDATE,UPGOFR
      SAVE FIRSTP,UPDATE
      CHARACTER*1 CHARD,CHARL,CHARS,CHART,CHARP,CHARW,CHARV
      SAVE CHARD,CHARL,CHARS,CHART,CHARP,CHARW,CHARV
      DATA CHARD,CHARL,CHARS,CHART,CHARP,CHARW,CHARV
     $/' ' , ' ' , ' ' , ' ' , ' ' , ' ' , ' '/
      DATA FIRSTP,UPDATE/.TRUE.,.TRUE./
      CPTIME=SECOND()
      IF(FIRSTP)THEN
      TEMP=CTRLCE
      TEMW=CTRLCE
      TEMN=CTRLCE
      TOTE=CTRLCE
      FIRSTP=.FALSE.
      ELSE
      IF((ITPART.EQ.2).OR.(ITWALL.EQ.2))THEN
      EKINC=0.D0
      IF(ITPART.EQ.2)EKINC=EKINC+EKINP
      IF(ITWALL.EQ.2)EKINC=EKINC+EKINW
      ENDIF
      ENDIF
      IF(NSCALE.GT.0)THEN
      CHARV=' '
      IF(MOD(NFI,NSCALE).EQ.0)THEN
      CHARV='V'
      PSCALE=VSCALE
      WSCALE=VSCALE
      ELSE
      PSCALE=1.D0
      WSCALE=1.D0
      ENDIF
      ELSE
      PSCALE=1.D0
      WSCALE=1.D0
      ENDIF
      IF(ITPART.NE.0)THEN
      CHARP=' '
      IF(ITPART.EQ.1)THEN
      IF((TEMN.LT.CTRLMI).OR.(TEMN.GT.CTRLMA))THEN
      IF(TEMN.GT.0.D0)THEN
      PSCALE=SQRT(CTRLCE/TEMN)
      CHARP='P'
      ENDIF
      ENDIF
      ELSE
      IF((TOTE.LT.CTRLMI).OR.(TOTE.GT.CTRLMA))THEN
      IF(EKINC.GT.0.D0)THEN
      PSCAL2=1.D0+(CTRLCE-TOTE)/EKINC
      IF(PSCAL2.GT.0.D0)THEN
      PSCALE=SQRT(PSCAL2)
      ELSE
      PSCALE=0.D0
      ENDIF
      CHARP='P'
      ENDIF
      ENDIF
      ENDIF
      ENDIF
      IF(ITWALL.NE.0)THEN
      CHARW=' '
      IF(ITWALL.EQ.1)THEN
      IF((TEMW.LT.CTRLMI).OR.(TEMW.GT.CTRLMA))THEN
      IF(TEMW.GT.0.D0)THEN
      WSCALE=SQRT(TEMP/TEMW)
      CHARW='W'
      ENDIF
      ENDIF
      ELSE
      IF((TOTE.LT.CTRLMI).OR.(TOTE.GT.CTRLMA))THEN
      IF(EKINC.GT.0.D0)THEN
      WSCAL2=1.D0+(CTRLCE-TOTE)/EKINC
      IF(WSCAL2.GT.0.D0)THEN
      WSCALE=SQRT(WSCAL2)
      ELSE
      WSCALE=0.D0
      ENDIF
      CHARW='W'
      ENDIF
      ENDIF
      ENDIF
      ENDIF
      TIMREF=SECOND()
      DO 102 K=1,3
      DO 101 I=-2,MOLSP
      SCAFAC=PSCALE
      IF(I.LE.0)SCAFAC=WSCALE
      X(K,I,1)=X(K,I,1)*SCAFAC
      X0(K,I)=X0(K,I)+X(K,I,1)+X(K,I,2)+X(K,I,3)
     $+X(K,I,4)+X(K,I,5)
      X(K,I,1)=X(K,I,1)+2.D0*X(K,I,2)+3.D0*X(K,I,3)
     $+4.D0*X(K,I,4)+5.D0*X(K,I,5)
      X(K,I,2)=X(K,I,2)+3.D0*X(K,I,3)+6.D0*X(K,I,4)
     $+10.D0*X(K,I,5)
      X(K,I,3)=X(K,I,3)+4.D0*X(K,I,4)+10.D0*X(K,I,5)
      X(K,I,4)=X(K,I,4)+5.D0*X(K,I,5)
      XNP(K,I)=0.D0
101   CONTINUE
102   CONTINUE
      TIMINT=TIMINT+(SECOND()-TIMREF)
      DO 1946 I=1,3
      DO 1947 J=1,3
      HT(I,J)=H(J,I)
      HDT(I,J)=HD(J,I)
1947  CONTINUE
1946  CONTINUE
      CALL MTXINV(H,HI,DH)
      CALL MTXMTP(HT,H,G)
      CALL MTXMTP(HDT,HD,HHH)
      URANS=0.5D0*(HHH(1,1)+HHH(2,2)+HHH(3,3))/(DELTA2*VMAS)
      DO 2400 K=1,3
      IF(LPBC(K).EQ.0)THEN
      SCM(K)=0.D0
      SMIN(K)=0.D0
      SMAX(K)=0.D0
      DO 2500 I=1,MOLSA
      SSKI=X0(K,I)
      SCM(K)=SCM(K)+SSKI
      IF(SSKI.GT.SMAX(K))SMAX(K)=SSKI
      IF(SSKI.LT.SMIN(K))SMIN(K)=SSKI
2500  CONTINUE
      SCM(K)=SCM(K)/MOLSA
      IF(LOCKCM.NE.0)THEN
      DO 2600 I=1,MOLSA
      X0(K,I)=X0(K,I)-SCM(K)
2600  CONTINUE
      SCM(K)=0.D0
      ENDIF
      ELSE
      SCM(K)=0.D0
      SMIN(K)=-0.5D0
      SMAX(K)=+0.5D0
      ENDIF
2400  CONTINUE
      VOLUME=DH*(SMAX(1)-SMIN(1))*
     $(SMAX(2)-SMIN(2))*
     $(SMAX(3)-SMIN(3))
      AREA=HI(3,3)*DH*(SMAX(1)-SMIN(1))*(SMAX(2)-SMIN(2))
      HEIG=VOLUME/AREA
      CALL RESET(VIRKIN,9)
      DO 270 I=1,MOLSP
      VELOCX=H(1,1)*X(1,I,1)+H(1,2)*X(2,I,1)+H(1,3)*X(3,I,1)
      VELOCY=H(2,1)*X(1,I,1)+H(2,2)*X(2,I,1)+H(2,3)*X(3,I,1)
      VELOCZ=H(3,1)*X(1,I,1)+H(3,2)*X(2,I,1)+H(3,3)*X(3,I,1)
      VIRKIN(1,1)=VIRKIN(1,1)+VELOCX*X(1,I,1)
      VIRKIN(2,1)=VIRKIN(2,1)+VELOCY*X(1,I,1)
      VIRKIN(3,1)=VIRKIN(3,1)+VELOCZ*X(1,I,1)
      VIRKIN(1,2)=VIRKIN(1,2)+VELOCX*X(2,I,1)
      VIRKIN(2,2)=VIRKIN(2,2)+VELOCY*X(2,I,1)
      VIRKIN(3,2)=VIRKIN(3,2)+VELOCZ*X(2,I,1)
      VIRKIN(1,3)=VIRKIN(1,3)+VELOCX*X(3,I,1)
      VIRKIN(2,3)=VIRKIN(2,3)+VELOCY*X(3,I,1)
      VIRKIN(3,3)=VIRKIN(3,3)+VELOCZ*X(3,I,1)
270   CONTINUE
      UPGOFR=.FALSE.
      IF(NGOFR.GT.0)THEN
      IF(MOD(NFI-1,NGOFR).EQ.0)THEN
      UPGOFR=.TRUE.
      UPDATE=.TRUE.
      ENDIF
      ENDIF
      IF(.NOT.UPDATE)THEN
      IF(LISTEM.GT.0)THEN
      UPDATE=(LISTER.EQ.LISTEM)
      ELSE
      DO 304 I=1,MOLSA
      RXI=H(1,1)*X0(1,I)+H(1,2)*X0(2,I)+H(1,3)*X0(3,I)
      RYI=H(2,1)*X0(1,I)+H(2,2)*X0(2,I)+H(2,3)*X0(3,I)
      RZI=H(3,1)*X0(1,I)+H(3,2)*X0(2,I)+H(3,3)*X0(3,I)
      RSQ=(RXI-RLIST(1,I))**2+
     $(RYI-RLIST(2,I))**2+
     $(RZI-RLIST(3,I))**2
      IF(RSQ.GT.RSQUPD)THEN
      UPDATE=.TRUE.
      GOTO 306
      ENDIF
304   CONTINUE
306   CONTINUE
      ENDIF
      ENDIF
      IF(UPDATE)THEN
      TIMREF=SECOND()
      LISTER=1
      IF(UPGOFR)THEN
      CHARL='G'
      CALL MLIST(-1)
      KNTGOR=KNTGOR+1
      ELSE
      CHARL='L'
      CALL MLIST(LMETHD)
      ENDIF
      TIMBLD=TIMBLD+(SECOND()-TIMREF)
      LCOUNT=LCOUNT+1
      IF(LISTEM.LE.0)THEN
      DO 302 I=1,MOLSA
      RLIST(1,I)=H(1,1)*X0(1,I)+H(1,2)*X0(2,I)+
     $H(1,3)*X0(3,I)
      RLIST(2,I)=H(2,1)*X0(1,I)+H(2,2)*X0(2,I)+
     $H(2,3)*X0(3,I)
      RLIST(3,I)=H(3,1)*X0(1,I)+H(3,2)*X0(2,I)+
     $H(3,3)*X0(3,I)
302   CONTINUE
      ENDIF
      UPDATE=.FALSE.
      ELSE
      CHARL=' '
      LISTER=LISTER+1
      ENDIF
      TIMREF=SECOND()
      CALL MFORCE(POT2,PGLU)
      POST=0.D0
      DO 747 I=1,3
      DO 748 J=1,3
      VM(I,J)=VMAS*(VIRKIN(I,J)/DELTA2+VIRPOT(I,J)
     $-(PAI(J,1)*H(I,1)+PAI(J,2)*H(I,2)
     $+PAI(J,3)*H(I,3))*(1-ILIN)
     $-PAI(I,J)*ILIN-PEXT*DH*HI(J,I))
     $-HD(I,J)*GAMMA/DELTA
      IF(LOCK(I,J).LT.0)VM(I,J)=0.D0
      POST=POST+(1-ILIN)*0.5D0*PAI(I,J)*G(I,J)
     $+ILIN*PAI(I,J)*H(I,J)
748   CONTINUE
747   CONTINUE
      CALL SYMM(VM)
      TIMFRC=TIMFRC+(SECOND()-TIMREF)
      DIFFUS=0.D0
      DO 107 I=1,MOLSP
      POSXI=H(1,1)*X0(1,I)+H(1,2)*X0(2,I)+H(1,3)*X0(3,I)
      POSYI=H(2,1)*X0(1,I)+H(2,2)*X0(2,I)+H(2,3)*X0(3,I)
      POSZI=H(3,1)*X0(1,I)+H(3,2)*X0(2,I)+H(3,3)*X0(3,I)
      DISXI=POSXI-(HIN(1,1)*XIN(1,I)
     $+HIN(1,2)*XIN(2,I)
     $+HIN(1,3)*XIN(3,I))
      DISYI=POSYI-(HIN(2,1)*XIN(1,I)
     $+HIN(2,2)*XIN(2,I)
     $+HIN(2,3)*XIN(3,I))
      DISZI=POSZI-(HIN(3,1)*XIN(1,I)
     $+HIN(3,2)*XIN(2,I)
     $+HIN(3,3)*XIN(3,I))
      DIFFUS=DIFFUS+DISXI**2+DISYI**2+DISZI**2
107   CONTINUE
      IF(MOLSP.GT.0)DIFFUS=DIFFUS/MOLSP
      TRANS=0.D0
      DO 450 I=1,MOLSP
      VELOCX=H(1,1)*X(1,I,1)+H(1,2)*X(2,I,1)+H(1,3)*X(3,I,1)
      VELOCY=H(2,1)*X(1,I,1)+H(2,2)*X(2,I,1)+H(2,3)*X(3,I,1)
      VELOCZ=H(3,1)*X(1,I,1)+H(3,2)*X(2,I,1)+H(3,3)*X(3,I,1)
      TRANS=TRANS+VELOCX**2+VELOCY**2+VELOCZ**2
450   CONTINUE
      TRANS=0.5D0*TRANS/DELTA2
      TIMREF=SECOND()
      CALL MTXMTP(HI,HD,HIHD)
      CALL MTXMTP(HDT,H,HIHDT)
      CALL MTXINV(G,GI,DG)
      CALL MTXMTP(GI,HIHDT,HCOM)
      DO 279 I=1,3
      HSAV(I,1)=HCOM(I,1)+HIHD(I,1)
      HSAV(I,2)=HCOM(I,2)+HIHD(I,2)
      HSAV(I,3)=HCOM(I,3)+HIHD(I,3)
279   CONTINUE
      DO 7110 I=1,MOLSP
      VELOCX=HSAV(1,1)*X(1,I,1)+HSAV(1,2)*X(2,I,1)+
     $HSAV(1,3)*X(3,I,1)
      VELOCY=HSAV(2,1)*X(1,I,1)+HSAV(2,2)*X(2,I,1)+
     $HSAV(2,3)*X(3,I,1)
      VELOCZ=HSAV(3,1)*X(1,I,1)+HSAV(3,2)*X(2,I,1)+
     $HSAV(3,3)*X(3,I,1)
      XNP(1,I)=XNP(1,I)-VELOCX/DELTA2
      XNP(2,I)=XNP(2,I)-VELOCY/DELTA2
      XNP(3,I)=XNP(3,I)-VELOCZ/DELTA2
7110  CONTINUE
      DO 110 K=1,3
      DO 111 I=-2,MOLSP
      XNP(K,I)=X(K,I,2)-0.5D0*DELTA2*XNP(K,I)
      X0(K,I)=X0(K,I)-XNP(K,I)*F02
      X(K,I,1)=X(K,I,1)-XNP(K,I)*F12
      X(K,I,2)=X(K,I,2)-XNP(K,I)
      X(K,I,3)=X(K,I,3)-XNP(K,I)*F32
      X(K,I,4)=X(K,I,4)-XNP(K,I)*F42
      X(K,I,5)=X(K,I,5)-XNP(K,I)*F52
111   CONTINUE
110   CONTINUE
      TIMINT=TIMINT+(SECOND()-TIMREF)
      SUMPX=0.D0
      SUMPY=0.D0
      SUMPZ=0.D0
      DO 1000 I=1,MOLSA
      VELOCX=H(1,1)*X(1,I,1)+H(1,2)*X(2,I,1)+H(1,3)*X(3,I,1)
      VELOCY=H(2,1)*X(1,I,1)+H(2,2)*X(2,I,1)+H(2,3)*X(3,I,1)
      VELOCZ=H(3,1)*X(1,I,1)+H(3,2)*X(2,I,1)+H(3,3)*X(3,I,1)
      SUMPX=SUMPX+VELOCX
      SUMPY=SUMPY+VELOCY
      SUMPZ=SUMPZ+VELOCZ
1000  CONTINUE
      SUMPX=SUMPX/DELTA
      SUMPY=SUMPY/DELTA
      SUMPZ=SUMPZ/DELTA
      AL=(VOLUME/(NBX*NBY*NBZ))**(1.D0/3.D0)
      DENS=MOLSA/VOLUME
      EKINP=TRANS/MOLSA
      EKINW=URANS/MOLSA
      EKIN=EKINP+EKINW
      POT2=POT2/MOLSA
      PGLU=PGLU/MOLSA
      POTE=POT2+PGLU
      POST=(POST+PEXT*DH)/MOLSA
      TOTE=EKIN+POTE+POST
      IF(NFREED.GT.0)THEN
      TEMW=BOLTZ*2.D0*URANS/NFREED
      ELSE
      TEMW=0.D0
      ENDIF
      TEMP=BOLTZ*2.D0*TRANS/(3*MOLSP)
      TEMN=BOLTZ*2.D0*(URANS+TRANS)/(3*MOLSP+NFREED)
      TEMPSM=TEMPSM+TEMP
      TEMWSM=TEMWSM+TEMW
      EKINSM=EKINSM+EKIN
      POT2SM=POT2SM+POT2
      PGLUSM=PGLUSM+PGLU
      POSTSM=POSTSM+POST
      TOTESM=TOTESM+TOTE
      DENSSM=DENSSM+DENS
      ALSM=ALSM+AL
      VOLUSM=VOLUSM+VOLUME
      AREASM=AREASM+AREA
      HEIGSM=HEIGSM+HEIG
      TIMMST=SECOND()-CPTIME
      IF(LNGPRT.GT.0)THEN
      PRINT '(1X,I5,1X,2A1,3F8.4,F9.4,3E9.1)',
     $NFI,CHARL,CHARP,
     $EKIN,POTE,TOTE,DIFFUS,SUMPX,SUMPY,SUMPZ
      ENDIF
      RETURN
      END
      SUBROUTINE MLIST(LMETHD)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/COUNT/NFI,LCOUNT,LISTER,KNTSTA,KNTGOR,LEP,MANYON
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/LISCOM/LIST(LL),MRKR1(NM),MRKR2(NM),LISLEN
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/PARAM/DELTA,DELTA2,GAMMA,VSCALE,CTRLCE
     $,CTRLMI,CTRLMA,RSQUPD,RANSQ,VMAS,BOX(3,3)
      COMMON/PBCS/HALF,PBCX,PBCY,PBCZ
      COMMON/PRINT/LNGPRT,IPRIND
      COMMON/SCRATC/DUMMY1(NM),DUMMY2(NM),DUMMY3(NM),DUMMY4(NM)
      COMMON/STATIS/FGS(NG),GRANG,FACNG,SCABY2,RESZ,DONTR,FONTR,SIG2,
     $NGS(NG),NGMAX,NZHIGH,NZLOW,MULTIP
      COMMON/WALLS/HI(3,3),G(3,3),DH,AREA,VOLUME,SCM(3)
      DIMENSION H(3,3),HIN(3,3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      EQUIVALENCE(XIN(1,-2),HIN(1,1))
2     FORMAT(1X,2A)
3     FORMAT(/,1X,2A)
      IF(LPBCSM.GT.0)THEN
      DO 999 I=1,MOLSP
      BOXJMP=PBCX*INT(X0(1,I)+SIGN(HALF,X0(1,I)))
      X0(1,I)=X0(1,I)-BOXJMP
      XIN(1,I)=XIN(1,I)-BOXJMP
      BOXJMP=PBCY*INT(X0(2,I)+SIGN(HALF,X0(2,I)))
      X0(2,I)=X0(2,I)-BOXJMP
      XIN(2,I)=XIN(2,I)-BOXJMP
      BOXJMP=PBCZ*INT(X0(3,I)+SIGN(HALF,X0(3,I)))
      X0(3,I)=X0(3,I)-BOXJMP
      XIN(3,I)=XIN(3,I)-BOXJMP
999   CONTINUE
      ENDIF
      IF(LMETHD.EQ.1)THEN
      CALL CBUILD(RANSQ,ICODE)
      IF(ICODE.EQ.2)THEN
      CALL FBUILD(RANSQ,ICODE)
      ENDIF
      ELSE IF(LMETHD.EQ.-1)THEN
      CALL GBUILD(RANSQ,ICODE)
      ELSE
      CALL FBUILD(RANSQ,ICODE)
      ENDIF
      IF(ICODE.NE.0)THEN
      IF(ICODE.EQ.1)THEN
      PRINT 3,'''LIST'' ARRAY OVERFLOW IN CBUILD/FBUILD',
     $'(''LL'' TOO SMALL)'
      ELSE IF(ICODE.EQ.3)THEN
      PRINT 2,'''NPC'' ARRAY OVERFLOW IN CBUILD,',
     $'(''NCEMAX'' TOO SMALL).'
      ELSE IF(ICODE.EQ.4)THEN
      PRINT 2,'''KNTNTS'' ARRAY OVERFLOW IN CBUILD',
     $'(''KNTSIZ'' TOO SMALL).'
      ELSE IF(ICODE.EQ.5)THEN
      PRINT 2,'''NEIGH'' ARRAY OVERFLOW IN CBUILD',
     $'(''NNEMAX'' TOO SMALL).'
      ENDIF
      PRINT 2,'INCREASE DIMENSIONS, RECOMPILE AND RERUN.'
      STOP
      ENDIF
      RETURN
      END
      SUBROUTINE GBUILD(RANSQ,ICODE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/COUNT/NFI,LCOUNT,LISTER,KNTSTA,KNTGOR,LEP,MANYON
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/LISCOM/LIST(LL),MRKR1(NM),MRKR2(NM),LISLEN
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/PBCS/HALF,PBCX,PBCY,PBCZ
      COMMON/PRINT/LNGPRT,IPRIND
      COMMON/SCRATC/DUMMY1(NM),DUMMY2(NM),DUMMY3(NM),DUMMY4(NM)
      DIMENSION SLICE(NM)
      EQUIVALENCE(SLICE,DUMMY1)
      DIMENSION RRPOS(3,NM)
      EQUIVALENCE(RRPOS,DUMMY2)
      COMMON/STATIS/FGS(NG),GRANG,FACNG,SCABY2,RESZ,DONTR,FONTR,SIG2,
     $NGS(NG),NGMAX,NZHIGH,NZLOW,MULTIP
      COMMON/WALLS/HI(3,3),G(3,3),DH,AREA,VOLUME,SCM(3)
      DIMENSION H(3,3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      G11=G(1,1)
      G22=G(2,2)
      G33=G(3,3)
      G12D=G(1,2)+G(2,1)
      G13D=G(1,3)+G(3,1)
      G23D=G(2,3)+G(3,2)
      CONTR=VOLUME/DONTR/MOLSA**2
      IF(LPBCSM.GT.1)THEN
      GONTR=AREA/FONTR
      GSLSC=FACNG
      ELSE
      GONTR=1.D0/FONTR
      GSLSC=NGMAX/2.D0
      ENDIF
      DO 200 I=1,MOLSP
      RRPOS(1,I)=H(1,1)*X0(1,I)+H(1,2)*X0(2,I)+H(1,3)*X0(3,I)
      RRPOS(2,I)=H(2,1)*X0(1,I)+H(2,2)*X0(2,I)+H(2,3)*X0(3,I)
      RRPOS(3,I)=H(3,1)*X0(1,I)+H(3,2)*X0(2,I)+H(3,3)*X0(3,I)
      IF(LPBCSM.GT.1)THEN
      SLICE(I)=RRPOS(3,I)
      ELSE
      SLICE(I)=RRPOS(1,I)**2+RRPOS(2,I)**2+RRPOS(3,I)**2
      ENDIF
200   CONTINUE
      NLL=1
      DO 300 I=1,MOLSA
      MRKR1(I)=NLL
      DO 280 J=I+1,MOLSA
      XIJ=X0(1,I)-X0(1,J)
      IF(XIJ.GT.+HALF)XIJ=XIJ-PBCX
      IF(XIJ.LT.-HALF)XIJ=XIJ+PBCX
      YIJ=X0(2,I)-X0(2,J)
      IF(YIJ.GT.+HALF)YIJ=YIJ-PBCY
      IF(YIJ.LT.-HALF)YIJ=YIJ+PBCY
      ZIJ=X0(3,I)-X0(3,J)
      IF(ZIJ.GT.+HALF)ZIJ=ZIJ-PBCZ
      IF(ZIJ.LT.-HALF)ZIJ=ZIJ+PBCZ
      RSQ=XIJ*(G11*XIJ+G12D*YIJ+G13D*ZIJ)+
     $YIJ*(G22*YIJ+G23D*ZIJ)+G33*ZIJ*ZIJ
      NBET=INT(FACNG*RSQ)+1
      IF(NBET.LE.NGMAX)THEN
      FGS(NBET)=FGS(NBET)+CONTR
      NGS(NBET)=NGS(NBET)+1
      ENDIF
      IF(RSQ.LT.RANSQ)THEN
      LIST(NLL)=J
      NLL=NLL+1
      ENDIF
280   CONTINUE
      MRKR2(I)=NLL-1
300   CONTINUE
      LISLEN=NLL-1
      IF(LISLEN.GT.LL)THEN
      ICODE=1
      ELSE
      ICODE=0
      ENDIF
      IF(LNGPRT.GT.0)THEN
      PRINT '(1X,A8,I8,A1,I8)',
     $'LENGTH =',LISLEN,'/',LL
      ENDIF
      RETURN
      END
      SUBROUTINE FBUILD(RANSQ,ICODE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/LISCOM/LIST(LL),MRKR1(NM),MRKR2(NM),LISLEN
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/PBCS/HALF,PBCX,PBCY,PBCZ
      COMMON/PRINT/LNGPRT,IPRIND
      COMMON/SCRATC/DUMMY1(NM),DUMMY2(NM),DUMMY3(NM),DUMMY4(NM)
      INTEGER MARK(NM)
      EQUIVALENCE(MARK,DUMMY1)
      COMMON/WALLS/HI(3,3),G(3,3),DH,AREA,VOLUME,SCM(3)
      DIMENSION H(3,3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      G11=G(1,1)
      G22=G(2,2)
      G33=G(3,3)
      G12D=G(1,2)+G(2,1)
      G13D=G(1,3)+G(3,1)
      G23D=G(2,3)+G(3,2)
      NLL=1
      DO 300 I=1,MOLSA
      DO 280 J=I+1,MOLSA
      XIJ=X0(1,I)-X0(1,J)
      IF(XIJ.GT.+HALF)XIJ=XIJ-PBCX
      IF(XIJ.LT.-HALF)XIJ=XIJ+PBCX
      YIJ=X0(2,I)-X0(2,J)
      IF(YIJ.GT.+HALF)YIJ=YIJ-PBCY
      IF(YIJ.LT.-HALF)YIJ=YIJ+PBCY
      ZIJ=X0(3,I)-X0(3,J)
      IF(ZIJ.GT.+HALF)ZIJ=ZIJ-PBCZ
      IF(ZIJ.LT.-HALF)ZIJ=ZIJ+PBCZ
      RSQ=XIJ*(G11*XIJ+G12D*YIJ+G13D*ZIJ)+
     $YIJ*(G22*YIJ+G23D*ZIJ)+G33*ZIJ*ZIJ
      MARK(J)=0
      IF(RSQ.LT.RANSQ)MARK(J)=1
280   CONTINUE
      MRKR1(I)=NLL
      DO 290 J=I+1,MOLSA
      LIST(NLL)=J
      NLL=NLL+MARK(J)
290   CONTINUE
      MRKR2(I)=NLL-1
300   CONTINUE
      LISLEN=NLL-1
      IF(LISLEN.GT.LL)THEN
      ICODE=1
      ELSE
      ICODE=0
      ENDIF
      IF(LNGPRT.GT.0)THEN
      PRINT '(1X,A8,I8,A1,I8)',
     $'LENGTH =',LISLEN,'/',LL
      ENDIF
      RETURN
      END
      SUBROUTINE CBUILD(RANSQ,ICODE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/LISCOM/LIST(LL),MRKR1(NM),MRKR2(NM),LISLEN
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/PBCS/HALF,PBCX,PBCY,PBCZ
      COMMON/PRINT/LNGPRT,IPRIND
      COMMON/SCRATC/DUMMY1(NM),DUMMY2(NM),DUMMY3(NM),DUMMY4(NM)
      INTEGER LCELL(NM)
      EQUIVALENCE(LCELL,DUMMY1)
      INTEGER MARK(NM)
      EQUIVALENCE(MARK,DUMMY1)
      EQUIVALENCE(KNTNTS,DUMMY2)
      COMMON/WALLS/HI(3,3),G(3,3),DH,AREA,VOLUME,SCM(3)
      DIMENSION H(3,3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      PARAMETER(NNEMAX=512)
      INTEGER NEIGH(NNEMAX)
      PARAMETER(NMHALF=NM/2,NMHAL1=NMHALF+1)
      PARAMETER(NCEMAX=NMHALF)
      PARAMETER(NCEMA1=NCEMAX-1)
      INTEGER NPC(0:NCEMA1)
      PARAMETER(KNTSIZ=3*NM)
      INTEGER KNTNTS(KNTSIZ)
      LOGICAL NOFLDX,NOFLDY,NOFLDZ
      INTEGER NIX(13),NIY(13),NIZ(13)
      DATA NIX/-1,-1,-1, 0, 0,-1, 1,-1, 0, 1,-1, 0, 1/
      DATA NIY/0,-1, 1, 1, 0, 0, 0,-1,-1,-1, 1, 1, 1/
      DATA NIZ/0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1/
      RANGE=SQRT(RANSQ)
      G11=G(1,1)
      G22=G(2,2)
      G33=G(3,3)
      G12D=G(1,2)+G(2,1)
      G13D=G(1,3)+G(3,1)
      G23D=G(2,3)+G(3,2)
      NLX=INT(1.D0/(RANGE*HI(1,1)))
      NLY=INT(1.D0/(RANGE*HI(2,2)))
      NLZ=INT(1.D0/(RANGE*HI(3,3)))
      IF((NLX.LT.3).AND.(LPBC(1).NE.0))NLX=1
      IF((NLY.LT.3).AND.(LPBC(2).NE.0))NLY=1
      IF((NLZ.LT.3).AND.(LPBC(3).NE.0))NLZ=1
      NOFLDX=(LPBC(1).EQ.0).OR.(NLX.EQ.1)
      NOFLDY=(LPBC(2).EQ.0).OR.(NLY.EQ.1)
      NOFLDZ=(LPBC(3).EQ.0).OR.(NLZ.EQ.1)
      FNLX=NLX
      FNLY=NLY
      FNLZ=NLZ
      NCELLS=NLX*NLY*NLZ
      IF(NCELLS.GT.NCEMAX)THEN
      ICODE=3
      RETURN
      ENDIF
      NPCMAX=KNTSIZ/NCELLS
      DO 1100 I=1,MOLSA
      IX=INT((X0(1,I)+HALF)*FNLX)
      IF(IX.GE.NLX)IX=NLX-1
      IF(IX.LT.0)IX=0
      IY=INT((X0(2,I)+HALF)*FNLY)
      IF(IY.GE.NLY)IY=NLY-1
      IF(IY.LT.0)IY=0
      IZ=INT((X0(3,I)+HALF)*FNLZ)
      IF(IZ.GE.NLZ)IZ=NLZ-1
      IF(IZ.LT.0)IZ=0
      LCELL(I)=IX+NLX*(IY+NLY*IZ)
1100  CONTINUE
      DO 1200 ICELL=0,NCELLS-1
      NPC(ICELL)=0
1200  CONTINUE
      DO 1300 I=1,MOLSA
      ICELL=LCELL(I)
      NPC(ICELL)=NPC(ICELL)+1
      KNTNTS(NPC(ICELL)+NPCMAX*ICELL)=I
1300  CONTINUE
      NPCUSD=0
      DO 1400 ICELL=0,NCELLS-1
      IF(NPC(ICELL).GT.NPCUSD)NPCUSD=NPC(ICELL)
1400  CONTINUE
      IF(NPCUSD.GT.NPCMAX)THEN
      ICODE=4
      RETURN
      ENDIF
      NLL=1
      NNEUSD=0
      DO 2100 ICELL=0,NCELLS-1
      IF(NPC(ICELL).NE.0)THEN
      ICOFFS=NPCMAX*ICELL
      IZ=ICELL/(NLX*NLY)
      ICELXY=ICELL-(NLX*NLY)*IZ
      IY=ICELXY/NLX
      IX=ICELXY-NLX*IY
      NNEIGC=0
      DO 2200 KC=1,13
      JZ=IZ+NIZ(KC)
      IF(JZ.GE.NLZ)THEN
      IF(NOFLDZ)GOTO 2200
      JZ=0
      ELSE IF(JZ.LT.0)THEN
      IF(NOFLDZ)GOTO 2200
      JZ=NLZ-1
      ENDIF
      JY=IY+NIY(KC)
      IF(JY.GE.NLY)THEN
      IF(NOFLDY)GOTO 2200
      JY=0
      ELSE IF(JY.LT.0)THEN
      IF(NOFLDY)GOTO 2200
      JY=NLY-1
      ENDIF
      JX=IX+NIX(KC)
      IF(JX.GE.NLX)THEN
      IF(NOFLDX)GOTO 2200
      JX=0
      ELSE IF(JX.LT.0)THEN
      IF(NOFLDX)GOTO 2200
      JX=NLX-1
      ENDIF
      JCELL=JX+NLX*(JY+NLY*JZ)
      JCOFFS=NPCMAX*JCELL
      DO 2300 JPC=1,NPC(JCELL)
      NEIGH(NNEIGC+JPC)=KNTNTS(JPC+JCOFFS)
2300  CONTINUE
      NNEIGC=NNEIGC+NPC(JCELL)
2200  CONTINUE
      DO 2500 IPC=1,NPC(ICELL)
      I=KNTNTS(IPC+ICOFFS)
      DO 2700 JPC=IPC+1,NPC(ICELL)
      NEIGH(NNEIGC+JPC-IPC)=KNTNTS(JPC+ICOFFS)
2700  CONTINUE
      NNEIGI=NNEIGC+NPC(ICELL)-IPC
      IF(NNEIGI.GT.NNEUSD)NNEUSD=NNEIGI
      DO 280 ICAND=1,NNEIGI
      J=NEIGH(ICAND)
      XIJ=X0(1,I)-X0(1,J)
      IF(XIJ.GT.+HALF)XIJ=XIJ-PBCX
      IF(XIJ.LT.-HALF)XIJ=XIJ+PBCX
      YIJ=X0(2,I)-X0(2,J)
      IF(YIJ.GT.+HALF)YIJ=YIJ-PBCY
      IF(YIJ.LT.-HALF)YIJ=YIJ+PBCY
      ZIJ=X0(3,I)-X0(3,J)
      IF(ZIJ.GT.+HALF)ZIJ=ZIJ-PBCZ
      IF(ZIJ.LT.-HALF)ZIJ=ZIJ+PBCZ
      RSQ=XIJ*(G11*XIJ+G12D*YIJ+G13D*ZIJ)+
     $YIJ*(G22*YIJ+G23D*ZIJ)+G33*ZIJ*ZIJ
      MARK(ICAND)=0
      IF(RSQ.LT.RANSQ)MARK(ICAND)=1
280   CONTINUE
      MRKR1(I)=NLL
      DO 290 ICAND=1,NNEIGI
      LIST(NLL)=NEIGH(ICAND)
      NLL=NLL+MARK(ICAND)
290   CONTINUE
      MRKR2(I)=NLL-1
2500  CONTINUE
      ENDIF
2100  CONTINUE
      LISLEN=NLL-1
      IF(LISLEN.GT.LL)THEN
      ICODE=1
      ELSE IF(NNEUSD.GT.NNEMAX)THEN
      ICODE=5
      ELSE
      ICODE=0
      ENDIF
      IF(LNGPRT.GT.0)THEN
      PRINT '(1X,A7,I7,A1,I7,
     $A8,I5,A1,I5,
     $A8,I3,A1,I3,
     $A8,I5,A1,I5)',
     $'LENGTH=',LISLEN,'/',LL,
     $', CELLS=',NCELLS,'/',NCEMAX,
     $', PA/CL=',NPCUSD,'/',NPCMAX,
     $', CA/PA=',NNEUSD,'/',NNEMAX
      ENDIF
      RETURN
      END
      SUBROUTINE MFORCE(POT2,PGLU)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/COUNT/NFI,LCOUNT,LISTER,KNTSTA,KNTGOR,LEP,MANYON
      COMMON/DEN/RNRHO,RCRHO,RNRHO2,RCRHO2,DR2RHO,RHO(KR),DRHO(KR)
      COMMON/GLUE/DMIN,DMAX,DD,UJ(KG),DUJ(KG)
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/LISCOM/LIST(LL),MRKR1(NM),MRKR2(NM),LISLEN
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/MOTION/VIRKIN(3,3),VIRPOT(3,3),XNP(3,-2:NM)
      COMMON/PBCS/HALF,PBCX,PBCY,PBCZ
      COMMON/POT/RN,RC,RN2,RC2,DR2,VJ(KP),FJ(KP)
      COMMON/SCRATC/DUMMY1(NM),DUMMY2(NM),DUMMY3(NM),DUMMY4(NM)
      DIMENSION DERU(NM)
      EQUIVALENCE(DERU,DUMMY1)
      DIMENSION DNSTY(NM)
      EQUIVALENCE(DNSTY,DUMMY2)
      DIMENSION PP(NM)
      EQUIVALENCE(PP,DUMMY3)
      COMMON/WALLS/HI(3,3),G(3,3),DH,AREA,VOLUME,SCM(3)
      DIMENSION H(3,3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      CALL RESET(VIRPOT,9)
      LEP=0
      PGLU=0.D0
      G11=G(1,1)
      G22=G(2,2)
      G33=G(3,3)
      G12D=G(1,2)+G(2,1)
      G13D=G(1,3)+G(3,1)
      G23D=G(2,3)+G(3,2)
      IF(MANYON.NE.0)THEN
      DO 190 I=1,MOLSA
      DNSTY(I)=0.D0
190   CONTINUE
      DO 100 I=1,MOLSA
      DO 110 NLL=MRKR1(I),MRKR2(I)
      J=LIST(NLL)
      XIJ=X0(1,I)-X0(1,J)
      IF(XIJ.GT.+HALF)XIJ=XIJ-PBCX
      IF(XIJ.LT.-HALF)XIJ=XIJ+PBCX
      YIJ=X0(2,I)-X0(2,J)
      IF(YIJ.GT.+HALF)YIJ=YIJ-PBCY
      IF(YIJ.LT.-HALF)YIJ=YIJ+PBCY
      ZIJ=X0(3,I)-X0(3,J)
      IF(ZIJ.GT.+HALF)ZIJ=ZIJ-PBCZ
      IF(ZIJ.LT.-HALF)ZIJ=ZIJ+PBCZ
      RSQ=XIJ*(G11*XIJ+G12D*YIJ+G13D*ZIJ)+
     $YIJ*(G22*YIJ+G23D*ZIJ)+G33*ZIJ*ZIJ
      IF(RSQ.LT.RCRHO2)THEN
      SQD=(RSQ-RNRHO2)/DR2RHO+1.D0
      M=INT(SQD)
      M=MAX(M,1)
      DELRSQ=SQD-M
      RTERM=RHO(M+1)*DELRSQ+RHO(M)*(1.D0-DELRSQ)
      DNSTY(I)=DNSTY(I)+RTERM
      DNSTY(J)=DNSTY(J)+RTERM
      ENDIF
110   CONTINUE
100   CONTINUE
      DO 200 I=1,MOLSA
      SQD=(DNSTY(I)-DMIN)/DD+1.D0
      M=INT(SQD)
      M=MAX(M,1)
      M=MIN(M,KG-1)
      DELD=SQD-M
      EELD=1.D0-DELD
      DERU(I)=DUJ(M+1)*DELD+DUJ(M)*EELD
      PTERM=UJ(M+1)*DELD+UJ(M)*EELD
      PP(I)=PTERM
      PGLU=PGLU+PTERM
200   CONTINUE
      DO 300 I=1,MOLSA
      DO 310 NLL=MRKR1(I),MRKR2(I)
      J=LIST(NLL)
      XIJ=X0(1,I)-X0(1,J)
      IF(XIJ.GT.+HALF)XIJ=XIJ-PBCX
      IF(XIJ.LT.-HALF)XIJ=XIJ+PBCX
      YIJ=X0(2,I)-X0(2,J)
      IF(YIJ.GT.+HALF)YIJ=YIJ-PBCY
      IF(YIJ.LT.-HALF)YIJ=YIJ+PBCY
      ZIJ=X0(3,I)-X0(3,J)
      IF(ZIJ.GT.+HALF)ZIJ=ZIJ-PBCZ
      IF(ZIJ.LT.-HALF)ZIJ=ZIJ+PBCZ
      RSQ=XIJ*(G11*XIJ+G12D*YIJ+G13D*ZIJ)+
     $YIJ*(G22*YIJ+G23D*ZIJ)+G33*ZIJ*ZIJ
      CCELT=0.D0
      IF(RSQ.LT.RC2)THEN
      SQD=(RSQ-RN2)/DR2+1.D0
      M=INT(SQD)
      IF(M.LE.0)LEP=LEP+1
      M=MAX(M,1)
      DELRSQ=SQD-M
      EELRSQ=1.D0-DELRSQ
      PTERM=0.5D0*(VJ(M+1)*DELRSQ+VJ(M)*EELRSQ)
      CCELT=FJ(M+1)*DELRSQ+FJ(M)*EELRSQ
      PP(I)=PP(I)+PTERM
      PP(J)=PP(J)+PTERM
      ENDIF
      CCELG=0.D0
      IF(RSQ.LT.RCRHO2)THEN
      SQD=(RSQ-RNRHO2)/DR2RHO+1.D0
      M=INT(SQD)
      M=MAX(M,1)
      DELRSQ=SQD-M
      RTERM=DRHO(M+1)*DELRSQ+DRHO(M)*(1.D0-DELRSQ)
      CCELG=(DERU(I)+DERU(J))*RTERM
      ENDIF
      CCEL=CCELT+CCELG
      IF(CCEL.NE.0.D0)THEN
      RXIJ=H(1,1)*XIJ+H(1,2)*YIJ+H(1,3)*ZIJ
      RYIJ=H(2,1)*XIJ+H(2,2)*YIJ+H(2,3)*ZIJ
      RZIJ=H(3,1)*XIJ+H(3,2)*YIJ+H(3,3)*ZIJ
      XIJ=XIJ*CCEL
      YIJ=YIJ*CCEL
      ZIJ=ZIJ*CCEL
      VIRPOT(1,1)=VIRPOT(1,1)+RXIJ*XIJ
      VIRPOT(2,1)=VIRPOT(2,1)+RYIJ*XIJ
      VIRPOT(3,1)=VIRPOT(3,1)+RZIJ*XIJ
      VIRPOT(1,2)=VIRPOT(1,2)+RXIJ*YIJ
      VIRPOT(2,2)=VIRPOT(2,2)+RYIJ*YIJ
      VIRPOT(3,2)=VIRPOT(3,2)+RZIJ*YIJ
      VIRPOT(1,3)=VIRPOT(1,3)+RXIJ*ZIJ
      VIRPOT(2,3)=VIRPOT(2,3)+RYIJ*ZIJ
      VIRPOT(3,3)=VIRPOT(3,3)+RZIJ*ZIJ
      XNP(1,I)=XNP(1,I)+XIJ
      XNP(2,I)=XNP(2,I)+YIJ
      XNP(3,I)=XNP(3,I)+ZIJ
      XNP(1,J)=XNP(1,J)-XIJ
      XNP(2,J)=XNP(2,J)-YIJ
      XNP(3,J)=XNP(3,J)-ZIJ
      ENDIF
310   CONTINUE
300   CONTINUE
      ELSE
      DO 1090 I=1,MOLSA
      PP(I)=0.D0
1090  CONTINUE
      DO 1300 I=1,MOLSA
      DO 1310 NLL=MRKR1(I),MRKR2(I)
      J=LIST(NLL)
      XIJ=X0(1,I)-X0(1,J)
      IF(XIJ.GT.+HALF)XIJ=XIJ-PBCX
      IF(XIJ.LT.-HALF)XIJ=XIJ+PBCX
      YIJ=X0(2,I)-X0(2,J)
      IF(YIJ.GT.+HALF)YIJ=YIJ-PBCY
      IF(YIJ.LT.-HALF)YIJ=YIJ+PBCY
      ZIJ=X0(3,I)-X0(3,J)
      IF(ZIJ.GT.+HALF)ZIJ=ZIJ-PBCZ
      IF(ZIJ.LT.-HALF)ZIJ=ZIJ+PBCZ
      RSQ=XIJ*(G11*XIJ+G12D*YIJ+G13D*ZIJ)+
     $YIJ*(G22*YIJ+G23D*ZIJ)+G33*ZIJ*ZIJ
      IF(RSQ.LT.RC2)THEN
      SQD=(RSQ-RN2)/DR2+1.D0
      M=INT(SQD)
      IF(M.LE.0)LEP=LEP+1
      M=MAX(M,1)
      DELRSQ=SQD-M
      EELRSQ=1.D0-DELRSQ
      PTERM=0.5D0*(VJ(M+1)*DELRSQ+VJ(M)*EELRSQ)
      CCEL=FJ(M+1)*DELRSQ+FJ(M)*EELRSQ
      PP(I)=PP(I)+PTERM
      PP(J)=PP(J)+PTERM
      RXIJ=H(1,1)*XIJ+H(1,2)*YIJ+H(1,3)*ZIJ
      RYIJ=H(2,1)*XIJ+H(2,2)*YIJ+H(2,3)*ZIJ
      RZIJ=H(3,1)*XIJ+H(3,2)*YIJ+H(3,3)*ZIJ
      XIJ=XIJ*CCEL
      YIJ=YIJ*CCEL
      ZIJ=ZIJ*CCEL
      VIRPOT(1,1)=VIRPOT(1,1)+RXIJ*XIJ
      VIRPOT(2,1)=VIRPOT(2,1)+RYIJ*XIJ
      VIRPOT(3,1)=VIRPOT(3,1)+RZIJ*XIJ
      VIRPOT(1,2)=VIRPOT(1,2)+RXIJ*YIJ
      VIRPOT(2,2)=VIRPOT(2,2)+RYIJ*YIJ
      VIRPOT(3,2)=VIRPOT(3,2)+RZIJ*YIJ
      VIRPOT(1,3)=VIRPOT(1,3)+RXIJ*ZIJ
      VIRPOT(2,3)=VIRPOT(2,3)+RYIJ*ZIJ
      VIRPOT(3,3)=VIRPOT(3,3)+RZIJ*ZIJ
      XNP(1,I)=XNP(1,I)+XIJ
      XNP(2,I)=XNP(2,I)+YIJ
      XNP(3,I)=XNP(3,I)+ZIJ
      XNP(1,J)=XNP(1,J)-XIJ
      XNP(2,J)=XNP(2,J)-YIJ
      XNP(3,J)=XNP(3,J)-ZIJ
      ENDIF
1310  CONTINUE
1300  CONTINUE
      ENDIF
      PTOT=0.D0
      DO 2100 I=1,MOLSA
      PTOT=PTOT+PP(I)
2100  CONTINUE
      POT2=PTOT-PGLU
      RETURN
      END
      SUBROUTINE IRESET(IARRAY,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION IARRAY(N)
      DO 1 I=1,N
1     IARRAY(I)=0
      RETURN
      END
      SUBROUTINE RESET(ARRAY,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION ARRAY(N)
      DO 1 I=1,N
1     ARRAY(I)=0.D0
      RETURN
      END
      SUBROUTINE MTXMTP(A,B,C)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(3,3),B(3,3),C(3,3)
      DO 1 J=1,3
      C(1,J)=A(1,1)*B(1,J)+A(1,2)*B(2,J)+A(1,3)*B(3,J)
      C(2,J)=A(2,1)*B(1,J)+A(2,2)*B(2,J)+A(2,3)*B(3,J)
      C(3,J)=A(3,1)*B(1,J)+A(3,2)*B(2,J)+A(3,3)*B(3,J)
1     CONTINUE
      RETURN
      END
      SUBROUTINE MTXINV(HM,HI,DH)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION HM(3,3),HI(3,3)
      D11=HM(2,2)*HM(3,3)-HM(2,3)*HM(3,2)
      D12=HM(2,3)*HM(3,1)-HM(2,1)*HM(3,3)
      D21=HM(3,2)*HM(1,3)-HM(1,2)*HM(3,3)
      D31=HM(1,2)*HM(2,3)-HM(2,2)*HM(1,3)
      D32=HM(1,3)*HM(2,1)-HM(1,1)*HM(2,3)
      D13=HM(2,1)*HM(3,2)-HM(3,1)*HM(2,2)
      D22=HM(1,1)*HM(3,3)-HM(1,3)*HM(3,1)
      D23=HM(3,1)*HM(1,2)-HM(1,1)*HM(3,2)
      D33=HM(1,1)*HM(2,2)-HM(1,2)*HM(2,1)
      DH=HM(1,1)*D11+HM(1,2)*D12+HM(1,3)*D13
      IF(DH.LE.0.D0)THEN
      IF(DH.EQ.0.D0)THEN
      PRINT '(''0MTXINV ERROR: DH=0'')'
      STOP
      ELSE
      PRINT '(''0MTXINV WARNING: DH<0'')'
      ENDIF
      ENDIF
      HI(1,1)=D11/DH
      HI(2,2)=D22/DH
      HI(3,3)=D33/DH
      HI(1,2)=D21/DH
      HI(1,3)=D31/DH
      HI(2,3)=D32/DH
      HI(2,1)=D12/DH
      HI(3,1)=D13/DH
      HI(3,2)=D23/DH
      RETURN
      END
      SUBROUTINE MTE(NBSIZE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/COUNT/NFI,LCOUNT,LISTER,KNTSTA,KNTGOR,LEP,MANYON
      COMMON/IDENT/ELEMEN,REF,TODAY,NOW
      CHARACTER ELEMEN*10,REF*72,TODAY*10,NOW*10
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/PARAM/DELTA,DELTA2,GAMMA,VSCALE,CTRLCE
     $,CTRLMI,CTRLMA,RSQUPD,RANSQ,VMAS,BOX(3,3)
      COMMON/STATIS/FGS(NG),GRANG,FACNG,SCABY2,RESZ,DONTR,FONTR,SIG2,
     $NGS(NG),NGMAX,NZHIGH,NZLOW,MULTIP
      COMMON/SUMS/TEMPSM,TEMWSM,EKINSM,POT2SM,PGLUSM,POSTSM
     $,TOTESM,DENSSM,ALSM,VOLUSM,AREASM,HEIGSM
      COMMON/THRU/ATMASS,ECOH,R0,SPAREF
      DIMENSION H(3,3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      PARAMETER(HALF=0.5D0 , TWO=2.D0)
      PARAMETER(PI=3.141592653589793D0)
      NFI=0
      KNTSTA=0
      KNTGOR=0
      DELTA=0.D0
      TEMPSM=0.D0
      TEMWSM=0.D0
      EKINSM=0.D0
      POT2SM=0.D0
      PGLUSM=0.D0
      POSTSM=0.D0
      TOTESM=0.D0
      DENSSM=0.D0
      ALSM=0.D0
      VOLUSM=0.D0
      AREASM=0.D0
      HEIGSM=0.D0
      CALL RESET(BOX,9)
      CALL RESET(X0,3*(NM+3))
      CALL RESET(X,3*(NM+3)*5)
      CALL RESET(XIN,3*(NM+3))
      ELEMEN='GOLD'
      ALAT=4.0704D0
      ATMASS=196.967D0
      ECOH=3.78D0
      SPAREF=ALAT
      R0=SPAREF/SQRT(TWO)
      LPBC(1)=1
      LPBC(2)=1
      LPBC(3)=1
      LPBCSM=LPBC(1)+LPBC(2)+LPBC(3)
      CALL CRYSTL(R0,NBSIZE)
      NGMAX=NG
      NZHIGH=NH
      NZLOW=NL
      CALL DEFLTS(LPBC,H,R0,MOLSA,SCADEF,GRAMAX)
      SCALE=SCADEF
      GRANG=GRAMAX
      SCABY2=SCALE/TWO
      RESZ=NZHIGH/SCALE
      MULTIP=NZHIGH/NZLOW
      FACNG=NGMAX/(GRANG**2)
      DONTR=PI/(FACNG*SQRT(FACNG))
      FONTR=HALF*PI/FACNG
      REF='IN-MEMORY GENERATED SAMPLE FOR BENCHMARKING'
      TODAY='*****NEW '
      NOW='SAMPLE***'
      AMP=0.5D0
      CALL RANPOS(AMP)
      CALL COPYIN
      CALL CENTCM
      RETURN
      END
      SUBROUTINE CRYSTL(R0,NBSIZE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      COMMON/PARAM/DELTA,DELTA2,GAMMA,VSCALE,CTRLCE
     $,CTRLMI,CTRLMA,RSQUPD,RANSQ,VMAS,BOX(3,3)
      DIMENSION H(3,3)
      EQUIVALENCE(X0(1,-2),H(1,1))
      PARAMETER(NREC=7)
      PARAMETER(NBASE=4)
      DIMENSION MM(3,NBASE,NREC),CMUL(3,NREC),NPLANE(NREC)
      CHARACTER*72 NAME(0:NREC)
      PARAMETER(HALF=0.5D0 , ONE=1.D0 , TWO=2.D0 ,
     $FOUR=4.D0  , TWELVE=12.D0)
      PARAMETER(SQRT2=1.41421356237310D0 ,
     $SQRT3=1.73205080756888D0 ,
     $SQRT8=TWO*SQRT2 ,
     $SQ8BY3=SQRT8/SQRT3 ,
     $SQ4BY3=TWO/SQRT3 ,
     $SQ16B3=FOUR/SQRT3)
      DATA MM/0,0,0,   6,6,0,   0,6,6,   6,0,6,
     $0,0,0,   6,6,3,   0,0,6,   6,6,9,
     $0,0,0,   6,6,3,   0,0,6,   6,6,9,
     $0,0,0,   6,6,0,   0,4,6,   6,10,6,
     $0,0,0,   6,6,0,   0,4,6,   6,10,6,
     $0,0,0,   6,6,3,   0,0,6,   6,6,9,
     $0,0,0,   6,6,3,   0,0,6,   6,6,9/
      DATA CMUL/SQRT2 , SQRT2 ,SQRT2 ,
     $ONE , ONE , SQRT8 ,
     $ONE , SQRT2 , TWO ,
     $ONE , SQRT3 , SQ8BY3 ,
     $ONE , SQRT3 , SQ8BY3 ,
     $ONE , ONE , SQRT8 ,
     $SQ4BY3 , SQ4BY3 , SQ16B3/
      DATA NPLANE/2 , 4 , 4 , 2 , 2 , 4 , 4/
      DATA NAME/
     $'READ COORDINATES FROM EXTERNAL FILE' ,
     $'FCC 100 , LATERAL FACES 100' ,
     $'FCC 100 , LATERAL FACES 110' ,
     $'FCC 110 , CLEAN OR RECONSTRUCTED' ,
     $'FCC 111 , CLEAN OR RECONSTRUCTED' ,
     $'HCP HEXAGONAL TOP FACE' ,
     $'FCC 100 , LATERAL FACES 110 , TOP LAYER RECONSTRUCTED' ,
     $'BCC 100 , LATERAL FACES 100'/
2     FORMAT(1X,8A)
3     FORMAT(/,1X,8A)
4     FORMAT(1X,8A)
5     FORMAT(/,1X,8A)
6     FORMAT(4X,8A)
      NBX=NBSIZE
      NBY=NBSIZE
      NBZ=NBSIZE
      NSTR=1
      BOX(1,1)=NBX*CMUL(1,NSTR)
      BOX(2,2)=NBY*CMUL(2,NSTR)
      BOX(3,3)=NBZ*CMUL(3,NSTR)
      NPLA=NBZ*NPLANE(NSTR)
      M=0
      DO 50 K=1,NBZ
      DO 60 L=1,NBASE
      DO 70 J=1,NBY
      DO 80 I=1,NBX
      M=M+1
      IF(M.GT.NM)GOTO 9950
      X0(1,M)=((I-1)+MM(1,L,NSTR)/TWELVE)/NBX
      X0(2,M)=((J-1)+MM(2,L,NSTR)/TWELVE)/NBY
      X0(3,M)=((K-1)+MM(3,L,NSTR)/TWELVE)/NBZ
80    CONTINUE
70    CONTINUE
60    CONTINUE
50    CONTINUE
      GOTO 9951
9950  CONTINUE
      PRINT '(/,1X,A,I7,A,/)',
     $'***ROOM FOR',NM,
     $' PARTICLES ONLY: CRYSTAL TRUNCATED.***'
      M=M-1
9951  CONTINUE
      MOLSA=M
      MOLSP=MOLSA
      H(1,1)=R0*BOX(1,1)
      H(2,2)=R0*BOX(2,2)
      H(3,3)=R0*BOX(3,3)
      CALL COPYIN
      CALL CENTCM
      RETURN
      END
      SUBROUTINE CENTCM
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      CM1=0.D0
      CM2=0.D0
      CM3=0.D0
      DO 1 I=1,MOLSA
      CM1=CM1+X0(1,I)
      CM2=CM2+X0(2,I)
      CM3=CM3+X0(3,I)
1     CONTINUE
      CM1=CM1/MOLSA
      CM2=CM2/MOLSA
      CM3=CM3/MOLSA
      IF((CM1.EQ.0.D0).AND.(CM2.EQ.0.D0).AND.(CM3.EQ.0.D0))
     $RETURN
      DO 2 I=1,MOLSA
      X0(1,I)=X0(1,I)-CM1
      X0(2,I)=X0(2,I)-CM2
      X0(3,I)=X0(3,I)-CM3
      XIN(1,I)=XIN(1,I)-CM1
      XIN(2,I)=XIN(2,I)-CM2
      XIN(3,I)=XIN(3,I)-CM3
2     CONTINUE
      RETURN
      END
      SUBROUTINE RANPOS(AMP)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      EQUIVALENCE(X0(1,-2),H(1,1))
      DIMENSION H(3,3),HI(3,3),DELTA(3)
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      EXTERNAL RANFM
      DOUBLE PRECISION RANFM
      CALL MTXINV(H,HI,DH)
      IDUM=-1
      DO 10 I=1,MOLSP
      DO 20 K=1,3
      DELTA(K)=2.D0*AMP*(RANFM(IDUM)-0.5D0)
20    CONTINUE
      DO 30 K=1,3
      X0(K,I)=X0(K,I)+
     $HI(K,1)*DELTA(1)+HI(K,2)*DELTA(2)+HI(K,3)*DELTA(3)
30    CONTINUE
10    CONTINUE
      RETURN
      END
      SUBROUTINE COPYIN
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NM=16384)
      PARAMETER(NG=100)
      PARAMETER(NH=100)
      PARAMETER(MU=20)
      PARAMETER(NL=1)
      PARAMETER(LL=10*NM)
      PARAMETER(KP=2001,KR=2001,KG=2001)
      COMMON/LCS/X0(3,-2:NM),X(3,-2:NM,5),XIN(3,-2:NM)
      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM
      DO 1 I=-2,MOLSA
      XIN(1,I)=X0(1,I)
      XIN(2,I)=X0(2,I)
      XIN(3,I)=X0(3,I)
1     CONTINUE
      RETURN
      END
      SUBROUTINE DEFLTS(LPBC,H,R0,MOLSA,SCADEF,GRAMAX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION LPBC(3),H(3,3)
      PARAMETER(HALF=0.5D0 , TWO=2.D0)
      PARAMETER(THIRD=1.D0/3.D0)
      PARAMETER(PI=3.141592653589793D0)
      PARAMETER(SQRT2=1.41421356237310D0)
      SCADEF=1.25D0*H(3,3)
      GRAMAX=MAX(HALF*H(1,1), HALF*H(2,2), HALF*H(3,3))
      DO 100 J=1,3
      IF(LPBC(J).NE.0)THEN
      GRAMAJ=HALF*H(J,J)
      IF(GRAMAX.GT.GRAMAJ)GRAMAX=GRAMAJ
      ENDIF
100   CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION RANFM(IDUM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(M=714025,IA=1366,IC=150889,RM=1.4005112D-6)
      DIMENSION IR(97)
      SAVE IFF,IR,IY
      DATA IFF/0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
      IFF=1
      IDUM=MOD(IC-IDUM,M)
      DO 11 J=1,97
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
11    CONTINUE
      IDUM=MOD(IA*IDUM+IC,M)
      IY=IDUM
      ENDIF
      J=1+(97*IY)/M
      IY=IR(J)
      RANFM=IY*RM
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
      RETURN
      END
.
