periodic boundary condition

ffs.toml

@@ -3,14 +3,14 @@
 dt = 0.001
 # maximum time
 t = 100.0
-# gravitational acceleration
-g = 0.000001
 # lennard jones parameters
-lennardE = 15
-lennardS = 0.5
+lennardE = 100
+lennardS = 1.0
+# neighbor distance, in lennardS
+neighborCutoff = 3.0
 # boundary details
-boundary = 30.0
+boundary = 10.0
 # save to disk every n time steps
-saveInterval = 500
+saveInterval = 200
 # initial velocity randomizer - adds extraVel at the start of calculation in a random direction to every particle
 extraVel = 1.0

forces.nim

@@ -4,6 +4,10 @@ import linalg
 import math
 import strformat
 import random
+import parsexyz
+
+import threadpool
+const nThreads = 8
 
 type
     ForceField = ref object
@@ -11,7 +15,7 @@ type
         # should be split in two - the actual ff and calculation parameters
         lennardE*: float
         lennardS*: float
-        gravAcc*: float
+        neighborCutoff*: float
         boundary*: float
         dt*: float
         maxt*: float
@@ -31,12 +35,12 @@ proc newForceField*(path: string, name: string): ForceField =
     let config = input[name]
     result.dt = config["dt"].getFloat()
     result.maxt = config["t"].getFloat()
-    result.gravAcc = config["g"].getFloat()
     result.lennardE = config["lennardE"].getFloat()
     result.lennardS = config["lennardS"].getFloat()
     result.boundary = config["boundary"].getFloat()
     result.saveInterval = config["saveInterval"].getInt()
     result.extraVel = config["extraVel"].getFloat()
+    result.neighborCutoff = config["neighborCutoff"].getFloat()
 
 proc applyExtraVel*(particles: var seq[Particle], ff: ForceField) =
     for i in 0..particles.high():
@@ -46,44 +50,66 @@ proc getAcc(particles: seq[Particle], ff: ForceField, index: int): Vector =
     result = vector0()
     let posThis = particles[index].pos
     let mThis = particles[index].mass
+
+    let s6 = ff.lennardS ^ 6
+    let lennard_pre = 24f * ff.lennardE * s6
+    let boundary = ff.boundary
+
     for i in 0..particles.high:
         if i == index:
             continue
-        let diff = particles[i].pos - posThis
+        var diff = particles[i].pos - posThis
+        # periodic condition
+        if diff.x > boundary:
+            diff.x -= 2f * boundary
+        elif diff.x < -boundary:
+            diff.x += 2f * boundary
+        if diff.y > boundary:
+            diff.y -= 2f * boundary
+        elif diff.y < -boundary:
+            diff.y += 2f * boundary
+        if diff.z > boundary:
+            diff.z -= 2f * boundary
+        elif diff.z < -boundary:
+            diff.z += 2f * boundary
+
         if diff.len() == 0.0:
             raise newException(Defect, &"overlap between {i} and {index}")
-        let direction = diff.normalize()
         let distance = diff.len()
 
+        if distance > ff.neighborCutoff * ff.lennardS:
+            continue
+
+        let direction = diff / distance
+
         # Lennard-Jones
         # F = 4e (6s^6/r^7 - 12s^12/r^13)
-        result -= direction * (4f * ff.lennardE * (6.0 * (ff.lennardS^6)/(distance^7) - 12.0 * (ff.lennardS^12)/(distance^13))) / mThis
+        let rm1 = 1/distance
+        let r6 = rm1 ^ 6
+        result -= direction * (lennard_pre * rm1 * (r6) * (1 - 2f * s6 * r6) / mThis)
         
-        # Gravity downwards
-        result += vector(0, 0, ff.gravAcc)
-
 proc keepInside(particles: var seq[Particle], ff: ForceField, index: int) =
-    # TODO: thermostatic walls and periodic boundary condition
+    # periodic boundary condition
     let posThis = particles[index].pos
     let velThis = particles[index].vel
 
     if posThis.x > ff.boundary and velThis.x > 0:
-        particles[index].vel.x = -velThis.x
+        particles[index].pos.x = -ff.boundary
 
     if posThis.x < -ff.boundary and velThis.x < 0:
-        particles[index].vel.x = -velThis.x
+        particles[index].pos.x = ff.boundary
 
     if posThis.y > ff.boundary and velThis.y > 0:
-        particles[index].vel.y = -velThis.y
+        particles[index].pos.y = -ff.boundary
 
     if posThis.y < -ff.boundary and velThis.y < 0:
-        particles[index].vel.y = -velThis.y
+        particles[index].pos.y = ff.boundary
     
     if posThis.z > ff.boundary and velThis.z > 0:
-        particles[index].vel.z = -velThis.z
+        particles[index].pos.z = -ff.boundary
     
     if posThis.z < -ff.boundary and velThis.z < 0:
-        particles[index].vel.z = -velThis.z
+        particles[index].pos.z = ff.boundary
 
 proc iterate*(particles: var seq[Particle], ff: ForceField, t: float) =
     let dt = ff.dt
@@ -95,9 +121,15 @@ proc iterate*(particles: var seq[Particle], ff: ForceField, t: float) =
             particles[i].pos = part.pos + part.vel * dt + dt * dt * 0.5 * part.acc
 
         # update acceleration and velocity
+        var newAccs = newSeq[Vector](particles.len)
+        #for i in 0..nThreads:
+            # spawn a thread
+        for i in 0..particles.high:
+            newAccs[i] = getAcc(particles, ff, i)
+
         for i in 0..particles.high:
             let part = particles[i]
-            let newAcc = getAcc(particles, ff, i)
+            let newAcc = newAccs[i]
             let newVel = part.vel + (part.acc + newAcc) * (dt * 0.5)
             particles[i].vel = newVel
             particles[i].acc = newAcc
@@ -106,3 +138,26 @@ proc iterate*(particles: var seq[Particle], ff: ForceField, t: float) =
         for i in 0..particles.high:
             keepInside(particles, ff, i)
         
+
+proc simulate*(ffName: string, xyzPath: string) =
+    let ff = newForceField("ffs.toml", ffName)
+
+    var particles: seq[Particle] = parseXyz(xyzPath)
+
+    applyExtraVel(particles, ff) # TODO replace this
+
+    # open the output file
+    let f = open(xyzPath, fmAppend)
+    var saveCounter = 0
+
+    var t = 0.0
+    while t < ff.maxt:
+        iterate(particles, ff, t)
+        inc saveCounter
+        if saveCounter >= ff.saveInterval:
+            saveCounter = 0
+            appendXyz(f, particles, t)
+        t += ff.dt
+
+    # cleanup
+    f.close()

genrandom.nim

@@ -46,6 +46,7 @@ for i in 0..atomCount:
             particles.add(newParticle(name, pos, 1.0))
             break thisPart
 
+echo &"{particles.len} particles generated."
 
 # and write them to xyz
 var f = open(filename, fmWrite)

main.nim

@@ -16,24 +16,4 @@ let xyzPath = paramStr(1)
 let ffName = paramStr(2)
 echo &"Using {xyzPath} as input, {ffName} as calculation settings"
 
-let ff = newForceField("ffs.toml", ffName)
-
-var particles: seq[Particle] = parseXyz(xyzPath)
-
-applyExtraVel(particles, ff) # TODO replace this
-
-# open the output file
-let f = open(xyzPath, fmAppend)
-var saveCounter = 0
-
-var t = 0.0
-while t < ff.maxt:
-    iterate(particles, ff, t)
-    inc saveCounter
-    if saveCounter >= ff.saveInterval:
-        saveCounter = 0
-        appendXyz(f, particles, t)
-    t += ff.dt
-
-# cleanup
-f.close()
+simulate(ffName, xyzPath)

rebuild.sh

@@ -0,0 +1,4 @@
+rm main
+rm genrandom
+nim c --threads:on --gc:arc -d:release main
+nim c --threads:on --d:release genrandom