"""A collection of classes and functions related to defining potentials"""
from . import _tablereaders
from ._potential import Potential
from ._eam_potential import EAMPotential
from .spline import SplinePotential # noqa
from ._util import gradient, num_deriv, deriv # noqa
from ._dlpoly_writeTABEAM import writeTABEAM # noqa
from ._dlpoly_writeTABEAM import writeTABEAMFinnisSinclair # noqa
from ._lammpsWriteEAM import writeFuncFL # noqa
from ._lammpsWriteEAM import writeSetFL # noqa
from ._lammpsWriteEAM import writeSetFLFinnisSinclair # noqa
from .potentialforms import * # noqa
from . import potentialfunctions
from . import spline
from . import pair_tabulation
from . import eam_tabulation
from . import tableforms
from ._multi_range_potential_form import create_Multi_Range_Potential_Form, Multi_Range_Defn
import sys
import math
[docs]def plus(a,b):
"""Takes two functions and returns a third which when evaluated returns the result of a(r) + b(r)
This function is useful for combining existing potentials.
**Derivatives:**
If either of the potential callables (`a` and `b`) provide a .deriv() method the function returned by
`plus()` will also have a `.deriv()` method. This allows analytical derivatives to be specified. If
only one of `a` or `b` provide `.deriv()` then the derivative of the other callable will be evaluated
numerically.
If neither function has a .deriv() method then the function returned here will also *not* have a .deriv()
method.
**Example:**
To combine :func:`~atsim.potentials.potentialforms.buck` and :func:`~atsim.potentials.potentialforms.hbnd` functions from the :mod:`atsim.potentials.potentialforms` module to give:
.. code:: python
A*(-r/rho) + C/r**6 + D/r**12 - E/r**10
this function can then be used as follows:
.. code:: python
plus(buck(A,rho,C), hbnd(D,E))
:param a: First callable
:param b: Second callable
:return: Function that when evaulated returns ``a(r) + b(r)``"""
def potential(r):
return a(r) + b(r)
# Set derivatives
if hasattr(a, 'deriv') or hasattr(b, 'deriv'):
deriv_a = gradient(a)
deriv_b = gradient(b)
def deriv(r):
return deriv_a(r) + deriv_b(r)
potential.deriv = deriv
if hasattr(deriv_a, 'deriv') or hasattr(deriv_b, 'deriv'):
deriv2_a = gradient(deriv_a)
deriv2_b = gradient(deriv_b)
def deriv2(r):
return deriv2_a(r) + deriv2_b(r)
potential.deriv2 = deriv2
return potential
[docs]def product(a,b):
"""Takes two callables and returns a third which when evaluated returns the result of a(r) * b(r)
This function is useful for combining existing potentials.
**Derivatives:**
If either of the potential callables (`a` and `b`) provide a .deriv() method the function returned by
`product()` will also have a `.deriv()` method. This allows analytical derivatives to be specified. If
only one of `a` or `b` provide `.deriv()` then the derivative of the other callable will be evaluated
numerically.
If neither function has a .deriv() method then the function returned here will also *not* have a .deriv()
method.
:param a: First callable
:param b: Second callable
:return: Function that when evaulated returns ``a(r) * b(r)``"""
def potential(r):
return a(r) * b(r)
# Set derivatives
if hasattr(a, 'deriv') or hasattr(b, 'deriv'):
deriv_a = gradient(a)
deriv_b = gradient(b)
def deriv(r):
return a(r) * deriv_b(r) + b(r) * deriv_a(r)
potential.deriv = deriv
if hasattr(deriv_a, 'deriv') or hasattr(deriv_b, 'deriv'):
deriv2_a = gradient(deriv_a)
deriv2_b = gradient(deriv_b)
def deriv2(r):
return deriv2_a(r) * b(r) + 2.0*deriv_a(r)*deriv_b(r) + a(r)*deriv2_b(r)
potential.deriv2 = deriv2
return potential
[docs]def pow(a,b):
"""Takes two callables and returns a third which when evaluated returns the result of a(r)**b(r)
This function is useful for combining existing potentials.
**Derivatives:**
If either of the potential callables (`a` and `b`) provide a .deriv() method the function returned by
`pow()` will also have a `.deriv()` method. This allows analytical derivatives to be specified. If
only one of `a` or `b` provide `.deriv()` then the derivative of the other callable will be evaluated
numerically.
If neither function has a .deriv() method then the function returned here will also *not* have a .deriv()
method.
:param a: First callable
:param b: Second callable
:return: Function that when evaulated returns ``a(r)**b(r)`` (a to the power of b)"""
def potential(r):
return a(r)**b(r)
# Set derivatives
if hasattr(a, 'deriv') or hasattr(b, 'deriv'):
deriv_a = gradient(a)
deriv_b = gradient(b)
def deriv(r):
ar = a(r)
return potential(r) * (deriv_b(r) * math.log(ar) + b(r) * deriv_a(r)/ar)
potential.deriv = deriv
if hasattr(deriv_a, 'deriv') or hasattr(deriv_b, 'deriv'):
deriv2_a = gradient(deriv_a)
deriv2_b = gradient(deriv_b)
def deriv2(r):
ar = a(r)
br = b(r)
dr = deriv(r)
p = potential(r)
da = deriv_a(r)
db = deriv_b(r)
d2a = deriv2_a(r)
d2b = deriv2_b(r)
# value = (deriv_b(r)*log(a(r)) + b(r)*deriv_a(r)/a(r))*deriv(r) + (math.log(a(r))*deriv2_b(r) + b(r)*deriv2_a(r)/a(r) + deriv_a(r)*deriv2_b(r)/a(r) + deriv_b(r)*deriv2_a(r)/a(r) - b(r)*deriv_a(r)*deriv2_a(r)/a(r)**2)*potential(r)
value = (db*math.log(ar) + (br*da)/ar)*dr + (math.log(ar)*d2b + (br*d2a)/ar + (da*db)/ar + (db*da)/ar - (br*da*da)/(ar**2))*p
return value
potential.deriv2 = deriv2
return potential
[docs]class TableReader(object):
"""Callable that allows pretabulated data to be used with a Potential object."""
def __init__(self, fileobject):
"""The file passed to tablereader() is assumed to have one separation, energy pair per line, separated by
a space.
:param fileobject: Python file object from which separation, energy pairs should be read"""
self._tablereader = _tablereaders.DatReader(fileobject)
@property
[docs] def datReader(self):
""":return: _tablereaders.DatReader associated with this callable"""
return self._tablereader
[docs] def __call__(self, separation):
return self._tablereader.getValue(separation)
[docs]def plotToFile(fileobj,lowx, highx, func, steps=10000):
"""Convenience function for plotting the potential functions contained herein.
Data is written to a text file as two columns (r and E) separated by spaces
with no header.
:param fileobj: Python file object into which data should be plotted
:param lowx: X-axis lower value
:param highx: X-axis upper value
:param func: Function to be plotted
:param steps: Number of data points to be plotted"""
step = (highx - lowx) / float(steps)
for i in range(steps):
v = lowx + float(i)*step
y = func(v)
fileobj.write("{0} {1}\n".format(v,y))
[docs]def plot(filename, lowx, highx, func, steps=10000):
"""Convenience function for plotting the potential functions contained herein.
Data is written to a text file as two columns (r and E) separated by spaces
with no header.
:param filename: File into which data should be plotted
:param lowx: X-axis lower value
:param highx: X-axis upper value
:param func: Function to be plotted
:param steps: Number of data points to be plotted"""
with open(filename, 'w') as outfile:
plotToFile(outfile, lowx, highx, func, steps)
[docs]def plotPotentialObject(filename, lowx, highx, potentialObject, steps=10000):
"""Convenience function for plotting energy of pair interactions
given by instances of :class:`atsim.potentials.Potential` obtained by calling
`potential` `.energy()` method.
Data is written to a text file as two columns (r and E) separated by spaces
with no header.
:param filename: File into which data should be plotted
:param lowx: X-axis lower value
:param highx: X-axis upper value
:param func: :class:`atsim.potentials.Potential` object.
:param steps: Number of data points to be plotted"""
with open(filename, 'w') as outfile:
plotPotentialObjectToFile(outfile, lowx, highx, potentialObject, steps)
[docs]def plotPotentialObjectToFile(fileobj, lowx, highx, potentialObject, steps=10000):
"""Convenience function for plotting energy of pair interactions
given by instances of :class:`atsim.potentials.Potential` obtained by calling
`potential` `.energy()` method.
Data is written to a text file as two columns (r and E) separated by spaces
with no header.
:param fileobj: Python file object into which data should be plotted
:param lowx: X-axis lower value
:param highx: X-axis upper value
:param func: :class:`atsim.potentials.Potential` object.
:param steps: Number of data points to be plotted"""
def f(r):
return potentialObject.energy(r)
plotToFile(fileobj, lowx, highx, f, steps)
def _LAMMPS_writePotentials(potentialList, cutoff, gridPoints, out):
"""Wrapper function that adapts lammps.writeTABLE.writePotentials() to the API
expected by potentials.writePotentials"""
minr = cutoff / float(gridPoints)
from ._lammps_writeTABLE import writePotentials
writePotentials(potentialList,minr, cutoff, gridPoints, out)
[docs]class UnsupportedTabulationType(Exception):
"""Exception thrown by writePotentials() when unknown tabulation type specified"""
pass
[docs]def writePotentials(outputType, potentialList, cutoff, gridPoints, out = sys.stdout):
"""Tabulates pair-potentials in formats suitable for multiple simulation codes.
* The ``outputType`` parameter can be one of the following:
* ``DL_POLY``:
* This function creates output that can be written to a ``TABLE`` and used
within DL_POLY.
* for a working example see :ref:`Quick-Start: DL_POLY <quick-start>`.
* ``GULP``:
* Creates output for the `GULP code <https://nanochemistry.curtin.edu.au/gulp/>`_
* Output is in the form of a series of `spline potential forms <https://nanochemistry.curtin.edu.au/gulp/help/new_help_40_txt.html#spline>`_
* The generated file can be loaded into GULP using the `library command <https://nanochemistry.curtin.edu.au/gulp/help/new_help_40_txt.html#library>`_
* ``LAMMPS``:
* Creates files readable by LAMMPS `pair_style table <http://lammps.sandia.gov/doc/pair_table.html>`_
* Each file can contain multiple potentials:
* the block representing each potential has a title formed from the ``speciesA``
and ``speciesB`` attributes of the :class:`~atsim.potentials.Potential`
instance represented by the block. These are sorted into their natural
order and separated by a hyphen to form the title.
* **Example:**
* For a :class:`~atsim.potentials.Potential` where ``speciesA`` = Xe
and ``speciesB`` = O the block title would be: ``O-Xe``.
* If ``speciesA`` = B
and ``speciesB`` = O the block title would be: ``B-O``.
* within LAMMPSthe block title is used as the ``keyword`` argument to the
`pair_style table <http://lammps.sandia.gov/doc/pair_table.html>`_
``pair_coeff`` directive.
:param outputType: The type of output that should be created can be one of: ``DL_POLY`` or ``LAMMPS``
:type outputType: str
:param potentialList: List of Potential objects to be tabulated.
:type potentialList: list
:param cutoff: Largest separation to be tabulated.
:type cutoff: float
:param gridPoints: Number of rows in tabulation.
:type gridPoints: int
:param out: Python file like object to which tabulation should be written
:type out: file"""
from .pair_tabulation import DLPoly_PairTabulation, LAMMPS_PairTabulation, GULP_PairTabulation
supportedTabulations = {
'DL_POLY' : DLPoly_PairTabulation,
'LAMMPS' : LAMMPS_PairTabulation,
'GULP' : GULP_PairTabulation
}
if not outputType in supportedTabulations:
raise UnsupportedTabulationType("Unsupported tabulation type: '{}' should be one of: {}".format(
outputType, ",".join(sorted(supportedTabulations.keys()))
))
tabulation = supportedTabulations[outputType](potentialList,cutoff, gridPoints)
tabulation.write(out)