Initial public commit of this basic MD code

Simulation of interacting particles

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