From ddbe08ad050c8c687d15b86871f2b3b4aeb61750 Mon Sep 17 00:00:00 2001 From: prod2 Date: Sat, 15 Apr 2023 15:59:24 +0200 Subject: [PATCH] periodic boundary condition --- forces.nim | 87 +++++++++++++++++++++++++++++++++++++++++---------- genrandom.nim | 1 + main.nim | 22 +------------ rebuild.sh | 4 +++ 4 files changed, 77 insertions(+), 37 deletions(-) create mode 100755 rebuild.sh diff --git a/forces.nim b/forces.nim index 1108100..092e782 100644 --- a/forces.nim +++ b/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() \ No newline at end of file diff --git a/genrandom.nim b/genrandom.nim index f99ec55..2f4f8a2 100644 --- a/genrandom.nim +++ b/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) diff --git a/main.nim b/main.nim index f0c5f57..0fa6feb 100644 --- a/main.nim +++ b/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() \ No newline at end of file +simulate(ffName, xyzPath) \ No newline at end of file diff --git a/rebuild.sh b/rebuild.sh new file mode 100755 index 0000000..961ce59 --- /dev/null +++ b/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 \ No newline at end of file