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
