Add harmonic potential to constraint (#445)

* Add harmonic potential to constraint

* Add include

* Adjust json output

* Formatting

* Change optional dereference style

docs/_docs/energy.md

@@ -1067,6 +1067,16 @@ energy:
     - constrain: {type: molecule, index: 0, property: end2end, range: [0,200]}
 ~~~
 
+Instead of a range, the constraint can be in the form of a harmonic potential
+where `req` is the value of the reaction coordinate at the minimum and `k` is
+the force constant. If `harmonic` is given, the `range` is ignored.
+The following constrains one box side length to 100 angstroms:
+
+~~~ yaml
+energy:
+    - constrain: {type: system, property: Lx, harmonic: {req: 100.0, k: 500.0}}
+~~~
+
 Tip: placing `constrain` at the _top_ of the energy list is more efficient as the remaining
 energy terms are skipped should an infinite energy arise.
 

src/energy.cpp

@@ -805,24 +805,37 @@ TEST_CASE("Energy::Isobaric") {
     }
 }
 
-Constrain::Constrain(const json &j, Space &spc) {
+Constrain::Constrain(const json& j, Space& spc) {
     name = "constrain";
     type = j.at("type").get<std::string>();
+    if (const auto it = j.find("harmonic"); it != j.end()) {
+        harmonic = std::make_optional<Faunus::pairpotential::HarmonicBond>();
+        harmonic->from_json(*it);
+    }
     coordinate = ReactionCoordinate::createReactionCoordinate({{type, j}}, spc);
 }
 
 double Constrain::energy(const Change& change) {
     if (change) {
         const auto value = (*coordinate)(); // calculate reaction coordinate
-        if (not coordinate->inRange(value)) // is it within allowed range?
-            return pc::infty;               // if not, return infinite energy
+        if (harmonic) {
+            return harmonic->half_force_constant * std::pow(harmonic->equilibrium_distance - value, 2);
+        }
+        if (not coordinate->inRange(value)) { // if outside allowed range ...
+            return pc::infty;                 // ... return infinite energy
+        }
     }
     return 0.0;
 }
-void Constrain::to_json(json &j) const {
+void Constrain::to_json(json& j) const {
     j = json(*coordinate).at(type);
     j.erase("resolution");
     j["type"] = type;
+    if (harmonic) {
+        harmonic->to_json(j["harmonic"]);
+        j["harmonic"].erase("index");
+        j.erase("range");
+    }
 }
 
 void Bonded::updateGroupBonds(const Space::GroupType& group) {

src/energy.h

@@ -17,6 +17,7 @@
 #include <numeric>
 #include <algorithm>
 #include <concepts>
+#include <optional>
 
 struct freesasa_parameters_fwd; // workaround for freesasa unnamed struct that cannot be forward declared
 
@@ -278,6 +279,7 @@ class Constrain : public EnergyTerm {
   private:
     std::string type;
     std::unique_ptr<ReactionCoordinate::ReactionCoordinateBase> coordinate;
+    std::optional<Faunus::pairpotential::HarmonicBond> harmonic;
 
   public:
     Constrain(const json& j, Space& space);