summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAdrian Kummerlaender2021-03-27 22:35:43 +0100
committerAdrian Kummerlaender2021-03-27 22:35:43 +0100
commit75ef8db0d058158c10951a4184186f8c6cc27acc (patch)
tree03a33785b65472fc5eff9cc5d2c03379747b83d4
downloadinteracticle-75ef8db0d058158c10951a4184186f8c6cc27acc.tar
interacticle-75ef8db0d058158c10951a4184186f8c6cc27acc.tar.gz
interacticle-75ef8db0d058158c10951a4184186f8c6cc27acc.tar.bz2
interacticle-75ef8db0d058158c10951a4184186f8c6cc27acc.tar.lz
interacticle-75ef8db0d058158c10951a4184186f8c6cc27acc.tar.xz
interacticle-75ef8db0d058158c10951a4184186f8c6cc27acc.tar.zst
interacticle-75ef8db0d058158c10951a4184186f8c6cc27acc.zip
Initial public commit of this basic MD code
Simulation of interacting particles
-rw-r--r--benchmark.py46
-rw-r--r--examples/argon.py39
-rw-r--r--examples/dimer_only.py36
-rw-r--r--examples/water.py33
-rw-r--r--interacticle/__init__.py4
-rw-r--r--interacticle/kernel.cl460
-rw-r--r--interacticle/simulation.py394
-rw-r--r--interacticle/visual/__init__.py5
-rw-r--r--interacticle/visual/box.py67
-rw-r--r--interacticle/visual/camera.py106
-rw-r--r--interacticle/visual/links.py48
-rw-r--r--interacticle/visual/plot.py248
-rw-r--r--interacticle/visual/shader.py28
-rw-r--r--interacticle/visual/view.py242
-rw-r--r--interacticle/visualizer.py99
-rw-r--r--library.py84
-rw-r--r--shell.nix46
17 files changed, 1985 insertions, 0 deletions
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())
+