index | modules | gitlab

# Source code for ase.calculators.morse

    import numpy as np

    from ase.calculators.calculator import Calculator, all_changes
    from ase.neighborlist import neighbor_list
    from ase.stress import full_3x3_to_voigt_6_stress

    def fcut(r: np.ndarray, r0: float, r1: float) -> np.ndarray:
        """Piecewise quintic C^{2,1} regular polynomial for use as a smooth cutoff.

        Ported from JuLIP.jl.

        https://github.com/JuliaMolSim/JuLIP.jl/blob/master/src/cutoffs.jl

        https://en.wikipedia.org/wiki/Smoothstep

        Parameters
        ----------
        r : np.ndarray
            Distances between atoms.
        r0 : float
            Inner cutoff radius.
        r1 : float
            Outder cutoff radius.

        Returns
        -------
        np.ndarray
            Sigmoid-like function smoothly interpolating (r0, 1) and (r1, 0).

        """""
        s = 1.0 - (r - r0) / (r1 - r0)
        return (s >= 1.0) + ((s > 0.0) & (s < 1.0)) * (
            6.0 * s**5 - 15.0 * s**4 + 10.0 * s**3
        )

    def fcut_d(r: np.ndarray, r0: float, r1: float) -> np.ndarray:
        """Derivative of fcut() function defined above."""
        s = 1.0 - (r - r0) / (r1 - r0)
        return -(
            ((s > 0.0) & (s < 1.0))
            * (30.0 * s**4 - 60.0 * s**3 + 30.0 * s**2)
            / (r1 - r0)
        )

    class MorsePotential(Calculator):
        """Morse potential."""

        implemented_properties = [
            'energy', 'energies', 'free_energy', 'forces', 'stress',
        ]
        default_parameters = {'epsilon': 1.0,
                              'rho0': 6.0,
                              'r0': 1.0,
                              'rcut1': 1.9,
                              'rcut2': 2.7}
        nolabel = True

        def __init__(self, **kwargs):
            r"""

            The pairwise energy between atoms *i* and *j* is given by

            .. math::

                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

              .. math::

                  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.

            """
            Calculator.__init__(self, **kwargs)

        def calculate(self, atoms=None, properties=['energy'],
                      system_changes=all_changes):
            Calculator.calculate(self, atoms, properties, system_changes)
            epsilon = self.parameters['epsilon']
            rho0 = self.parameters['rho0']
            r0 = self.parameters['r0']
            rcut1 = self.parameters['rcut1'] * r0
            rcut2 = self.parameters['rcut2'] * r0

            number_of_atoms = len(self.atoms)

            forces = np.zeros((number_of_atoms, 3))

            i, _j, d, D = neighbor_list('ijdD', atoms, rcut2)
            dhat = (D / d[:, None]).T

            expf = np.exp(rho0 * (1.0 - d / r0))

            cutoff_fn = fcut(d, rcut1, rcut2)
            d_cutoff_fn = fcut_d(d, rcut1, rcut2)

            pairwise_energies = epsilon * expf * (expf - 2.0)
            self.results['energies'] = np.bincount(
                i,
                weights=0.5 * (pairwise_energies * cutoff_fn),
                minlength=number_of_atoms,
            )
            self.results['energy'] = self.results['energies'].sum()
            self.results['free_energy'] = self.results['energy']

            # derivatives of `pair_energies` with respect to `d`
            de = (-2.0 * epsilon * rho0 / r0) * expf * (expf - 1.0)

            # smoothened `de`
            de = de * cutoff_fn + pairwise_energies * d_cutoff_fn

            de_vec = (de * dhat).T
            for dim in range(3):
                forces[:, dim] = np.bincount(
                    i,
                    weights=de_vec[:, dim],
                    minlength=number_of_atoms,
                )
            self.results['forces'] = forces

            if self.atoms.cell.rank == 3:
                stress = 0.5 * (D.T @ de_vec) / self.atoms.get_volume()
                self.results['stress'] = full_3x3_to_voigt_6_stress(stress)