From 371cb7cd7ba0350c8c38e1015563b6066e80778b Mon Sep 17 00:00:00 2001 From: Mikael Lund Date: Tue, 2 Jul 2024 11:18:57 +0200 Subject: [PATCH] Make atomic displacement parameters optional (#451) * Add get_optional json function * Make dp and dprot optional * Add default dp and dprot to atomtransrot move * Clang format * Update schema.yml with dp/dprot --- docs/_docs/moves.md | 6 +++++- docs/_docs/topology.md | 4 ++-- docs/schema.yml | 2 ++ src/atomdata.cpp | 29 ++++++++++++++++++++++------- src/atomdata.h | 31 ++++++++++++++++--------------- src/core.h | 13 +++++++++++++ src/move.cpp | 8 ++++++-- src/move.h | 4 +++- 8 files changed, 69 insertions(+), 28 deletions(-) diff --git a/docs/_docs/moves.md b/docs/_docs/moves.md index 647d4d204..bc29d6d6d 100644 --- a/docs/_docs/moves.md +++ b/docs/_docs/moves.md @@ -77,11 +77,15 @@ Upon MC movement, the mean squared displacement will be tracked. `molecule` | Molecule name to operate on `dir=[1,1,1]` | Translational directions `energy_resolution` | If set to a non-zero value (kT), an energy histogram will be generated. +`dp=0` | Default translational displacement parameter (Å) +`dprot=0` | Default rotational displacement parameter (radians) As `moltransrot` but instead of operating on the molecular mass center, this translates -and rotates individual atoms in the group. The repeat is set to the number of atoms in the specified group and the +and rotates individual atoms in the group. +The repeat is set to the number of atoms in the specified group and the displacement parameters `dp` and `dprot` for the individual atoms are taken from the atom properties defined in the [topology](topology). +If `dp` and `dprot` are not defined for an atom, the default values for the move are used. Atomic _rotation_ affects only anisotropic particles such as dipoles, spherocylinders, quadrupoles etc. An energy histogram of each participating species will be written to disk if the `energy_resolution` diff --git a/docs/_docs/topology.md b/docs/_docs/topology.md index d7de2b34d..2b27b93e5 100644 --- a/docs/_docs/topology.md +++ b/docs/_docs/topology.md @@ -61,8 +61,8 @@ Atoms are the smallest possible particle entities with properties defined below. `activity=0` | Chemical activity for grand canonical MC [mol/l] `pactivity` | −log10 of chemical activity (will be converted to activity) `alphax=0` | Excess polarizability (unit-less) -`dp=0` | Translational displacement parameter [Å] -`dprot=0` | Rotational displacement parameter [radians] +`dp` | Translational displacement parameter [Å] (optional) +`dprot` | Rotational displacement parameter [radians] (optional) `eps=0` | Lennard-Jones/WCA energy parameter [kJ/mol] `mu=[0,0,0]` | Dipole moment vector [eÅ] `mulen=|mu|` | Dipole moment scalar [eÅ] diff --git a/docs/schema.yml b/docs/schema.yml index e984457b6..b51230367 100644 --- a/docs/schema.yml +++ b/docs/schema.yml @@ -775,6 +775,8 @@ properties: minItems: 3 maxItems: 3 default: [1,1,1] + dp: {type: number, minimum: 0.0, description: "default translational displacement", default: 0.0} + dprot: {type: number, minimum: 0.0, description: "default rotational displacement", default: 0.0} required: [molecule] additionalProperties: false type: object diff --git a/src/atomdata.cpp b/src/atomdata.cpp index c7c01cb9e..1e036e3b8 100644 --- a/src/atomdata.cpp +++ b/src/atomdata.cpp @@ -93,14 +93,18 @@ void to_json(json& j, const AtomData& a) {"alphax", a.alphax}, {"mw", a.mw}, {"q", a.charge}, - {"dp", a.dp / 1.0_angstrom}, - {"dprot", a.dprot / 1.0_rad}, {"tension", a.tension * 1.0_angstrom * 1.0_angstrom / 1.0_kJmol}, {"tfe", a.tfe * 1.0_angstrom * 1.0_angstrom * 1.0_molar / 1.0_kJmol}, {"mu", a.mu}, {"mulen", a.mulen}, {"psc", a.sphero_cylinder}, {"id", a.id()}}; + if (a.dp.has_value()) { + _j["dp"] = a.dp.value() / 1.0_angstrom; + } + if (a.dprot.has_value()) { + _j["dprot"] = a.dprot.value() / 1.0_rad; + } to_json(_j, a.interaction); // append other interactions if (a.hydrophobic) { _j["hydrophobic"] = a.hydrophobic; @@ -110,6 +114,21 @@ void to_json(json& j, const AtomData& a) } } +/// Handles optional translational and rotational displacement +void set_dp_and_dprot(const json& j, AtomData& atomdata) +{ + // todo: use `std::optional::and_then()` when C++23 is available + if (const auto dp = get_optional(j, "dp")) { + atomdata.dp = dp.value() * 1.0_angstrom; + } + if (const auto dprot = get_optional(j, "dprot")) { + atomdata.dprot = dprot.value() * 1.0_rad; + if (std::fabs(atomdata.dprot.value()) > 2.0 * pc::pi) { + faunus_logger->warn("rotational displacement should be between [0:2π]"); + } + } +} + void from_json(const json& j, AtomData& a) { if (!j.is_object() || j.size() != 1) { @@ -120,11 +139,6 @@ void from_json(const json& j, AtomData& a) auto val = SingleUseJSON(atom_iter.value()); a.alphax = val.value("alphax", a.alphax); a.charge = val.value("q", a.charge); - a.dp = val.value("dp", a.dp) * 1.0_angstrom; - a.dprot = val.value("dprot", a.dprot) * 1.0_rad; - if (std::fabs(a.dprot) > 2.0 * pc::pi) { - faunus_logger->warn("rotational displacement should be between [0:2π]"); - } a.id() = val.value("id", a.id()); a.mu = val.value("mu", a.mu); a.mulen = val.value("mulen", a.mulen); @@ -140,6 +154,7 @@ void from_json(const json& j, AtomData& a) a.tfe = val.value("tfe", a.tfe) * 1.0_kJmol / (1.0_angstrom * 1.0_angstrom * 1.0_molar); a.hydrophobic = val.value("hydrophobic", false); a.implicit = val.value("implicit", false); + set_dp_and_dprot(val, a); if (val.contains("activity")) { a.activity = val.at("activity").get() * 1.0_molar; } diff --git a/src/atomdata.h b/src/atomdata.h index edc2fd5f7..63159143c 100644 --- a/src/atomdata.h +++ b/src/atomdata.h @@ -1,5 +1,6 @@ #pragma once #include "core.h" +#include #include #include #include @@ -79,21 +80,21 @@ class AtomData friend void from_json(const json&, AtomData&); public: - std::string name; //!< Name - double charge = 0; //!< Particle charge [e] - double mw = 1; //!< Weight - double sigma = 0; //!< Diameter for e.g Lennard-Jones etc. [angstrom] - //!< Do not set! Only a temporal class member during the refactorization - double activity = 0; //!< Chemical activity [mol/l] - double alphax = 0; //!< Excess polarisability (unit-less) - double dp = 0; //!< Translational displacement parameter [angstrom] - double dprot = 0; //!< Rotational displacement parameter [degrees] - double tension = 0; //!< Surface tension [kT/Å^2] - double tfe = 0; //!< Transfer free energy [J/mol/angstrom^2/M] - Point mu = {0, 0, 0}; //!< Dipole moment unit vector - double mulen = 0; //!< Dipole moment length - bool hydrophobic = false; //!< Is the particle hydrophobic? - bool implicit = false; //!< Is the particle implicit (e.g. proton)? + std::string name; //!< Name + double charge = 0; //!< Particle charge [e] + double mw = 1; //!< Weight + double sigma = 0; //!< Diameter for e.g Lennard-Jones etc. [angstrom] + //!< Do not set! Only a temporal class member during the refactorization + double activity = 0; //!< Chemical activity [mol/l] + double alphax = 0; //!< Excess polarisability (unit-less) + std::optional dp = std::nullopt; //!< Translational displacement parameter [angstrom] + std::optional dprot = std::nullopt; //!< Rotational displacement parameter [degrees] + double tension = 0; //!< Surface tension [kT/Å^2] + double tfe = 0; //!< Transfer free energy [J/mol/angstrom^2/M] + Point mu = {0, 0, 0}; //!< Dipole moment unit vector + double mulen = 0; //!< Dipole moment length + bool hydrophobic = false; //!< Is the particle hydrophobic? + bool implicit = false; //!< Is the particle implicit (e.g. proton)? InteractionData interaction; //!< Arbitrary interaction parameters, e.g., epsilons in various potentials SpheroCylinderData sphero_cylinder; //!< Data for patchy sphero cylinders (PSCs) diff --git a/src/core.h b/src/core.h index da63ba660..cdb664c11 100644 --- a/src/core.h +++ b/src/core.h @@ -254,4 +254,17 @@ class Electrolyte void to_json(json& j, const Electrolyte& electrolyte); std::optional makeElectrolyte(const json& j); //!< Create ionic salt object from json +/// Extract JSON value associated with `key` into `std::optional` +/// +/// If the value does not exist, return `std::nullopt` +/// +/// @throws If value exists but cannot be extracted as `T`. +template std::optional get_optional(const json& j, std::string_view key) +{ + if (const auto it = j.find(key); it != j.end()) { + return it->get(); // may throw exception + } + return std::nullopt; +} + } // namespace Faunus diff --git a/src/move.cpp b/src/move.cpp index 1276c4c96..86a48bfea 100644 --- a/src/move.cpp +++ b/src/move.cpp @@ -189,6 +189,8 @@ void AtomicTranslateRotate::_to_json(json& j) const { j = {{"dir", directions}, {"molid", molid}, + {"dp", default_dp}, + {"dprot", default_dprot}, {unicode::rootof + unicode::bracket("r" + unicode::squared), std::sqrt(mean_square_displacement.avg())}, {"molecule", molecule_name}}; @@ -214,6 +216,8 @@ void AtomicTranslateRotate::_from_json(const json& j) } } energy_resolution = j.value("energy_resolution", 0.0); + default_dp = j.value("dp", 0.0); + default_dprot = j.value("dprot", 0.0); } void AtomicTranslateRotate::translateParticle(ParticleVector::iterator particle, @@ -262,8 +266,8 @@ void AtomicTranslateRotate::_move(Change& change) { if (auto particle = randomAtom(); particle != spc.particles.end()) { latest_particle = particle; - const auto translational_displacement = particle->traits().dp; - const auto rotational_displacement = particle->traits().dprot; + const double translational_displacement = particle->traits().dp.value_or(default_dp); + const double rotational_displacement = particle->traits().dprot.value_or(default_dprot); if (translational_displacement > 0.0) { // translate translateParticle(particle, translational_displacement); diff --git a/src/move.h b/src/move.h index 8cb8607e1..280071cf0 100644 --- a/src/move.h +++ b/src/move.h @@ -139,7 +139,7 @@ class [[maybe_unused]] AtomicSwapCharge : public Move }; /** - * @brief Translate and rotate a molecular group + * @brief Translate and rotate individual atoms */ class AtomicTranslateRotate : public Move { @@ -149,6 +149,8 @@ class AtomicTranslateRotate : public Move energy_histogram; //!< Energy histogram (value) for each particle type (key) double energy_resolution = 0.0; //!< Resolution of sampled energy histogram double latest_displacement_squared; //!< temporary squared displacement + double default_dp = 0.0; //!< Default translational displacement (Å) + double default_dprot = 0.0; //!< Default rotational displacement (rad) void sampleEnergyHistogram(); //!< Update energy histogram based on latest move void saveHistograms(); //!< Write histograms for file void checkMassCenter(