Other built-in calculators

TIP3P

class ase.calculators.tip3p.TIP3P(rc=5.0, width=1.0)[source]

TIP3P potential.

rc: float

Cutoff radius for Coulomb part.

width: float

Width for cutoff function for Coulomb part.

TIP4P

class ase.calculators.tip4p.TIP4P(rc=7.0, width=1.0)[source]

TIP4P potential for water.

doi: 10.1063/1.445869

Requires an atoms object of OHH,OHH, … sequence Correct TIP4P charges and LJ parameters set automatically.

Virtual interaction sites implemented in the following scheme: Original atoms object has no virtual sites. When energy/forces are requested:

  • virtual sites added to temporary xatoms object

  • energy / forces calculated

  • forces redistributed from virtual sites to actual atoms object

This means you do not get into trouble when propagating your system with MD while having to skip / account for massless virtual sites.

This also means that if using for QM/MM MD with GPAW, the EmbedTIP4P class must be used.

Lennard-Jones

class ase.calculators.lj.LennardJones(**kwargs)[source]

Lennard Jones potential calculator

see https://en.wikipedia.org/wiki/Lennard-Jones_potential

The fundamental definition of this potential is a pairwise energy:

u_ij = 4 epsilon ( sigma^12/r_ij^12 - sigma^6/r_ij^6 )

For convenience, we’ll use d_ij to refer to “distance vector” and r_ij to refer to “scalar distance”. So, with position vectors \(r_i\):

r_ij = | r_j - r_i | = | d_ij |

Therefore:

d r_ij / d d_ij = + d_ij / r_ij d r_ij / d d_i  = - d_ij / r_ij

The derivative of u_ij is:

d u_ij / d r_ij
= (-24 epsilon / r_ij) ( 2 sigma^12/r_ij^12 - sigma^6/r_ij^6 )

We can define a “pairwise force”

f_ij = d u_ij / d d_ij = d u_ij / d r_ij * d_ij / r_ij

The terms in front of d_ij are combined into a “general derivative”.

du_ij = (d u_ij / d d_ij) / r_ij

We do this for convenience: \(du_ij\) is purely scalar The pairwise force is:

f_ij = du_ij * d_ij

The total force on an atom is:

f_i = sum_(j != i) f_ij

There is some freedom of choice in assigning atomic energies, i.e. choosing a way to partition the total energy into atomic contributions.

We choose a symmetric approach (\(bothways=True\) in the neighbor list):

u_i = 1/2 sum_(j != i) u_ij

The total energy of a system of atoms is then:

u = sum_i u_i = 1/2 sum_(i, j != i) u_ij

Differentiating \(u\) with respect to \(r_i\) yields the force, independent of the choice of partitioning.

f_i = - d u / d r_i = - sum_ij d u_ij / d r_i
    = - sum_ij d u_ij / d r_ij * d r_ij / d r_i
    = sum_ij du_ij d_ij = sum_ij f_ij

This justifies calling \(f_ij\) pairwise forces.

The stress can be written as ( \((x)\) denoting outer product):

sigma = 1/2 sum_(i, j != i) f_ij (x) d_ij = sum_i sigma_i , with atomic contributions

sigma_i  = 1/2 sum_(j != i) f_ij (x) d_ij

Another consideration is the cutoff. We have to ensure that the potential goes to zero smoothly as an atom moves across the cutoff threshold, otherwise the potential is not continuous. In cases where the cutoff is so large that u_ij is very small at the cutoff this is automatically ensured, but in general, \(u_ij(rc) != 0\).

This implementation offers two ways to deal with this:

Either, we shift the pairwise energy

u'_ij = u_ij - u_ij(rc)

which ensures that it is precisely zero at the cutoff. However, this means that the energy effectively depends on the cutoff, which might lead to unexpected results! If this option is chosen, the forces discontinuously jump to zero at the cutoff.

An alternative is to modify the pairwise potential by multiplying it with a cutoff function that goes from 1 to 0 between an onset radius ro and the cutoff rc. If the function is chosen suitably, it can also smoothly push the forces down to zero, ensuring continuous forces as well. In order for this to work well, the onset radius has to be set suitably, typically around 2*sigma.

In this case, we introduce a modified pairwise potential:

u'_ij = fc * u_ij

The pairwise forces have to be modified accordingly:

f'_ij = fc * f_ij + fc' * u_ij

Where \(fc' = d fc / d d_ij\).

This approach is taken from Jax-MD (https://github.com/google/jax-md), which in turn is inspired by HOOMD Blue (https://glotzerlab.engin.umich.edu/hoomd-blue/).

Parameters:
  • sigma (float) – The potential minimum is at 2**(1/6) * sigma, default 1.0

  • epsilon (float) – The potential depth, default 1.0

  • rc (float, None) – Cut-off for the NeighborList is set to 3 * sigma if None. The energy is upshifted to be continuous at rc. Default None

  • ro (float, None) – Onset of cutoff function in ‘smooth’ mode. Defaults to 0.66 * rc.

  • smooth (bool, False) – Cutoff mode. False means that the pairwise energy is simply shifted to be 0 at r = rc, leading to the energy going to 0 continuously, but the forces jumping to zero discontinuously at the cutoff. True means that a smooth cutoff function is multiplied to the pairwise energy that smoothly goes to 0 between ro and rc. Both energy and forces are continuous in that case. If smooth=True, make sure to check the tail of the forces for kinks, ro might have to be adjusted to avoid distorting the potential too much.

Morse

class ase.calculators.morse.MorsePotential(**kwargs)[source]

Morse potential.

The pairwise energy between atoms i and j is given by

\[V_{ij} = \epsilon \left( \mathrm{e}^{-2 \rho_0 (r_{ij} / r_0 - 1)} - 2 \mathrm{e}^{- \rho_0 (r_{ij} / r_0 - 1)} \right)\]
Parameters:
  • epsilon (float, default 1.0) – Absolute minimum depth.

  • r0 (float, default 1.0) – Minimum distance.

  • rho0 (float, default 6.0) –

    Exponential prefactor. The force constant in the potential minimum is given by

    \[k = 2 \epsilon \left(\frac{\rho_0}{r_0}\right)^2.\]

  • rcut1 (float, default 1.9) – Distance starting a smooth cutoff normalized by r0.

  • rcut2 (float, default 2.7) – Distance ending a smooth cutoff normalized by r0.

Notes

The default values are chosen to be similar as Lennard-Jones.