adjusting energy/force methods in classes

pylj/forcefields.py

@@ -8,47 +8,51 @@ class lennard_jones(object):
 
     Parameters
     ----------
-    dr: float, array_like
-        The distances between the all pairs of particles.
     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, dr, constants):
-        self.dr = dr
+    def __init__(self, constants):
         self.a = constants[0]
         self.b = constants[1]
-        self.energy = self.update_energy()
-        self.force = self.update_force()
 
-    def update_energy(self):
+    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(self.dr, -12) - (self.b * np.power(self.dr, -6))
+        self.energy = self.a * np.power(dr, -12) - (self.b * np.power(dr, -6))
         return self.energy
 
-    def update_force(self):
+    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(self.dr, -13) - (6 * self.b * np.power(self.dr, -7))
+        self.force = 12 * self.a * np.power(dr, -13) - (6 * self.b * np.power(dr, -7))
         return self.force
 
 
@@ -59,49 +63,53 @@ class lennard_jones_sigma_epsilon(object):
 
     Parameters
     ----------
-    dr: float, array_like
-        The distances between the all pairs of particles.
     constants: float, array_like
         An array of length two consisting of the sigma (a) and epsilon (e)
         parameters for the 12-6 Lennard-Jones function
+
     """
-    def __init__(self, dr, constants):
-        self.dr = dr
+    def __init__(self, constants):
         self.sigma = constants[0]
         self.epsilon = constants[1]
-        self.energy = self.update_energy()
-        self.force = self.update_force()
 
-    def update_energy(self):
+    def energy(self, dr):
         r"""Calculate the energy for a pair of particles using the
         Lennard-Jones (sigma/epsilon variant) forcefield.
 
         .. math::
             E = \frac{4e*a^{12}}{dr^{12}} - \frac{4e*a^{6}}{dr^6}
 
+        Attributes:
+        ----------
+        dr (float): The distance between particles.
+ 
         Returns
         -------
         float: array_like
         The potential energy between the particles.
         """
-        self.energy = 4 * self.epsilon * np.power(self.sigma, 12) * np.power(self.dr, -12) - (
-                        4 * self.epsilon * np.power(self.sigma, 6) * np.power(self.dr, -6)) 
+        self.energy = 4 * self.epsilon * np.power(self.sigma, 12) * np.power(dr, -12) - (
+                        4 * self.epsilon * np.power(self.sigma, 6) * np.power(dr, -6)) 
         return self.energy 
 
-    def update_force(self):
+    def force(self, dr):
         r"""Calculate the force for a pair of particles using the
         Lennard-Jones (sigma/epsilon variant) forcefield.
 
         .. math::
             f = \frac{48e*a^{12}}{dr^{13}} - \frac{24e*a^{6}}{dr^7}
 
+        Attributes:
+        ----------
+        dr (float): The distance between particles.
+ 
         Returns
         -------
         float: array_like
         The force between the particles.
         """
         self.force = 48 * self.epsilon * np.power(self.sigma, 12) * np.power(
-            self.dr, -13) - (24 * self.epsilon * np.power(self.sigma, 6) * np.power(self.dr, -7))
+            dr, -13) - (24 * self.epsilon * np.power(self.sigma, 6) * np.power(dr, -7))
         return self.force
 
 
@@ -112,48 +120,52 @@ class buckingham(object):
 
     Parameters
     ----------
-    dr: float, array_like
-        The distances between all the pairs of particles.
     constants: float, array_like
         An array of length three consisting of the A, B and C parameters for
         the Buckingham function.
+
     """
-    def __init__(self, dr, constants):
-        self.dr = dr
+    def __init__(self, constants):
         self.a = constants[0]
         self.b = constants[1]
         self.c = constants[2]
-        self.energy = self.update_energy()
-        self.force = self.update_force()
 
-    def update_energy(self):
+    def energy(self, dr):
         r"""Calculate the energy for a pair of particles using the
         Buckingham forcefield.
 
         .. math::
             E = Ae^{(-Bdr)} - \frac{C}{dr^6}
 
+        Attributes:
+        ----------
+        dr (float): The distance between particles.
+ 
         Returns
         -------
         float: array_like
         The potential energy between the particles.
         """
-        self.energy = self.a * np.exp(- np.multiply(self.b, self.dr)) - self.c / np.power(self.dr, 6)
+        self.energy = self.a * np.exp(- np.multiply(self.b, dr)) - self.c / np.power(dr, 6)
         return self.energy
 
-    def update_force(self):
+    def force(self, dr):
         r"""Calculate the force for a pair of particles using the
         Buckingham forcefield.
 
         .. math::
             f = ABe^{(-Bdr)} - \frac{6C}{dr^7}
 
+        Attributes:
+        ----------
+        dr (float): The distance between particles.
+ 
         Returns
         -------
         float: array_like
         The force between the particles.
         """
-        self.force = self.a * self.b * np.exp(- np.multiply(self.b, self.dr)) - 6 * self.c / np.power(self.dr, 7)
+        self.force = self.a * self.b * np.exp(- np.multiply(self.b, dr)) - 6 * self.c / np.power(dr, 7)
         return self.force
 
 
@@ -164,31 +176,21 @@ class square_well(object):
 
     Parameters
     ----------
-    dr: float, array_like
-        The distances between all the pairs of particles.
     constants: float, array_like
         An array of length three consisting of the epsilon, sigma, and lambda
         parameters for the square well model.
     max_val: int (optional)
         Upper bound for values in square well - replaces usual infinite values
+        
     '''
-    def __init__(self, dr, constants, max_val=np.inf):
+    def __init__(self, constants, max_val=np.inf):
 
-        if not isinstance(dr, np.ndarray):
-            if isinstance(dr, list):
-                dr = np.array(dr, dtype='float')
-            elif isinstance(dr, float):
-                dr = np.array([dr], dtype='float')
-
-        self.dr = dr
         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.energy = self.update_energy()
-        self.force = 0
 
-    def update_energy(self):
+    def energy(self, dr):
         r'''Calculate the energy for a pair of particles using a
         square well model.
 
@@ -202,19 +204,30 @@ def update_energy(self):
                 E = 0
             }
 
+        Attributes:
+        ----------
+        dr (float): The distance between particles.
+ 
         Returns
         -------
         float: array_like
             The potential energy between the particles.
         '''
-        E = np.zeros_like(self.dr)
-        E = np.zeros_like(self.dr)
-        E[np.where(self.dr < self.epsilon)] = self.max_val
-        E[np.where(self.dr >= self.lamda * self.sigma)] = 0
-
-        # apply mask for sigma <= self.dr < lambda * sigma
-        a = self.sigma <= self.dr
-        b = self.dr < self.lamda * self.sigma
+
+        if not isinstance(dr, np.ndarray):
+            if isinstance(dr, list):
+                dr = np.array(dr, dtype='float')
+            elif isinstance(dr, float):
+                dr = np.array([dr], dtype='float')
+
+        E = np.zeros_like(dr)
+        E = np.zeros_like(dr)
+        E[np.where(dr < self.epsilon)] = self.max_val
+        E[np.where(dr >= self.lamda * self.sigma)] = 0
+
+        # apply mask for sigma <= dr < lambda * sigma
+        a = self.sigma <= dr
+        b = dr < self.lamda * self.sigma
         E[np.where(a & b)] = -self.epsilon
 
         if len(E) == 1:
@@ -224,7 +237,7 @@ def update_energy(self):
 
         return self.energy
 
-    def update_force(self):
+    def force(self):
         r'''The force of a pair of particles using a square well model is given by:
 
         .. math::

pylj/pairwise.py

@@ -48,8 +48,9 @@ def compute_force(particles, box_length, cut_off, constants, forcefield, mass):
     distances, dx, dy = heavy.dist(
         particles["xposition"], particles["yposition"], box_length
     )
-    forces = forcefield(distances, constants).force
-    energies = forcefield(distances, constants).energy
+    ff = forcefield(constants)
+    forces = ff.force(distances)
+    energies = ff.energy(distances)
     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)
@@ -211,7 +212,8 @@ def compute_energy(particles, box_length, cut_off, constants, forcefield):
     distances, dx, dy = heavy.dist(
         particles["xposition"], particles["yposition"], box_length
     )
-    energies = forcefield(distances, constants).energy
+    ff = forcefield(constants)
+    energies = ff.energy(distances)
     energies[np.where(distances > cut_off)] = 0.0
     return distances, energies
 
@@ -248,7 +250,8 @@ def calculate_pressure(
     distances, dx, dy = heavy.dist(
         particles["xposition"], particles["yposition"], box_length
     )
-    forces = forcefield(distances, constants).force
+    ff = forcefield(constants)
+    forces = ff.force(distances)
     forces[np.where(distances > cut_off)] = 0.0
     pres = np.sum(forces * distances)
     boltzmann_constant = 1.3806e-23  # joules / kelvin