From 0cbb7c2cafa1f2d0f3e54228a56d7b2ea0664064 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Tue, 28 May 2024 14:12:52 -0500 Subject: [PATCH] Remove pair.slj --- src/PairEvaluatorShiftedLJ.h | 132 --------------------- src/PairPotentialShiftedLJ.cu | 18 --- src/PairPotentials.h | 1 - src/module.cc | 1 - src/pair.py | 93 --------------- src/pytest/test_pair_slj.py | 212 ---------------------------------- 6 files changed, 457 deletions(-) delete mode 100644 src/PairEvaluatorShiftedLJ.h delete mode 100644 src/PairPotentialShiftedLJ.cu delete mode 100644 src/pytest/test_pair_slj.py diff --git a/src/PairEvaluatorShiftedLJ.h b/src/PairEvaluatorShiftedLJ.h deleted file mode 100644 index 96b0abc8..00000000 --- a/src/PairEvaluatorShiftedLJ.h +++ /dev/null @@ -1,132 +0,0 @@ -// Copyright (c) 2018-2020, Michael P. Howard -// Copyright (c) 2021-2022, Auburn University -// This file is part of the azplugins project, released under the Modified BSD License. - -/*! - * \file PairEvaluatorShiftedLJ.h - * \brief Defines the pair force evaluator class for core-shifted Lennard-Jones potential - */ - -#ifndef AZPLUGINS_PAIR_EVALUATOR_SHIFTED_LJ_H_ -#define AZPLUGINS_PAIR_EVALUATOR_SHIFTED_LJ_H_ - -#include "PairEvaluator.h" - -#ifdef NVCC -#define DEVICE __device__ -#else -#define DEVICE -#endif - -namespace azplugins -{ - -namespace detail -{ -//! Class for evaluating the core-shifted Lennard-Jones pair potential -/*! - * PairEvaluatorShiftedLJ evaluates the function: - * \f{eqnarray*}{ - * V(r) = & 4 \varepsilon \left[ \left(\frac{\sigma}{r-\Delta} \right)^12 - - \alpha \left(\frac{\sigma}{r-\Delta} \right)^6 \right] & r < r_{\rm cut} \\ - * = & 0 & r \ge r_{\rm cut} - * \f} - * - * where \f$\varepsilon\f$, \f$\sigma\f$, and \f$\alpha\f$ (default: 1.0) are the standard Lennard-Jones parameters. - * \f$\Delta\f$ is a shift factor for the minimum of the potential, which moves to \f$r_{\rm min} = 2^{1/6} \sigma + \Delta \f$. - * This potential is analogous to the HOOMD shifted Lennard-Jones potential (see EvaluatorPairSLJ), but does - * not rely on reading the diameter to set the shift factor. This can be more convenient for certain systems where - * the shift factor is a parameterization, but not necessarily connected to the diameter of the particle. - * - * The core-shifted Lennard-Jones potential does not need diameter or charge. Three parameters are specified and stored in a - * Scalar3. These are related to the standard Lennard-Jones parameters by: - * - \a lj1 = 4.0 * epsilon * pow(sigma,12.0) - * - \a lj2 = alpha * 4.0 * epsilon * pow(sigma,6.0); - * - \a Delta is the amount to shift the potential minimum by. - */ -class PairEvaluatorShiftedLJ : public PairEvaluator - { - public: - //! Define the parameter type used by this pair potential evaluator - typedef Scalar3 param_type; - - //! Constructor - /*! - * \param _rsq Squared distance between particles - * \param _rcutsq Cutoff radius squared - * \param _params Pair potential parameters, given by typedef above - * - * The functor initializes its members from \a _params. - */ - DEVICE PairEvaluatorShiftedLJ(Scalar _rsq, Scalar _rcutsq, const param_type& _params) - : PairEvaluator(_rsq,_rcutsq), lj1(_params.x), lj2(_params.y), delta(_params.z) - {} - - //! Evaluate the force and energy - /*! - * \param force_divr Holds the computed force divided by r - * \param pair_eng Holds the computed pair energy - * \param energy_shift If true, the potential is shifted to zero at the cutoff - * - * \returns True if the energy calculation occurs - * - * The calculation does not occur if the pair distance is greater than the cutoff - * or if the potential is scaled to zero. - */ - DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& pair_eng, bool energy_shift) - { - if (rsq < rcutsq && lj1 != 0) - { - Scalar rinv = fast::rsqrt(rsq); - Scalar r = Scalar(1.0) / rinv; - - Scalar rmd = r - delta; - Scalar rmdinv = Scalar(1.0) / rmd; - Scalar rmd2inv = rmdinv * rmdinv; - Scalar rmd6inv = rmd2inv * rmd2inv * rmd2inv; - force_divr= rinv * rmdinv * rmd6inv * (Scalar(12.0)*lj1*rmd6inv - Scalar(6.0)*lj2); - - pair_eng = rmd6inv * (lj1*rmd6inv - lj2); - - /* - * Energy shift requires another sqrt call. If this pair becomes - * performance limiting, it could be replaced by caching the shift - * into a Scalar4. - */ - if (energy_shift) - { - Scalar rcutinv = fast::rsqrt(rcutsq); - Scalar rcut = Scalar(1.0) / rcutinv; - - Scalar rcutmd = rcut - delta; - Scalar rcutmdinv = Scalar(1.0) / rcutmd; - Scalar rcutmd2inv = rcutmdinv * rcutmdinv; - Scalar rcutmd6inv = rcutmd2inv * rcutmd2inv * rcutmd2inv; - pair_eng -= rcutmd6inv * (lj1*rcutmd6inv - lj2); - } - return true; - } - else - return false; - } - - #ifndef NVCC - //! Return the name of this potential - static std::string getName() - { - return std::string("slj"); - } - #endif - - protected: - Scalar lj1; //!< lj1 parameter extracted from the params passed to the constructor - Scalar lj2; //!< lj2 parameter extracted from the params passed to the constructor - Scalar delta; //!< shift parameter - }; - -} // end namespace detail -} // end namespace azplugins - -#undef DEVICE - -#endif // AZPLUGINS_PAIR_EVALUATOR_SHIFTED_LJ_H_ diff --git a/src/PairPotentialShiftedLJ.cu b/src/PairPotentialShiftedLJ.cu deleted file mode 100644 index 2b3db64b..00000000 --- a/src/PairPotentialShiftedLJ.cu +++ /dev/null @@ -1,18 +0,0 @@ -// Copyright (c) 2018-2020, Michael P. Howard -// Copyright (c) 2021-2022, Auburn University -// This file is part of the azplugins project, released under the Modified BSD License. - -#include "PairPotentials.cuh" - -namespace azplugins -{ -namespace gpu -{ - -//! Kernel driver for core-shifted Lennard-Jones pair potential -template cudaError_t compute_pair_potential - (const pair_args_t& pair_args, - const typename azplugins::detail::PairEvaluatorShiftedLJ::param_type *d_params); - -} // end namespace gpu -} // end namespace azplugins diff --git a/src/PairPotentials.h b/src/PairPotentials.h index 6658e808..8d506a57 100644 --- a/src/PairPotentials.h +++ b/src/PairPotentials.h @@ -33,7 +33,6 @@ #include "PairEvaluatorAshbaugh.h" #include "PairEvaluatorColloid.h" #include "PairEvaluatorHertz.h" -#include "PairEvaluatorShiftedLJ.h" #include "PairEvaluatorSpline.h" /* diff --git a/src/module.cc b/src/module.cc index ddf2000a..48f174d8 100644 --- a/src/module.cc +++ b/src/module.cc @@ -192,7 +192,6 @@ PYBIND11_MODULE(_azplugins, m) azplugins::detail::export_ashbaugh_params(m); azplugins::detail::export_pair_potential(m, "PairPotentialColloid"); azplugins::detail::export_pair_potential(m, "PairPotentialHertz"); - azplugins::detail::export_pair_potential(m, "PairPotentialShiftedLJ"); azplugins::detail::export_pair_potential(m, "PairPotentialSpline"); /* Anisotropic pair potentials */ diff --git a/src/pair.py b/src/pair.py index cc8a3916..9765d7e9 100644 --- a/src/pair.py +++ b/src/pair.py @@ -11,14 +11,12 @@ ashbaugh colloid hertz - slj spline two_patch_morse .. autoclass:: ashbaugh .. autoclass:: colloid .. autoclass:: hertz -.. autoclass:: slj .. autoclass:: spline .. autoclass:: two_patch_morse @@ -279,97 +277,6 @@ def __init__(self, r_cut, nlist, name=None): def process_coeff(self, coeff): return coeff['epsilon'] -class slj(hoomd.md.pair.pair): - R""" Core-shifted Lennard-Jones potential - - Args: - r_cut (float): Default cutoff radius (in distance units). - nlist (:py:mod:`hoomd.md.nlist`): Neighbor list - name (str): Name of the force instance. - - :py:class:`slj` is Lennard-Jones potential with the core (minimum) shifted by - an amount :math:`\Delta`. The form of the potential is similar to the standard - Lennard-Jones potential - - .. math:: - :nowrap: - - \begin{eqnarray*} - V(r) = & 4 \varepsilon \left[ \left(\frac{\sigma}{r-\Delta} \right)^12 - - \alpha \left(\frac{\sigma}{r-\Delta} \right)^6 \right] & r < r_{\rm cut} \\ - = & 0 & r \ge r_{\rm cut} - \end{eqnarray*} - - Here, :math:`\varepsilon`, :math:`\sigma`, and :math:`\alpha` are the stanard - Lennard-Jones potential parameters, and :math:`\Delta` is the amount the potential - is shifted by. The minimum of the potential :math:`r_{\rm min}` is shifted to - - .. math:: - - r_{\rm min} = 2^{1/6} \sigma + \Delta - - Setting :math:`\Delta = 0` recovers the standard Lennard-Jones potential. - See :py:class:`hoomd.md.pair.pair` for details on how forces are calculated - and the available energy shifting and smoothing modes. - Use :py:meth:`pair_coeff.set ` to set potential coefficients. - - The following coefficients must be set per unique pair of particle types: - - - :math:`\varepsilon` - *epsilon* (in energy units) - - :math:`\sigma` - *sigma* (in distance units) - - :math:`\Delta` - *delta* (in distance units) - - :math:`\alpha` - *alpha* (unitless) - *optional*: defaults to 1.0 - - :math:`r_{\mathrm{cut}}` - *r_cut* (in distance units) - - *optional*: defaults to the global r_cut specified in the pair command - - :math:`r_{\mathrm{on}}`- *r_on* (in distance units) - - *optional*: defaults to the global r_cut specified in the pair command - - Example:: - - nl = hoomd.md.nlist.cell() - slj = azplugins.pair.slj(r_cut=3.0, nlist=nl) - slj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0, delta=1.0) - slj.pair_coeff.set(['A','B'], 'B', epsilon=2.0, sigma=1.0, alpha=0.5, delta=0.0, r_cut=3.0, r_on=2.0) - - .. note:: - Because of the form of the potential, square-root calls are necessary - to evaluate the potential and also to perform energy shifting. This will - incur a corresponding performance hit compared to the standard - Lennard-Jones potential even when :math:`\Delta=0`. - - """ - def __init__(self, r_cut, nlist, name=None): - hoomd.util.print_status_line() - - # initialize the base class - hoomd.md.pair.pair.__init__(self, r_cut, nlist, name) - - # create the c++ mirror class - if not hoomd.context.exec_conf.isCUDAEnabled(): - self.cpp_class = _azplugins.PairPotentialShiftedLJ - else: - self.cpp_class = _azplugins.PairPotentialShiftedLJGPU - self.nlist.cpp_nlist.setStorageMode(_md.NeighborList.storageMode.full) - self.cpp_force = self.cpp_class(hoomd.context.current.system_definition, self.nlist.cpp_nlist, self.name) - - hoomd.context.current.system.addCompute(self.cpp_force, self.force_name) - - # setup the coefficent options - self.required_coeffs = ['epsilon', 'sigma', 'delta', 'alpha'] - self.pair_coeff.set_default_coeff('alpha', 1.0) - - def process_coeff(self, coeff): - epsilon = coeff['epsilon'] - sigma = coeff['sigma'] - delta = coeff['delta'] - alpha = coeff['alpha'] - - lj1 = 4.0 * epsilon * math.pow(sigma, 12.0) - lj2 = alpha * 4.0 * epsilon * math.pow(sigma, 6.0) - return _hoomd.make_scalar3(lj1, lj2, delta) - - - class spline(hoomd.md.pair.pair): R""" Spline potential diff --git a/src/pytest/test_pair_slj.py b/src/pytest/test_pair_slj.py deleted file mode 100644 index 5476b174..00000000 --- a/src/pytest/test_pair_slj.py +++ /dev/null @@ -1,212 +0,0 @@ -# Copyright (c) 2018-2020, Michael P. Howard -# Copyright (c) 2021-2022, Auburn University -# This file is part of the azplugins project, released under the Modified BSD License. - -import hoomd -from hoomd import md -hoomd.context.initialize() -try: - from hoomd import azplugins -except ImportError: - import azplugins -import unittest - -# azplugins.pair.slj -class pair_slj_tests(unittest.TestCase): - def setUp(self): - # raw snapshot is fine, just needs to have the types - snap = hoomd.data.make_snapshot(N=100, box=hoomd.data.boxdim(L=20), particle_types=['A']) - self.s = hoomd.init.read_snapshot(snap) - self.nl = md.nlist.cell() - hoomd.context.current.sorter.set_params(grid=8) - - # basic test of creation - def test(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - slj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0, alpha=1.0, delta=1.0, r_cut=2.5, r_on=2.0) - slj.update_coeffs() - - # test missing epsilon coefficient - def test_set_missing_epsilon(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - slj.pair_coeff.set('A', 'A', sigma=1.0, delta=1.0) - self.assertRaises(RuntimeError, slj.update_coeffs) - - # test missing sigma coefficient - def test_set_missing_sigma(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - slj.pair_coeff.set('A', 'A', epsilon=1.0, delta=1.0) - self.assertRaises(RuntimeError, slj.update_coeffs) - - # test missing delta coefficient - def test_set_missing_delta(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - slj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0) - self.assertRaises(RuntimeError, slj.update_coeffs) - - # test missing type coefficients - def test_missing_AA(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - self.assertRaises(RuntimeError, slj.update_coeffs) - - # test set params - def test_set_params(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - slj.set_params(mode="no_shift") - slj.set_params(mode="shift") - slj.set_params(mode="xplor") - self.assertRaises(RuntimeError, slj.set_params, mode="blah") - - # test default coefficients - def test_default_coeff(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - # (r_cut, and r_on are default) - slj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0, delta=1.0) - slj.update_coeffs() - - # test max rcut - def test_max_rcut(self): - slj = azplugins.pair.slj(r_cut=2.5, nlist = self.nl) - slj.pair_coeff.set('A', 'A', sigma=1.0, epsilon=1.0, delta=1.0) - self.assertAlmostEqual(2.5, slj.get_max_rcut()) - slj.pair_coeff.set('A', 'A', r_cut = 2.0) - self.assertAlmostEqual(2.0, slj.get_max_rcut()) - - # test specific nlist subscription - def test_nlist_subscribe(self): - slj = azplugins.pair.slj(r_cut=2.5, nlist = self.nl) - - slj.pair_coeff.set('A', 'A', sigma = 1.0, epsilon=1.0, delta=1.0) - self.nl.update_rcut() - self.assertAlmostEqual(2.5, self.nl.r_cut.get_pair('A','A')) - - slj.pair_coeff.set('A', 'A', r_cut = 2.0) - self.nl.update_rcut() - self.assertAlmostEqual(2.0, self.nl.r_cut.get_pair('A','A')) - - # test coeff list - def test_coeff_list(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - slj.pair_coeff.set(['A', 'B'], ['A', 'C'], epsilon=1.0, sigma=1.0, delta=1.0, r_cut=2.5, r_on=2.0) - slj.update_coeffs() - - # test adding types - def test_type_add(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - slj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0, delta=1.0) - self.s.particles.types.add('B') - self.assertRaises(RuntimeError, slj.update_coeffs) - slj.pair_coeff.set('A', 'B', epsilon=1.0, sigma=1.0, delta=1.0) - slj.pair_coeff.set('B', 'B', epsilon=1.0, sigma=1.0, delta=1.0) - slj.update_coeffs() - - def tearDown(self): - del self.s, self.nl - hoomd.context.initialize() - -# test the validity of the pair potential -class potential_slj_tests(unittest.TestCase): - def setUp(self): - snap = hoomd.data.make_snapshot(N=2, box=hoomd.data.boxdim(L=20),particle_types=['A']) - if hoomd.comm.get_rank() == 0: - snap.particles.position[0] = (0,0,0) - snap.particles.position[1] = (1.05,0,0) - hoomd.init.read_snapshot(snap) - self.nl = md.nlist.cell() - - # test the calculation of force and potential - def test_potential(self): - slj = azplugins.pair.slj(r_cut=2.0, nlist = self.nl) - - # some set of parameters which will give a potential of zero by shifting - slj.pair_coeff.set('A','A', epsilon=2.0, sigma=0.95, delta=0.1) - slj.set_params(mode="no_shift") - - md.integrate.mode_standard(dt=0) - nve = md.integrate.nve(group = hoomd.group.all()) - hoomd.run(1) - U = 0.0 - F = -50.5263 - f0 = slj.forces[0].force - f1 = slj.forces[1].force - e0 = slj.forces[0].energy - e1 = slj.forces[1].energy - - self.assertAlmostEqual(e0,0.5*U,3) - self.assertAlmostEqual(e1,0.5*U,3) - - self.assertAlmostEqual(f0[0],F,3) - self.assertAlmostEqual(f0[1],0) - self.assertAlmostEqual(f0[2],0) - - self.assertAlmostEqual(f1[0],-F,3) - self.assertAlmostEqual(f1[1],0) - self.assertAlmostEqual(f1[2],0) - - # shift the potential and check the energy (force stays the same) - slj.set_params(mode='shift') - hoomd.run(1) - U -= -0.123046875 - f0 = slj.forces[0].force - f1 = slj.forces[1].force - e0 = slj.forces[0].energy - e1 = slj.forces[1].energy - - self.assertAlmostEqual(e0,0.5*U,3) - self.assertAlmostEqual(e1,0.5*U,3) - - self.assertAlmostEqual(f0[0],F,3) - self.assertAlmostEqual(f0[1],0) - self.assertAlmostEqual(f0[2],0) - - self.assertAlmostEqual(f1[0],-F,3) - self.assertAlmostEqual(f1[1],0) - self.assertAlmostEqual(f1[2],0) - - # test alpha parameter in potential - def test_alpha(self): - slj = azplugins.pair.slj(r_cut=3.0, nlist = self.nl) - slj.pair_coeff.set('A','A', epsilon=2.0, sigma=0.95, delta=0.1, alpha=0.5) - slj.set_params(mode="no_shift") - - md.integrate.mode_standard(dt=0) - nve = md.integrate.nve(group = hoomd.group.all()) - hoomd.run(1) - - U = 4.0 - F = 75.78947368421039 - self.assertAlmostEqual(slj.forces[0].energy,0.5*U,3) - self.assertAlmostEqual(slj.forces[1].energy,0.5*U,3) - self.assertAlmostEqual(slj.forces[0].force[0], -F,3) - self.assertAlmostEqual(slj.forces[1].force[0], F,3) - - # test the cases where the potential should be zero - def test_noninteract(self): - slj = azplugins.pair.slj(r_cut=1.0, nlist = self.nl) - - # outside cutoff - slj.pair_coeff.set('A','A', epsilon=1.0, sigma=0.95, delta=0.5) - slj.set_params(mode="no_shift") - - md.integrate.mode_standard(dt=0) - nve = md.integrate.nve(group = hoomd.group.all()) - hoomd.run(1) - self.assertAlmostEqual(slj.forces[0].energy, 0) - self.assertAlmostEqual(slj.forces[1].energy, 0) - self.assertAlmostEqual(slj.forces[0].force[0], 0) - self.assertAlmostEqual(slj.forces[1].force[0], 0) - - # inside cutoff but epsilon = 0 - slj.pair_coeff.set('A','A', epsilon=0.0, r_cut=3.0) - hoomd.run(1) - self.assertAlmostEqual(slj.forces[0].energy, 0) - self.assertAlmostEqual(slj.forces[1].energy, 0) - self.assertAlmostEqual(slj.forces[0].force[0], 0) - self.assertAlmostEqual(slj.forces[1].force[0], 0) - - def tearDown(self): - del self.nl - hoomd.context.initialize() - -if __name__ == '__main__': - unittest.main(argv = ['test.py', '-v'])