From 75ef8db0d058158c10951a4184186f8c6cc27acc Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Sat, 27 Mar 2021 22:35:43 +0100 Subject: Initial public commit of this basic MD code Simulation of interacting particles --- benchmark.py | 46 ++++ examples/argon.py | 39 ++++ examples/dimer_only.py | 36 ++++ examples/water.py | 33 +++ interacticle/__init__.py | 4 + interacticle/kernel.cl | 460 ++++++++++++++++++++++++++++++++++++++++ interacticle/simulation.py | 394 ++++++++++++++++++++++++++++++++++ interacticle/visual/__init__.py | 5 + interacticle/visual/box.py | 67 ++++++ interacticle/visual/camera.py | 106 +++++++++ interacticle/visual/links.py | 48 +++++ interacticle/visual/plot.py | 248 ++++++++++++++++++++++ interacticle/visual/shader.py | 28 +++ interacticle/visual/view.py | 242 +++++++++++++++++++++ interacticle/visualizer.py | 99 +++++++++ library.py | 84 ++++++++ shell.nix | 46 ++++ 17 files changed, 1985 insertions(+) create mode 100644 benchmark.py create mode 100644 examples/argon.py create mode 100644 examples/dimer_only.py create mode 100644 examples/water.py create mode 100644 interacticle/__init__.py create mode 100644 interacticle/kernel.cl create mode 100644 interacticle/simulation.py create mode 100644 interacticle/visual/__init__.py create mode 100644 interacticle/visual/box.py create mode 100644 interacticle/visual/camera.py create mode 100644 interacticle/visual/links.py create mode 100644 interacticle/visual/plot.py create mode 100644 interacticle/visual/shader.py create mode 100644 interacticle/visual/view.py create mode 100644 interacticle/visualizer.py create mode 100644 library.py create mode 100644 shell.nix diff --git a/benchmark.py b/benchmark.py new file mode 100644 index 0000000..26bde57 --- /dev/null +++ b/benchmark.py @@ -0,0 +1,46 @@ +import sys +import numpy as np +from timeit import default_timer as timer + +from interacticle import MoleculeCollection, LennardJones, Coulomb, Simulation + +from library import Argon + +setup = MoleculeCollection() + +setup.potential(LennardJones(39.948, 0.3395, 0.971)) + +for x in range(10): + for y in range(10): + for z in range(10): + setup.add(Argon, (0.2+0.5*x,0.2+0.5*y,0.2+0.5*z)) + +setup.domain_size = 5 +setup.tau = 0.0005 +setup.cutoff = 0.3395*2.5 +setup.target_temperature = 300 + +setup.max_lennard_jones = 2000 +setup.max_coulomb = 2000 + +setup.statistics_step = 100 +setup.neighborhood_step = int(sys.argv[1]) if len(sys.argv) == 2 else 100 + +simulation = Simulation(setup, opengl = False) +simulation.setup() +simulation.verbose = True + +neighborhood_sizes = [ ] + +start = timer() + +for i in range(10000): + simulation.evolve() + if i % setup.statistics_step == 0: + neighborhood_sizes.append(simulation.get_neighborhood_characteristics()[1]) + +duration = timer() - start + +updates_per_second = 10000 / duration + +print(f"{setup.neighborhood_step}, {updates_per_second:.2f}, {np.max(neighborhood_sizes):.2f}") diff --git a/examples/argon.py b/examples/argon.py new file mode 100644 index 0000000..6eead38 --- /dev/null +++ b/examples/argon.py @@ -0,0 +1,39 @@ +from interacticle import MoleculeCollection, LennardJones, Coulomb, Simulation + +import interacticle.visualizer +from interacticle.visual import WireBox, MolecularLinks, VelocityHistogram, TemperaturePlot + +from library import Argon + +setup = MoleculeCollection() + +setup.potential(LennardJones(39.948, 0.3395, 0.971)) + +for x in range(10): + for y in range(10): + for z in range(10): + setup.add(Argon, (2+0.4*x,2+0.4*y,2+0.4*z)) + +def target_temperature(step, temperature): + if step < 4e4: # equilibration for 20 ps + return None + elif step >= 4e4 and step <= 2e5: # Heat up to 300 Kelvin in 20 ps and keep there until t=100ps + return min(max(temperature, (step-4e4) / 4e4 * 300), 300) + else: # Cool down to 10 Kelvin within 20 ps + return max(min(temperature, (4e4 - (step-2e5)) / 4e4 * 300), 10) + +setup.domain_size = 8 +setup.tau = 0.0005 +setup.cutoff = 0.3395*2.5 +setup.target_temperature = target_temperature + +setup.max_lennard_jones = 200 +setup.max_coulomb = 200 + +simulation = Simulation(setup, opengl = True) +simulation.verbose = False + +histogram = VelocityHistogram(simulation, [1.1,0.5], [1,0.5]) +temperature = TemperaturePlot(simulation, [1.1,0], [1,0.5]) + +interacticle.visualizer.simulate(setup, simulation, [ WireBox(setup.domain_size), MolecularLinks(simulation), histogram, temperature ], 100) diff --git a/examples/dimer_only.py b/examples/dimer_only.py new file mode 100644 index 0000000..c58eb39 --- /dev/null +++ b/examples/dimer_only.py @@ -0,0 +1,36 @@ +from interacticle import MoleculeCollection, LennardJones, Coulomb, Simulation + +import interacticle.visualizer +from interacticle.visual import WireBox, MolecularLinks + +from library import Benzene + +setup = MoleculeCollection() + +setup.potential(LennardJones(12.011, 0.3550, 0.070)) +setup.potential(LennardJones( 1.008, 0.2420, 0.030)) + +setup.potential(Coulomb(12.011, -0.115)) +setup.potential(Coulomb( 1.008, 0.115)) + +# T-shaped dimer +setup.add(Benzene, (1,1,1.519), rotations=[([0,1,0], 15), ([1,0,0], -9), ([1,0,0], 90), ([0,1,0], 90)]) +setup.add(Benzene, (1,1,1.0), rotations=[([0,1,0], 15), ([1,0,0], -9), ([0,0,1], 30)]) + +# Sandwich-shaped dimer +#setup.add(Benzene, (1,1,1), rotations=[([0,1,0], 15), ([1,0,0], -9), ([1,0,0], 90), ([0,1,0], 90)]) +#setup.add(Benzene, (1,1.377,1), rotations=[([0,1,0], 15), ([1,0,0], -9), ([1,0,0], 90), ([0,1,0], 60)]) + +# Parallel displaced dimer +#setup.add(Benzene, (1,1,1), rotations=[([0,1,0], 15), ([1,0,0], -9), ([1,0,0], 90), ([0,1,0], 90)]) +#setup.add(Benzene, (1.5,1.3,1), rotations=[([0,1,0], 15), ([1,0,0], -9), ([1,0,0], 90), ([0,1,0], 90)]) + +setup.domain_size = 2 +setup.tau = 0.0005 +setup.cutoff = 0.242*2.5 +setup.statistics_step = 100 +setup.neighborhood_step = 100 + +simulation = Simulation(setup, opengl = True) + +interacticle.visualizer.simulate(setup, simulation, [ WireBox(setup.domain_size), MolecularLinks(simulation) ], 100) diff --git a/examples/water.py b/examples/water.py new file mode 100644 index 0000000..64d62d3 --- /dev/null +++ b/examples/water.py @@ -0,0 +1,33 @@ +from interacticle import MoleculeCollection, LennardJones, Coulomb, Simulation + +import interacticle.visualizer +from interacticle.visual import WireBox, MolecularLinks, VelocityHistogram, TemperaturePlot + +from library import Water + +setup = MoleculeCollection() + +setup.potential(LennardJones(15.9992, 0.31776, 0.792324)) +setup.potential(Coulomb(15.9992, -0.8450)) +setup.potential(Coulomb( 1.0079, 0.4225)) + +for x in range(4): + for y in range(4): + for z in range(4): + setup.add(Water, (0.2+0.7*x,0.2+0.7*y,0.2+0.7*z)) + +setup.domain_size = 5 +setup.tau = 0.0001 +setup.cutoff = 2.6 + +setup.max_lennard_jones = 500 +setup.max_coulomb = 500 + +setup.statistics_step = 100 +setup.neighborhood_step = 500 + +simulation = Simulation(setup, opengl = True) + +histogram = VelocityHistogram(simulation, [1.1,0.5], [1,0.5]) + +interacticle.visualizer.simulate(setup, simulation, [ WireBox(setup.domain_size), MolecularLinks(simulation), histogram ], 100) diff --git a/interacticle/__init__.py b/interacticle/__init__.py new file mode 100644 index 0000000..444c908 --- /dev/null +++ b/interacticle/__init__.py @@ -0,0 +1,4 @@ +import interacticle.visual +import interacticle.visualizer + +from interacticle.simulation import Atom, Bond, Angle, Torsion, LennardJones, Coulomb, Molecule, MoleculeCollection, Simulation diff --git a/interacticle/kernel.cl b/interacticle/kernel.cl new file mode 100644 index 0000000..5574efd --- /dev/null +++ b/interacticle/kernel.cl @@ -0,0 +1,460 @@ +typedef float scalar_t; +typedef float3 vec_t; +typedef float4 data_vec_t; + +vec_t get(__global data_vec_t* data, + unsigned int iParticle) +{ + return data[iParticle].xyz; +} + +scalar_t sq(scalar_t x) { + return pown(x,2); +} + +vec_t lennard_jones(__global data_vec_t* pos, unsigned int i, unsigned int j, scalar_t sigma, scalar_t epsilon, vec_t shift, scalar_t cutoff) { + vec_t iPos = get(pos, i); + vec_t jPos = get(pos, j) + shift; + + vec_t rij = jPos - iPos; + + scalar_t r = length(rij); + scalar_t s = pown(sigma / r, 6); + scalar_t f = 24 * epsilon * s / sq(r) * (1 - 2*s); + + return (r < cutoff) * f * rij; +} + +vec_t coulomb(__global data_vec_t* pos, unsigned int i, unsigned int j, scalar_t charge, vec_t shift, scalar_t cutoff) { + vec_t iPos = get(pos, i); + vec_t jPos = get(pos, j) + shift; + + vec_t rij = jPos - iPos; + scalar_t r = length(rij); + // Coulomb's constant is 138.935458 kJ mol^-1 nm e^-2 in MD units + scalar_t f = -138.935458 * charge * (1 / pown(r,3)); + + return (r < cutoff) * f * rij; +} + +vec_t bond(__global data_vec_t* pos, unsigned int i, unsigned int j, scalar_t kb, scalar_t r0) { + vec_t iPos = get(pos, i); + vec_t jPos = get(pos, j); + + vec_t rij = jPos - iPos; + scalar_t r = length(rij); + + return kb * (1 - r0 / r) * rij; +} + +vec_t angle(__global data_vec_t* pos, unsigned int iParticle, unsigned int i, unsigned int j, unsigned int k, scalar_t theta0, scalar_t k0) { + vec_t iPos = get(pos, i); + vec_t jPos = get(pos, j); + vec_t kPos = get(pos, k); + + vec_t rij = jPos - iPos; + vec_t rjk = kPos - jPos; + vec_t rkj = jPos - kPos; + + scalar_t S = dot(rij, rkj); + scalar_t D = length(rij) * length(rkj); + + scalar_t theta = acos(S/D); + + scalar_t sd2 = S / pown(D,2); + + scalar_t factor = k0 * sin(theta - theta0) / sin(theta) * (1/D); + vec_t fi = factor * (sd2 * sq(length(rkj)) * rij - rkj); + vec_t fk = factor * (sd2 * sq(length(rij)) * rkj - rij); + + return (iParticle == i) * fi + + (iParticle == k) * fk + + (iParticle == j) * -(fi+fk); +} + +vec_t torsion(__global data_vec_t* pos, unsigned int iParticle, unsigned int i, unsigned int j, unsigned int k, unsigned int l, scalar_t phi0, scalar_t k0) { + vec_t iPos = get(pos, i); + vec_t jPos = get(pos, j); + vec_t kPos = get(pos, k); + vec_t lPos = get(pos, l); + + vec_t rij = jPos - iPos; + vec_t rjk = kPos - jPos; + vec_t rik = kPos - iPos; + vec_t rkl = lPos - kPos; + + vec_t m = cross(rij,rjk); + vec_t n = cross(rjk,rkl); + + scalar_t phi = M_PI + sign(dot(rij,n)) * acos(clamp(dot(m,n)/(length(m)*length(n)),-1.,1.)); + scalar_t pot = -k0 * sin(phi - phi0); + + vec_t fi = - pot * length(rjk) / sq(length(m)) * m; + vec_t fl = pot * length(rjk) / sq(length(n)) * n; + vec_t fj = -dot(rij,rjk)/sq(length(rjk)) * fi - fi + dot(rkl,rjk)/sq(length(rjk)) * fl; + vec_t fk = -dot(rkl,rjk)/sq(length(rjk)) * fl - fl + dot(rij,rjk)/sq(length(rjk)) * fi; + + return (iParticle == i) * fi + + (iParticle == l) * fl + + (iParticle == j) * fj + + (iParticle == k) * fk; +} + +__kernel void compute_lennard_jones(__global data_vec_t* pos, + __global data_vec_t* force, + __global unsigned int* lj_count, + __global unsigned int* lj_indices, + __global scalar_t* lj_sigma, + __global scalar_t* lj_epsilon, + __global data_vec_t* lj_shift, + scalar_t cutoff) +{ + unsigned int iParticle = get_global_id(0); + unsigned int nNeighbors = lj_count[iParticle]; + + vec_t f = 0; + + for (unsigned int iNeighbor=iParticle*$max_lj; iNeighbor < iParticle*$max_lj + nNeighbors; ++iNeighbor) { + f += lennard_jones(pos, iParticle, + lj_indices[iNeighbor], + lj_sigma[iNeighbor], + lj_epsilon[iNeighbor], + lj_shift[iNeighbor].xyz, + cutoff); + } + + force[iParticle].xyz += f; +} + +__kernel void compute_coulomb(__global data_vec_t* pos, + __global data_vec_t* force, + __global unsigned int* coulomb_count, + __global unsigned int* coulomb_indices, + __global scalar_t* coulomb_charge, + __global data_vec_t* coulomb_shift, + scalar_t cutoff) +{ + unsigned int iParticle = get_global_id(0); + unsigned int nNeighbors = coulomb_count[iParticle]; + + vec_t f = 0; + + for (unsigned int iNeighbor=iParticle*$max_coulomb; iNeighbor < iParticle*$max_coulomb + nNeighbors; ++iNeighbor) { + f += coulomb(pos, iParticle, + coulomb_indices[iNeighbor], + coulomb_charge[iNeighbor], + coulomb_shift[iNeighbor].xyz, + cutoff); + } + + force[iParticle].xyz += f; +} + +__kernel void compute_bonds(__global data_vec_t* pos, + __global data_vec_t* force, + + __global unsigned int* bond_count, + __global unsigned int* bond_indices, + __global scalar_t* bond_kb, + __global scalar_t* bond_k0) +{ + unsigned int iParticle = get_global_id(0); + + vec_t f = 0; + + unsigned int nBonds = bond_count[iParticle]; + for (unsigned int iBond=iParticle*$max_bonds; iBond < iParticle*$max_bonds + nBonds; ++iBond) { + f += bond(pos, + iParticle, + bond_indices[iBond], + bond_kb[iBond], + bond_k0[iBond]); + } + + force[iParticle].xyz += f; +} + +__kernel void compute_angles(__global data_vec_t* pos, + __global data_vec_t* force, + + __global unsigned int* angle_count, + __global unsigned int* angle_indices, + __global scalar_t* angle_theta0, + __global scalar_t* angle_r0) +{ + unsigned int iParticle = get_global_id(0); + + vec_t f = 0; + + unsigned int nAngles = angle_count[iParticle]; + for (unsigned int iAngle=iParticle*$max_angles; iAngle < iParticle*$max_angles + nAngles; ++iAngle) { + f += angle(pos, iParticle, + angle_indices[0*$n_angles+iAngle], + angle_indices[1*$n_angles+iAngle], + angle_indices[2*$n_angles+iAngle], + angle_theta0[iAngle], + angle_r0[iAngle]); + } + + force[iParticle].xyz += f; +} + +__kernel void compute_torsions(__global data_vec_t* pos, + __global data_vec_t* force, + + __global unsigned int* torsion_count, + __global unsigned int* torsion_indices, + __global scalar_t* torsion_phi0, + __global scalar_t* torsion_r0) +{ + unsigned int iParticle = get_global_id(0); + + vec_t f = 0; + + unsigned int nTorsions = torsion_count[iParticle]; + for (unsigned int iTorsion=iParticle*$max_torsions; iTorsion < iParticle*$max_torsions + nTorsions; ++iTorsion) { + f += torsion(pos, iParticle, + torsion_indices[0*$n_torsions+iTorsion], + torsion_indices[1*$n_torsions+iTorsion], + torsion_indices[2*$n_torsions+iTorsion], + torsion_indices[3*$n_torsions+iTorsion], + torsion_phi0[iTorsion], + torsion_r0[iTorsion]); + } + + force[iParticle].xyz += f; +} + +__kernel void compute_intramolecular(__global data_vec_t* pos, + __global data_vec_t* force, + + __global unsigned int* bond_count, + __global unsigned int* bond_indices, + __global scalar_t* bond_kb, + __global scalar_t* bond_k0, + + __global unsigned int* angle_count, + __global unsigned int* angle_indices, + __global scalar_t* angle_theta0, + __global scalar_t* angle_r0, + + __global unsigned int* torsion_count, + __global unsigned int* torsion_indices, + __global scalar_t* torsion_phi0, + __global scalar_t* torsion_r0) +{ + unsigned int iParticle = get_global_id(0); + + vec_t f = 0; + + unsigned int nBonds = bond_count[iParticle]; + for (unsigned int iBond=iParticle*$max_bonds; iBond < iParticle*$max_bonds + nBonds; ++iBond) { + f += bond(pos, + iParticle, + bond_indices[iBond], + bond_kb[iBond], + bond_k0[iBond]); + } + + unsigned int nAngles = angle_count[iParticle]; + for (unsigned int iAngle=iParticle*$max_angles; iAngle < iParticle*$max_angles + nAngles; ++iAngle) { + f += angle(pos, iParticle, + angle_indices[0*$n_angles+iAngle], + angle_indices[1*$n_angles+iAngle], + angle_indices[2*$n_angles+iAngle], + angle_theta0[iAngle], + angle_r0[iAngle]); + } + + unsigned int nTorsions = torsion_count[iParticle]; + for (unsigned int iTorsion=iParticle*$max_torsions; iTorsion < iParticle*$max_torsions + nTorsions; ++iTorsion) { + f += torsion(pos, iParticle, + torsion_indices[0*$n_torsions+iTorsion], + torsion_indices[1*$n_torsions+iTorsion], + torsion_indices[2*$n_torsions+iTorsion], + torsion_indices[3*$n_torsions+iTorsion], + torsion_phi0[iTorsion], + torsion_r0[iTorsion]); + } + + force[iParticle].xyz += f; +} + +__kernel void evolve_x(__global data_vec_t* pos, + __global data_vec_t* vel, + __global data_vec_t* force_prev, + __global data_vec_t* force_curr, + scalar_t tau) +{ + unsigned int iParticle = get_global_id(0); + + vec_t f_prev = get(force_prev, iParticle); + vec_t f_curr = get(force_curr, iParticle); + + vec_t p = get(pos, iParticle); + vec_t v = get(vel, iParticle); + + scalar_t a = tau * 0.5 / pos[iParticle].w; + p += tau * (v + a * f_curr); + + pos[iParticle].xyz = p; + force_prev[iParticle].xyz = f_curr; + force_curr[iParticle].xyz = 0; +} + +__kernel void evolve_v(__global data_vec_t* pos, + __global data_vec_t* vel, + __global data_vec_t* force_prev, + __global data_vec_t* force_curr, + scalar_t tau) +{ + unsigned int iParticle = get_global_id(0); + + vec_t f_prev = get(force_prev, iParticle); + vec_t f_curr = get(force_curr, iParticle); + + vec_t v = get(vel, iParticle); + + scalar_t a = tau * 0.5 / pos[iParticle].w; + v += a * (f_prev + f_curr); + + vel[iParticle].xyz = v; +} + +scalar_t lennard_jones_sigma(scalar_t iMass) +{ + return $lj_sigma_expr; +} + +scalar_t lennard_jones_epsilon(scalar_t iMass) +{ + return $lj_epsilon_expr; +} + +scalar_t coulomb_q(scalar_t iMass) +{ + return $coulomb_charge_expr; +} + +__kernel void update_potential_neighborhoods(__global data_vec_t* pos, + __global unsigned int* molecules, + + __global unsigned int* lj_count, + __global unsigned int* lj_indices, + __global scalar_t* lj_sigma, + __global scalar_t* lj_epsilon, + __global data_vec_t* lj_shift, + + __global unsigned int* coulomb_count, + __global unsigned int* coulomb_indices, + __global scalar_t* coulomb_charge, + __global data_vec_t* coulomb_shift, + + scalar_t cutoff, + scalar_t skin) +{ + unsigned int iParticle = get_global_id(0); + + __global unsigned int* lj_i_count = lj_count + $max_lj * iParticle; + __global unsigned int* lj_i_indices = lj_indices + $max_lj * iParticle; + __global scalar_t* lj_i_sigma = lj_sigma + $max_lj * iParticle; + __global scalar_t* lj_i_epsilon = lj_epsilon + $max_lj * iParticle; + __global data_vec_t* lj_i_shift = lj_shift + $max_lj * iParticle; + + __global unsigned int* coulomb_i_count = coulomb_count + $max_coulomb * iParticle; + __global unsigned int* coulomb_i_indices = coulomb_indices + $max_coulomb * iParticle; + __global scalar_t* coulomb_i_charge = coulomb_charge + $max_coulomb * iParticle; + __global data_vec_t* coulomb_i_shift = coulomb_shift + $max_coulomb * iParticle; + + for (unsigned int i=0; i < $max_lj; ++i) { + lj_i_indices[i] = 0; + lj_i_sigma[i] = 0; + lj_i_epsilon[i] = 0; + lj_i_shift[i] = 0; + + coulomb_i_indices[i] = 0; + coulomb_i_charge[i] = 0; + coulomb_i_shift[i] = 0; + } + + unsigned int idx = 0; + + vec_t iPos = get(pos, iParticle); + + for (unsigned int jParticle=0; jParticle < $n_atoms; ++jParticle) { + for (int sX=-1; sX <= 1; ++sX) { + for (int sY=-1; sY <= 1; ++sY) { + for (int sZ=-1; sZ <= 1; ++sZ) { + vec_t shift = $domain_size * (vec_t)(sX,sY,sZ); + vec_t jPos = get(pos, jParticle) + shift; + + scalar_t r = length(jPos - iPos); + + unsigned int iMolecule = molecules[iParticle]; + unsigned int jMolecule = molecules[jParticle]; + + scalar_t iMass = pos[iParticle].w; + scalar_t jMass = pos[jParticle].w; + + scalar_t mask = (r < cutoff + skin) * (iParticle != jParticle) * (iMolecule != jMolecule); + + lj_i_indices[idx] += mask * jParticle; + lj_i_sigma[idx] += mask * 0.5*(lennard_jones_sigma(iMass) + lennard_jones_sigma(jMass)); + lj_i_epsilon[idx] += mask * sqrt(lennard_jones_epsilon(iMass) * lennard_jones_epsilon(jMass)); + lj_i_shift[idx].xyz += mask * shift; + + coulomb_i_indices[idx] += mask * jParticle; + coulomb_i_charge[idx] += mask * coulomb_q(iMass) * coulomb_q(jMass); + coulomb_i_shift[idx].xyz += mask * shift; + + idx += mask; + } + } + } + } + + lj_count[iParticle] = idx; + coulomb_count[iParticle] = idx; +} + +__kernel void wrap_molecules(__global data_vec_t* pos, + __global unsigned int* molecules) +{ + unsigned int iMolecule = get_global_id(0); + + int nAtoms = 0; + int xShift = 0; + int yShift = 0; + int zShift = 0; + + for (unsigned int jParticle=0; jParticle < $n_atoms; ++jParticle) { + nAtoms += (iMolecule == molecules[jParticle]); + xShift -= (iMolecule == molecules[jParticle]) * (pos[jParticle].x > $domain_size); + xShift += (iMolecule == molecules[jParticle]) * (pos[jParticle].x < 0.); + yShift -= (iMolecule == molecules[jParticle]) * (pos[jParticle].y > $domain_size); + yShift += (iMolecule == molecules[jParticle]) * (pos[jParticle].y < 0.); + zShift -= (iMolecule == molecules[jParticle]) * (pos[jParticle].z > $domain_size); + zShift += (iMolecule == molecules[jParticle]) * (pos[jParticle].z < 0.); + } + + for (unsigned int jParticle=0; jParticle < $n_atoms; ++jParticle) { + if (molecules[jParticle] == iMolecule) { + pos[jParticle].x += (xShift == nAtoms) * $domain_size; + pos[jParticle].x -= (xShift == -nAtoms) * $domain_size; + pos[jParticle].y += (yShift == nAtoms) * $domain_size; + pos[jParticle].y -= (yShift == -nAtoms) * $domain_size; + pos[jParticle].z += (zShift == nAtoms) * $domain_size; + pos[jParticle].z -= (zShift == -nAtoms) * $domain_size; + } + } +} + +__kernel void update_view_links(__global data_vec_t* pos, + __global unsigned int* link_indices, + __global data_vec_t* link_pos) +{ + unsigned int iLink = get_global_id(0); + + link_pos[2*iLink+0].xyz = pos[link_indices[2*iLink+0]].xyz; + link_pos[2*iLink+1].xyz = pos[link_indices[2*iLink+1]].xyz; +} diff --git a/interacticle/simulation.py b/interacticle/simulation.py new file mode 100644 index 0000000..6da1dd1 --- /dev/null +++ b/interacticle/simulation.py @@ -0,0 +1,394 @@ +import numpy as np +import scipy.linalg as sp + +import pyopencl as cl +mf = cl.mem_flags +from pyopencl.tools import get_gl_sharing_context_properties + +import OpenGL.GL as gl +from OpenGL.arrays import vbo + +from string import Template + +from collections import namedtuple, defaultdict +from itertools import product + +Atom = namedtuple('Atom', 'x y z vx vy vz mass') + +Bond = namedtuple('Bond', 'i j kb k0') +Angle = namedtuple('Angle', 'i j k theta0 r0') +Torsion = namedtuple('Torsion', 'i j k l phi0 r0') + +LennardJones = namedtuple('LennardJones', 'mass sigma epsilon') +Coulomb = namedtuple('Coulomb', 'mass charge') + +def rotationMatrix(axis, theta): + return sp.expm(np.cross(np.eye(3), axis/sp.norm(axis)*np.radians(theta))) + +class Molecule: + def __init__(self): + self.atoms = [ ] + self.connections = [ ] + +class MoleculeCollection: + def __init__(self): + self.index = 0 + self.current_molecule = 0 + self.indices_of_molecule = defaultdict(lambda: [ ]) + + self.atoms = [ ] + + self.molecules = [ ] + + self.neighborhood_step = 100 + self.statistics_step = 100 + + self.max_lennard_jones = 100 + self.max_coulomb = 100 + + self.lennard_jones = defaultdict(lambda: None) + self.coulomb = defaultdict(lambda: None) + + self.bonds = defaultdict(lambda: [ ]) + self.angles = defaultdict(lambda: [ ]) + self.torsions = defaultdict(lambda: [ ]) + + self.cutoff = 0.242*2.5 + self.skin = 0.1 + + self.tau = 0.0005 + self.domain_size = 2 + self.target_temperature = None + + def add(self, molecule, origin=(0,0,0), rotations=None): + for i, atom in enumerate(molecule.atoms): + x, y, z = atom.x, atom.y, atom.z + if rotations: + for rotation in rotations: + x, y, z = np.dot(rotationMatrix(*rotation), [x, y, z]) + self.atoms.append(Atom(origin[0]+x, origin[1]+y, origin[2]+z, atom.vy, atom.vy, atom.vz, atom.mass)) + self.indices_of_molecule[self.current_molecule].append(self.index+i) + self.molecules.append(self.current_molecule) + + for c in molecule.connections: + if isinstance(c, Bond): + self.bonds[self.index+c.i].append(Bond(self.index+c.i, self.index+c.j, c.kb, c.k0)) + self.bonds[self.index+c.j].append(Bond(self.index+c.j, self.index+c.i, c.kb, c.k0)) + elif isinstance(c, Angle): + tmp = Angle(self.index+c.i, self.index+c.j, self.index+c.k, + c.theta0, c.r0) + self.angles[self.index+c.i].append(tmp) + self.angles[self.index+c.j].append(tmp) + self.angles[self.index+c.k].append(tmp) + elif isinstance(c, Torsion): + tmp = Torsion(self.index+c.i, self.index+c.j, self.index+c.k, self.index+c.l, + c.phi0, c.r0) + self.torsions[self.index+c.i].append(tmp) + self.torsions[self.index+c.j].append(tmp) + self.torsions[self.index+c.k].append(tmp) + self.torsions[self.index+c.l].append(tmp) + else: + raise Exception('Unknown connection type') + + self.index = len(self.atoms) + self.current_molecule = self.current_molecule + 1 + + def potential(self, config): + if isinstance(config, LennardJones): + self.lennard_jones[config.mass] = config + elif isinstance(config, Coulomb): + self.coulomb[config.mass] = config + + def serialize(self): + n_atoms = len(self.atoms) + atom_pos, atom_vel = np.zeros((n_atoms,4), dtype=np.float32), np.zeros((n_atoms,4), dtype=np.float32) + for i, atom in enumerate(self.atoms): + atom_pos[i,:] = [atom.x, atom.y, atom.z, atom.mass] + atom_vel[i,:] = [atom.vx, atom.vy, atom.vz, 0] + + molecules = np.array(self.molecules, dtype=np.uint32) + + max_bonds = max(map(lambda x: len(x[1]), self.bonds.items()), default=0) + n_bonds = n_atoms * max_bonds + bond_count, bond_indices, bond_kb, bond_k0 = np.zeros((n_atoms,1), dtype=np.uint32), np.zeros((n_bonds,1), dtype=np.uint32), np.zeros((n_bonds,1), dtype=np.float32), np.zeros((n_bonds,1), dtype=np.float32) + for i, _ in enumerate(self.atoms): + bond_count[i] = len(self.bonds[i]) + curr = i * max_bonds + for j, bond in enumerate(self.bonds[i]): + bond_indices[curr+j] = bond.j + bond_kb[curr+j] = bond.kb + bond_k0[curr+j] = bond.k0 + + max_angles = max(map(lambda x: len(x[1]), self.angles.items()), default=0) + n_angles = n_atoms * max_angles + angle_count, angle_indices, angle_theta0, angle_r0 = np.zeros((n_atoms,1), dtype=np.uint32), np.zeros((3*n_angles,1), dtype=np.uint32), np.zeros((n_angles,1), dtype=np.float32), np.zeros((n_angles,1), dtype=np.float32) + for i, _ in enumerate(self.atoms): + angle_count[i] = len(self.angles[i]) + curr = i * max_angles + for j, angle in enumerate(self.angles[i]): + angle_indices[0*n_angles+curr+j] = angle.i + angle_indices[1*n_angles+curr+j] = angle.j + angle_indices[2*n_angles+curr+j] = angle.k + angle_theta0[curr+j] = np.radians(angle.theta0) + angle_r0[curr+j] = angle.r0 + + max_torsions = max(map(lambda x: len(x[1]), self.torsions.items()), default=0) + n_torsions = n_atoms * max_torsions + torsion_count, torsion_indices, torsion_phi0, torsion_r0 = np.zeros((n_atoms,1), dtype=np.uint32), np.zeros((4*n_torsions,1), dtype=np.uint32), np.zeros((n_torsions,1), dtype=np.float32), np.zeros((n_torsions,1), dtype=np.float32) + curr = 0 + for i, _ in enumerate(self.atoms): + torsion_count[i] = len(self.torsions[i]) + curr = i * max_torsions + for j, torsion in enumerate(self.torsions[i]): + torsion_indices[0*n_torsions+curr+j] = torsion.i + torsion_indices[1*n_torsions+curr+j] = torsion.j + torsion_indices[2*n_torsions+curr+j] = torsion.k + torsion_indices[3*n_torsions+curr+j] = torsion.l + torsion_phi0[curr+j] = np.radians(torsion.phi0) + torsion_r0[curr+j] = torsion.r0 + + return atom_pos, atom_vel, molecules, bond_count, bond_indices, bond_kb, bond_k0, angle_count, angle_indices, angle_theta0, angle_r0, torsion_count, torsion_indices, torsion_phi0, torsion_r0 + +def build_kernel(domain_size, n_atoms, n_bonds, n_angles, n_torsions, max_bonds, max_angles, max_torsions, max_lj, max_coulomb, lennard_jones, coulomb): + lj_sigma_expr = " + ".join(map(lambda lj: f"(iMass == {lj.mass}) * {lj.sigma}", lennard_jones)) + lj_epsilon_expr = " + ".join(map(lambda lj: f"(iMass == {lj.mass}) * {lj.epsilon}", lennard_jones)) + coulomb_charge_expr = " + ".join(map(lambda c: f"(iMass == {c.mass}) * {c.charge}", coulomb)) + + with open('interacticle/kernel.cl', 'r') as kernel_src: + return Template(kernel_src.read()).substitute( + domain_size = domain_size, + n_atoms = n_atoms, + n_bonds = n_bonds, + n_angles = n_angles, + n_torsions = n_torsions, + max_lj = max_lj, + max_coulomb = max_coulomb, + max_bonds = max_bonds, + max_angles = max_angles, + max_torsions = max_torsions, + lj_sigma_expr = lj_sigma_expr if lj_sigma_expr else "0", + lj_epsilon_expr = lj_epsilon_expr if lj_epsilon_expr else "0", + coulomb_charge_expr = coulomb_charge_expr if coulomb_charge_expr else "0") + +class Simulation: + def __init__(self, setup, opengl = False): + self.domain_size = np.float32(setup.domain_size) + self.cutoff = np.float32(setup.cutoff) + self.skin = np.float32(setup.skin) + self.adaptive_skin = setup.skin + self.neighborhood_step = setup.neighborhood_step + self.statistics_step = setup.statistics_step + self.tau = np.float32(setup.tau) + self.target_temperature = setup.target_temperature + self.verbose = True + + self.np_atom_pos, self.np_atom_vel, self.np_molecules, self.np_bond_count, self.np_bond_indices, self.np_bond_kb, self.np_bond_k0, self.np_angle_count, self.np_angle_indices, self.np_angle_theta0, self.np_angle_r0, self.np_torsion_count, self.np_torsion_indices, self.np_torsion_phi0, self.np_torsion_r0 = setup.serialize() + + self.n_molecules = np.max(self.np_molecules)+1 + self.n_atoms = self.np_atom_pos.shape[0] + self.n_bonds = self.np_bond_kb.shape[0] + self.n_angles = self.np_angle_r0.shape[0] + self.n_torsions = self.np_torsion_r0.shape[0] + + self.max_lj = setup.max_lennard_jones + self.max_coulomb = setup.max_coulomb + self.max_bonds = max(self.np_bond_count)[0] + self.max_angles = max(self.np_angle_count)[0] + self.max_torsions = max(self.np_torsion_count)[0] + + self.n_lj = self.n_atoms * setup.max_lennard_jones + self.np_lj_count, self.np_lj_indices, self.np_lj_sigma, self.np_lj_epsilon = np.zeros((self.n_atoms,1), dtype=np.uint32), np.zeros((self.n_lj,1), dtype=np.uint32), np.zeros((self.n_lj,1), dtype=np.float32), np.zeros((self.n_lj,1), dtype=np.float32) + + self.n_coulomb = self.n_atoms * setup.max_coulomb + self.np_coulomb_count, self.np_coulomb_indices, self.np_coulomb_charge = np.zeros((self.n_atoms,1), dtype=np.uint32), np.zeros((self.n_coulomb,1), dtype=np.uint32), np.zeros((self.n_coulomb,1), dtype=np.float32) + + self.step = 0 + self.opengl = opengl + + self.kernel_src = build_kernel(setup.domain_size, self.n_atoms, self.n_bonds, self.n_angles, self.n_torsions, self.max_bonds, self.max_angles, self.max_torsions, self.max_lj, self.max_coulomb, setup.lennard_jones.values(), setup.coulomb.values()) + + def setup(self): + self.platform = cl.get_platforms()[0] + if self.opengl: + self.context = cl.Context( + properties=[(cl.context_properties.PLATFORM, self.platform)] + get_gl_sharing_context_properties()) + else: + self.context = cl.Context( + properties=[(cl.context_properties.PLATFORM, self.platform)]) + self.queue = cl.CommandQueue(self.context) + self.program = cl.Program(self.context, self.kernel_src).build( + '-cl-single-precision-constant') + + self.cl_molecules = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_molecules) + + self.np_force = np.zeros((self.n_atoms,4), dtype=np.float32) + + self.cl_force_prev = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_force) + self.cl_force_curr = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_force) + + self.cl_velocity = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_atom_vel) + + self.np_lj_shift = np.zeros((self.n_lj,4), dtype=np.float32) + + self.cl_lj_count = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_lj_count) + self.cl_lj_indices = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_lj_indices) + self.cl_lj_sigma = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_lj_sigma) + self.cl_lj_epsilon = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_lj_epsilon) + self.cl_lj_shift = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_lj_shift) + + self.np_coulomb_shift = np.zeros((self.n_coulomb,4), dtype=np.float32) + + self.cl_coulomb_count = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_coulomb_count) + self.cl_coulomb_indices = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_coulomb_indices) + self.cl_coulomb_charge = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_coulomb_charge) + self.cl_coulomb_shift = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_coulomb_shift) + + if self.n_bonds > 0: + self.cl_bond_count = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_bond_count) + self.cl_bond_indices = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_bond_indices) + self.cl_bond_kb = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_bond_kb) + self.cl_bond_k0 = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_bond_k0) + + if self.n_angles > 0: + self.cl_angle_count = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_angle_count) + self.cl_angle_indices = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_angle_indices) + self.cl_angle_theta0 = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_angle_theta0) + self.cl_angle_r0 = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_angle_r0) + + if self.n_torsions > 0: + self.cl_torsion_count = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_torsion_count) + self.cl_torsion_indices = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_torsion_indices) + self.cl_torsion_phi0 = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_torsion_phi0) + self.cl_torsion_r0 = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_torsion_r0) + + if self.opengl: + self.gl_position = vbo.VBO(data=self.np_atom_pos, usage=gl.GL_DYNAMIC_DRAW, target=gl.GL_ARRAY_BUFFER) + self.gl_position.bind() + self.cl_position = cl.GLBuffer(self.context, mf.READ_WRITE, int(self.gl_position)) + cl.enqueue_acquire_gl_objects(self.queue, [self.cl_position]) + else: + self.cl_position = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_atom_pos) + + self.kernel_args = (self.queue, (self.n_atoms,), None, self.cl_position) + + if self.target_temperature and not callable(self.target_temperature): + self.randomize_velocities(self.target_temperature) + + def compute_intramolecular_forces(self): + if self.n_bonds > 0 and self.n_angles > 0 and self.n_torsions > 0: + self.program.compute_intramolecular( + *self.kernel_args, + self.cl_force_curr, + self.cl_bond_count, self.cl_bond_indices, self.cl_bond_kb, self.cl_bond_k0, + self.cl_angle_count, self.cl_angle_indices, self.cl_angle_theta0, self.cl_angle_r0, + self.cl_torsion_count, self.cl_torsion_indices, self.cl_torsion_phi0, self.cl_torsion_r0) + else: + if self.n_bonds > 0: + self.program.compute_bonds( + *self.kernel_args, + self.cl_force_curr, + self.cl_bond_count, self.cl_bond_indices, self.cl_bond_kb, self.cl_bond_k0) + if self.n_angles > 0: + self.program.compute_angles( + *self.kernel_args, + self.cl_force_curr, + self.cl_angle_count, self.cl_angle_indices, self.cl_angle_theta0, self.cl_angle_r0) + if self.n_torsions > 0: + self.program.compute_torsions( + *self.kernel_args, + self.cl_force_curr, + self.cl_torsion_count, self.cl_torsion_indices, self.cl_torsion_phi0, self.cl_torsion_r0) + + def compute_intermolecular_forces(self): + self.program.compute_lennard_jones( + *self.kernel_args, + self.cl_force_curr, + self.cl_lj_count, self.cl_lj_indices, self.cl_lj_sigma, self.cl_lj_epsilon, self.cl_lj_shift, + self.cutoff) + self.program.compute_coulomb( + *self.kernel_args, + self.cl_force_curr, + self.cl_coulomb_count, self.cl_coulomb_indices, self.cl_coulomb_charge, self.cl_coulomb_shift, + self.cutoff) + + def update_neighborhoods(self): + self.program.wrap_molecules(*self.kernel_args, self.cl_molecules) + self.program.update_potential_neighborhoods( + *self.kernel_args, + self.cl_molecules, + self.cl_lj_count, self.cl_lj_indices, self.cl_lj_sigma, self.cl_lj_epsilon, self.cl_lj_shift, + self.cl_coulomb_count, self.cl_coulomb_indices, self.cl_coulomb_charge, self.cl_coulomb_shift, + self.cutoff, self.skin).wait() + + def evolve(self): + if self.opengl: + cl.enqueue_acquire_gl_objects(self.queue, [self.cl_position]) + + if self.step % self.statistics_step == 0: + self.read_velocities() + temperature, rmsv = self.get_temperature(), self.get_rms_velocity() + if self.target_temperature and temperature > 0: + if callable(self.target_temperature): + target = self.target_temperature(self.step, temperature) + if target: + self.scale_velocities(np.sqrt(target/temperature)) + else: + self.scale_velocities(np.sqrt(self.target_temperature/temperature)) + if self.verbose: + print(f"t={self.step*self.tau:.2f} T={temperature:.0f}, V_rms={rmsv:.0f}") + + if self.step % self.neighborhood_step == 0: + if self.adaptive_skin: + skin_bound = np.max(np.linalg.norm(self.np_atom_vel[:,0:3], axis=1)) * self.tau * self.neighborhood_step + self.skin = np.float32(max(2*skin_bound, self.adaptive_skin)) + self.update_neighborhoods() + if self.verbose: + min_n, max_n, mean_n = self.get_neighborhood_characteristics() + print(f"t={self.step*self.tau:.2f} minN={min_n}, maxN={max_n}, meanN={mean_n}, skin={self.skin:.3f}") + + self.program.evolve_x(*self.kernel_args, self.cl_velocity, self.cl_force_prev, self.cl_force_curr, self.tau) + if self.n_bonds > 0 or self.n_angles > 0 or self.n_torsions > 0: + self.compute_intramolecular_forces() + self.compute_intermolecular_forces() + self.program.evolve_v(*self.kernel_args, self.cl_velocity, self.cl_force_prev, self.cl_force_curr, self.tau).wait() + self.step = self.step + 1 + + def read_velocities(self): + cl.enqueue_copy(self.queue, self.np_atom_vel, self.cl_velocity).wait() + + def get_neighborhood_characteristics(self): + cl.enqueue_copy(self.queue, self.np_lj_count, self.cl_lj_count).wait() + return np.min(self.np_lj_count), np.max(self.np_lj_count), int(np.mean(self.np_lj_count)) + + def get_velocity_norms(self): + return np.linalg.norm(self.np_atom_vel[:,0:3], axis=1) + + def get_temperature(self): + velocity_norms = np.square(np.linalg.norm(self.np_atom_vel[:,0:3], axis=1)) # nm^2 ps^-2 + masses = self.np_atom_pos[:,3:4].flatten() # u + energy = np.sum(velocity_norms * masses * 0.5) # u nm^2 ps^-2 + kb = 0.00831446262 # Boltzmann's constant [kJ mol^-1 K^-1] + return energy * (2 / (3*self.n_atoms*kb)) # K + + def get_rms_velocity(self): + velocity_norms = (np.linalg.norm(self.np_atom_vel[:,0:3], axis=1)) # nm^2 ps^-2 + return np.sqrt(np.mean(velocity_norms) * 1e6) + + def randomize_velocities(self, temperature): + kb = 1.38064852e-23 # Boltzmann's constant [J/K] + masses = self.np_atom_pos[:,3:4].flatten() * 1.66053907e-27 # Mass in [kg] + self.np_atom_vel[:,0] = np.random.normal(0, np.sqrt(kb*temperature/masses), self.n_atoms) * 0.001 + self.np_atom_vel[:,1] = np.random.normal(0, np.sqrt(kb*temperature/masses), self.n_atoms) * 0.001 + self.np_atom_vel[:,2] = np.random.normal(0, np.sqrt(kb*temperature/masses), self.n_atoms) * 0.001 + cl.enqueue_copy(self.queue, self.cl_velocity, self.np_atom_vel).wait() + + def scale_velocities(self, scale): + self.np_atom_vel = self.np_atom_vel * scale; + cl.enqueue_copy(self.queue, self.cl_velocity, self.np_atom_vel).wait() + + def gl_draw_particles(self): + gl.glEnableClientState(gl.GL_VERTEX_ARRAY) + self.gl_position.bind() + gl.glVertexPointer(4, gl.GL_FLOAT, 0, self.gl_position) + gl.glDrawArrays(gl.GL_POINTS, 0, self.n_atoms) + gl.glDisableClientState(gl.GL_VERTEX_ARRAY) diff --git a/interacticle/visual/__init__.py b/interacticle/visual/__init__.py new file mode 100644 index 0000000..a9f9870 --- /dev/null +++ b/interacticle/visual/__init__.py @@ -0,0 +1,5 @@ +from .box import WireBox +from .view import View +from .links import MolecularLinks +from .plot import VelocityHistogram +from .plot import TemperaturePlot diff --git a/interacticle/visual/box.py b/interacticle/visual/box.py new file mode 100644 index 0000000..fb18d84 --- /dev/null +++ b/interacticle/visual/box.py @@ -0,0 +1,67 @@ +from OpenGL.GL import * + +class WireBox: + def __init__(self, size): + self.x0 = 0 + self.x1 = size + self.y0 = 0 + self.y1 = size + self.z0 = 0 + self.z1 = size + + def setup(self): + pass + + def update(self): + pass + + def shutdown(self): + pass + + def display_decoration(self, uniform): + glUniform3fv(uniform['color'], 1, [0.6,0.6,0.6]) + glBegin(GL_LINE_STRIP) + glVertex(self.x0, self.y0, self.z0) + glVertex(self.x0, self.y1, self.z0) + glVertex(self.x0, self.y1, self.z1) + glVertex(self.x0, self.y0, self.z1) + glVertex(self.x0, self.y0, self.z0) + glEnd() + glBegin(GL_LINE_STRIP) + glVertex(self.x1, self.y0, self.z0) + glVertex(self.x1, self.y1, self.z0) + glVertex(self.x1, self.y1, self.z1) + glVertex(self.x1, self.y0, self.z1) + glVertex(self.x1, self.y0, self.z0) + glEnd() + glBegin(GL_LINE_STRIP) + glVertex(self.x0, self.y0, self.z1) + glVertex(self.x1, self.y0, self.z1) + glVertex(self.x1, self.y1, self.z1) + glVertex(self.x0, self.y1, self.z1) + glVertex(self.x0, self.y0, self.z1) + glEnd() + glBegin(GL_LINE_STRIP) + glVertex(self.x0, self.y0, self.z0) + glVertex(self.x1, self.y0, self.z0) + glVertex(self.x1, self.y1, self.z0) + glVertex(self.x0, self.y1, self.z0) + glVertex(self.x0, self.y0, self.z0) + glEnd() + glBegin(GL_LINE_STRIP) + glVertex(self.x0, self.y0, self.z0) + glVertex(self.x1, self.y0, self.z0) + glVertex(self.x1, self.y0, self.z1) + glVertex(self.x0, self.y0, self.z1) + glVertex(self.x0, self.y0, self.z0) + glEnd() + glBegin(GL_LINE_STRIP) + glVertex(self.x0,self.y1,self.z0) + glVertex(self.x1,self.y1,self.z0) + glVertex(self.x1,self.y1,self.z1) + glVertex(self.x0,self.y1,self.z1) + glVertex(self.x0,self.y1,self.z0) + glEnd() + + def display_window(self, uniform): + pass diff --git a/interacticle/visual/camera.py b/interacticle/visual/camera.py new file mode 100644 index 0000000..b1e9650 --- /dev/null +++ b/interacticle/visual/camera.py @@ -0,0 +1,106 @@ +from OpenGL.GL import * +from OpenGL.GLUT import * + +from pyrr import matrix44, quaternion + +import numpy as np + +class Projection: + def __init__(self, distance): + self.distance = distance + self.ratio = 4./3. + self.x = 0 + self.y = 0 + self.update() + + def update(self): + projection = matrix44.create_perspective_projection(20.0, self.ratio, 0.1, 5000.0) + look = matrix44.create_look_at( + eye = [self.x, -self.distance, self.y], + target = [self.x, 0, self.y], + up = [ 0, 0,-1]) + + self.matrix = np.matmul(look, projection) + + def update_ratio(self, width, height, update_viewport = True): + if update_viewport: + glViewport(0,0,width,height) + + self.ratio = width/height + self.update() + + def update_distance(self, change): + self.distance += change + self.update() + + def shift(self, x, y): + self.x -= x + self.y -= y + self.update() + + def get(self): + return self.matrix + +class Rotation: + def __init__(self, shift, x = np.pi, y = np.pi): + self.matrix = matrix44.create_from_translation(shift), + self.rotation_x = quaternion.Quaternion() + self.update(x,y) + + def shift(self, x, y): + self.matrix = np.matmul( + self.matrix, + matrix44.create_from_translation([x,y,0]) + ) + self.inverse_matrix = np.linalg.inv(self.matrix) + + def update(self, x, y): + rotation_x = quaternion.Quaternion(quaternion.create_from_eulers([x,0,0])) + rotation_y = self.rotation_x.conjugate.cross( + quaternion.Quaternion(quaternion.create_from_eulers([0,y,0]))) + self.rotation_x = self.rotation_x.cross(rotation_x) + + self.matrix = np.matmul( + self.matrix, + matrix44.create_from_quaternion(rotation_y.cross(self.rotation_x)) + ) + self.inverse_matrix = np.linalg.inv(self.matrix) + + def get(self): + return self.matrix + + def get_inverse(self): + return self.inverse_matrix + +class MouseDragMonitor: + def __init__(self, button, callback): + self.button = button + self.active = False + self.callback = callback + + def on_mouse(self, button, state, x, y): + if button == self.button: + self.active = (state == GLUT_DOWN) + self.last_x = x + self.last_y = y + + def on_mouse_move(self, x, y): + if self.active: + delta_x = self.last_x - x + delta_y = y - self.last_y + self.last_x = x + self.last_y = y + self.callback(delta_x, delta_y) + +class MouseScrollMonitor: + def __init__(self, callback): + self.callback = callback + + def on_mouse(self, button, state, x, y): + if button == 3: + self.callback(-1.0) + elif button == 4: + self.callback(1.0) + + def on_mouse_move(self, x, y): + pass diff --git a/interacticle/visual/links.py b/interacticle/visual/links.py new file mode 100644 index 0000000..52e9d52 --- /dev/null +++ b/interacticle/visual/links.py @@ -0,0 +1,48 @@ +import numpy as np + +import pyopencl as cl +mf = cl.mem_flags + +import OpenGL.GL as gl +from OpenGL.arrays import vbo + +class MolecularLinks: + def __init__(self, simulation): + self.simulation = simulation + + def setup(self): + if self.simulation.n_bonds > 0: + self.np_view_link_indices = np.zeros((int(np.sum(self.simulation.np_bond_count)*2),1), dtype=np.uint32) + curr = 0 + for i in range(self.simulation.n_atoms): + for j in range(int(self.simulation.np_bond_count[i])): + self.np_view_link_indices[2*curr+0] = i + self.np_view_link_indices[2*curr+1] = self.simulation.np_bond_indices[self.simulation.max_bonds*i+j][0] + curr = curr + 1 + + self.cl_view_link_indices = cl.Buffer(self.simulation.context, mf.COPY_HOST_PTR, hostbuf=self.np_view_link_indices) + self.np_view_link_pos = np.zeros((self.np_view_link_indices.shape[0],4), dtype=np.float32) + self.gl_view_link_pos = vbo.VBO(data=self.np_view_link_pos, usage=gl.GL_DYNAMIC_DRAW, target=gl.GL_ARRAY_BUFFER) + self.gl_view_link_pos.bind() + self.cl_view_link_pos = cl.GLBuffer(self.simulation.context, mf.READ_WRITE, int(self.gl_view_link_pos)) + + def update(self): + pass + + def shutdown(self): + pass + + def display_decoration(self, uniform): + if self.simulation.n_bonds > 0: + cl.enqueue_acquire_gl_objects(self.simulation.queue, [self.cl_view_link_pos]) + self.simulation.program.update_view_links(self.simulation.queue, (int(self.np_view_link_indices.shape[0] / 2),), None, self.simulation.cl_position, self.cl_view_link_indices, self.cl_view_link_pos).wait() + + gl.glEnableClientState(gl.GL_VERTEX_ARRAY) + gl.glUniform3fv(uniform['color'], 1, [0.5,0.5,0.5]) + self.gl_view_link_pos.bind() + gl.glVertexPointer(4, gl.GL_FLOAT, 0, self.simulation.gl_position) + gl.glDrawArrays(gl.GL_LINES, 0, self.np_view_link_indices.shape[0]) + gl.glDisableClientState(gl.GL_VERTEX_ARRAY) + + def display_window(self, uniform): + pass diff --git a/interacticle/visual/plot.py b/interacticle/visual/plot.py new file mode 100644 index 0000000..60b4937 --- /dev/null +++ b/interacticle/visual/plot.py @@ -0,0 +1,248 @@ +from OpenGL.GL import * + +import numpy as np + +import scipy.stats as stats +import scipy.constants as const +from scipy.optimize import minimize + +import matplotlib +import matplotlib.pyplot as plt + +from concurrent.futures import Future, ProcessPoolExecutor + +def get_histogram(velocities): + upper = np.ceil(np.max(velocities) / 100) * 100 + maxwellian = stats.maxwell.fit(velocities) + + plt.style.use('dark_background') + + fig, ax = plt.subplots(1, figsize=(8,4)) + + #plt.ylim(0, 0.004) + plt.ylabel('Probability') + + plt.xlim(0, upper) + plt.xlabel('Velocity magnitude [m/s]') + + plt.hist(velocities, bins=50, density=True, alpha=0.5, label='Simulated velocities') + + xs = np.linspace(0, upper, 100) + ys = stats.maxwell.pdf(xs, *maxwellian) + plt.plot(xs, ys, label='Maxwell-Boltzmann distribution') + + plt.legend(loc='upper right') + + fig.canvas.draw() + buf = fig.canvas.tostring_rgb() + width, height = fig.canvas.get_width_height() + texture = np.frombuffer(buf, dtype=np.uint8).reshape(width, height, 3) + + plt.close(fig) + + return (texture, width, height) + + +class VelocityHistogram: + def __init__(self, gas, origin, extend): + self.gas = gas + self.origin = origin + self.extend = extend + self.steps = -1 + + self.pool = ProcessPoolExecutor(max_workers=1) + self.plotter = None + + self.tick = False + self.mixing = 0.0 + + self.vertices = np.array([ + self.origin[0] , self.origin[1] , 0., 1., + self.origin[0] + self.extend[0], self.origin[1] , 1., 1., + self.origin[0] , self.origin[1] + self.extend[1], 0., 0., + self.origin[0] + self.extend[0], self.origin[1] + self.extend[1], 1., 0. + ], dtype=np.float32) + + def setup(self): + self.vao = glGenVertexArrays(1) + self.vbo = glGenBuffers(1) + + glBindVertexArray(self.vao) + glBindBuffer(GL_ARRAY_BUFFER, self.vbo) + glBufferData(GL_ARRAY_BUFFER, self.vertices.tostring(), GL_STATIC_DRAW) + + glEnableVertexAttribArray(0) + glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(0)) + glEnableVertexAttribArray(1) + glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(2*np.dtype(np.float32).itemsize)) + + self.texture_id = glGenTextures(2) + + glBindTexture(GL_TEXTURE_2D, self.texture_id[0]) + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + + glBindTexture(GL_TEXTURE_2D, self.texture_id[1]) + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + + def shutdown(self): + self.pool.shutdown(wait=True) + + def update(self): + self.steps = self.steps + 1 + + if self.steps % 50 == 0 and self.plotter == None: + self.plotter = self.pool.submit(get_histogram, self.gas.get_velocity_norms() * 1000) + + else: + if self.plotter != None and self.plotter.done(): + texture, width, height = self.plotter.result() + if self.tick: + glBindTexture(GL_TEXTURE_2D, self.texture_id[0]) + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture) + self.tick = False + self.mixing = 1.0 + + else: + glBindTexture(GL_TEXTURE_2D, self.texture_id[1]) + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture) + self.tick = True + self.mixing = 0.0 + + self.plotter = None + + def display_decoration(self, uniform): + pass + + def display_window(self, uniform): + if self.tick: + self.mixing = min(self.mixing+0.05, 1.0); + else: + self.mixing = max(self.mixing-0.05, 0.0); + + glBindTextures(self.texture_id[0], 2, self.texture_id) + glUniform1iv(uniform['picture'], len(self.texture_id), self.texture_id) + glUniform1f(uniform['mixing'], self.mixing) + + glBindVertexArray(self.vao); + glDrawArrays(GL_TRIANGLE_STRIP, 0, 4) + glBindVertexArray(0) + +def get_plot(data): + xs = data[0] + ys = data[1] + + plt.style.use('dark_background') + + fig, ax = plt.subplots(1, figsize=(8,4)) + + plt.ylim(0, 350) + plt.ylabel('Temperature [K]') + + plt.xlim(0, 200) + plt.xlabel('Time [ps]') + + plt.plot(xs, ys, label='Temperature') + + plt.legend(loc='upper right') + + fig.canvas.draw() + buf = fig.canvas.tostring_rgb() + width, height = fig.canvas.get_width_height() + texture = np.frombuffer(buf, dtype=np.uint8).reshape(width, height, 3) + + plt.close(fig) + + return (texture, width, height) + +class TemperaturePlot: + def __init__(self, gas, origin, extend): + self.gas = gas + self.origin = origin + self.extend = extend + self.steps = -1 + + self.pool = ProcessPoolExecutor(max_workers=1) + self.plotter = None + + self.tick = False + self.mixing = 0.0 + + self.xs = [] + self.ys = [] + + self.vertices = np.array([ + self.origin[0] , self.origin[1] , 0., 1., + self.origin[0] + self.extend[0], self.origin[1] , 1., 1., + self.origin[0] , self.origin[1] + self.extend[1], 0., 0., + self.origin[0] + self.extend[0], self.origin[1] + self.extend[1], 1., 0. + ], dtype=np.float32) + + def setup(self): + self.vao = glGenVertexArrays(1) + self.vbo = glGenBuffers(1) + + glBindVertexArray(self.vao) + glBindBuffer(GL_ARRAY_BUFFER, self.vbo) + glBufferData(GL_ARRAY_BUFFER, self.vertices.tostring(), GL_STATIC_DRAW) + + glEnableVertexAttribArray(0) + glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(0)) + glEnableVertexAttribArray(1) + glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(2*np.dtype(np.float32).itemsize)) + + self.texture_id = glGenTextures(2) + + glBindTexture(GL_TEXTURE_2D, self.texture_id[0]) + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + + glBindTexture(GL_TEXTURE_2D, self.texture_id[1]) + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + + def shutdown(self): + self.pool.shutdown(wait=True) + + def update(self): + self.steps = self.steps + 1 + + if self.steps % 10 == 0 and self.plotter == None: + self.xs.append(self.gas.step * self.gas.tau) + self.ys.append(self.gas.get_temperature()) + self.plotter = self.pool.submit(get_plot, (self.xs, self.ys)) + + else: + if self.plotter != None and self.plotter.done(): + texture, width, height = self.plotter.result() + if self.tick: + glBindTexture(GL_TEXTURE_2D, self.texture_id[0]) + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture) + self.tick = False + self.mixing = 1.0 + + else: + glBindTexture(GL_TEXTURE_2D, self.texture_id[1]) + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture) + self.tick = True + self.mixing = 0.0 + + self.plotter = None + + def display_decoration(self, uniform): + pass + + def display_window(self, uniform): + if self.tick: + self.mixing = min(self.mixing+0.05, 1.0); + else: + self.mixing = max(self.mixing-0.05, 0.0); + + glBindTextures(self.texture_id[0], 2, self.texture_id) + glUniform1iv(uniform['picture'], len(self.texture_id), self.texture_id) + glUniform1f(uniform['mixing'], self.mixing) + + glBindVertexArray(self.vao); + glDrawArrays(GL_TRIANGLE_STRIP, 0, 4) + glBindVertexArray(0) diff --git a/interacticle/visual/shader.py b/interacticle/visual/shader.py new file mode 100644 index 0000000..c6c035b --- /dev/null +++ b/interacticle/visual/shader.py @@ -0,0 +1,28 @@ +from OpenGL.GL import GL_VERTEX_SHADER, GL_FRAGMENT_SHADER, GL_GEOMETRY_SHADER +from OpenGL.GL import shaders + +class Shader: + def __init__(self, fragment_src, vertex_src, uniform): + self.program = shaders.compileProgram( + shaders.compileShader(vertex_src, GL_VERTEX_SHADER), + shaders.compileShader(fragment_src, GL_FRAGMENT_SHADER)) + self.uniform = { } + for name in uniform: + self.uniform[name] = shaders.glGetUniformLocation(self.program, name) + + def use(self): + shaders.glUseProgram(self.program) + +class GeometryShader: + def __init__(self, fragment_src, geometry_src, vertex_src, uniform): + self.program = shaders.compileProgram( + shaders.compileShader(vertex_src, GL_VERTEX_SHADER), + shaders.compileShader(geometry_src, GL_GEOMETRY_SHADER), + shaders.compileShader(fragment_src, GL_FRAGMENT_SHADER)) + self.uniform = { } + for name in uniform: + self.uniform[name] = shaders.glGetUniformLocation(self.program, name) + + def use(self): + shaders.glUseProgram(self.program) + diff --git a/interacticle/visual/view.py b/interacticle/visual/view.py new file mode 100644 index 0000000..14ff487 --- /dev/null +++ b/interacticle/visual/view.py @@ -0,0 +1,242 @@ +from OpenGL.GL import * +from OpenGL.GLUT import * + +from pyrr import matrix44 + +import numpy as np + +from .shader import Shader, GeometryShader +from .camera import Projection, Rotation, MouseDragMonitor, MouseScrollMonitor + +particle_shader = ( + """ + #version 430 + + uniform vec4 camera_pos; + + in GS_OUT + { + vec3 particle_pos; + vec3 atom_color; + vec2 tex_coord; + } fs_in; + + void main(){ + if (length(fs_in.tex_coord - vec2(0.5)) > 0.5) { + discard; + } + + vec3 n = vec3(fs_in.tex_coord - vec2(0.5), 0.); + n.z = -sqrt(1 - length(n))*0.6; + n = normalize(n); + + vec3 dir = normalize(camera_pos.xyz - fs_in.particle_pos); + vec3 color = fs_in.atom_color * dot(dir, n); + + gl_FragColor = vec4(color.xyz, 1.0); +