C ******************************************************************* C ** FICHE F.3 ** C ** TWO SEPARATE PARTS: FORTRAN AND BASIC VERSIONS. ** C ******************************************************************* C ******************************************************************* C ** FICHE F.3 - PART A ** C ** FORTRAN PROGRAM USING THE LEAPFROG ALGORITHM. ** C ******************************************************************* PROGRAM FROGGY COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ C ******************************************************************* C ** FORTRAN PROGRAM TO CONDUCT MOLECULAR DYNAMICS OF ATOMS. ** C ** ** C ** A SPECIAL LOW-STORAGE FORM OF THE LEAPFROG VERLET ALGORITHM ** C ** IS USED AND WE TAKE THE LENNARD-JONES POTENTIAL AS AN EXAMPLE.** C ** ** C ** REFERENCE: ** C ** ** C ** FINCHAM AND HEYES, CCP5 QUARTERLY, 6, 4, 1982. ** C ** ** C ** ROUTINES REFERENCED: ** C ** ** C ** SUBROUTINE READCN ( CNFILE ) ** C ** READS IN CONFIGURATION ** C ** SUBROUTINE FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW ) ** C ** CALCULATES THE ACCELERATIONS AND ADVANCES THE VELOCITIES ** C ** FROM T - DT/2 TO T + DT/2. ALSO CALCULATES POTENTIAL ** C ** ENERGY AND VIRIAL AT TIMESTEP T. ** C ** SUBROUTINE MOVE ( DT ) ** C ** ADVANCES POSITIONS FROM T TO T + DT. ** C ** SUBROUTINE KINET ( NEWK ) ** C ** CALCULATES KINETIC ENERGY. ** C ** SUBROUTINE WRITCN ( CNFILE ) ** C ** WRITES OUT CONFIGURATION ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** INTEGER N NUMBER OF ATOMS ** C ** REAL DT TIMESTEP ** C ** REAL RX(N),RY(N),RZ(N) ATOMIC POSITIONS ** C ** REAL VX(N),VY(N),VZ(N) ATOMIC VELOCITIES ** C ** REAL ACV,ACK ETC. AVERAGE VALUE ACCUMULATORS ** C ** REAL AVV,AVK ETC. AVERAGE VALUES ** C ** REAL ACVSQ, ACKSQ ETC. AVERAGE SQUARED VALUE ACCUMULATORS ** C ** REAL FLV,FLK ETC. FLUCTUATION AVERAGES ** C ** ** C ** USAGE: ** C ** ** C ** THE LEAPFROG ALGORITHM IS OF THE FORM ** C ** VX(T + DT/2) = VX(T - DT/2) + DT * AX(T) (SIMILARLY FOR Y,Z) ** C ** RX(T + DT) = RX(T) + DT * VX(T + DT/2) (SIMILARLY FOR Y,Z) ** C ** TO SAVE STORAGE IN THIS PROGRAM WE ACCUMULATE VALUES AX(T) ** C ** DIRECTLY ONTO THE VELOCITIES VX(T - DT/2) IN THE FORCE LOOP. ** C ** THIS MEANS THAT AN APPROXIMATION MUST BE USED FOR THE KINETIC ** C ** ENERGY AT EACH TIME STEP: ** C ** K = ( OLDK + NEWK ) / 2 + ( OLDV - 2 * V + NEWV ) / 8 ** C ** WHERE K = K( T ), V = V( T ) ** C ** OLDK = K( T - DT/2 ), OLDV = V( T - DT ) ** C ** NEWK = K( T + DT/2 ), NEWV = V( T + DT ) ** C ** AT THE START OF A STEP THE FOLLOWING VARIABLES ARE STORED: ** C ** R : R(STEP) V : V(STEP+1/2) ** C ** OLDV : V(STEP-2) V : V(STEP-1) NEWV : V(STEP) ** C ** OLDK : K(STEP-1/2) NEWK : K(STEP+1/2) ** C ** THIS PROGRAM USES UNIT 10 FOR CONFIGURATION INPUT AND OUTPUT ** C ** ** C ** UNITS: ** C ** ** C ** THE PROGRAM USES LENNARD-JONES REDUCED UNITS FOR USER INPUT ** C ** AND OUTPUT BUT CONDUCTS SIMULATION IN A BOX OF UNIT LENGTH. ** C ** SUMMARY FOR BOX LENGTH L, ATOMIC MASS M, AND LENNARD-JONES ** C ** POTENTIAL PARAMETERS SIGMA AND EPSILON: ** C ** OUR PROGRAM LENNARD-JONES SYSTEM ** C ** LENGTH L SIGMA ** C ** MASS M M ** C ** ENERGY EPSILON EPSILON ** C ** TIME SQRT(M*L**2/EPSILON) SQRT(M*SIGMA**2/EPSILON) ** C ** VELOCITY SQRT(EPSILON/M) SQRT(EPSILON/M) ** C ** PRESSURE EPSILON/L**3 EPSILON/SIGMA**3 ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) INTEGER STEP, NSTEP, IPRINT REAL DENS, NORM, RCUT, DT, SIGMA REAL V, K, E, W REAL OLDK, NEWK, OLDV, NEWV, NEWW REAL OLDVC, NEWVC, VC, KC, EC REAL VN, KN, EN, ECN, PRES, TEMP REAL ACV, ACK, ACE, ACEC, ACP, ACT REAL AVV, AVK, AVE, AVEC, AVP, AVT REAL ACVSQ, ACKSQ, ACESQ, ACECSQ, ACPSQ, ACTSQ REAL FLV, FLK, FLE, FLEC, FLP, FLT REAL SR3, SR9, VLRC, WLRC, PI, SIGCUB CHARACTER CNFILE*30, TITLE*80 REAL FREE PARAMETER ( FREE = 3.0 ) PARAMETER ( PI = 3.1415927 ) C ******************************************************************* C ** READ IN INITIAL PARAMETERS ** WRITE(*,'(1H1,'' **** PROGRAM FROGGY **** '')') WRITE(*,'(//, '' MOLECULAR DYNAMICS OF LENNARD-JONES ATOMS '')') WRITE(*,'( '' LEAPFROG ALGORITHM WITH MINIMAL STORAGE '')') WRITE(*,'('' ENTER RUN TITLE '')') READ (*,'(A)') TITLE WRITE(*,'('' ENTER NUMBER OF STEPS '')') READ (*,*) NSTEP WRITE(*,'('' ENTER NUMBER OF STEPS BETWEEN OUTPUT LINES '')') READ (*,*) IPRINT WRITE(*,'('' ENTER CONFIGURATION FILENAME '')') READ (*,'(A)') CNFILE WRITE(*,'('' ENTER THE FOLLOWING IN LENNARD-JONES UNITS '')') WRITE(*,'('' ENTER DENSITY '')') READ (*,*) DENS WRITE(*,'('' ENTER POTENTIAL CUTOFF DISTANCE '')') READ (*,*) RCUT WRITE(*,'('' ENTER TIME STEP '')') READ (*,*) DT WRITE(*,'(//1X,A)') TITLE WRITE(*,'('' NUMBER OF ATOMS = '',I10 )') N WRITE(*,'('' NUMBER OF STEPS = '',I10 )') NSTEP WRITE(*,'('' OUTPUT FREQUENCY = '',I10 )') IPRINT WRITE(*,'('' POTENTIAL CUTOFF = '',F10.4)') RCUT WRITE(*,'('' DENSITY = '',F10.4)') DENS WRITE(*,'('' TIME STEP = '',F10.6)') DT C ** READ CONFIGURATION INTO COMMON / BLOCK1 / VARIABLES ** CALL READCN ( CNFILE ) C ** CONVERT INPUT DATA TO PROGRAM UNITS ** SIGMA = ( DENS / REAL ( N ) ) ** ( 1.0 / 3.0 ) RCUT = RCUT * SIGMA DT = DT * SIGMA DENS = DENS / ( SIGMA ** 3 ) C ** CALCULATE LONG-RANGE CORRECTIONS ** C ** NOTE: SPECIFIC TO LENNARD-JONES ** SR3 = ( SIGMA / RCUT ) ** 3 SR9 = SR3 ** 3 SIGCUB = SIGMA ** 3 VLRC = ( 8.0 /9.0 ) * PI * DENS * SIGCUB * REAL ( N ) : * ( SR9 - 3.0 * SR3 ) WLRC = ( 16.0 / 9.0 ) * PI * DENS * SIGCUB * REAL ( N ) : * ( 2.0 * SR9 - 3.0 * SR3 ) C ** ZERO ACCUMULATORS ** ACV = 0.0 ACK = 0.0 ACE = 0.0 ACEC = 0.0 ACP = 0.0 ACT = 0.0 ACVSQ = 0.0 ACKSQ = 0.0 ACESQ = 0.0 ACECSQ = 0.0 ACPSQ = 0.0 ACTSQ = 0.0 FLV = 0.0 FLK = 0.0 FLE = 0.0 FLEC = 0.0 FLP = 0.0 FLT = 0.0 C ** CALCULATE INITIAL VALUES ** CALL FORCE ( -DT, SIGMA, RCUT, NEWV, NEWVC, NEWW ) CALL MOVE ( -DT ) CALL FORCE ( -DT, SIGMA, RCUT, V, VC, W ) CALL FORCE ( DT, SIGMA, RCUT, V, VC, W ) CALL KINET ( OLDK ) CALL MOVE ( DT ) CALL FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW ) CALL KINET ( NEWK ) C ** INCLUDE LONG-RANGE CORRECTIONS ** V = V + VLRC W = W + WLRC NEWV = NEWV + VLRC NEWW = NEWW + WLRC IF ( IPRINT .LE. 0 ) IPRINT = NSTEP + 1 WRITE(*,'(//1X,''**** START OF DYNAMICS ****'')') WRITE(*,10001) C ******************************************************************* C ** MAIN LOOP BEGINS ** C ******************************************************************* DO 1000 STEP = 1, NSTEP C ** IMPLEMENT ALGORITHM ** CALL MOVE ( DT ) OLDV = V V = NEWV OLDVC = VC VC = NEWVC W = NEWW CALL FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW ) C ** INCLUDE LONG-RANGE CORRECTIONS ** NEWV = NEWV + VLRC NEWW = NEWW + WLRC C ** CALCULATE KINETIC ENERGY AT CURRENT STEP ** K = ( NEWK + OLDK ) / 2.0 : + ( NEWV - 2.0 * V + OLDV ) / 8.0 KC = ( NEWK + OLDK ) / 2.0 : + ( NEWVC - 2.0 * VC + OLDVC ) / 8.0 OLDK = NEWK CALL KINET ( NEWK ) C ** CALCULATE INSTANTANEOUS VALUES ** E = K + V EC = KC + VC VN = V / REAL ( N ) KN = K / REAL ( N ) EN = E / REAL ( N ) ECN = EC / REAL ( N ) TEMP = 2.0 * KN / FREE PRES = DENS * TEMP + W C ** CONVERT TO LENNARD-JONES UNITS ** PRES = PRES * SIGMA ** 3 C ** INCREMENT ACCUMULATORS ** ACE = ACE + EN ACEC = ACEC + ECN ACK = ACK + KN ACV = ACV + VN ACP = ACP + PRES ACESQ = ACESQ + EN ** 2 ACECSQ = ACECSQ + ECN ** 2 ACKSQ = ACKSQ + KN ** 2 ACVSQ = ACVSQ + VN ** 2 ACPSQ = ACPSQ + PRES ** 2 C ** OPTIONALLY PRINT INFORMATION ** IF ( MOD( STEP, IPRINT ) .EQ. 0 ) THEN WRITE(*,'(1X,I8,6(2X,F10.4))') : STEP, EN, ECN, KN, VN, PRES, TEMP ENDIF 1000 CONTINUE C ******************************************************************* C ** MAIN LOOP ENDS ** C ******************************************************************* WRITE(*,'(/1X,''**** END OF DYNAMICS **** ''//)') C ** OUTPUT RUN AVERAGES ** NORM = REAL ( NSTEP ) AVE = ACE / NORM AVEC = ACEC / NORM AVK = ACK / NORM AVV = ACV / NORM AVP = ACP / NORM ACESQ = ( ACESQ / NORM ) - AVE ** 2 ACECSQ = ( ACECSQ / NORM ) - AVEC ** 2 ACKSQ = ( ACKSQ / NORM ) - AVK ** 2 ACVSQ = ( ACVSQ / NORM ) - AVV ** 2 ACPSQ = ( ACPSQ / NORM ) - AVP ** 2 IF ( ACESQ .GT. 0.0 ) FLE = SQRT ( ACESQ ) IF ( ACECSQ .GT. 0.0 ) FLEC = SQRT ( ACECSQ ) IF ( ACKSQ .GT. 0.0 ) FLK = SQRT ( ACKSQ ) IF ( ACVSQ .GT. 0.0 ) FLV = SQRT ( ACVSQ ) IF ( ACPSQ .GT. 0.0 ) FLP = SQRT ( ACPSQ ) AVT = AVK * 2.0 / FREE FLT = FLK * 2.0 / FREE WRITE(*,'('' AVERAGES'',6(2X,F10.5))') : AVE, AVEC, AVK, AVV, AVP, AVT WRITE(*,'('' FLUCTS '',6(2X,F10.5))') : FLE, FLEC, FLK, FLV, FLP, FLT C ** WRITE OUT FINAL CONFIGURATION ** CALL WRITCN ( CNFILE ) STOP 10001 FORMAT(//1X,'TIMESTEP ..ENERGY.. CUTENERGY.' : ' ..KINETIC. ..POTENT..', : ' .PRESSURE. ..TEMPER..'/) END SUBROUTINE FORCE ( DT, SIGMA, RCUT, V, VC, W ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ C ******************************************************************* C ** ROUTINE TO COMPUTE FORCES WITH CUTOFF & MINIMUM IMAGING. ** C ** ** C ** THE ROUTINE ACTUALLY RETURNS TWO DIFFERENT POTENTIAL ENERGIES.** C ** V IS CALCULATED USING THE UNSHIFTED LENNARD-JONES POTENTIAL. ** C ** WHEN LONG-RANGE TAIL CORRECTIONS ARE ADDED, THIS MAY BE USED ** C ** TO CALCULATE THERMODYNAMIC INTERNAL ENERGY ETC. ** C ** VC IS CALCULATED USING THE SHIFTED LENNARD-JONES POTENTIAL ** C ** WITH NO DISCONTINUITY AT THE CUTOFF. THIS MAY BE USED TO ** C ** CHECK ENERGY CONSERVATION. ** C ** PROGRAM UNITS MAKE EPSILON = MASS = 1, BUT SIGMA IS NOT UNITY ** C ** SINCE WE TAKE A BOX OF UNIT LENGTH. ** C ** THE ROUTINE IS COMBINED WITH THE LEAPFROG ALGORITHM FOR ** C ** ADVANCING VELOCITIES, I.E. FORCES ARE ACCUMULATED DIRECTLY ** C ** ONTO VELOCITIES. ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) REAL V, RCUT, DT, SIGMA, W, VC INTEGER I, J, NCUT REAL RCUTSQ, VCUT, SIGSQ REAL RXI, RYI, RZI, VXI, VYI, VZI REAL RXIJ, RYIJ, RZIJ, RIJSQ REAL SR2, SR6, SR12 REAL VIJ, WIJ, VELIJ, DVX, DVY, DVZ C ******************************************************************* C ** CALCULATE SQUARED DISTANCES FOR INNER LOOP ** SIGSQ = SIGMA ** 2 RCUTSQ = RCUT ** 2 C ** CONDUCT OUTER LOOP OVER ATOMS ** V = 0.0 W = 0.0 NCUT = 0 DO 1000 I = 1, N - 1 RXI = RX(I) RYI = RY(I) RZI = RZ(I) VXI = VX(I) VYI = VY(I) VZI = VZ(I) C ** CONDUCT INNER LOOP OVER ATOMS ** DO 900 J = I + 1, N RXIJ = RXI - RX(J) RXIJ = RXIJ - ANINT ( RXIJ ) IF ( ABS ( RXIJ ) .LT. RCUT ) THEN RYIJ = RYI - RY(J) RYIJ = RYIJ - ANINT ( RYIJ ) RIJSQ = RXIJ ** 2 + RYIJ ** 2 IF ( RIJSQ .LT. RCUTSQ ) THEN RZIJ = RZI - RZ(J) RZIJ = RZIJ - ANINT ( RZIJ ) RIJSQ = RIJSQ + RZIJ ** 2 IF ( RIJSQ .LT. RCUTSQ ) THEN SR2 = SIGSQ / RIJSQ SR6 = SR2 ** 3 SR12 = SR6 ** 2 VIJ = 4.0 * ( SR12 - SR6 ) WIJ = 24.0 * ( 2.0 * SR12 - SR6 ) VELIJ = WIJ * DT / RIJSQ DVX = VELIJ * RXIJ DVY = VELIJ * RYIJ DVZ = VELIJ * RZIJ VXI = VXI + DVX VYI = VYI + DVY VZI = VZI + DVZ VX(J) = VX(J) - DVX VY(J) = VY(J) - DVY VZ(J) = VZ(J) - DVZ V = V + VIJ W = W + WIJ NCUT = NCUT + 1 ENDIF ENDIF ENDIF 900 CONTINUE C ** END OF INNER LOOP ** VX(I) = VXI VY(I) = VYI VZ(I) = VZI 1000 CONTINUE C ** END OF OUTER LOOP ** C ** CALCULATE POTENTIAL SHIFT ** SR2 = SIGSQ / RCUTSQ SR6 = SR2 * SR2 * SR2 SR12 = SR6 * SR6 VIJ = 4.0 * ( SR12 - SR6 ) VC = V - REAL ( NCUT ) * VIJ W = W / 3.0 RETURN END SUBROUTINE READCN ( CNFILE ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ C ******************************************************************* C ** SUBROUTINE TO READ IN INITIAL CONFIGURATION FROM UNIT 10 ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) CHARACTER CNFILE*(*) REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) INTEGER CNUNIT, NN PARAMETER ( CNUNIT = 10 ) C ****************************************************************** OPEN ( UNIT = CNUNIT, FILE = CNFILE, STATUS = 'OLD', : FORM = 'UNFORMATTED' ) READ ( CNUNIT ) NN IF ( NN .NE. N ) STOP ' INCORRECT NUMBER OF ATOMS ' READ ( CNUNIT ) RX, RY, RZ READ ( CNUNIT ) VX, VY, VZ CLOSE ( UNIT = CNUNIT ) RETURN END SUBROUTINE WRITCN ( CNFILE ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ C ******************************************************************* C ** ROUTINE TO WRITE OUT FINAL CONFIGURATION TO UNIT 10 ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) CHARACTER CNFILE*(*) REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) INTEGER CNUNIT PARAMETER ( CNUNIT = 10 ) C **************************************************************** OPEN ( UNIT = CNUNIT, FILE = CNFILE, STATUS = 'UNKNOWN', : FORM = 'UNFORMATTED' ) WRITE ( CNUNIT ) N WRITE ( CNUNIT ) RX, RY, RZ WRITE ( CNUNIT ) VX, VY, VZ CLOSE ( UNIT = CNUNIT ) RETURN END SUBROUTINE MOVE ( DT ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ C ******************************************************************* C ** LEAPFROG VERLET ALGORITHM FOR ADVANCING POSITIONS ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) REAL DT INTEGER I C ******************************************************************* DO 1000 I = 1, N RX(I) = RX(I) + VX(I) * DT RY(I) = RY(I) + VY(I) * DT RZ(I) = RZ(I) + VZ(I) * DT 1000 CONTINUE RETURN END SUBROUTINE KINET ( K ) COMMON / BLOCK1 / RX, RY, RZ, VX, VY, VZ C ******************************************************************* C ** COMPUTES KINETIC ENERGY ** C ******************************************************************* INTEGER N PARAMETER ( N = 108 ) REAL RX(N), RY(N), RZ(N) REAL VX(N), VY(N), VZ(N) REAL K INTEGER I C ******************************************************************* K = 0.0 DO 100 I = 1, N K = K + VX(I) ** 2 + VY(I) ** 2 + VZ(I) ** 2 100 CONTINUE K = K * 0.5 RETURN END C ******************************************************************* C ** FICHE F.3 - PART B ** C ** BASIC PROGRAM USING THE LEAPFROG VERLET ALGORITHM. ** C ******************************************************************* C ******************************************************************* C ** THE FOLLOWING PROGRAMS WERE WRITTEN ON AN ACORN MODEL B MICRO,** C ** USING BBC BASIC AND SOME MACHINE CODE FOR THE GRAPHICS. ** C ** CONSEQUENTLY MINOR ALTERATIONS MAY BE NEEDED FOR OTHER MICROS.** C ** THE FIRST PROGRAM RUNS MOLECULAR DYNAMICS FROM AN INITIAL ** C ** CONFIGURATION WHICH MAY BE GENERATED USING THE SECOND PROGRAM.** C ** DEFAULT INPUT PARAMETERS ARE PROVIDED IN THE MD PROGRAM: ** C ** THESE CAN BE ALTERED TO SUIT THE USER. ** C ******************************************************************* MOLECULAR DYNAMICS PROGRAM 10 IF PAGE<>&1700 THEN PAGE=&1700:CHAIN"2.FROG1" 20 MODE1 30 CLS:PRINT''''TAB(12);"Program Froggy" 40PRINT'TAB(10);"Molecular Dynamics"' 50PRINTTAB(7);" of Lennard-Jones atoms" 60 PRINT'TAB(1);"Leapfrog algorithm with minimal storage" 70 PRINT'TAB(7);"Repulsive WCA potential"' 80 PRINTTAB(9);"in two dimensions" 90 PRINT'TAB(15);"Fiche 3b" 100 INPUTTAB(10,31)"( press )"A$ 110 CLS 120*KEY0"100" 130*KEY1"2.DATA" 140*KEY2"0.3" 150*KEY3"0.02" 160 N=10:FREE=2:x1=200:y1=200:x2=512:y2=512 170 HIMEM=&2E00 180DIM X(12),Y(12),RX(N),RY(N),VX(N),VY(N) 190HIMEM=&2E00:GCOL 3,4:VDU 19,2,1,0,0,0 200 *FX138,0,128 210NSTEP%= FNgetnum( "Enter number of steps") 220 *FX138,0,129 230PRINT"Enter configuration filename ";:INPUTCNFILE$ 240PRINT'TAB(5);"IN Lennard-Jones units," 250 *FX138,0,130 260DENS= FNgetnum("Enter density") 270 *FX138,0,131 280DT= FNgetnum("Enter time step") 290PROCreadcn(CNFILE$):SIGMA=(DENS/N)^(1/2):RCUT=SIGMA*(2.0^(1.0/6.0)) 300DT=DT*SIGMA 310CLS 320R=(DENS*(x2-x1)*(y2-y1)/N)^0.5/2 330PROCcircle(R) 340PROCcode 350PROCforce( -DT,SIGMA,RCUT ) 360PROCmove( -DT ) 370PROCforce( -DT,SIGMA,RCUT ) 380PROCforce( DT,SIGMA,RCUT ) 390PROCmove( DT ) 400PROCforce(DT,SIGMA,RCUT ) 410GCOL 3,1 420PROCbox(x1-2*R,y1,x1+x2-2*R,y1+y2) 430GCOL 3,3 440PROCgraphics(&2F00,FALSE) 450PROCerase 460FOR STP%=1 TO NSTEP% 470PROCmove( DT ) 480PROCforce(DT,SIGMA,RCUT ) 490PROCgraphics(&2F04,TRUE) 500NEXT STP% 510END 520DEF FNgetnum(S$) 530PRINT S$;" ";:INPUT A 540=A 550DEF PROCreadcn( FILE$ ) 560LOCAL PORT 570PORT=OPENUP( FILE$ ) 580IF PORT=0 THEN PRINT''CHR$(34);FILE$;CHR$(34);" - File not found"':END 590FOR I%= 1 TO N:INPUT#PORT,RX(I%),RY(I%):NEXT I% 600FOR I%= 1 TO N:INPUT#PORT,VX(I%),VY(I%):NEXT I% 610CLOSE #PORT 620DEF PROCmove( DT ) 630FOR I%= 1 TO N 640RX(I%)= RX(I%)+VX(I%)*DT 650IF RX(I%)>0.5 REPEAT:RX(I%)=RX(I%)-1:UNTIL RX(I%)<0.5 660IF RX(I%)<-0.5 REPEAT:RX(I%)=RX(I%)+1:UNTIL RX(I%)>-0.5 670RY(I%)= RY(I%)+VY(I%)*DT 680IF RY(I%)>0.5 REPEAT:RY(I%)=RY(I%)-1:UNTIL RY(I%)<0.5 690IF RY(I%)<-0.5 REPEAT:RY(I%)=RY(I%)+1:UNTIL RY(I%)>-0.5 700NEXT I% 710ENDPROC 720DEF PROCforce( DT,SIGMA,RCUT ) 730A=SIGMA*SIGMA:B=RCUT*RCUT 740FOR I%=1 TO N-1 750C=RX(I%):D=RY(I%):H=VX(I%):I=VY(I%) 760FOR J%= I%+1 TO N 770E=C-RX(J%):E=E-INT(E+0.5) 780IF ABS(E)>=RCUT THEN 830 790F=D-RY(J%):F=F-INT(F+0.5):G=E*E+F*F 800IF G>=B THEN 830 810S=A/G:T=S*S*S:U=T*T:W=24*(U+U-S):J=W*DT/G:M=J*E:O=J*F:H=H+M:I=I+O 820VX(J%)=VX(J%)-M:VY(J%)=VY(J%)-O 830NEXT J% 840VX(I%)=H:VY(I%)=I 850NEXT I% 860ENDPROC 870DEF PROCcircle(R) 880D=100:T=0:X(1)=R*SIN(T):Y(1)=SQR(R^2-X(1)^2) 890FOR A=2 TO 12 900X(A)=X(A-1)*COS(D)+Y(A-1)*SIN(D):Y(A)=Y(A-1)*COS(D)-X(A-1)*SIN(D) 910NEXT 920B%=&2E00 930FOR A%=1 TO 12 940?B%=25:B%?1=1:B%?2=X(A%) MOD 256:B%?4=Y(A%) MOD 256 950IF X(A%)<0 B%?3=255-(X(A%) DIV 256) ELSE B%?3=X(A%) DIV 256 960IF Y(A%)<0 B%?5=255-(Y(A%) DIV 256) ELSE B%?5=Y(A%) DIV 256 970B%=B%+6 980NEXT 990?&2E17=0 1000ENDPROC 1010DEF PROCcode 1020P%=&1100 1030FOR I%=0 TO 2 STEP2 1040\OPT I% 1050 .circle LDY#0 1060.loop LDA &2E00,Y 1070JSR &FFEE 1080INY 1090CPY #72 1100BNE loop 1110RTS 1120.code LDX #20 1130 LDA #&00 1140 STA &70 1150 LDA #&2F 1160 STA &71 1170 .l2 JSR getnum 1180 JSR pl2 1190 JSR image 1200 INC &70 1210 INC &70 1220 INC &70 1230 INC &70 1240 DEX 1250 BNE l2 1260 RTS 1270.getnum LDY#3 1280.getnum1 LDA(&70),Y 1290 STA &72,Y 1300 DEY 1310 BPL getnum1 1320 RTS 1330.pl2 LDA #25 1340 JSR &FFEE 1350 LDA #4 1360 JSR &FFEE 1370 LDY #0 1380.pl21 LDA &72,Y 1390 JSR &FFEE 1400 INY 1410 CPY #4 1420 BNE pl21 1430 JSR circle 1440 RTS 1450 1460.image 1470 LDA &73 1480 ADC #1 1490 STA &73 1500 CMP #3 1510 BPL P%+5 1520 JSR pl2 1530 JSR getnum 1540 .n7 LDA &75 1550 ADC #1 1560 STA &75 1570 CMP #3 1580 BPL P%+5 1590 JSR pl2 1600JSR getnum 1610LDA &75 1620SEC 1630SBC #2 1640STA &75 1650CMP #0 1660BNE P%+11 1670 LDA &74 1680 CMP #&80 1690 BMI P%+5 1700JSR pl2 1710JSR getnum 1720LDA &73 1730SEC 1740SBC #2 1750STA &73 1760CMP #0 1770BNE P%+11 1780 LDA &72 1790 CMP #&80 1800 BMI P%+5 1810 JSR pl2 1820 RTS 1830.go JSR code 1840 LDY #0 1850 .l3 LDA &2F04,Y 1860 STA &2F00,Y 1870 INY 1880 CPY #81 1890 BNE l3 1900 RTS 1910| 1920NEXT 1930ENDPROC 1940DEF PROCgraphics(B%,run) 1950FOR A%=1 TO 10 1960 D%=x1+(RX(A%)+0.5)*x2:C%=y1+(RY(A%)+0.5)*y2:?(B%+(A%-1)*8)=D% MOD 256 1970?(B%+(A%-1)*8+1)=D% DIV 256:?(B%+(A%-1)*8+2)=C% MOD 256: 1980?(B%+(A%-1)*8+3)=C% DIV 256 1990 NEXT 2000 IF run CALL go 2010ENDPROC 2020DEF PROCerase 2030FOR A%=1 TO 10 2040 X%=x1+(RX(A%)+0.5)*x2 2050 Y%=y1+(RY(A%)+0.5)*y2 2060MOVE x1+(RX(A%)+0.5)*x2,y1+(RY(A%)+0.5)*y2 2070CALL circle 2080 IF X%+512<256*3 MOVEX%+512,Y%:CALL circle 2090 IF Y%+512<256*3 MOVEX%,Y%+512:CALL circle 2100 IF (Y%-512<256) AND (Y%-512)>128 MOVEX%,Y%-512:CALL circle 2110 IF X%-512<256 AND X%-512>128 MOVEX%-512,Y%:CALL circle 2120NEXT 2130ENDPROC 2140DEF PROCbox(A%,B%,C%,D%) 2150MOVE A%,B%:DRAW A%,D%:DRAW C%,D%:DRAW C%,B%:DRAW A%,B% 2160ENDPROC CONFIGURATION GENERATOR 10 MODE 1 20REM PROGRAM SETUP 30REM 40N= 10: REM CONSTANT 50DIM RX(N),RY(N) 60DIM VX(N),VY(N) 70 A1= 3.949846138: A3= 0.252408784: A5= 0.076542912 80 A7= 0.008355968: A9= 0.029899776 90PROCfcc 100 *KEY0"1.75" 110 *KEY1"2.DATA" 120 130 *FX138,0,128 140PRINT'" Enter Temperature in LJ units ";:INPUT TEMP 150PROCcomvel( TEMP ) 160 170 *FX138,0,129 180PRINT'"Enter Output Data File ";:INPUT CNFILE$ 190 200PORT= OPENOUT( CNFILE$ ) 210IF PORT=0:PRINT''"Error in opening ";CHR$(34);CNFILE$;CHR$(34)'':END 220 230FOR I%=1 TO N:PRINT #PORT,RX(I%),RY(I%):NEXT I% 240FOR I%=1 TO N:PRINT #PORT,VX(I%),VY(I%):NEXT I% 250CLOSE #PORT 260END 270 280 DEF PROCfcc 290 LOCAL X,Y,M,NC,KL,KX 300 NC=INT((N/4)^(1/2) +0.5) 310 X=1/NC:Y=1 320RX(1)=0:RY(1)=0 330RX(2)=2*Y/6:RY(2)=0 340RX(3)=4*Y/6:RY(3)=0 350RX(4)=Y/6:RY(4)=Y/5 360RX(5)=3*Y/6:RY(5)=Y/5 370RX(6)=0:RY(6)=2*Y/5 380RX(7)=2*Y/6:RY(7)=2*Y/5 390RX(8)=4*Y/6:RY(8)=2*Y/5 400RX(9)=Y/6:RY(9)=3*Y/5 410RX(10)=3*Y/6:RY(10)=3*Y/5 420 FOR A%=1 TO 10:RX(A%)=RX(A%)-0.5:RY(A%)=RY(A%)-0.5:NEXT 430 ENDPROC 440DEF PROCcomvel( TEMP ) 450 460LOCAL RTEMP,SUMX,SUMY,GAUSS,DUMMY 470 480RTEMP= SQR( TEMP ) 490FOR I%= 1 TO N 500VX(I%)= RTEMP*FNgauss 510VY(I%)= RTEMP*FNgauss 520NEXT I% 530SUX= 0:SUMY= 0 540FOR I%=1 TO N 550SUMX= SUMX+VX(I%) 560SUMY= SUMY+VY(I%) 570NEXT I% 580SUMX= SUMX/N: SUMY= SUMY/N 590FOR I%=1 TO N 600VX(I%)= VX(I%)-SUMX 610VY(I%)= VY(I%)-SUMY 620NEXT I% 630ENDPROC 640DEF FNgauss 650LOCAL SUM,R,R2 660FOR L%=1 TO 12 670SUM= SUM+RND(1) 680NEXT L% 690R= (SUM-6)/4: R2= R*R 700= ((((A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R 710END