diff --git a/hoomd/md/CMakeLists.txt b/hoomd/md/CMakeLists.txt index 95f24e2214..aefbd22e69 100644 --- a/hoomd/md/CMakeLists.txt +++ b/hoomd/md/CMakeLists.txt @@ -87,6 +87,7 @@ set(_md_headers ActiveForceComputeGPU.h EvaluatorSpecialPairLJ.h EvaluatorSpecialPairCoulomb.h EvaluatorExternalElectricField.h + EvaluatorExternalMagneticField.h EvaluatorExternalPeriodic.h EvaluatorPairALJ.h EvaluatorPairBuckingham.h @@ -426,7 +427,7 @@ foreach(_evaluator ${_triplets}) endif() endforeach() -set(_external_evaluators Periodic ElectricField) +set(_external_evaluators Periodic ElectricField MagneticField) foreach(_evaluator ${_external_evaluators}) configure_file(export_PotentialExternal.cc.inc diff --git a/hoomd/md/EvaluatorExternalElectricField.h b/hoomd/md/EvaluatorExternalElectricField.h index 46a649437d..c04688874c 100644 --- a/hoomd/md/EvaluatorExternalElectricField.h +++ b/hoomd/md/EvaluatorExternalElectricField.h @@ -10,6 +10,7 @@ #include "hoomd/BoxDim.h" #include "hoomd/HOOMDMath.h" +#include "hoomd/VectorMath.h" #include /*! \file EvaluatorExternalElectricField.h @@ -75,6 +76,7 @@ class EvaluatorExternalElectricField \param params per-type parameters of external potential */ DEVICE EvaluatorExternalElectricField(Scalar3 X, + quat q, const BoxDim& box, const param_type& params, const field_type& field) @@ -82,6 +84,11 @@ class EvaluatorExternalElectricField { } + DEVICE static bool isAnisotropic() + { + return false; + } + //! ExternalElectricField needs charges DEVICE static bool needsCharge() { @@ -106,10 +113,12 @@ class EvaluatorExternalElectricField //! Evaluate the force, energy and virial /*! \param F force vector + \param T torque vector \param energy value of the energy \param virial array of six scalars for the upper triangular virial tensor */ - DEVICE void evalForceEnergyAndVirial(Scalar3& F, Scalar& energy, Scalar* virial) + DEVICE void + evalForceTorqueEnergyAndVirial(Scalar3& F, Scalar3& T, Scalar& energy, Scalar* virial) { F = m_qi * m_E; energy = -m_qi * dot(m_E, m_pos); @@ -120,6 +129,10 @@ class EvaluatorExternalElectricField virial[3] = F.y * m_pos.y; virial[4] = F.y * m_pos.z; virial[5] = F.z * m_pos.z; + + T.x = Scalar(0.0); + T.y = Scalar(0.0); + T.z = Scalar(0.0); } #ifndef __HIPCC__ diff --git a/hoomd/md/EvaluatorExternalMagneticField.h b/hoomd/md/EvaluatorExternalMagneticField.h new file mode 100644 index 0000000000..d3fb39e211 --- /dev/null +++ b/hoomd/md/EvaluatorExternalMagneticField.h @@ -0,0 +1,164 @@ +// Copyright (c) 2009-2023 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#ifndef __EVALUATOR_EXTERNAL_MAGNETIC_FIELD_H__ +#define __EVALUATOR_EXTERNAL_MAGNETIC_FIELD_H__ + +#ifndef __HIPCC__ +#include +#endif + +#include "hoomd/BoxDim.h" +#include "hoomd/HOOMDMath.h" +#include "hoomd/VectorMath.h" +#include + +/*! \file EvaluatorExternalMagneticField.h + \brief Defines the external potential evaluator to induce a magnetic field +*/ + +// need to declare these class methods with __device__ qualifiers when building in nvcc +// DEVICE is __host__ __device__ when included in nvcc and blank when included into the host +// compiler +#ifdef __HIPCC__ +#define DEVICE __device__ +#else +#define DEVICE +#endif + +namespace hoomd + { +namespace md + { +//! Class for evaluating an magnetic field +/*! General Overview + The external potential \f$V(\theta) \f$ is implemented using the following formula: + + \f[ + V(\theta}) = - \vec{B} \cdot \vec{n}_i(\theta) + \f] + + where \f$B\f$ is the strength of the magnetic field and \f$\vec{n}_i\f$ is the magnetic moment + of particle i. +*/ +class EvaluatorExternalMagneticField + { + public: + //! type of parameters this external potential accepts + struct param_type + { + Scalar3 B; + Scalar3 mu; + +#ifndef __HIPCC__ + param_type() : B(make_scalar3(0, 0, 0)), mu(make_scalar3(0, 0, 0)) { } + + param_type(pybind11::dict params) + { + pybind11::tuple py_B = params["B"]; + B.x = pybind11::cast(py_B[0]); + B.y = pybind11::cast(py_B[1]); + B.z = pybind11::cast(py_B[2]); + pybind11::tuple py_mu = params["mu"]; + mu.x = pybind11::cast(py_mu[0]); + mu.y = pybind11::cast(py_mu[1]); + mu.z = pybind11::cast(py_mu[2]); + } + + pybind11::dict toPython() + { + pybind11::dict d; + d["B"] = pybind11::make_tuple(B.x, B.y, B.z); + d["mu"] = pybind11::make_tuple(mu.x, mu.y, mu.z); + return d; + } +#endif // ifndef __HIPCC__ + } __attribute__((aligned(16))); + + typedef void* field_type; + + //! Constructs the external field evaluator + /*! \param X position of particle + \param box box dimensions + \param params per-type parameters of external potential + */ + DEVICE EvaluatorExternalMagneticField(Scalar3 X, + quat q, + const BoxDim& box, + const param_type& params, + const field_type& field) + : m_q(q), m_B(params.B), m_mu(params.mu) + { + } + + DEVICE static bool isAnisotropic() + { + return true; + } + + //! ExternalMagneticField needs charges + DEVICE static bool needsCharge() + { + return false; + } + + //! Accept the optional charge value + /*! \param qi Charge of particle i + */ + DEVICE void setCharge(Scalar qi) { } + + //! Declares additional virial contributions are needed for the external field + /*! No contribution + */ + DEVICE static bool requestFieldVirialTerm() + { + return false; + } + + //! Evaluate the force, energy and virial + /*! \param F force vector + \param T torque vector + \param energy value of the energy + \param virial array of six scalars for the upper triangular virial tensor + */ + DEVICE void + evalForceTorqueEnergyAndVirial(Scalar3& F, Scalar3& T, Scalar& energy, Scalar* virial) + { + vec3 dir = rotate(m_q, m_mu); + + vec3 T_vec = cross(dir, m_B); + + T.x = T_vec.x; + T.y = T_vec.y; + T.z = T_vec.z; + + energy = -dot(dir, m_B); + + F.x = Scalar(0.0); + F.y = Scalar(0.0); + F.z = Scalar(0.0); + + for (unsigned int i = 0; i < 6; i++) + virial[i] = Scalar(0.0); + } + +#ifndef __HIPCC__ + //! Get the name of this potential + /*! \returns The potential name. + */ + static std::string getName() + { + return std::string("b_field"); + } +#endif + + protected: + quat m_q; //!< Particle orientation + vec3 m_B; //!< Magnetic field vector (box frame). + vec3 m_mu; //!< Magnetic dipole moment (particle frame). + }; + + } // end namespace md + } // end namespace hoomd + +#endif // __EVALUATOR_EXTERNAL_LAMELLAR_H__ diff --git a/hoomd/md/EvaluatorExternalPeriodic.h b/hoomd/md/EvaluatorExternalPeriodic.h index 1d3a08688b..ddc4b24bc6 100644 --- a/hoomd/md/EvaluatorExternalPeriodic.h +++ b/hoomd/md/EvaluatorExternalPeriodic.h @@ -10,6 +10,7 @@ #include "hoomd/BoxDim.h" #include "hoomd/HOOMDMath.h" +#include "hoomd/VectorMath.h" #include /*! \file EvaluatorExternalPeriodic.h @@ -89,6 +90,7 @@ class EvaluatorExternalPeriodic \param params per-type parameters of external potential */ DEVICE EvaluatorExternalPeriodic(Scalar3 X, + quat q, const BoxDim& box, const param_type& params, const field_type& field) @@ -97,6 +99,11 @@ class EvaluatorExternalPeriodic { } + DEVICE static bool isAnisotropic() + { + return false; + } + //! External Periodic doesn't need charges DEVICE static bool needsCharge() { @@ -117,10 +124,12 @@ class EvaluatorExternalPeriodic //! Evaluate the force, energy and virial /*! \param F force vector + \param T torque vector \param energy value of the energy \param virial array of six scalars for the upper triangular virial tensor */ - DEVICE void evalForceEnergyAndVirial(Scalar3& F, Scalar& energy, Scalar* virial) + DEVICE void + evalForceTorqueEnergyAndVirial(Scalar3& F, Scalar3& T, Scalar& energy, Scalar* virial) { Scalar3 a2 = make_scalar3(0, 0, 0); Scalar3 a3 = make_scalar3(0, 0, 0); @@ -128,6 +137,11 @@ class EvaluatorExternalPeriodic F.x = Scalar(0.0); F.y = Scalar(0.0); F.z = Scalar(0.0); + + T.x = Scalar(0.0); + T.y = Scalar(0.0); + T.z = Scalar(0.0); + energy = Scalar(0.0); // For this potential, since it uses scaled positions, the virial is always zero. diff --git a/hoomd/md/EvaluatorWalls.h b/hoomd/md/EvaluatorWalls.h index 95b390caf8..4535c54e92 100644 --- a/hoomd/md/EvaluatorWalls.h +++ b/hoomd/md/EvaluatorWalls.h @@ -134,11 +134,20 @@ template class EvaluatorWalls typedef wall_type field_type; //! Constructs the external wall potential evaluator - DEVICE EvaluatorWalls(Scalar3 pos, const BoxDim& box, const param_type& p, const field_type& f) + DEVICE EvaluatorWalls(Scalar3 pos, + quat q, + const BoxDim& box, + const param_type& p, + const field_type& f) : m_pos(pos), m_field(f), m_params(p) { } + DEVICE static bool isAnisotropic() + { + return false; + } + //! Charges not supported by walls evals DEVICE static bool needsCharge() { @@ -229,11 +238,17 @@ template class EvaluatorWalls } //! Generates force and energy from standard evaluators using wall geometry functions - DEVICE void evalForceEnergyAndVirial(Scalar3& F, Scalar& energy, Scalar* virial) + DEVICE void + evalForceTorqueEnergyAndVirial(Scalar3& F, Scalar3& T, Scalar& energy, Scalar* virial) { F.x = Scalar(0.0); F.y = Scalar(0.0); F.z = Scalar(0.0); + + T.x = Scalar(0.0); + T.y = Scalar(0.0); + T.z = Scalar(0.0); + energy = Scalar(0.0); // initialize virial for (unsigned int i = 0; i < 6; i++) diff --git a/hoomd/md/PotentialExternal.h b/hoomd/md/PotentialExternal.h index 82d181988c..aac631a0a3 100644 --- a/hoomd/md/PotentialExternal.h +++ b/hoomd/md/PotentialExternal.h @@ -4,6 +4,7 @@ #include "hoomd/ForceCompute.h" #include "hoomd/GPUArray.h" #include "hoomd/GlobalArray.h" +#include "hoomd/VectorMath.h" #include "hoomd/managed_allocator.h" #include "hoomd/md/EvaluatorExternalPeriodic.h" #include @@ -43,6 +44,8 @@ template class PotentialExternal : public ForceCompute typedef typename evaluator::param_type param_type; typedef typename evaluator::field_type field_type; + bool isAnisotropic(); + //! Sets parameters of the evaluator pybind11::object getParams(std::string type); @@ -99,8 +102,13 @@ template void PotentialExternal::computeForces(uint6 assert(m_pdata); // access the particle data arrays ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); + ArrayHandle h_orientation(m_pdata->getOrientationArray(), + access_location::host, + access_mode::read); ArrayHandle h_force(m_force, access_location::host, access_mode::overwrite); + ArrayHandle h_torque(m_torque, access_location::host, access_mode::overwrite); + ArrayHandle h_virial(m_virial, access_location::host, access_mode::overwrite); ArrayHandle h_charge(m_pdata->getCharges(), access_location::host, access_mode::read); @@ -112,10 +120,12 @@ template void PotentialExternal::computeForces(uint6 // Zero data for force calculation. memset((void*)h_force.data, 0, sizeof(Scalar4) * m_force.getNumElements()); + memset((void*)h_torque.data, 0, sizeof(Scalar4) * m_torque.getNumElements()); memset((void*)h_virial.data, 0, sizeof(Scalar) * m_virial.getNumElements()); // there are enough other checks on the input data: but it doesn't hurt to be safe assert(h_force.data); + assert(h_torque.data); assert(h_virial.data); // for each of the particles @@ -124,18 +134,19 @@ template void PotentialExternal::computeForces(uint6 // get the current particle properties Scalar3 X = make_scalar3(h_pos.data[idx].x, h_pos.data[idx].y, h_pos.data[idx].z); unsigned int type = __scalar_as_int(h_pos.data[idx].w); - Scalar3 F; + quat q(h_orientation.data[idx]); + Scalar3 F, T; Scalar energy; Scalar virial[6]; - evaluator eval(X, box, h_params.data[type], *m_field); + evaluator eval(X, q, box, h_params.data[type], *m_field); if (evaluator::needsCharge()) { Scalar qi = h_charge.data[idx]; eval.setCharge(qi); } - eval.evalForceEnergyAndVirial(F, energy, virial); + eval.evalForceTorqueEnergyAndVirial(F, T, energy, virial); // apply the constraint force h_force.data[idx].x = F.x; @@ -144,6 +155,10 @@ template void PotentialExternal::computeForces(uint6 h_force.data[idx].w = energy; for (int k = 0; k < 6; k++) h_virial.data[k * m_virial_pitch + idx] = virial[k]; + + h_torque.data[idx].x = T.x; + h_torque.data[idx].y = T.y; + h_torque.data[idx].z = T.z; } } @@ -156,6 +171,13 @@ void PotentialExternal::validateType(unsigned int type, std::string a } } +//! Returns true if this ForceCompute requires anisotropic integration +template bool PotentialExternal::isAnisotropic() + { + // by default, only translational degrees of freedom are integrated + return evaluator::isAnisotropic(); + } + //! Set the parameters for this potential /*! \param type type for which to set parameters \param params value of parameters diff --git a/hoomd/md/PotentialExternalGPU.cuh b/hoomd/md/PotentialExternalGPU.cuh index 4ee4a7a108..c33495c534 100644 --- a/hoomd/md/PotentialExternalGPU.cuh +++ b/hoomd/md/PotentialExternalGPU.cuh @@ -28,23 +28,28 @@ struct external_potential_args_t { //! Construct a external_potential_args_t external_potential_args_t(Scalar4* _d_force, + Scalar4* _d_torque, Scalar* _d_virial, const size_t _virial_pitch, const unsigned int _N, const Scalar4* _d_pos, + const Scalar4* _d_orientation, const Scalar* _d_charge, const BoxDim& _box, const unsigned int _block_size, const hipDeviceProp_t& _devprop) - : d_force(_d_force), d_virial(_d_virial), virial_pitch(_virial_pitch), box(_box), N(_N), - d_pos(_d_pos), d_charge(_d_charge), block_size(_block_size), devprop(_devprop) {}; + : d_force(_d_force), d_torque(_d_torque), d_virial(_d_virial), virial_pitch(_virial_pitch), + box(_box), N(_N), d_pos(_d_pos), d_orientation(_d_orientation), d_charge(_d_charge), + block_size(_block_size), devprop(_devprop) {}; Scalar4* d_force; //!< Force to write out + Scalar4* d_torque; //!< Torque to write out Scalar* d_virial; //!< Virial to write out const size_t virial_pitch; //!< The pitch of the 2D array of virial matrix elements const BoxDim box; //!< Simulation box in GPU format const unsigned int N; //!< Number of particles const Scalar4* d_pos; //!< Device array of particle positions + const Scalar4* d_orientation; //!< Device array of particle orientations const Scalar* d_charge; //!< particle charges const unsigned int block_size; //!< Block size to execute const hipDeviceProp_t& devprop; //!< Device properties @@ -69,20 +74,24 @@ hipError_t __attribute__((visibility("default"))) gpu_compute_potential_external the potentials and forces for each particle is handled via the template class \a evaluator. \param d_force Device memory to write computed forces + \param d_torque Device memory to write computed torques \param d_virial Device memory to write computed virials \param virial_pitch pitch of 2D virial array \param N number of particles \param d_pos device array of particle positions + \param d_orientation device array of particle orientations \param box Box dimensions used to implement periodic boundary conditions \param params per-type array of parameters for the potential */ template __global__ void gpu_compute_external_forces_kernel(Scalar4* d_force, + Scalar4* d_torque, Scalar* d_virial, const size_t virial_pitch, const unsigned int N, const Scalar4* d_pos, + const Scalar4* d_orientation, const Scalar* d_charge, const BoxDim box, const typename evaluator::param_type* params, @@ -126,6 +135,7 @@ __global__ void gpu_compute_external_forces_kernel(Scalar4* d_force, // initialize the force to 0 Scalar3 force = make_scalar3(Scalar(0.0), Scalar(0.0), Scalar(0.0)); + Scalar3 torque = make_scalar3(Scalar(0.0), Scalar(0.0), Scalar(0.0)); Scalar virial[6]; for (unsigned int k = 0; k < 6; k++) virial[k] = Scalar(0.0); @@ -133,12 +143,13 @@ __global__ void gpu_compute_external_forces_kernel(Scalar4* d_force, unsigned int typei = __scalar_as_int(posi.w); Scalar3 Xi = make_scalar3(posi.x, posi.y, posi.z); - evaluator eval(Xi, box, params[typei], field); + quat q(d_orientation[idx]); + evaluator eval(Xi, q, box, params[typei], field); if (evaluator::needsCharge()) eval.setCharge(qi); - eval.evalForceEnergyAndVirial(force, energy, virial); + eval.evalForceTorqueEnergyAndVirial(force, torque, energy, virial); // now that the force calculation is complete, write out the result) d_force[idx].x = force.x; @@ -148,6 +159,10 @@ __global__ void gpu_compute_external_forces_kernel(Scalar4* d_force, for (unsigned int k = 0; k < 6; k++) d_virial[k * virial_pitch + idx] = virial[k]; + + d_torque[idx].x = torque.x; + d_torque[idx].y = torque.y; + d_torque[idx].z = torque.z; } /*! @@ -187,10 +202,12 @@ hipError_t gpu_compute_potential_external_forces( bytes, 0, external_potential_args.d_force, + external_potential_args.d_torque, external_potential_args.d_virial, external_potential_args.virial_pitch, external_potential_args.N, external_potential_args.d_pos, + external_potential_args.d_orientation, external_potential_args.d_charge, external_potential_args.box, d_params, diff --git a/hoomd/md/PotentialExternalGPU.h b/hoomd/md/PotentialExternalGPU.h index fd7c36f331..8d4ecf4c3c 100644 --- a/hoomd/md/PotentialExternalGPU.h +++ b/hoomd/md/PotentialExternalGPU.h @@ -61,6 +61,9 @@ template void PotentialExternalGPU::computeForces(ui ArrayHandle d_pos(this->m_pdata->getPositions(), access_location::device, access_mode::read); + ArrayHandle d_orientation(this->m_pdata->getOrientationArray(), + access_location::device, + access_mode::read); ArrayHandle d_charge(this->m_pdata->getCharges(), access_location::device, access_mode::read); @@ -68,6 +71,7 @@ template void PotentialExternalGPU::computeForces(ui const BoxDim box = this->m_pdata->getGlobalBox(); ArrayHandle d_force(this->m_force, access_location::device, access_mode::overwrite); + ArrayHandle d_torque(this->m_torque, access_location::device, access_mode::overwrite); ArrayHandle d_virial(this->m_virial, access_location::device, access_mode::overwrite); ArrayHandle d_params(this->m_params, access_location::device, @@ -76,10 +80,12 @@ template void PotentialExternalGPU::computeForces(ui m_tuner->begin(); kernel::gpu_compute_potential_external_forces( kernel::external_potential_args_t(d_force.data, + d_torque.data, d_virial.data, this->m_virial.getPitch(), this->m_pdata->getN(), d_pos.data, + d_orientation.data, d_charge.data, box, m_tuner->getParam()[0], diff --git a/hoomd/md/external/field.py b/hoomd/md/external/field.py index a2f4a68459..b053588444 100644 --- a/hoomd/md/external/field.py +++ b/hoomd/md/external/field.py @@ -1,7 +1,13 @@ # Copyright (c) 2009-2023 The Regents of the University of Michigan. # Part of HOOMD-blue, released under the BSD 3-Clause License. -"""External field forces.""" +"""External field forces. + +.. invisible-code-block: python + + simulation = hoomd.util.make_example_simulation() + simulation.operations.integrator = hoomd.md.Integrator(dt=0.001) +""" import hoomd from hoomd.md import _md @@ -67,12 +73,14 @@ class Periodic(Field): Type: `TypeParameter` [``particle_type``, `dict`] - Example:: + .. rubric:: Example: - # Apply a periodic composition modulation along the first lattice vector - periodic = external.field.Periodic() + .. code-block:: python + + periodic = hoomd.md.external.field.Periodic() periodic.params['A'] = dict(A=1.0, i=0, w=0.02, p=3) periodic.params['B'] = dict(A=-1.0, i=0, w=0.02, p=3) + simulation.operations.integrator.forces = [periodic] """ _cpp_class_name = "PotentialExternalPeriodic" @@ -86,8 +94,8 @@ def __init__(self): class Electric(Field): """Electric field force. - `Electric` computes forces, and virials, and energies on all particles in - the in the simulation state with consistent with: + `Electric` computes forces, virials, and energies on all particles in + the simulation state which are consistent with: .. math:: @@ -107,11 +115,13 @@ class Electric(Field): Type: `TypeParameter` [``particle_type``, `tuple` [`float`, `float`, `float`]] - Example:: + .. rubric:: Example: + + .. code-block:: python - # Apply an electric field in the x-direction - e_field = external.field.Electric() - e_field.E['A'] = (1, 0, 0) + electric = hoomd.md.external.field.Electric() + electric.E['A'] = (1, 0, 0) + simulation.operations.integrator.forces = [electric] """ _cpp_class_name = "PotentialExternalElectricField" @@ -120,3 +130,52 @@ def __init__(self): 'E', 'particle_types', TypeParameterDict((float, float, float), len_keys=1)) self._add_typeparam(params) + + +class Magnetic(Field): + """Magnetic field torque on a magnetic dipole. + + `Magnetic` computes torques and energies on all particles in the simulation + state which are consistent with: + + .. math:: + + U_i = -\\vec{\\mu}_i \\cdot \\vec{B} + + + where :math:`\\vec{\\mu}_i` is the magnetic dipole moment of particle + :math:`i` and :math:`\\vec{B}` is the field vector. + + .. py:attribute:: params + + The `Magnetic` external potential parameters. The dictionary has the + following keys: + + * ``B`` (`tuple` [`float`, `float`, `float`], **required**) - The + magnetic field vector in the global reference frame + :math:`[\\mathrm{energy} \\cdot \\mathrm{time} \\cdot + \\mathrm{charge}^{-1} \\cdot \\mathrm{length}^{-2} ]`. + * ``mu`` (`tuple` [`float`, `float`, `float`], **required**) - The + magnetic dipole moment of the particles type in the particle + reference frame :math:`[\\mathrm{charge} \\cdot \\mathrm{length}^2 + \\cdot \\mathrm{time}^{-1}]`. + + Type: `TypeParameter` [``particle_type``, `dict`] + + .. rubric:: Example: + + .. code-block:: python + + magnetic = hoomd.md.external.field.Magnetic() + magnetic.params['A'] = dict(B=(1.0,0.0,0.0), mu=(1.0,0.0,0.0)) + simulation.operations.integrator.forces = [magnetic] + """ + _cpp_class_name = "PotentialExternalMagneticField" + + def __init__(self): + params = TypeParameter( + 'params', 'particle_types', + TypeParameterDict(B=(float, float, float), + mu=(float, float, float), + len_keys=1)) + self._add_typeparam(params) diff --git a/hoomd/md/module-md.cc b/hoomd/md/module-md.cc index f3e574ffed..ed9d3ea0d0 100644 --- a/hoomd/md/module-md.cc +++ b/hoomd/md/module-md.cc @@ -88,6 +88,7 @@ void export_PotentialRevCross(pybind11::module& m); void export_PotentialExternalPeriodic(pybind11::module& m); void export_PotentialExternalElectricField(pybind11::module& m); +void export_PotentialExternalMagneticField(pybind11::module& m); void export_PotentialExternalWallLJ(pybind11::module& m); void export_PotentialExternalWallYukawa(pybind11::module& m); @@ -230,6 +231,7 @@ void export_PotentialRevCrossGPU(pybind11::module& m); void export_PotentialExternalPeriodicGPU(pybind11::module& m); void export_PotentialExternalElectricFieldGPU(pybind11::module& m); +void export_PotentialExternalMagneticFieldGPU(pybind11::module& m); void export_PotentialExternalWallLJGPU(pybind11::module& m); void export_PotentialExternalWallYukawaGPU(pybind11::module& m); @@ -369,6 +371,7 @@ PYBIND11_MODULE(_md, m) export_PotentialExternalPeriodic(m); export_PotentialExternalElectricField(m); + export_PotentialExternalMagneticField(m); export_wall_data(m); export_wall_field(m); @@ -457,6 +460,7 @@ PYBIND11_MODULE(_md, m) export_ConstantForceComputeGPU(m); export_PotentialExternalPeriodicGPU(m); export_PotentialExternalElectricFieldGPU(m); + export_PotentialExternalMagneticFieldGPU(m); export_PotentialExternalWallLJGPU(m); export_PotentialExternalWallYukawaGPU(m); diff --git a/hoomd/md/pytest/test_external.py b/hoomd/md/pytest/test_external.py index e4a655b26e..cc980e7d58 100644 --- a/hoomd/md/pytest/test_external.py +++ b/hoomd/md/pytest/test_external.py @@ -35,7 +35,10 @@ def _evaluate_periodic(snapshot, params): forces *= 1 - (np.tanh( np.cos(p * np.dot(positions, b)) / (2 * np.pi * p * w)))**2 forces = np.outer(forces, b) - return forces, energies + torques = [ + [0, 0, 0], + ] * len(positions) + return forces, torques, energies def _evaluate_electric(snapshot, params): @@ -45,7 +48,24 @@ def _evaluate_electric(snapshot, params): E_field = params energies = -charges * np.dot(positions, E_field) forces = np.outer(charges, E_field) - return forces, energies + torques = [ + [0, 0, 0], + ] * len(positions) + return forces, torques, energies + + +def _evaluate_magnetic(snapshot, params): + """Evaluate force and energy in python for MagneticField.""" + positions = snapshot.particles.position + N = len(positions) + B_field = params['B'] + b_moment = params['mu'] + energies = np.repeat(-np.dot(b_moment, B_field), N) + torques = np.tile(np.cross(b_moment, B_field), (N, 1)) + forces = [ + [0, 0, 0], + ] * N + return forces, torques, energies def _external_params(): @@ -60,6 +80,11 @@ def _external_params(): (1, 0, 0), (0, 2, 0), ]), _evaluate_electric)) + list_ext_params.append((hoomd.md.external.field.Magnetic, "params", + list([ + dict(B=(0, 2, -11.5), mu=(1, 2, 3)), + dict(B=(1, 0, 1), mu=(1, 1, 1)) + ]), _evaluate_magnetic)) return list_ext_params @@ -131,19 +156,22 @@ def test_forces_and_energies(simulation_factory, lattice_snapshot_factory, # test energies new_snap = sim.state.get_snapshot() forces = sim.operations.integrator.forces[0].forces + torques = sim.operations.integrator.forces[0].torques energies = sim.operations.integrator.forces[0].energies if new_snap.communicator.rank == 0: - expected_forces, expected_energies = evaluator(new_snap, param) + expected_forces, expected_torques, expected_energies = evaluator( + new_snap, param) # Set atol as the energies and forces very close to 0. # It would be better to run a test that applies appreciable forces # and energies. np.testing.assert_allclose(expected_forces, forces, atol=1e-5) + np.testing.assert_allclose(expected_torques, torques, atol=1e-5) np.testing.assert_allclose(expected_energies, energies, atol=1e-5) # Test Logging _potential_cls = (md.external.field.Field, md.external.field.Periodic, - md.external.field.Electric) + md.external.field.Electric, md.external.field.Magnetic) @pytest.mark.parametrize( diff --git a/sphinx-doc/module-md-external-field.rst b/sphinx-doc/module-md-external-field.rst index ed27eeffbd..6d7ba6ff57 100644 --- a/sphinx-doc/module-md-external-field.rst +++ b/sphinx-doc/module-md-external-field.rst @@ -13,6 +13,7 @@ md.external.field Field Electric + Magnetic Periodic .. rubric:: Details @@ -21,4 +22,5 @@ md.external.field :synopsis: External field potentials. :members: Field, Electric, + Magnetic, Periodic