added different particle types exclusive interactions to md

pylj/md.py

@@ -10,8 +10,9 @@ def initialise(
     init_conf,
     timestep_length=1e-14,
     mass=39.948,
-    constants=[1.363e-134, 9.273e-78],
+    constants=[[1.363e-134, 9.273e-78]],
     forcefield=ff.lennard_jones,
+    mixing = False
 ):
     """Initialise the particle positions (this can be either as a square or
     random arrangement), velocities (based on the temperature defined, and
@@ -54,6 +55,7 @@ def initialise(
         mass,
         init_conf=init_conf,
         timestep_length=timestep_length,
+        mixing = mixing
     )
     v = np.random.rand(system.particles.size, 2, 12)
     v = np.sum(v, axis=2) - 6.0
@@ -78,7 +80,7 @@ def initialize(
 
 #Jit tag here had to be removed
 def velocity_verlet(
-    particles, timestep_length, box_length, cut_off, constants, forcefield, mass
+    particles, timestep_length, box_length, cut_off, constants, forcefield, mass, type_identifiers
 ):
     """Uses the Velocity-Verlet integrator to move forward in time. The
     Updates the particles positions and velocities in terms of the Velocity
@@ -112,7 +114,7 @@ def velocity_verlet(
     xacceleration_store = list(particles["xacceleration"])
     yacceleration_store = list(particles["yacceleration"])
     particles, distances, forces, energies = heavy.compute_force(
-        particles, box_length, cut_off, constants, forcefield, mass
+        particles, box_length, cut_off, constants, forcefield, mass, type_identifiers
     )
     [particles["xvelocity"], particles["yvelocity"]] = update_velocities(
         [particles["xvelocity"], particles["yvelocity"]],
@@ -153,6 +155,7 @@ def sample(particles, box_length, initial_particles, system):
         system.cut_off,
         system.constants,
         system.forcefield,
+        system.type_identifiers
     )
     msd_new = calculate_msd(particles, initial_particles, box_length)
     system.pressure_sample = np.append(system.pressure_sample, pressure_new)

pylj/pairwise.py

@@ -3,7 +3,7 @@
 from pylj import pairwise as heavy
 
 #Jit tag here had to be removed
-def compute_force(particles, box_length, cut_off, constants, forcefield, mass):
+def compute_force(particles, box_length, cut_off, constants, forcefield, mass, type_identifiers):
     r"""Calculates the forces and therefore the accelerations on each of the
     particles in the simulation.
     Parameters
@@ -39,17 +39,23 @@ def compute_force(particles, box_length, cut_off, constants, forcefield, mass):
         (particles["xacceleration"].size - 1) * particles["xacceleration"].size / 2
     )
     forces = np.zeros(pairs)
-    distances = np.zeros(pairs)
     energies = np.zeros(pairs)
     atomic_mass_unit = 1.660539e-27  # kilograms
     mass_amu = mass  # amu
     mass_kg = mass_amu * atomic_mass_unit  # kilograms
     distances, dx, dy = heavy.dist(
         particles["xposition"], particles["yposition"], box_length
     )
-    ff = forcefield(constants)
-    forces = ff.force(distances)
-    energies = ff.energy(distances)
+    for type, constants in enumerate(constants):
+        identifier = type_identifiers[type]
+        ff = forcefield(constants)
+        type_distances = distances * create_dist_identifiers(identifier)
+        type_forces = ff.force(type_distances)
+        type_energies = ff.energy(distances)
+        type_forces = np.nan_to_num(type_forces)
+        type_energies = np.nan_to_num(type_energies)
+        forces+=type_forces
+        energies+=type_energies
     forces[np.where(distances > cut_off)] = 0.0
     energies[np.where(distances > cut_off)] = 0.0
     particles = update_accelerations(particles, forces, mass_kg, dx, dy, distances)
@@ -95,11 +101,12 @@ def update_accelerations(particles, f, m, dx, dy, dr):
     k = 0
     for i in range(0, particles.size - 1):
         for j in range(i + 1, particles.size):
-            particles["xacceleration"][i] += second_law(f[k], m, dx[k], dr[k])
-            particles["yacceleration"][i] += second_law(f[k], m, dy[k], dr[k])
-            particles["xacceleration"][j] -= second_law(f[k], m, dx[k], dr[k])
-            particles["yacceleration"][j] -= second_law(f[k], m, dy[k], dr[k])
+            particles["xacceleration"][i] += second_law(f[k], m, dx[k], dr[k]) if f[k]!=0 else 0
+            particles["yacceleration"][i] += second_law(f[k], m, dy[k], dr[k]) if f[k]!=0 else 0
+            particles["xacceleration"][j] -= second_law(f[k], m, dx[k], dr[k]) if f[k]!=0 else 0
+            particles["yacceleration"][j] -= second_law(f[k], m, dy[k], dr[k]) if f[k]!=0 else 0
             k += 1
+
     return particles
 
 
@@ -171,7 +178,7 @@ def lennard_jones_force(A, B, dr):
     float:
         The force between the two particles.
     """
-    print(
+    xacceleration(
         "pairwise.lennard_jones_energy has been deprecated, please use "
         "forcefields.lennard_jones with force=True instead"
     )
@@ -218,7 +225,7 @@ def compute_energy(particles, box_length, cut_off, constants, forcefield):
 
 
 def calculate_pressure(
-    particles, box_length, temperature, cut_off, constants, forcefield
+    particles, box_length, temperature, cut_off, constants, forcefield, type_identifiers
 ):
     r"""Calculates the instantaneous pressure of the simulation cell, found
     with the following relationship:
@@ -248,9 +255,19 @@ def calculate_pressure(
     """
     distances, dx, dy = heavy.dist(
         particles["xposition"], particles["yposition"], box_length
-    )
-    ff = forcefield(constants)
-    forces = ff.force(distances)
+    ) 
+    forces = np.zeros(len(distances))
+    energies = np.zeros(len(distances))
+    for type, constants in enumerate(constants):
+        identifier = type_identifiers[type]
+        ff = forcefield(constants)
+        type_distances = distances * create_dist_identifiers(identifier)
+        type_forces = ff.force(type_distances)
+        type_energies = ff.energy(distances)
+        type_forces = np.nan_to_num(type_forces)
+        type_energies = np.nan_to_num(type_energies)
+        forces+=type_forces
+        energies+=type_energies
     forces[np.where(distances > cut_off)] = 0.0
     pres = np.sum(forces * distances)
     boltzmann_constant = 1.3806e-23  # joules / kelvin
@@ -345,3 +362,28 @@ def pbc_correction(position, cell):
     if np.abs(position) > 0.5 * cell:
         position *= 1 - cell / np.abs(position)
     return position
+
+def create_dist_identifiers(type_identifier):
+    '''
+    Creates correct distance identifier matrix for particular type
+    of particle
+    Parameters
+    -------
+    type identifiers:
+        the identifier array listing 1 for particles of that type
+        or 0 for particles of a different type
+    Returns
+    -------
+    distances: float, array_like 
+        the distance identifier for interactions between each particle 
+        of that type, or 0 for interactions involving particles of a
+        different type
+
+    '''
+    distance_type_identifier = np.array([])
+    for index in range(len(type_identifier)):
+        if type_identifier[index]:
+            distance_type_identifier = np.append(distance_type_identifier,type_identifier[index+1:])
+        else:
+            distance_type_identifier = np.append(distance_type_identifier,np.zeros(len(type_identifier[index+1:])))
+    return distance_type_identifier

pylj/sample.py

@@ -494,12 +494,13 @@ def setup_cellview(ax, system, scale=1):  # pragma: no cover
     scale: float (optional)
         The amount by which the particle size should be scaled down.
     """
-    xpos = system.particles["xposition"]
-    ypos = system.particles["yposition"]
     # These numbers are chosen to scale the size of the particles nicely to
     # allow the particles to appear to interact appropriately
     mk = 6.00555e-8 / (system.box_length - 2.2727e-10) - 1e-10 / scale
-    ax.plot(xpos, ypos, "o", markersize=mk, markeredgecolor="black", color="#34a5daff")
+    for type in system.type_identifiers:
+        xpos = system.particles["xposition"] * type
+        ypos = system.particles["yposition"] * type   
+        ax.plot(xpos, ypos, "o", markersize=mk, markeredgecolor="black")
     ax.set_xlim([0, system.box_length])
     ax.set_ylim([0, system.box_length])
     ax.set_xticks([])
@@ -612,11 +613,12 @@ def update_cellview(ax, system):  # pragma: no cover
     system: System
         The whole system information.
     """
-    x3 = system.particles["xposition"]
-    y3 = system.particles["yposition"]
-    line = ax.lines[0]
-    line.set_ydata(y3)
-    line.set_xdata(x3)
+    for index, type in enumerate(system.type_identifiers):
+        x3 = system.particles["xposition"] * type
+        y3 = system.particles["yposition"] * type
+        line = ax.lines[index]
+        line.set_ydata(y3)
+        line.set_xdata(x3)
 
 
 def update_rdfview(ax, system, average_rdf, r):  # pragma: no cover

pylj/tests/test_md.py

@@ -15,7 +15,7 @@ def test_initialise_square(self):
     def test_velocity_verlet(self):
         a = md.initialise(2, 300, 8, "square")
         a.particles, a.distances, a.forces, a.energies = md.velocity_verlet(
-            a.particles, 1, a.box_length, a.cut_off, a.constants, a.forcefield, a.mass
+            a.particles, 1, a.box_length, a.cut_off, a.constants, a.forcefield, a.mass, a.type_identifiers
         )
         assert_almost_equal(a.particles["xprevious_position"] * 1e10, [2, 2])
         assert_almost_equal(a.particles["yprevious_position"] * 1e10, [2, 6])

pylj/tests/test_pairwise.py

@@ -34,9 +34,10 @@ def test_compute_forces(self):
             particles,
             30,
             15,
-            constants=[1.363e-134, 9.273e-78],
+            constants=[[1.363e-134, 9.273e-78]],
             forcefield=ff.lennard_jones,
             mass=39.948,
+            type_identifiers = [[1,1]]
         )
         assert_almost_equal(distances, [4e-10])
         assert_almost_equal(energies, [-1.4515047e-21])
@@ -70,8 +71,9 @@ def test_calculate_pressure(self):
             30,
             300,
             15,
-            constants=[1.363e-134, 9.273e-78],
+            constants=[[1.363e-134, 9.273e-78]],
             forcefield=ff.lennard_jones,
+            type_identifiers = [[1,1]]
         )
         assert_almost_equal(p * 1e24, 7.07368869)
 

pylj/util.py

@@ -46,10 +46,23 @@ def __init__(
         init_conf="square",
         timestep_length=1e-14,
         cut_off=15,
+        mixing = False
     ):
         self.number_of_particles = number_of_particles
         self.init_temp = temperature
-        self.constants = constants
+        self.constants = constants #if mixing else constants[0]
+
+        # Creates arrays to identify which particle is in which type
+        number_of_types = len(constants)
+        type_identifiers = np.zeros((number_of_types,number_of_particles))
+        i = 0
+        while i < number_of_particles:
+            for k in range(number_of_types):
+                if i < number_of_particles:
+                    type_identifiers[k][i] = 1
+                    i+=1
+        self.type_identifiers = type_identifiers
+
         self.forcefield = forcefield
         self.mass = mass
         if box_length <= 600:
@@ -180,6 +193,7 @@ def integrate(self, method):
             self.constants,
             self.forcefield,
             self.mass,
+            self.type_identifiers
         )
 
     def md_sample(self):