Compare commits
2 Commits
16a4a8f2d9
...
ddbe08ad05
Author | SHA1 | Date |
---|---|---|
Art | ddbe08ad05 | |
Art | 8bc009246d |
12
ffs.toml
12
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
|
87
forces.nim
87
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
|
||||
|
||||
# Gravity downwards
|
||||
result += vector(0, 0, ff.gravAcc)
|
||||
let rm1 = 1/distance
|
||||
let r6 = rm1 ^ 6
|
||||
result -= direction * (lennard_pre * rm1 * (r6) * (1 - 2f * s6 * r6) / mThis)
|
||||
|
||||
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()
|
|
@ -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)
|
||||
|
|
22
main.nim
22
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)
|
|
@ -0,0 +1,4 @@
|
|||
rm main
|
||||
rm genrandom
|
||||
nim c --threads:on --gc:arc -d:release main
|
||||
nim c --threads:on --d:release genrandom
|
Loading…
Reference in New Issue