diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..2ef716b --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +main +genrandom +*.xyz \ No newline at end of file diff --git a/FreeMono.ttf b/FreeMono.ttf deleted file mode 100644 index 96b830e..0000000 Binary files a/FreeMono.ttf and /dev/null differ diff --git a/ffs.toml b/ffs.toml new file mode 100644 index 0000000..3bd7835 --- /dev/null +++ b/ffs.toml @@ -0,0 +1,16 @@ +[test] +# time step +dt = 0.001 +# maximum time +t = 100.0 +# gravitational acceleration +g = 0.000001 +# lennard jones parameters +lennardE = 15 +lennardS = 0.5 +# boundary details +boundary = 30.0 +# save to disk every n time steps +saveInterval = 500 +# initial velocity randomizer - adds extraVel at the start of calculation in a random direction to every particle +extraVel = 1.0 \ No newline at end of file diff --git a/forces.nim b/forces.nim new file mode 100644 index 0000000..1108100 --- /dev/null +++ b/forces.nim @@ -0,0 +1,108 @@ +import parsetoml +import particle +import linalg +import math +import strformat +import random + +type + ForceField = ref object + # this is stretching the definition of an ff beyond what's reasonable + # should be split in two - the actual ff and calculation parameters + lennardE*: float + lennardS*: float + gravAcc*: float + boundary*: float + dt*: float + maxt*: float + saveInterval*: int + extraVel*: float # TODO remove this once diff velocity is implemented + # TODO: new options - thermostat and periodic boundaries + # cutoff distance for calculations + # load velocity as diff from last two frames options and adapt genrandom for that option + +proc newForceField*(path: string, name: string): ForceField = + new(result) + + # read config file + let input = parsetoml.parseString(path.readFile()) + + # initialize variables based on input file + 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() + +proc applyExtraVel*(particles: var seq[Particle], ff: ForceField) = + for i in 0..particles.high(): + particles[i].vel = vector(rand(ff.extraVel * 2) - ff.extraVel, rand(ff.extraVel * 2) - ff.extraVel, rand(ff.extraVel * 2) - ff.extraVel) + +proc getAcc(particles: seq[Particle], ff: ForceField, index: int): Vector = + result = vector0() + let posThis = particles[index].pos + let mThis = particles[index].mass + for i in 0..particles.high: + if i == index: + continue + let diff = particles[i].pos - posThis + if diff.len() == 0.0: + raise newException(Defect, &"overlap between {i} and {index}") + let direction = diff.normalize() + let distance = diff.len() + + # 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) + +proc keepInside(particles: var seq[Particle], ff: ForceField, index: int) = + # TODO: thermostatic walls and 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 + + if posThis.x < -ff.boundary and velThis.x < 0: + particles[index].vel.x = -velThis.x + + if posThis.y > ff.boundary and velThis.y > 0: + particles[index].vel.y = -velThis.y + + if posThis.y < -ff.boundary and velThis.y < 0: + particles[index].vel.y = -velThis.y + + if posThis.z > ff.boundary and velThis.z > 0: + particles[index].vel.z = -velThis.z + + if posThis.z < -ff.boundary and velThis.z < 0: + particles[index].vel.z = -velThis.z + +proc iterate*(particles: var seq[Particle], ff: ForceField, t: float) = + let dt = ff.dt + block step: + # iterate through particles + # update positions + for i in 0..particles.high: + let part = particles[i] + particles[i].pos = part.pos + part.vel * dt + dt * dt * 0.5 * part.acc + + # update acceleration and velocity + for i in 0..particles.high: + let part = particles[i] + let newAcc = getAcc(particles, ff, i) + let newVel = part.vel + (part.acc + newAcc) * (dt * 0.5) + particles[i].vel = newVel + particles[i].acc = newAcc + + # enforce boundary + for i in 0..particles.high: + keepInside(particles, ff, i) + diff --git a/genrandom.nim b/genrandom.nim new file mode 100644 index 0000000..f99ec55 --- /dev/null +++ b/genrandom.nim @@ -0,0 +1,53 @@ +# generate input xyz based on some parameters + +import random +import os +import strutils +import strformat +import particle +import linalg +import parsexyz +import forces + +if paramCount() != 3: + echo "usage ./genrandom n filename.xyz settingsName" + quit 1 + +let atomCount = parseInt(paramStr(1)) +let filename = paramStr(2) +let settingsName = paramStr(3) + +let ff = newForceField("ffs.toml", settingsName) +let boundary = ff.boundary +let lennardS = ff.lennardS + +var particles: seq[Particle] + +# get particles +for i in 0..atomCount: + block thisPart: + let name = &"c{i}" + var attempt = 0 + while true: + # too many attempts: stop trying to spawn this particle + inc attempt + if attempt > 10: + break thisPart + + # random position + let pos = vector(rand(2 * boundary)-boundary, rand(2*boundary)-boundary, rand(2*boundary) - boundary) + + # check for overlap, go back to start of while if there is one + for p in particles: + if (p.pos - pos).len() < lennardS * 2.5: + break + + # this reached? particle successfully spawned + particles.add(newParticle(name, pos, 1.0)) + break thisPart + + +# and write them to xyz +var f = open(filename, fmWrite) +appendXyz(f, particles, 0.0) +f.close() \ No newline at end of file diff --git a/input.toml b/input.toml deleted file mode 100644 index 9cb7742..0000000 --- a/input.toml +++ /dev/null @@ -1,35 +0,0 @@ -[config] -# time step -dt = 0.001 -# maximum time -t = 1000.0 -# gravitational constant -G = 0.0 -# gravitational acceleration -g = 0.000001 -# stabilize centre of mass? -stab_mass = false -# how often reset centre of mass -stab_interval = 100 -# lennard jones parameters -lennardE = 15 -lennardS = 0.5 -# boundary details -boundary = 30.0 -# save to disk every n time steps -saveInterval = 200 - -# format: after an initial config block, -# every block is one object, with the name being in square brackets -# mandatory fields: - # pos (vector) - initial position - # vel (vector) - initial velocity - # mass - its mass - -[random] -count = 80 -xmin = -29.0 -xmax = 29.0 -ymin = -29.0 -ymax = 29.0 -vmax = 0.0 diff --git a/linalg.nim b/linalg.nim index fd82ea1..0e0235f 100644 --- a/linalg.nim +++ b/linalg.nim @@ -5,27 +5,31 @@ type Vector* = object x*: float y*: float + z*: float -func vector*(x, y: float): Vector = - Vector(x: x, y: y) +func vector*(x, y, z: float): Vector = + Vector(x: x, y: y, z: z) + +func vector0*: Vector = + vector(0.0, 0.0, 0.0) func `+`*(a, b: Vector): Vector = - Vector(x: a.x + b.x, y: a.y + b.y) + Vector(x: a.x + b.x, y: a.y + b.y, z: a.z + b.z) func `-`*(a: Vector): Vector = - Vector(x: -a.x, y: -a.y) + Vector(x: -a.x, y: -a.y, z: -a.z) func `-`*(a, b: Vector): Vector = a + (-b) func `*`*(a: Vector, f: float): Vector = - Vector(x: a.x * f, y: a.y * f) + Vector(x: a.x * f, y: a.y * f, z: a.z * f) func `*`*(f: float, a: Vector): Vector = a * f func dot*(a, b: Vector): float = - a.x * b.x + a.y * b.y + a.x * b.x + a.y * b.y + a.z * b.z func len*(a: Vector): float = sqrt(dot(a, a)) @@ -43,4 +47,4 @@ func `-=`*(a: var Vector, b: Vector) = a = a - b func `$`*(a: Vector): string = - &"[{a.x}, {a.y}]" + &"{a.x} {a.y} {a.z}" diff --git a/main.nim b/main.nim index 9900d79..f0c5f57 100644 --- a/main.nim +++ b/main.nim @@ -1,178 +1,39 @@ -import linalg -import parsetoml +import particle +import parsexyz +import forces -import math import strformat -import random +import os -# define types -type - Particle = object - name: string - pos: Vector - vel: Vector - acc: Vector - mass: float +# what's the input file name? +if paramCount() != 2: + echo "Usage: ./main inputfile.xyz settingName" + quit 1 +if not fileExists("ffs.toml"): + echo "Error: input.toml doesn't exist" + quit 2 +let xyzPath = paramStr(1) +let ffName = paramStr(2) +echo &"Using {xyzPath} as input, {ffName} as calculation settings" -# read input file -let input = parsetoml.parseString("input.toml".readFile()) +let ff = newForceField("ffs.toml", ffName) -# initialize variables based on input file -let config = input["config"] -var t = 0.0 # current time -let dt = config["dt"].getFloat() -let maxt = config["t"].getFloat() -let gravConst = config["G"].getFloat() -let gravAcc = config["g"].getFloat() -let stabilizeMassCentre = config["stab_mass"].getBool() -let stabInterval = config["stab_interval"].getInt() -var stabCounter = 0 -let lennardE = config["lennardE"].getFloat() -let lennardS = config["lennardS"].getFloat() -let boundary = config["boundary"].getFloat() -let saveInterval = config["saveInterval"].getInt() -# a thing that runs every 100 frames to move the centre of mass back to the original one -var particles: seq[Particle] = @[] +var particles: seq[Particle] = parseXyz(xyzPath) -# process the rest of the input file -for key, val in input.tableVal.pairs: - if key == "config": - continue - elif key == "random": - let count = val["count"].getInt() - let xmin = val["xmin"].getFloat() - let xmax = val["xmax"].getFloat() - let ymin = val["ymin"].getFloat() - let ymax = val["ymax"].getFloat() - let vmax = val["vmax"].getFloat() - for i in 0..count: - block thisPart: - var particle = Particle(name: "random") - var pass = false - var attempt = 0 - while not pass: - inc attempt - if attempt > 10: - break thisPart - particle.pos = vector(rand(xmax-xmin)+xmin, rand(ymax-ymin)+ymin) - pass = true - for p in particles: - if (p.pos - particle.pos).len() < lennardS * 2.5: - pass = false - break - particle.vel = vector(rand(vmax), rand(vmax)) - particle.acc = vector(0.0, 0.0) - particle.mass = 1.0 - particles.add(particle) - else: - let x: float = val["pos"][0].getFloat() - let y: float = val["pos"][1].getFloat() - let vx: float = val["vel"][0].getFloat() - let vy: float = val["vel"][1].getFloat() - let pos = vector(x, y) - let vel = vector(vx, vy) - let mass = val["mass"].getFloat() - let particle = Particle(name: key, pos: pos, vel: vel, acc: vector(0,0), mass: mass) - particles.add(particle) +applyExtraVel(particles, ff) # TODO replace this # open the output file -let f = open("output.txt", fmWrite) -# do the simulation, writing results to a file (streaming) -f.writeLine("mdsim input file: input.toml") -proc oupStatus = - # output current status - var status: string = "" - for i in 0..particles.high: - status &= &"{particles[i].name} {$particles[i].pos} {$particles[i].vel} " - f.writeLine(&"{t} {status}") - -oupStatus() - -proc getForces(index: int): Vector = - result = vector(0, 0) - let posThis = particles[index].pos - let mThis = particles[index].mass - for i in 0..particles.high: - if i == index: - continue - let diff = particles[i].pos - posThis - let direction = diff.normalize() - let distance = diff.len() - let mThat = particles[i].mass - - # gravity between particles - # m2 * C / r^2 (so it's acceleartion, mass of this doesn't apply) - result += direction * (mThat * gravConst / distance / distance) - - # Lennard-Jones - # F = 4e (6s^6/r^7 - 12s^12/r^13) - result -= direction * (4f * lennardE * (6.0 * (lennardS^6)/(distance^7) - 12.0 * (lennardS^12)/(distance^13))) / mThis - - result += vector(0, gravAcc) - -proc keepInside(index: int) = - let posThis = particles[index].pos - let velThis = particles[index].vel - - if posThis.x > boundary and velThis.x > 0: - particles[index].vel.x = -velThis.x - - if posThis.x < -boundary and velThis.x < 0: - particles[index].vel.x = -velThis.x - - if posThis.y > boundary and velThis.y > 0: - particles[index].vel.y = -velThis.y - - if posThis.y < -boundary and velThis.y < 0: - particles[index].vel.y = -velThis.y - - -proc getCentreOfMass(): Vector = - var totMass = 0.0 - result = vector(0.0, 0.0) - for i in 0..particles.high: - let part = particles[i] - totMass += part.mass - result += part.pos * part.mass - result = result / totMass - -let initCentreOfMass = getCentreOfMass() - -proc resetCentreOfMass() = - let cCentre = getCentreOfMass() - let delta = cCentre - initCentreOfMass - for i in 0..particles.high: - particles[i].pos -= delta - +let f = open(xyzPath, fmAppend) var saveCounter = 0 -while t < maxt: - block step: - # iterate through particles - for i in 0..particles.high: - let part = particles[i] - particles[i].pos = part.pos + part.vel * dt + dt * dt * 0.5 * part.acc - - for i in 0..particles.high: - let part = particles[i] - let newAcc = getForces(i) - let newVel = part.vel + (part.acc + newAcc) * (dt * 0.5) - particles[i].vel = newVel - particles[i].acc = newAcc - - # enforce boundary - for i in 0..particles.high: - keepInside(i) +var t = 0.0 +while t < ff.maxt: + iterate(particles, ff, t) inc saveCounter - if saveCounter >= saveInterval: + if saveCounter >= ff.saveInterval: saveCounter = 0 - oupStatus() - t += dt - inc stabCounter - if stabCounter > stabInterval: - stabCounter = 0 - if stabilizeMassCentre: - resetCentreOfMass() + appendXyz(f, particles, t) + t += ff.dt # cleanup f.close() \ No newline at end of file diff --git a/parsexyz.nim b/parsexyz.nim new file mode 100644 index 0000000..dd08b99 --- /dev/null +++ b/parsexyz.nim @@ -0,0 +1,39 @@ +# Reader/writer for xyz files +import particle +import strutils +import strscans +import strformat +import linalg + +proc parseXyz*(path: string): seq[Particle] = + let f = readFile(path) + let lines = f.split("\n") + var i = 0 + while i <= lines.high(): + # read a single frame + # read atom count + let atomCount = parseInt(lines[i]) + inc i + + # the comment + let comment = lines[i] + inc i + + var x, y, z: float + var name: string + # the atoms + if i + atomCount + 1 >= lines.high(): # tolerate empty line at end + for j in 0..atomCount-1: + if not lines[i + j].scanf("$s$w$s$f$s$f$s$f$s", name, x, y, z): + raise newException(Defect, "Invalid format on line " & $i) + result.add(newParticle(name, vector(x, y, z), massFromName(name))) + break + i += atomCount + +proc appendXyz*(f: File, parts: seq[Particle], t: float) = + # f should be open in APPEND mode + f.writeLine(parts.len()) # atom count + f.writeLine("t = " & $t) # comment line containing the time + for atom in parts: + let pos = atom.pos + f.writeLine(&"{atom.name} {pos.x} {pos.y} {pos.z}") diff --git a/particle.nim b/particle.nim new file mode 100644 index 0000000..e747bde --- /dev/null +++ b/particle.nim @@ -0,0 +1,20 @@ +import linalg + +# define types +type + Particle* = object + name*: string + pos*: Vector + vel*: Vector + acc*: Vector + mass*: float + +func `$`*(part: Particle): string = + # Conversion to xyz + part.name & " " & $part.pos + +func newParticle*(name: string, pos: Vector, mass: float): Particle = + Particle(name: name, pos: pos, vel: vector0(), acc: vector0(), mass: mass) + +func massFromName*(name: string): float = + 1.0 # TODO implement \ No newline at end of file diff --git a/vis.nim b/vis.nim deleted file mode 100644 index ec21997..0000000 --- a/vis.nim +++ /dev/null @@ -1,254 +0,0 @@ -import sdl2 -import sdl2/ttf - -import strutils -import linalg -import sequtils -import strformat -import os - -var filename = "output.txt" -if paramCount() > 0: - filename = paramStr(1) -let lines = filename.readFile().split('\n').filterIt(it.len() > 0) - -const - windowWidth = 1000 - windowHeight = 1000 - - ballRadius = 8 - textHeight = 35 - - # assumed to be in ratio with the win width/height - xMin = -50.0 - xMax = 50.0 - yMin = -50.0 - yMax = 50.0 - -type - Particle = object - pos: Vector - vel: Vector - name: string - -var lineN = 1 -proc parseLine(line: string, t: var float): seq[Particle] = - type - TokenType = enum - ttFloat, ttString - Token = object - lineN: int - case kind: TokenType: - of ttFloat: - floatVal: float - of ttString: - stringVal: string - - var - c = 1 - tokens: seq[Token] - - proc advance: char = - c.inc - line[c-1] - - proc peek: char = - if c > line.high(): - '\0' - else: - line[c] - - proc scanNum: float = - let start = c - 1 - while peek() in {'0'..'9'}: - discard advance() - if peek() == '.': - discard advance() - while peek() in {'0'..'9'}: - discard advance() - if peek() == 'e': - discard advance() - if peek() == '-' or peek() == '+': - discard advance() - while peek() in {'0'..'9'}: - discard advance() - return parseFloat(line[start..c-1]) - - proc scanIdent: string = - let start = c - 1 - while peek() in {'a' .. 'z'}: - discard advance() - return line[start..c-1] - - t = scanNum() - while c < line.len(): - case advance(): - of '-': - tokens.add(Token(lineN: lineN, kind: ttFloat, floatVal: scanNum())) - of '0'..'9': - tokens.add(Token(lineN: lineN, kind: ttFloat, floatVal: scanNum())) - of 'a' .. 'z': - tokens.add(Token(lineN: lineN, kind: ttString, stringVal: scanIdent())) - else: - discard #whitespace or [ or ] - - var i = 0 - while i < tokens.len(): - try: - let name = tokens[i].stringVal - let x = tokens[i+1].floatVal - let y = tokens[i+2].floatVal - let vx = tokens[i+3].floatVal - let vy = tokens[i+4].floatVal - i += 5 - result.add(Particle(name: name, pos: vector(x, y), vel: vector(vx, vy))) - except: - echo "parsing crashed on line" & $tokens[i].lineN - raise newException(Defect, "Parser crash") - - inc lineN - - -proc vec2scr(v: Vector, x, y: var cint) = - let xper = (v.x - xMin) / (xMax - xMin) - let yper = (v.y - yMin) / (yMax - yMin) - x = cint(xper * windowWidth) - y = cint(yper * windowHeight) - -# traces -var trace: array[2000, Vector] -var traceIndex = 0 - -proc drawRectAt(v: Vector, renderer: RendererPtr) = - trace[traceIndex] = v - traceIndex.inc() - if traceIndex > trace.high(): - traceIndex = 0 - # how many % of the screen it is at - var x, y: cint - vec2scr(v, x, y) - if x < ballRadius or y < ballRadius: - return - if x + ballRadius > windowWidth or y + ballRadius > windowHeight: - return - x -= ballRadius - y -= ballRadius - let size = cint(2 * ballRadius) - var r = rect(x, y, size, size) - renderer.fillRect(r) - -proc drawArrow(start: Vector, delta: Vector, renderer: RendererPtr) = - var x, y: cint - var ex, ey: cint - vec2scr(start, x, y) - vec2scr(start+delta, ex, ey) - renderer.drawLine(x, y, ex, ey) - -proc drawTrace(renderer: RendererPtr) = - var x, y: cint - for v in trace: - vec2scr(v, x, y) - renderer.drawPoint(x, y) - -proc draw(renderer: RendererPtr, font: FontPtr, frame: int, speed: int) = - let line = lines[frame] - var t = 0.0 - let parts = parseLine(line, t) - - renderer.setDrawColor(0, 0, 0, 255) - renderer.clear() - - renderer.setDrawColor(255, 255, 255, 255) - for i in 0..parts.high: - let part = parts[i] - part.pos.drawRectAt(renderer) - #part.pos.drawArrow(part.vel, renderer) - - renderer.setDrawColor(0, 0, 255, 255) - #renderer.drawTrace() - - let color = color(255, 255, 255, 0) - let text = &"t = {t:>09.2f} frame = {frame} speed = {speed}" - let surface = ttf.renderTextSolid(font, text, color) - let texture = renderer.createTextureFromSurface(surface) - surface.freeSurface() - defer: texture.destroy() - let textWidth: int = text.len() * int(float(textHeight) * 0.6) - var r = rect( - cint((windowWidth - textWidth) div 2), - 0, - cint(textWidth), - textHeight - ) - renderer.copy(texture, nil, addr r) - - renderer.present() - -type SDLException = object of Defect -template sdlFailIf(cond: bool, reason: string) = - if cond: - raise newException(SDLException, reason & ", SDL error " & $getError()) - -proc main = - - sdlFailIf(not sdl2.init(INIT_VIDEO or INIT_TIMER or INIT_EVENTS)): - "SDL2 init required" - defer: sdl2.quit() - - let win = createWindow( - title = "vis", - x = SDL_WINDOWPOS_CENTERED, - y = SDL_WINDOWPOS_CENTERED, - w = windowWidth, - h = windowHeight, - flags = SDL_WINDOW_SHOWN - ) - - sdlFailIf (win.isNil): "window couldn't be created" - defer: win.destroy() - - let renderer = createRenderer( - window = win, - index = -1, - flags = Renderer_Accelerated or Renderer_PresentVsync or Renderer_TargetTexture - ) - - sdlFailIf renderer.isNil: "renderer couldn't be created" - defer: renderer.destroy() - - sdlFailIf(not ttfInit()): "TTF init failed" - defer: ttfQuit() - - let font = ttf.openFont("FreeMono.ttf", textHeight) - sdlFailIf(font.isNil): "font couldn't be created" - - var - running = true - frame = 0 - speed = 1 - - while running: - var event = defaultEvent - while pollEvent(event): - case event.kind: - of QuitEvent: - running = false - break - of KeyDown: - if event.key.keysym.scancode == SDL_SCANCODE_UP: - speed *= 2 - elif event.key.keysym.scancode == SDL_SCANCODE_DOWN: - if speed > 1 or speed < -1: - speed = speed div 2 - elif event.key.keysym.scancode == SDL_SCANCODE_LEFT: - speed = -speed - else: - discard - frame += speed - if frame < 1: - frame = 1 - if frame > lines.high(): - frame = lines.high() - draw(renderer, font, frame, speed) - -main() \ No newline at end of file