diff --git a/pylj/forcefields.py b/pylj/forcefields.py index 1221906..a644ea0 100644 --- a/pylj/forcefields.py +++ b/pylj/forcefields.py @@ -2,8 +2,7 @@ from numba import njit -#Jit tag here had to be removed -def lennard_jones(dr, constants, force=False): +class lennard_jones(object): r"""Calculate the energy or force for a pair of particles using the Lennard-Jones (A/B variant) forcefield. @@ -28,18 +27,22 @@ def lennard_jones(dr, constants, force=False): float: array_like The potential energy or force between the particles. """ - - if force: - return 12 * constants[0] * np.power(dr, -13) - ( - 6 * constants[1] * np.power(dr, -7)) - - else: - return constants[0] * np.power(dr, -12) - ( - constants[1] * np.power(dr, -6)) - - -#Jit tag here had to be removed -def lennard_jones_sigma_epsilon(dr, constants, force=False): + def __init__(self, dr, constants): + self.dr = dr + self.a = constants[0] + self.b = constants[1] + self.energy = self.update_energy() + self.force = self.update_force() + + def update_energy(self): + self.energy = self.a * np.power(self.dr, -12) - (self.b * np.power(self.dr, -6)) + return self.energy + + def update_force(self): + self.force = 12 * self.a * np.power(self.dr, -13) - (6 * self.b * np.power(self.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. @@ -64,18 +67,24 @@ def lennard_jones_sigma_epsilon(dr, constants, force=False): float: array_like The potential energy or force between the particles. """ - - if force: - return 48 * constants[1] * np.power(constants[0], 12) * np.power( - dr, -13) - (24 * constants[1] * np.power( - constants[0], 6) * np.power(dr, -7)) - else: - return 4 * constants[1] * np.power(constants[0], 12) * np.power(dr, -12) - ( - 4 * constants[1] * np.power(constants[0], 6) * np.power(dr, -6)) - - -#Jit tag here had to be removed -def buckingham(dr, constants, force=False): + def __init__(self, dr, constants): + self.dr = dr + self.sigma = constants[0] + self.epsilon = constants[1] + self.energy = self.update_energy() + self.force = self.update_force() + + def update_energy(self): + 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)) + return self.energy + + def update_force(self): + 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)) + return self.force + +class buckingham(object): r""" Calculate the energy or force for a pair of particles using the Buckingham forcefield. @@ -100,16 +109,23 @@ def buckingham(dr, constants, force=False): float: array_like The potential energy or force between the particles. """ - if force: - return constants[0] * constants[1] * np.exp( - - np.multiply(constants[1], dr)) - 6 * constants[2] / np.power( - dr, 7) - else: - return constants[0] * np.exp( - - np.multiply(constants[1], dr)) - constants[2] / np.power(dr, 6) - - -def square_well(dr, constants, max_val=np.inf, force=False): + def __init__(self, dr, constants): + self.dr = dr + 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): + self.energy = self.a * np.exp(- np.multiply(self.b, self.dr)) - self.c / np.power(self.dr, 6) + return self.energy + + def update_force(self): + self.force = self.a * self.b * np.exp(- np.multiply(self.b, self.dr)) - 6 * self.c / np.power(self.dr, 7) + return self.force + +class square_well(object): r'''Calculate the energy or force for a pair of particles using a square well model. @@ -146,26 +162,39 @@ def square_well(dr, constants, max_val=np.inf, force=False): float: array_like The potential energy between the particles. ''' - 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') - - if force: - raise ValueError("Force is infinite at sigma <= dr < lambda * sigma") - - else: - E = np.zeros_like(dr) - E[np.where(dr < constants[0])] = max_val - E[np.where(dr >= constants[2] * constants[1])] = 0 - - # apply mask for sigma <= dr < lambda * sigma - a = constants[1] <= dr - b = dr < constants[2] * constants[1] - E[np.where(a & b)] = -constants[0] - - if len(E) == 1: - return float(E[0]) - else: - return np.array(E, dtype='float') + def __init__(self, dr, 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): + 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 + E[np.where(a & b)] = -self.epsilon + + if len(E) == 1: + self.energy = float(E[0]) + else: + self.energy = np.array(E, dtype='float') + + return self.energy + + def update_force(self): + raise ValueError("Force is infinite at sigma <= dr < lambda * sigma") \ No newline at end of file diff --git a/pylj/pairwise.py b/pylj/pairwise.py index 8280a80..ce549f8 100644 --- a/pylj/pairwise.py +++ b/pylj/pairwise.py @@ -48,8 +48,8 @@ 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=True) - energies = forcefield(distances, constants, force=False) + forces = forcefield(distances, constants).force + energies = forcefield(distances, constants).energy 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 +211,7 @@ 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, force=False) + energies = forcefield(distances, constants).energy energies[np.where(distances > cut_off)] = 0.0 return distances, energies @@ -248,7 +248,7 @@ def calculate_pressure( distances, dx, dy = heavy.dist( particles["xposition"], particles["yposition"], box_length ) - forces = forcefield(distances, constants, force=True) + forces = forcefield(distances, constants).force forces[np.where(distances > cut_off)] = 0.0 pres = np.sum(forces * distances) boltzmann_constant = 1.3806e-23 # joules / kelvin diff --git a/pylj/tests/test_forcefields.py b/pylj/tests/test_forcefields.py index d343495..59bc629 100644 --- a/pylj/tests/test_forcefields.py +++ b/pylj/tests/test_forcefields.py @@ -6,93 +6,85 @@ class TestForcefields(unittest.TestCase): def test_lennard_jones_energy(self): a = forcefields.lennard_jones(2.0, [1.0, 1.0]) - assert_almost_equal(a, -0.015380859) + assert_almost_equal(a.energy, -0.015380859) b = forcefields.lennard_jones([2.0, 4.5], [1.0, 1.0]) - assert_almost_equal(b, [-0.015380859, -0.00012041]) + assert_almost_equal(b.energy, [-0.015380859, -0.00012041]) c = forcefields.lennard_jones([2.0, 4.5], [0.5, 3.5]) - assert_almost_equal(c, [-0.05456543, -0.00042149]) + assert_almost_equal(c.energy, [-0.05456543, -0.00042149]) d = forcefields.lennard_jones([1.0, 1.5, 20.0], [5.0, 3.5]) - assert_almost_equal( - d, [1.50000000, -0.268733500, -5.46874988e-08]) + assert_almost_equal(d.energy, [1.50000000, -0.268733500, -5.46874988e-08]) e = forcefields.lennard_jones([100.0, 200.0, 500.0], [100.0, 300.0]) - assert_almost_equal(e, [0, 0, 0]) + assert_almost_equal(e.energy, [0, 0, 0]) def test_lennard_jones_force(self): - a = forcefields.lennard_jones(2.0, [1.0, 1.0], force=True) - assert_almost_equal(a, -0.045410156) - b = forcefields.lennard_jones([2.0, 4.0, 6.0], [1.0, 1.0], force=True) - assert_almost_equal(b, [-0.045410156, -3.66032124e-04, -2.14325517e-05]) - c = forcefields.lennard_jones([2.0, 4.0, 6.0], [1.5, 4.0], force=True) - assert_almost_equal(c, [-0.185302734, -1.46457553e-03, -8.57325038e-05]) - d = forcefields.lennard_jones( - [150.0, 300.0, 500.0], [200.0, 500.0], force=True) - assert_almost_equal(d, [-1.7558299e-12, -1.3717421e-14, -3.8400000e-16]) + a = forcefields.lennard_jones(2.0, [1.0, 1.0]) + assert_almost_equal(a.force, -0.045410156) + b = forcefields.lennard_jones([2.0, 4.0, 6.0], [1.0, 1.0]) + assert_almost_equal(b.force, [-0.045410156, -3.66032124e-04, -2.14325517e-05]) + c = forcefields.lennard_jones([2.0, 4.0, 6.0], [1.5, 4.0]) + assert_almost_equal(c.force, [-0.185302734, -1.46457553e-03, -8.57325038e-05]) + d = forcefields.lennard_jones([150.0, 300.0, 500.0], [200.0, 500.0]) + assert_almost_equal(d.force, [-1.7558299e-12, -1.3717421e-14, -3.8400000e-16]) def test_lennard_jones_sigma_epsilon_energy(self): a = forcefields.lennard_jones_sigma_epsilon(2.0, [1.0, 0.25]) - assert_almost_equal(a, -0.015380859) + assert_almost_equal(a.energy, -0.015380859) b = forcefields.lennard_jones_sigma_epsilon([2.0, 1.0], [1.0, 0.25]) - assert_almost_equal(b, [-0.015380859, 0]) + assert_almost_equal(b.energy, [-0.015380859, 0]) c = forcefields.lennard_jones_sigma_epsilon( [2.0, 1.0, 1.5], [0.5, 0.75]) - assert_almost_equal(c, [-0.0007322, -0.0461425, -0.0041095]) + assert_almost_equal(c.energy, [-0.0007322, -0.0461425, -0.0041095]) d = forcefields.lennard_jones_sigma_epsilon( [400.0, 500.0, 600.0], [5e-10, 9e-9]) - assert_almost_equal(d, [0, 0, 0]) + assert_almost_equal(d.energy, [0, 0, 0]) def test_lennard_jones_sigma_epsilon_force(self): - a = forcefields.lennard_jones_sigma_epsilon( - 2.0, [1.0, 0.25], force=True) - assert_almost_equal(a, -0.045410156) - b = forcefields.lennard_jones_sigma_epsilon( - [2.0, 4.0], [1.0, 0.25], force=True) - assert_almost_equal( - b, [-0.0454102, -0.000366]) - c = forcefields.lennard_jones_sigma_epsilon( - [3.0, 4.0], [3.0, 1.0], force=True) - assert_almost_equal(c, [8.0, -0.6877549]) + a = forcefields.lennard_jones_sigma_epsilon(2.0, [1.0, 0.25]) + assert_almost_equal(a.force, -0.045410156) + b = forcefields.lennard_jones_sigma_epsilon([2.0, 4.0], [1.0, 0.25]) + assert_almost_equal(b.force, [-0.0454102, -0.000366]) + c = forcefields.lennard_jones_sigma_epsilon([3.0, 4.0], [3.0, 1.0]) + assert_almost_equal(c.force, [8.0, -0.6877549]) def test_buckingham_energy(self): a = forcefields.buckingham(2.0, [1.0, 1.0, 1.0]) - assert_almost_equal(a, 0.1197103832) + assert_almost_equal(a.energy, 0.1197103832) b = forcefields.buckingham([2.0], [1.0, 1.0, 1.0]) - assert_almost_equal(b, 0.1197103832) + assert_almost_equal(b.energy, 0.1197103832) c = forcefields.buckingham([2.0, 4.0], [1.0, 1.0, 1.0]) - assert_almost_equal(c, [0.1197103832, 0.0180715]) + assert_almost_equal(c.energy, [0.1197103832, 0.0180715]) d = forcefields.buckingham([2.0, 4.0, 5.0], [0.01, 0.01, 0.01]) - assert_almost_equal(d, [0.0096457, 0.0096055, 0.0095117]) + assert_almost_equal(d.energy, [0.0096457, 0.0096055, 0.0095117]) def test_buckingham_force(self): - a = forcefields.buckingham(2.0, [1.0, 1.0, 1.0], force=True) - assert_almost_equal(a, 0.08846028324) - b = forcefields.buckingham([2.0], [1.0, 1.0, 1.0], force=True) - assert_almost_equal(b, 0.08846028324) - c = forcefields.buckingham( - [2.0, 1.0, 4.0], [1.5, 0.1, 2.0], force=True) - assert_almost_equal(c, [0.0290596, -11.8642744, 0.0998156]) + a = forcefields.buckingham(2.0, [1.0, 1.0, 1.0]) + assert_almost_equal(a.force, 0.08846028324) + b = forcefields.buckingham([2.0], [1.0, 1.0, 1.0]) + assert_almost_equal(b.force, 0.08846028324) + c = forcefields.buckingham([2.0, 1.0, 4.0], [1.5, 0.1, 2.0]) + assert_almost_equal(c.force, [0.0290596, -11.8642744, 0.0998156]) def test_square_well_energy(self): a = forcefields.square_well(2.0, [1.0, 1.5, 2.0]) - assert_equal(a, -1.0) + assert_equal(a.energy, -1.0) b = forcefields.square_well(0.5, [1.0, 2.0, 1.25]) - assert_equal(b, float('inf')) + assert_equal(b.energy, float('inf')) c = forcefields.square_well(3.0, [0.5, 1.5, 1.25]) - assert_equal(c, 0) + assert_equal(c.energy, 0) d = forcefields.square_well([2.0, 0.5], [1.0, 1.5, 2.0]) - assert_equal(d, [-1.0, float('inf')]) + assert_equal(d.energy, [-1.0, float('inf')]) e = forcefields.square_well([3.0, 3.0, 0.25], [1.0, 1.5, 1.25]) - assert_equal(e, [0, 0, float('inf')]) - f = forcefields.square_well( - [3.0, 3.0, 0.25], [1.0, 1.5, 1.25], max_val=5000) - assert_equal(f, [0, 0, 5000]) + assert_equal(e.energy, [0, 0, float('inf')]) + f = forcefields.square_well([3.0, 3.0, 0.25], [1.0, 1.5, 1.25], max_val=5000) + assert_equal(f.energy, [0, 0, 5000]) def test_square_well_force(self): + a = forcefields.square_well(2.0, [1.0, 1.5, 2.0]) with self.assertRaises(ValueError): - forcefields.square_well( - 2.0, [1.0, 1.5, 2.0], force=True) + a.update_force() + b = forcefields.square_well([2.0], [1.0, 1.5, 2.0]) with self.assertRaises(ValueError): - forcefields.square_well( - [2.0], [1.0, 1.5, 2.0], force=True) + b.update_force() if __name__ == '__main__':