Skip to content

Commit

Permalink
updating forcefields docstrings
Browse files Browse the repository at this point in the history
  • Loading branch information
maxdolan committed Jun 17, 2024
1 parent 877d58d commit 449abc2
Showing 1 changed file with 118 additions and 78 deletions.
196 changes: 118 additions & 78 deletions pylj/forcefields.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,26 +6,13 @@ class lennard_jones(object):
r"""Calculate the energy or force for a pair of particles using the
Lennard-Jones (A/B variant) forcefield.
.. math::
E = \frac{A}{dr^{12}} - \frac{B}{dr^6}
.. math::
f = \frac{12A}{dr^{13}} - \frac{6B}{dr^7}
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
force: bool (optional)
If true, the negative first derivative will be found.
Returns
-------
float: array_like
The potential energy or force between the particles.
"""
def __init__(self, dr, constants):
self.dr = dr
Expand All @@ -35,37 +22,48 @@ def __init__(self, dr, constants):
self.force = self.update_force()

def update_energy(self):
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}
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))
return self.energy

def update_force(self):
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}
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))
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.
.. math::
E = \frac{4e*a^{12}}{dr^{12}} - \frac{4e*a^{6}}{dr^6}
.. math::
f = \frac{48e*a^{12}}{dr^{13}} - \frac{24e*a^{6}}{dr^7}
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
force: bool (optional)
If true, the negative first derivative will be found.
Returns
-------
float: array_like
The potential energy or force between the particles.
"""
def __init__(self, dr, constants):
self.dr = dr
Expand All @@ -75,39 +73,50 @@ def __init__(self, dr, constants):
self.force = self.update_force()

def update_energy(self):
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}
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))
return self.energy

def update_force(self):
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}
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))
return self.force



class buckingham(object):
r""" Calculate the energy or force for a pair of particles using the
Buckingham forcefield.
.. math::
E = Ae^{(-Bdr)} - \frac{C}{dr^6}
.. math::
f = ABe^{(-Bdr)} - \frac{6C}{dr^7}
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.
force: bool (optional)
If true, the negative first derivative will be found.
Returns
-------
float: array_like
The potential energy or force between the particles.
"""
def __init__(self, dr, constants):
self.dr = dr
Expand All @@ -118,33 +127,41 @@ def __init__(self, dr, constants):
self.force = self.update_force()

def update_energy(self):
r"""Calculate the energy for a pair of particles using the
Buckingham forcefield.
.. math::
E = Ae^{(-Bdr)} - \frac{C}{dr^6}
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)
return self.energy

def update_force(self):
r"""Calculate the force for a pair of particles using the
Buckingham forcefield.
.. math::
f = ABe^{(-Bdr)} - \frac{6C}{dr^7}
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)
return self.force



class square_well(object):
r'''Calculate the energy or force for a pair of particles using a
square well model.
.. math::
E = {
if dr < sigma:
E = max_val
elif sigma <= dr < lambda * sigma:
E = -epsilon
elif r >= lambda * sigma:
E = 0
}
.. math::
f = {
if sigma <= dr < lambda * sigma:
f = inf
else:
f = 0
}
Parameters
----------
dr: float, array_like
Expand All @@ -154,13 +171,6 @@ class square_well(object):
parameters for the square well model.
max_val: int (optional)
Upper bound for values in square well - replaces usual infinite values
force: bool (optional)
If true, the negative first derivative will be found.
Returns
-------
float: array_like
The potential energy between the particles.
'''
def __init__(self, dr, constants, max_val=np.inf):

Expand All @@ -179,22 +189,52 @@ def __init__(self, dr, constants, max_val=np.inf):
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')
r'''Calculate the energy for a pair of particles using a
square well model.
.. math::
E = {
if dr < sigma:
E = max_val
elif sigma <= dr < lambda * sigma:
E = -epsilon
elif r >= lambda * sigma:
E = 0
}
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
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
return self.energy

def update_force(self):
raise ValueError("Force is infinite at sigma <= dr < lambda * sigma")
r'''The force of a pair of particles using a square well model is given by:
.. math::
f = {
if sigma <= dr < lambda * sigma:
f = inf
else:
f = 0
}
Therefore the force here will always be infinite, and therefore not possible to simulate
'''
raise ValueError("Force is infinite at sigma <= dr < lambda * sigma")

0 comments on commit 449abc2

Please sign in to comment.