api_potential.rst

Potentials

A Potential object is a compiled interpolation of a given function. The Universe applies potentials to particles to calculate the net force on them.

For performance reasons, we found that implementing potentials as interpolations can be much faster than evaluating the function directly.

A potential can be treated just like any python callable object to evaluate it:

>>> pot = m.Potential.lennard_jones_12_6(0.1, 5, 9.5075e-06 , 6.1545e-03, 1.0e-3 )
>>> x = np.linspace(0.1,1,100)
>>> y=[pot(j) for j in x]
>>> plt.plot(x,y, 'r')

Potential is a callable object, we can invoke it like any Python function.

power([k], [r0], [min], max=[10], [tol=0.001])

param k

interaction strength, represents the potential energy peak value. Defaults to 1

param r0

potential rest length, zero of the potential, defaults to 0.

param alpha

Exponent, defaults to 1.

param min

minimal value potential is computed for, defaults to r0 / 2.

param max

cutoff distance, defaults to 3 * r0.

param tol

Tolerance, defaults to 0.01.

some stuff, r₀ − r₀/2.

The power potential the general form of many of the potential functions, such as linear, etc... power has the form:

U(r, k) = k(r − r₀)^α

Some examples of power potentials are:

>>> import mechanica as m
>>> p = m.Potential.power(k=1, r0=1, alpha=1.1)
>>> p.plot(potential=True)

power potential function of U(r) = 1.0(r − 1)^1.1

>>> import mechanica as m >>> p = m.Potential.power(k=1,alpha=0.5, r0=1) >>> p.plot(potential=True)

power potential function of U(r) = 1.0(r − 1)^0.5

>>> import mechanica as m >>> p = m.Potential.power(k=1,alpha=2.5, r0=1) >>> p.plot(potential=True)

power potential function of U(r) = 1.0(r − 1)^2.5

lennard_jones_12_6(min, max, a, b, )

Creates a Potential representing a 12-6 Lennard-Jones potential

param min

The smallest radius for which the potential will be constructed.

param max

The largest radius for which the potential will be constructed.

param A

The first parameter of the Lennard-Jones potential.

param B

The second parameter of the Lennard-Jones potential.

param tol

The tolerance to which the interpolation should match the exact potential., optional

The Lennard Jones potential has the form:

$$\left( \frac{A}{r^{12}} - \frac{B}{r^6} \right)$$

lennard_jones_12_6_coulomb(min, max, a, b, tol )

Creates a Potential representing the sum of a

12-6 Lennard-Jones potential and a shifted Coulomb potential.

param min

The smallest radius for which the potential will be constructed.

param max

The largest radius for which the potential will be constructed.

param A

The first parameter of the Lennard-Jones potential.

param B

The second parameter of the Lennard-Jones potential.

param q

The charge scaling of the potential.

param tol

The tolerance to which the interpolation should match the exact potential. (optional)

The 12-6 Lennard Jones - Coulomb potential has the form:

$$\left( \frac{A}{r^{12}} - \frac{B}{r^6} \right) + q \left(\frac{1}{r} - \frac{1}{max} \right)$$

ewald(min, max, q, kappa, tol)

Creates a potential representing the real-space part of an Ewald potential.

param min

The smallest radius for which the potential will be constructed.

param max

The largest radius for which the potential will be constructed.

param q

The charge scaling of the potential.

param kappa

The screening distance of the Ewald potential.

param tol

The tolerance to which the interpolation should match the exact potential.

The Ewald potential has the form:

$$q \frac{\mbox{erfc}( \kappa r)}{r}$$

well(k, n, r0, [min], [max], [tol])

Creates a continuous square well potential. Usefull for binding a particle to a region.

param float k

potential prefactor constant, should be decreased for larger n.

param float n

exponent of the potential, larger n makes a sharper potential.

param float r0

The extents of the potential, length units. Represents the maximum extents that a two objects connected with this potential should come appart.

param float min

[optional] The smallest radius for which the potential will be constructed, defaults to zero.

param float max

[optional] The largest radius for which the potential will be constructed, defaults to r0.

param float tol

[optional[ The tolerance to which the interpolation should match the exact potential, defaults to 0.01 * abs(min-max).

$$\frac{k}{\left(r_0 - r\right)^{n}}$$

As with all potentials, we can create one, and plot it like so:

>>> p = m.Potential.well(0.01, 2, 1)
>>> x=n.arange(0, 1, 0.0001)
>>> y = [p(xx) for xx in x]
>>> plt.plot(x, y)
>>> plt.title(r"Continuous Square Well Potential $\frac{0.01}{(1 - r)^{2}}$ \n",
...           fontsize=16, color='black')

A continuous square well potential.

harmonic_angle(k, theta0, [min], max, [tol])

Creates a harmonic angle potential

param k

The energy of the angle.

param theta0

The minimum energy angle.

param min

The smallest angle for which the potential will be constructed.

param max

The largest angle for which the potential will be constructed.

param tol

The tolerance to which the interpolation should match the exact potential.

returns A newly-allocated potential representing the potential

k(θ − θ₀)²

Note, for computational effeciency, this actually generates a function of r, where r is the cosine of the angle (calculated from the dot product of the two vectors. So, this actually evaluates internally,

k(arccos (r) − θ₀)²

harmonic(k, r0, [min], [max], [tol])

Creates a harmonic bond potential

param k

The energy of the bond.

param r0

The bond rest length

param min

[optional] The smallest radius for which the potential will be constructed. Defaults to r₀ − r₀/2.

param max

[optional] The largest radius for which the potential will be constructed. Defaults to r₀ + r₀/2.

param tol

[optional] The tolerance to which the interpolation should match the exact potential. Defaults to $0.01 \abs(max-min)$

return A newly-allocated potential

U(r) = k(r − r₀)²

We can plot the harmonic function to get an idea of it's behavior:

>>> import mechanica as m
>>> p = p = m.Potential.harmonic(k=10, r0=1, min=0.001, max=10)
>>> p.plot(s=2, potential=True)

soft_sphere(kappa, epsilon, r0, eta, min, max, tol)

Creates a soft sphere interaction potential. The soft sphere is a generalized Lennard Jones type potentail, but you can varry the exponents to create a softer interaction.

param kappa

param epsilon

param r0

param eta

param min

param max

param tol

glj(e, [m], [n], [r0], [min], [max], [tol], [shifted])

param e

effective energy of the potential

param m

order of potential, defaults to 3

param n

order of potential, defaults to 2*m

param r0

mimumum of the potential, defaults to 1

param min

minimum distance, defaults to 0.05 * r0

param max

max distance, defaults to 5 * r0

param tol

tolerance, defaults to 0.01

param shifted

is this shifted potential, defaults to true

Generalized Lennard-Jones potential.

$$V^{GLJ}_{m,n}(r) = \frac{\epsilon}{n-m} \left[ m \left( \frac{r_0}{r} \right)^n - n \left( \frac{r_0}{r} \right) ^ m \right]$$

where r_e is the effective radius, which is automatically computed as the sum of the interacting particle radii.

The Generalized Lennard-Jones potential for different exponents (m, n) with fixed n = 2m. As the exponents grow smaller, the potential flattens out and becomes softer, but as the exponents grow larger the potential becomes narrower and sharper, and approaches the hard sphere potential.

overlapping_sphere(mu=1, [kc=1], [kh=0], [r0=0], [min=0.001], max=[10], [tol=0.001])

param mu

interaction strength, represents the potential energy peak value.

param kc

decay strength of long range attraction. Larger values make a shorter ranged function.

param kh

Optionally add a harmonic long-range attraction, same as glj function.

param r0

Optional harmonic rest length, only used if kh is non-zero.

param min

Minimum value potential is computed for.

param max

Potential cutoff values.

param tol

Tolerance, defaults to 0.001.

The overlapping_sphere function implements the Overlapping Sphere, from Osborne:2017hk. This is a soft sphere, from our testing, it appears too soft, probably better suited for 2D models. This potential appears to allow particles to collapse too closely, probably needs more paramater fiddling.

Note

From the equation below, we can see that there is a log  term as the short range repulsion term. The logarithm is the radial Green's function for cylindrical (2D) geometry, however the Green's function for 3D is the 1/r function. This is possibly why Osborne has success in 2D but it's unclear if this was used in 3D geometries.

$$\begin{aligned} \mathbf{F}_{ij}= \begin{cases} \mu_{ij} s_{ij}(t) \hat{\mathbf{r}}_{ij} \log \left( 1 + \frac{||\mathbf{r}_{ij}|| - s_{ij}(t)}{s_{ij}(t)} \right) ,& \text{if } ||\mathbf{r}_{ij}|| < s_{ij}(t) \\\ \mu_{ij}\left(||\mathbf{r}_{ij}|| - s_{ij}(t)\right) \hat{\mathbf{r}}_{ij} \exp \left( -k_c \frac{||\mathbf{r}_{ij}|| - s_{ij}(t)}{s_{ij}(t)} \right) ,& \text{if } s_{ij}(t) \leq ||\mathbf{r}_{ij}|| \leq r_{max} \\\ 0, & \text{otherwise} \\\ \end{cases} \end{aligned}$$

Osborne refers to μ_ij as a "spring constant", this controls the size of the force, and is the potential energy peak value. $\hat{\mathbf{r}}_{ij}$ is the unit vector from particle i center to particle j center, k_C is a parameter that defines decay of the attractive force. Larger values of k_C result in a shaper peaked attraction, and thus a shorter ranged force. s_ij(t) is the is the sum of the radii of the two particles.

We can plot the overlapping sphere function to get an idea of it's behavior:

>>> import mechanica as m
>>> p = m.Potential.overlapping_sphere(mu=10, max=20)
>>> p.plot(s=2, force=True, potential=True, ymin=-10, ymax=8)

We can plot the overlapping sphere function like any other function.

linear(k, [min], max=[10], [tol=0.001])

param k

interaction strength, represents the potential energy peak value.

param min

Minimum value potential is computed for.

param max

Potential cutoff values.

param tol

Tolerance, defaults to 0.001.

The linear potential is the simplest one, it's simply a potential of the form

U(r) = kr

>>> import mechanica as m >>> p = p = m.Potential.linear(k=5, max=10) >>> p.plot(potential=True)

Linear potential function

Dissapative Particle Dynamics

Potential.dpd(alpha=1, gamma=1, sigma=1, cutoff=1)

param alpha

interaction strength of the conservative force

param gamma

interaction strength of dissapative force

param sigma

strength of random force.

The dissapative particle dynamics force is the sum of the conservative, dissapative and random forces:

F_ij = F_ij^C + F_ij^D + F_ij^R

The conservative force is:

$$\mathbf{F}^C_{ij} = \alpha \left(1 - \frac{r_{ij}}{r_c}\right) \mathbf{e}_{ij}$$

The dissapative force is:

$$\mathbf{F}^D_{ij} = -\gamma \left(1 - \frac{r_{ij}}{r_c}\right)^{2}(\mathbf{e}_{ij} \cdot \mathbf{v}_{ij}) \mathbf{e}_{ij}$$

And the random force is:

$$\mathbf{F}^R_{ij} = \sigma \left(1 - \frac{r_{ij}}{r_c}\right) \xi_{ij}\Delta t^{-1/2}\mathbf{e}_{ij}$$

A subtype of the Potential class that adds features specific to Dissapative Particle Dynamics.