Merge pull request #68 from maximillian-dolan/master

added different particle types exclusive interactions to md

examples/ideal_gas_law/ideal_gas_law.ipynb

@@ -450,7 +450,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.2"
+   "version": "3.10.12"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,

examples/simple_examples/molecular_dynamics.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -23,7 +23,8 @@
    "source": [
     "def md_simulation(number_of_particles, temperature, box_length, number_of_steps, sample_frequency):\n",
     "    # Initialise the system\n",
-    "    system = md.initialise(number_of_particles, temperature, box_length, 'square')\n",
+    "    constants=[[1.363e-134, 9.273e-78],[1.365e-130, 9.278e-77],[1.368e-130, 9.278e-77]]\n",
+    "    system = md.initialise(number_of_particles, temperature, box_length, 'square', constants=constants)\n",
     "    # This sets the sampling class\n",
     "    sample_system = sample.MaxBolt(system)# Start at time 0\n",
     "    system.time = 0\n",
@@ -63,7 +64,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "md_simulation(20, 1000, 100, 5000, 100)"
+    "system = md_simulation(20, 1000, 100, 5000, 25)"
    ]
   },
   {
@@ -90,7 +91,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.2"
+   "version": "3.12.4"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,

examples/simple_examples/monte_carlo.ipynb

@@ -22,7 +22,8 @@
    "source": [
     "def mc_simulation(number_of_particles, temperature, box_length, number_of_steps, sample_frequency):\n",
     "    # Initialise the system placing the particles on a square lattice\n",
-    "    system = mc.initialise(number_of_particles, temperature, box_length, 'square')\n",
+    "    constants=[[1.363e-134, 9.273e-78],[1.363e-134, 9.273e-78]]\n",
+    "    system = mc.initialise(number_of_particles, temperature, box_length, 'square', constants=constants, mixing = False)\n",
     "    # This sets the sampling class as Energy, which shows the energy of the system\n",
     "    sample_system = sample.Energy(system)\n",
     "    # Compute the energy of the system\n",
@@ -71,7 +72,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "system = mc_simulation(100, 273.15, 45, 5000, 25)"
+    "system = mc_simulation(50, 273.15, 45, 1000, 25)"
    ]
   },
   {
@@ -105,7 +106,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.2"
+   "version": "3.12.4"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,

pylj/forcefields.py

@@ -1,61 +1,6 @@
 import numpy as np
 
 
-class lennard_jones(object):
-    r"""Calculate the energy or force for a pair of particles using the
-    Lennard-Jones (A/B variant) forcefield.
-
-    Parameters
-    ----------
-    constants: float, array_like
-        An array of length two consisting of the A and B parameters for the
-        12-6 Lennard-Jones function
-
-    """
-    def __init__(self, constants):
-        self.a = constants[0]
-        self.b = constants[1]
-
-    def energy(self, dr ):
-        r"""Calculate the energy for a pair of particles using the
-        Lennard-Jones (A/B variant) forcefield.
-
-        .. math::
-            E = \frac{A}{dr^{12}} - \frac{B}{dr^6}
-
-        Attributes:
-        ----------
-        dr (float): The distance between particles.
-
-        Returns
-        -------
-        float: array_like
-        The potential energy between the particles.
-        """
-        self.energy = self.a * np.power(dr, -12) - (self.b * np.power(dr, -6))
-        return self.energy
-
-    def force(self, dr):
-        r"""Calculate the force for a pair of particles using the
-        Lennard-Jones (A/B variant) forcefield.
-
-        .. math::
-            f = \frac{12A}{dr^{13}} - \frac{6B}{dr^7}
-
-        Attributes:
-        ----------
-        dr (float): The distance between particles.
-        
-        Returns
-        -------
-        float: array_like
-        The force between the particles.
-        """
-        self.force = 12 * self.a * np.power(dr, -13) - (6 * self.b * np.power(dr, -7))
-        return self.force
-
-
-
 class lennard_jones_sigma_epsilon(object):
     r"""Calculate the energy or force for a pair of particles using the
     Lennard-Jones (sigma/epsilon variant) forcefield.
@@ -68,8 +13,12 @@ class lennard_jones_sigma_epsilon(object):
 
     """
     def __init__(self, constants):
+        if len(constants) != 2:
+            raise IndexError(f'There should be two constants per set, not {len(constants)}')
+
         self.sigma = constants[0]
         self.epsilon = constants[1]
+        self.point_size = 1.3e10 * (self.sigma*(2**(1/6)))
 
     def energy(self, dr):
         r"""Calculate the energy for a pair of particles using the
@@ -110,8 +59,70 @@ def force(self, dr):
         self.force = 48 * self.epsilon * np.power(self.sigma, 12) * np.power(
             dr, -13) - (24 * self.epsilon * np.power(self.sigma, 6) * np.power(dr, -7))
         return self.force
+
+    def mixing(self, constants_2):
+        r""" Calculates mixing for two sets of constants
+        
+        ..math::
+            \sigma_{12} = \frac{\sigma_1 + \sigma_2}{2}
+            \epsilon{12} = \sqrt{\epsilon_1 * \epsilon_2}
+        
+        Parameters:
+        ----------
+        constants_2: float, array_like
+            The second set of constantss
+        """
+        sigma2 = constants_2[0]
+        epsilon2 = constants_2[1]
+        self.sigma = (self.sigma+sigma2)/2
+        self.epsilon = np.sqrt(self.epsilon * epsilon2)
 
 
+class lennard_jones(lennard_jones_sigma_epsilon):
+    r"""Converts a/b variant values to sigma/epsilon variant
+    then maps to lennard_jones_sigma_epsilon class
+
+    ..math::
+        \sigma = \frac{a}{b}^(\frac{1}{6})
+        \sigma = \frace{b^2}{4*a}
+
+    Parameters
+    ----------
+    constants: float, array_like
+        An array of length two consisting of the A and B
+        parameters for the 12-6 Lennard-Jones function
+    """   
+    def __init__(self, constants):
+        if len(constants) != 2:
+            raise IndexError(f'There should be two constants per set, not {len(constants)}')
+        self.a = constants[0]
+        self.b = constants[1]
+        sigma = (self.a / self.b)**(1/6)
+        epsilon = (self.b**2)/(4*self.a)
+        super().__init__([sigma, epsilon])
+
+    def mixing(self, constants_2):
+        r"""Converts second set of a/b constants into sigma/epsilon
+        for use in mixing method. Then converts changed self sigma/epsilon
+        values back to a/b
+
+        ..math::
+            a = 4*\epsilon*(\sigma^12)
+            b = 4*\epsilon*(\sigma^6)
+
+        Parameters:
+        ----------
+        constants_2: float, array_like
+            The second set of constantss
+        """
+        a2 = constants_2[0]
+        b2 = constants_2[1]
+        sigma2 = (a2 / b2)**(1/6)
+        epsilon2 = (b2**2)/(4*a2)
+        super().mixing([sigma2,epsilon2])
+        self.a = 4 * self.epsilon * (self.sigma**12)
+        self.b = 4 * self.epsilon * (self.sigma**6)
+
 
 class buckingham(object):
     r""" Calculate the energy or force for a pair of particles using the
@@ -125,9 +136,12 @@ class buckingham(object):
 
     """
     def __init__(self, constants):
+        if len(constants) != 3:
+            raise IndexError(f'There should be three constants per set, not {len(constants)}')
         self.a = constants[0]
         self.b = constants[1]
         self.c = constants[2]
+        self.point_size = 8 # Needs better solution relevant to constants
 
     def energy(self, dr):
         r"""Calculate the energy for a pair of particles using the
@@ -145,7 +159,13 @@ def energy(self, dr):
         float: array_like
         The potential energy between the particles.
         """
-        self.energy = self.a * np.exp(- np.multiply(self.b, dr)) - self.c / np.power(dr, 6)
+        energy = self.a * np.exp(- np.multiply(self.b, dr)) - self.c / np.power(dr, 6)
+        # Cut out infinite values where r = 0
+        if type(dr) != float:
+            energy = np.array(energy)
+            energy[np.where(energy > 10e300)] = 0
+            energy[np.where(energy < -10e300)] = 0
+        self.energy = energy
         return self.energy
 
     def force(self, dr):
@@ -164,9 +184,31 @@ def force(self, dr):
         float: array_like
         The force between the particles.
         """
-        self.force = self.a * self.b * np.exp(- np.multiply(self.b, dr)) - 6 * self.c / np.power(dr, 7)
+        force = self.a * self.b * np.exp(- np.multiply(self.b, dr)) - 6 * self.c / np.power(dr, 7)
+        # Cut out infinite values where r = 0
+        if type(dr) != float:
+            force = np.array(force)
+            force[np.where(force > 10e300)] = 0
+            force[np.where(force < -10e300)] = 0
+        self.force = force
         return self.force
 
+    def mixing(self, constants2):
+        r""" Calculates mixing for two sets of constants
+        
+        ..math::
+            a_{12} = \sqrt{a_1 * a_2}
+            b_{12} = \sqrt{b_1 * b_2}
+            c_{12} = \sqrt{c_1 * c_2}
+        
+        Attributes:
+        ----------
+        constants_2: float, array_like
+            The second set of constantss
+        """
+        self.a = np.sqrt(self.a*constants2[0])
+        self.b = np.sqrt(self.b*constants2[1])
+        self.c = np.sqrt(self.c*constants2[2])
 
 
 class square_well(object):
@@ -183,11 +225,13 @@ class square_well(object):
         
     '''
     def __init__(self, constants, max_val=np.inf):
-
+        if len(constants) != 3:
+            raise IndexError(f'There should be three constants per set, not {len(constants)}')
         self.epsilon = constants[0]
         self.sigma = constants[1]
         self.lamda = constants[2] #Spelling as lamda not lambda to avoid calling python lambda function
         self.max_val = max_val
+        self.point_size = 10
 
     def energy(self, dr):
         r'''Calculate the energy for a pair of particles using a

pylj/mc.py

@@ -8,8 +8,8 @@ def initialise(
     box_length,
     init_conf,
     mass=39.948,
-    constants=[1.363e-134, 9.273e-78],
-    forcefield=ff.lennard_jones,
+    constants=[[1.363e-134, 9.273e-78]],
+    forcefield=ff.lennard_jones
 ):
     """Initialise the particle positions (this can be either as a square or
     random arrangement), velocities (based on the temperature defined, and #
@@ -48,7 +48,7 @@ def initialise(
         constants,
         forcefield,
         mass,
-        init_conf=init_conf,
+        init_conf=init_conf
     )
     system.particles["xvelocity"] = 0
     system.particles["yvelocity"] = 0

pylj/md.py

@@ -10,8 +10,8 @@ def initialise(
     init_conf,
     timestep_length=1e-14,
     mass=39.948,
-    constants=[1.363e-134, 9.273e-78],
-    forcefield=ff.lennard_jones,
+    constants=[[1.363e-134, 9.273e-78]],
+    forcefield=ff.lennard_jones
 ):
     """Initialise the particle positions (this can be either as a square or
     random arrangement), velocities (based on the temperature defined, and
@@ -53,7 +53,7 @@ def initialise(
         forcefield,
         mass,
         init_conf=init_conf,
-        timestep_length=timestep_length,
+        timestep_length=timestep_length
     )
     v = np.random.rand(system.particles.size, 2, 12)
     v = np.sum(v, axis=2) - 6.0
@@ -153,6 +153,7 @@ def sample(particles, box_length, initial_particles, system):
         system.cut_off,
         system.constants,
         system.forcefield,
+        system.mass
     )
     msd_new = calculate_msd(particles, initial_particles, box_length)
     system.pressure_sample = np.append(system.pressure_sample, pressure_new)
@@ -329,37 +330,6 @@ def compute_force(particles, box_length, cut_off, constants, forcefield, mass):
     return part, dist, forces, energies
 
 
-def compute_energy(particles, box_length, cut_off, constants, forcefield):
-    r"""Calculates the total energy of the simulation.
-    Parameters
-    ----------
-    particles: util.particle_dt, array_like
-        Information about the particles.
-    box_length: float
-        Length of a single dimension of the simulation square, in Angstrom.
-    cut_off: float
-        The distance greater than which the energies between particles is
-        taken as zero.
-    constants: float, array_like (optional)
-        The constants associated with the particular forcefield used, e.g. for
-        the function forcefields.lennard_jones, theses are [A, B]
-    forcefield: function (optional)
-        The particular forcefield to be used to find the energy and forces.
-    Returns
-    -------
-    util.particle_dt, array_like
-        Information about particles, with updated accelerations and forces.
-    float, array_like
-        Current distances between pairs of particles in the simulation.
-    float, array_like
-        Current energies between pairs of particles in the simulation.
-    """
-    dist, energies = heavy.compute_energy(
-        particles, box_length, cut_off, constants, forcefield
-    )
-    return dist, energies
-
-
 def heat_bath(particles, temperature_sample, bath_temperature):
     r"""Rescales the velocities of the particles in the system to control the
     temperature of the simulation. Thereby allowing for an NVT ensemble. The