* 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