I was reading an old book about computer simulation of liquids and wrote this little script that uses the Lennard-Jones potential to calculate the random distribution of atoms in a gas lattice.
OPTION NOLET
! This code illustrates the calculation of the potential energy in a system of
! Lennard-Jones atoms.
! I create a set of atoms in random locations, calculate the Lennard-Jones potential,
! and then move the atoms in a way that would minimize this energy calculation.
N = 20 ! Number of particles
V = 0.0
EPSILON = 164
SIGSQ = .383
STEPX = 0.1
STEPY = STEPX
STEPZ = STEPY
STEP = 1
NOMOVE = 0
InitialPlotString$ = "IP = point(["
FinalPlotString$ = "FP = point(["
filename$ = "output.txt"
OPEN #1: NAME filename$, CREATE NEWOLD
ERASE #1
! Coordinate array for the particles
DIM RX(20)
DIM RY(20)
DIM RZ(20)
! Read coordinate data into arrays
PRINT "Initial atomic positions:"
PRINT "X","Y","Z"
PRINT #1: "Initial atomic positions:"
PRINT #1: "X","Y","Z"
FOR I = 1 TO N
RX(I) = RND*10
RY(I) = RND*10
RZ(I) = RND*10
PRINT RX(I),RY(I),RZ(I)
PRINT #1:RX(I),RY(I),RZ(I)
IF I <> N THEN
InitialPlotString$ = InitialPlotString$ & "(" & STR$(RX(I)) & "," & STR$(RY(I)) & "," & STR$(RZ(I)) & "),"
ELSE
InitialPlotString$ = InitialPlotString$ & "(" & STR$(RX(I)) & "," & STR$(RY(I)) & "," & STR$(RZ(I)) & ")],color=red,size=20,opacity=1)"
END IF
NEXT I
CALL LennardJones
PRINT "Initial V = "; V
PRINT #1: "Initial V = "; V
OLDV = V
DO WHILE (NOMOVE < 3*N)
FOR I = 1 TO N
! Wiggle atom I in X and see if it reduces the potential.
RX(I) = RX(I) + STEPX
! PRINT "STEP ";STEP
STEP = STEP + 1
CALL LennardJones
IF ABS(V) > ABS(OLDV) THEN
RX(I) = RX(I) - 2*STEPX
! PRINT "STEP ";STEP
STEP = STEP + 1
CALL LennardJones
IF ABS(V) > ABS(OLDV) THEN
NOMOVE = NOMOVE + 1
RX(I) = RX(I) + STEPX
ELSE
OLDV = V
END IF
ELSE
OLDV = V
END IF
! Wiggle atom I in Y and see if it reduces the potential.
RY(I) = RY(I) + STEPY
! PRINT "STEP ";STEP
STEP = STEP + 1
CALL LennardJones
IF ABS(V) > ABS(OLDV) THEN
RY(I) = RY(I) - 2*STEPY
! PRINT "STEP ";STEP
STEP = STEP + 1
CALL LennardJones
IF ABS(V) > ABS(OLDV) THEN
NOMOVE = NOMOVE + 1
RY(I) = RY(I) + STEPY
ELSE
OLDV = V
END IF
ELSE
OLDV = V
END IF
! Wiggle atom I in Z and see if it reduces the potential.
RZ(I) = RZ(I) + STEPZ
!PRINT "STEP ";STEP
STEP = STEP + 1
CALL LennardJones
IF ABS(V) > ABS(OLDV) THEN
RZ(I) = RZ(I) - 2*STEPZ
!PRINT "STEP ";STEP
STEP = STEP + 1
CALL LennardJones
IF ABS(V) > ABS(OLDV) THEN
NOMOVE = NOMOVE + 1
RZ(I) = RZ(I) + STEPZ
ELSE
OLDV = V
END IF
ELSE
OLDV = V
END IF
NEXT I
PRINT "STEP: ";STEP
LOOP
PRINT "Minimum V: ";OLDV
PRINT "Atomic Coordinates:"
PRINT "X","Y","Z"
PRINT #1: "Minimum V: ";OLDV
PRINT #1: "Atomic Coordinates:"
PRINT #1: "X","Y","Z"
FOR I = 1 TO N
PRINT RX(I),RY(I),RZ(I)
PRINT #1: RX(I),RY(I),RZ(I)
IF I <> N THEN
FinalPlotString$ = FinalPlotString$ & "(" & STR$(RX(I)) & "," & STR$(RY(I)) & "," & STR$(RZ(I)) & "),"
ELSE
FinalPlotString$ = FinalPlotString$ & "(" & STR$(RX(I)) & "," & STR$(RY(I)) & "," & STR$(RZ(I)) & ")], color=blue,size=20,opacity=1)"
END IF
NEXT I
PRINT #1: "Atomic Positions for Sage plotting:"
PRINT #1: InitialPlotString$
PRINT #1: FinalPlotString$
PRINT #1: "Show(IP+FP)"
PRINT "Done."
SUB LennardJones
! This subroutine calculates the Lennard-Jones potential for the system
LOCAL I
V = 0.0
FOR I = 1 TO N - 1
RXI = RX(I)
RYI = RY(I)
RZI = RZ(I)
FOR J = I + 1 TO N
RXIJ = RXI - RX(J)
RYIJ = RYI - RY(J)
RZIJ = RZI - RZ(J)
RIJSQ = RXIJ^2 + RYIJ^2 + RZIJ^2
SR2 = SIGSQ / RIJSQ
SR6 = SR2 * SR2 * SR2
SR12 = SR6^2
V = V + SR12 - SR6
NEXT J
NEXT I
V = 4.0 * EPSILON * V
PRINT #1: "STEP "; STEP; ":"
FOR I = 1 TO N
PRINT #1: RX(I),RY(I),RZ(I)
NEXT I
!PRINT "V = "; V
PRINT #1:"V =";V
END SUB
END
Tuesday, September 26, 2017
Friday, May 12, 2017
Learning about CBM-BASIC
I've been solving Brilliant math problems with the simplest of BASIC programming languages lately. I do this because it's super convenient for me to write a quick little program to do some equation solving. CBM-BASIC is pretty nice for this because it works with Windows and Linux and my iPhone (by installing Hand-BASIC).
For a nice little example I solved this problem:
With the help of Mathematica and CBM-BASIC. I used Mathematica to solve for z and plot the resulting equation to help me understand if there were a finite number of solutions. I was being lazy, because I could have figured this out just by plugging a few exploratory numbers in for X,Y, and Z to help me understand the bounds of the problem. Once I understood the bounds of the problem I was able to find all the solutions with CBM-BASIC in a relatively short time period.
Here is the resulting BASIC program:
And there you have it. I'll let the reader build CBM-BASIC and find the solution themselves. I'm having a lot of fun using this little Commordore 64 BASIC 2 implementation for my mental exercise. Who said BASIC causes brain damage?
For a nice little example I solved this problem:
With the help of Mathematica and CBM-BASIC. I used Mathematica to solve for z and plot the resulting equation to help me understand if there were a finite number of solutions. I was being lazy, because I could have figured this out just by plugging a few exploratory numbers in for X,Y, and Z to help me understand the bounds of the problem. Once I understood the bounds of the problem I was able to find all the solutions with CBM-BASIC in a relatively short time period.
Here is the resulting BASIC program:
10 REM Solve Brilliant Problem 2 of May 8, 2017
20 S = 0
30 E1 = 0
40 MX = 150
50 MY = 50
60 MZ = 40
65 PRINT "FINDING SOLUTIONS..."
70 FOR X = 4 TO MX
80 FOR Y = 1 TO MY
90 FOR Z = 1 TO MZ
100 GOSUB 500
110 IF ABS(E1) < 0.1 THEN S=S+X+Y+Z : PRINT "S = ";S
120 NEXT Z : NEXT Y : NEXT X
130 PRINT "DONE." : END
500 REM EVALUATE FUNCTION E1
510 E1 = X^2*(Y^3+Z^3)-315*(X*Y*Z+7)
520 RETURNAnd there you have it. I'll let the reader build CBM-BASIC and find the solution themselves. I'm having a lot of fun using this little Commordore 64 BASIC 2 implementation for my mental exercise. Who said BASIC causes brain damage?
Friday, April 21, 2017
Unicon benchmark results on Raspberry Pi Zero W
This was from a version of Unicon compiled from unicon-code-5082-trunk. A recent development snapshot.
Times reported below reflect averages over three executions.
Expect 2-20 minutes for suite to run to completion.
Word Size Main Memory C Compiler clock OS
32 bit 424 MB gcc 4.9.2 1.0 GHz UNIX
CPU
1x ARMv6-compatible processor rev 7 (v6l)
Elapsed time h:m:s |Concurrent |
benchmark |Sequential| |Concurrent| |Performance|
concord concord.dat 01:0.685 N/A
deal 50000 01:0.504 N/A
ipxref ipxref.dat 35.286 N/A
queens 12 01:13.523 N/A
rsg rsg.dat 01:3.325 N/A
binary-trees 14 01:36.725 04:55.517 0.327x
chameneos-redux 65000 N/A 01:8.534
fannkuch 9 01:8.319 N/A
fasta 250000 01:22.580 N/A
k-nucleotide 150-thou.dat 03:19.441 N/A
mandelbrot 750 08:10.612 13:54.547 0.587x
meteor-contest 600 03:30.708 N/A
n-body 100000 03:4.327 N/A
pidigits 7000 14:37.462 N/A
regex-dna 700-thou.dat 01:52.136 03:33.752 0.524x
reverse-complement 15-mil.dat
Run-time error 306
File reverse-complement.icn; Line 35
inadequate space in string region
Traceback:
main(list_1 = [])
getavgtimes(procedure run_reversecomplement,"15-mil.dat") from line 186 in run-benchmark.icn
run_reversecomplement(list_5068944 = ["15-mil.dat"]) from line 129 in auxiliary.icn
mapseq(">ONE Homo sapien...") from line 46 in reverse-complement.icn
"
Run-time error 302
File reverse-complement.icn; Line 35
memory violation
Traceback:
main(list_1 = [])
getavgtimes(procedure run_reversecomplement,"15-mil.dat") from line 186 in run-benchmark.icn
run_reversecomplement(list_5068944 = ["15-mil.dat"]) from line 129 in auxiliary.icn
mapseq(">ONE Homo sapien...") from line 46 in reverse-complement.icn
Times reported below reflect averages over three executions.
Expect 2-20 minutes for suite to run to completion.
Word Size Main Memory C Compiler clock OS
32 bit 424 MB gcc 4.9.2 1.0 GHz UNIX
CPU
1x ARMv6-compatible processor rev 7 (v6l)
Elapsed time h:m:s |Concurrent |
benchmark |Sequential| |Concurrent| |Performance|
concord concord.dat 01:0.685 N/A
deal 50000 01:0.504 N/A
ipxref ipxref.dat 35.286 N/A
queens 12 01:13.523 N/A
rsg rsg.dat 01:3.325 N/A
binary-trees 14 01:36.725 04:55.517 0.327x
chameneos-redux 65000 N/A 01:8.534
fannkuch 9 01:8.319 N/A
fasta 250000 01:22.580 N/A
k-nucleotide 150-thou.dat 03:19.441 N/A
mandelbrot 750 08:10.612 13:54.547 0.587x
meteor-contest 600 03:30.708 N/A
n-body 100000 03:4.327 N/A
pidigits 7000 14:37.462 N/A
regex-dna 700-thou.dat 01:52.136 03:33.752 0.524x
reverse-complement 15-mil.dat
Run-time error 306
File reverse-complement.icn; Line 35
inadequate space in string region
Traceback:
main(list_1 = [])
getavgtimes(procedure run_reversecomplement,"15-mil.dat") from line 186 in run-benchmark.icn
run_reversecomplement(list_5068944 = ["15-mil.dat"]) from line 129 in auxiliary.icn
mapseq(">ONE Homo sapien...") from line 46 in reverse-complement.icn
"
Run-time error 302
File reverse-complement.icn; Line 35
memory violation
Traceback:
main(list_1 = [])
getavgtimes(procedure run_reversecomplement,"15-mil.dat") from line 186 in run-benchmark.icn
run_reversecomplement(list_5068944 = ["15-mil.dat"]) from line 129 in auxiliary.icn
mapseq(">ONE Homo sapien...") from line 46 in reverse-complement.icn
Subscribe to:
Posts (Atom)