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

Author: mlund

Diff for: 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 for: 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 for: 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 for: 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<double>(j, "dp")) {
+        atomdata.dp = dp.value() * 1.0_angstrom;
+    }
+    if (const auto dprot = get_optional<double>(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<double>() * 1.0_molar;
         }

Diff for: src/atomdata.h

@@ -1,5 +1,6 @@
 #pragma once
 #include "core.h"
+#include <optional>
 #include <range/v3/range/concepts.hpp>
 #include <range/v3/view/transform.hpp>
 #include <range/v3/algorithm/any_of.hpp>
@@ -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<double> dp = std::nullopt;    //!< Translational displacement parameter [angstrom]
+    std::optional<double> 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 for: src/core.h

@@ -254,4 +254,17 @@ class Electrolyte
 void to_json(json& j, const Electrolyte& electrolyte);
 std::optional<Electrolyte> 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 <typename T> std::optional<T> get_optional(const json& j, std::string_view key)
+{
+    if (const auto it = j.find(key); it != j.end()) {
+        return it->get<T>(); // may throw exception
+    }
+    return std::nullopt;
+}
+
 } // namespace Faunus

Diff for: 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 for: 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(