change to xyz usage, use jmol to visualize

This commit is contained in:
Art 2023-04-05 10:57:56 +02:00
parent fd610548c9
commit 16a4a8f2d9
Signed by: prod2
GPG Key ID: F3BB5A97A70A8DDE
11 changed files with 274 additions and 459 deletions

3
.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
main
genrandom
*.xyz

Binary file not shown.

16
ffs.toml Normal file
View File

@ -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

108
forces.nim Normal file
View File

@ -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)

53
genrandom.nim Normal file
View File

@ -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()

View File

@ -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

View File

@ -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}"

187
main.nim
View File

@ -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()

39
parsexyz.nim Normal file
View File

@ -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}")

20
particle.nim Normal file
View File

@ -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

254
vis.nim
View File

@ -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()