GlowScript 2.1 VPython

num = 50
atoms = []
stiffness = 1000     # N/m - stiffness of elastic force
speed = 2    # initial speeds in m/s
freq = 200 
kb = 1380.649

scale = 0.5   # size scale of system in meters
atomR = scale / 20     # radius of atoms

for i in arange(num):
  atoms[i] = sphere( radius = atomR)
  atoms[i].pos = 0.5 * scale * vector.random()
  atoms[i].mass = 0.1    # in kg
  atoms[i].vel = speed * vector.random()
  atoms[i].ke = 0.5 * atoms[i].mass * atoms[i].vel.mag * atoms[i].vel.mag

dt = 1/freq
time = 0

side = scale
thk = 2 * atomR
s2 = 2*side + thk
s3 = 2*side + thk

wallR = box (pos=vector( scale + thk, 0, 0), size=vector(thk, s2, s3),  color = color.red)
wallL = box (pos=vector(-scale - thk, 0, 0), size=vector(thk, s2, s3),  color = color.red)
wallB = box (pos=vector(0, -scale-thk, 0), size=vector(s3, thk, s3),  color = color.blue)
wallT = box (pos=vector(0,  scale+thk, 0), size=vector(s3, thk, s3),  color = color.blue)
wallBK = box(pos=vector(0, 0, -scale-thk), size=vector(s2, s2, thk), color = vector(0.7,0.7,0.7))

# choose some random atoms to color for visibility 
for i in arange(num/20):
  atoms[i].color = vec(0,1,1)

kinetic_energy = 0
potential_energy = 0
rms_speed = 0

while time < 50:  
  rate(freq)
  for i in arange(num):
    
    f_normal = vec(0, 0, 0)
    
    # check for collisions w/ walls and apply elastic force
    if (atoms[i].pos.x < -scale):
      f_normal = f_normal + vec(stiffness * (-scale - atoms[i].pos.x) , 0, 0)
    if (atoms[i].pos.x > scale):
      f_normal = f_normal + vec(stiffness * (scale - atoms[i].pos.x), 0, 0)

    if (atoms[i].pos.y < -scale):
      f_normal = f_normal + vec(0, stiffness * (-scale - atoms[i].pos.y) , 0)
    if (atoms[i].pos.y > scale):
      f_normal = f_normal + vec(0, stiffness * (scale - atoms[i].pos.y), 0)

    if (atoms[i].pos.z < -scale):
      f_normal = f_normal + vec(0, 0, stiffness * (-scale - atoms[i].pos.z))
    if (atoms[i].pos.z > scale):
      f_normal = f_normal + vec(0, 0, stiffness * (scale - atoms[i].pos.z))

    # check for collisions with other atoms 
    for j in arange(num):
      if (j == i): 
        continue
      sep = atoms[i].pos - atoms[j].pos 
      if (sep.mag < 2.0 * atomR):   # two atoms overlapping
        
        f_normal = f_normal - stiffness * (sep - 2.0 * atomR * sep.hat)

    # calculate change in motion
    atoms[i].vel = atoms[i].vel + f_normal * dt / atoms[i].mass
    
    kinetic_energy += 0.5 * atoms[i].mass * atoms[i].vel.mag * atoms[i].vel.mag
    #potential_energy += stiffness * (scale - atom)
    rms_speed += (atoms[i].vel.mag * atoms[i].vel.mag) ** 0.5
  
  for i in arange(num):
    # update positions
    atoms[i].pos = atoms[i].pos + atoms[i].vel * dt
  
  scene.caption = "Total Kinetic Energy: {0:.3f}, Average Kinetic energy: {1:.3f}, Average rms speed: {2:.3f}, Average Temperature: {3:.3e}".format(kinetic_energy, kinetic_energy/num, rms_speed/num, 2*kinetic_energy/(num*kb*3))
  kinetic_energy = 0
  rms_speed = 0
  time = time + dt