ewald initial energy assertion fix; added bulk example (albeit still no

saved output)

CMakeLists.txt

@@ -311,14 +311,19 @@ add_test(
     ; ${PYTHON_EXECUTABLE} ../scripts/jsoncompare.py minimal.out.json out.json --tol 0.02"
     WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/examples)
 
+add_test(
+    NAME bulk
+    COMMAND sh -c "${PYTHON_EXECUTABLE} ../scripts/yason.py bulk.yml\
+    | $<TARGET_FILE:faunus> --quiet"
+    WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/examples)
+
 add_test(
     NAME membrane
     COMMAND sh -c "${PYTHON_EXECUTABLE} ../scripts/yason.py membrane.yml\
     | $<TARGET_FILE:faunus> --quiet --state membrane.state.json\
     ; ${PYTHON_EXECUTABLE} ../scripts/jsoncompare.py membrane.out.json out.json --tol 0.065"
     WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/examples)
 
-
 # INSTALL TARGETS
 
 install(TARGETS faunus DESTINATION bin)

docs/_docs/energy.md

@@ -304,7 +304,7 @@ and hydrophobic/hydrophilic interactions.
 `sasa`       | Description
 ------------ | ----------------------------------------------------------
 `radius=1.4` | Probe radius for SASA calculation (angstrom)
-`molarity=0` | Molar concentration of co-solute, $c_s$.
+`molarity`   | Molar concentration of co-solute, $c_s$.
 
 
 ## Bonded Interactions

examples/bulk.yml

@@ -1,24 +1,25 @@
 temperature: 1100
-mcloop: {macro: 10, micro: 100}
-geometry: {length: 19.73}
+mcloop: {macro: 2, micro: 25}
+geometry: {length: 42.5} # 19.73
 energy:
+    - isobaric: {P/mM: 11}
     - nonbonded_coulomblj:
         lennardjones: {mixing: LB}
-        coulomb: {type: ewald, epsr: 1, cutoff: 9, kcutoff: 7, alpha: 0.2}
+        coulomb: {type: ewald, epsr: 1, cutoff: 14, kcutoff: 7, alpha: 0.2}
 atomlist:
-    - Na: {q:  1.0, sigma: 3.33, eps: 0.01158968, dp: 1} 
-    - Cl: {q: -1.0, sigma: 4.4, eps: 0.4184, dp: 1} 
+    - Na: {q:  1.0, sigma: 3.33, eps: 0.01158968, dp: 0.5} 
+    - Cl: {q: -1.0, sigma: 4.4, eps: 0.4184, dp: 0.5} 
 moleculelist:
     - salt: {atoms: [Na,Cl], atomic: true}
 insertmolecules:
-    - salt: {N: 152}
+    - salt: {N: 1152} # 152
 moves:
     - transrot: {molecule: salt}
+    - volume: {dV: 0.03}
 analysis:
-    - atomrdf: {file: rdf.dat, nstep: 10, dr: 0.1, name1: Na, name2: Cl}
-    - virtualvolume: {dV: 0.2, nstep: 5}
-    - systemenergy: {file: energy.dat, nstep: 10}
-    - xtcfile: {file: traj.xtc, nstep: 10}
+    - atomrdf: {file: rdf.dat, nstep: 20, dr: 0.1, name1: Na, name2: Cl}
+    - systemenergy: {file: energy.dat, nstep: 100}
+    - xtcfile: {file: traj.xtc, nstep: 50}
     - savestate: {file: confout.pqr}
-    - savestate: {file: state}
-    - density: {nstep: 10}
+    - savestate: {file: state.json}
+    - density: {nstep: 50}

src/analysis.cpp

@@ -129,18 +129,22 @@ void Faunus::Analysis::PairFunctionBase::normalize() {
 
 void Faunus::Analysis::PairFunctionBase::_to_json(Faunus::json &j) const {
     j = {
-        {"dr", dr}, {"name1", name1}, {"name2", name2}, {"file", file},
-        {"file2", file2}, {"dim", dim}, {"Rhyper", Rhypersphere}
+        {"dr", dr/1.0_angstrom},
+        {"name1", name1},
+        {"name2", name2},
+        {"file", file},
+        {"dim", dim}
     };
+    if (Rhypersphere>0)
+        j["Rhyper"] = Rhypersphere;
 }
 
 void Faunus::Analysis::PairFunctionBase::_from_json(const Faunus::json &j) {
     file = j.at("file");
-    file2 = j.value("file2", file+".avg");
     name1 = j.at("name1");
     name2 = j.at("name2");
     dim = j.value("dim", 3);
-    dr = j.value("dr", 0.1);
+    dr = j.value("dr", 0.1) * 1.0_angstrom;
     hist.setResolution(dr);
     hist2.setResolution(dr);
     Rhypersphere = j.value("Rhyper", -1.0);
@@ -160,7 +164,7 @@ void Faunus::Analysis::VirtualVolume::_sample() {
 void Faunus::Analysis::VirtualVolume::_from_json(const Faunus::json &j) { dV = j.at("dV"); }
 
 void Faunus::Analysis::VirtualVolume::_to_json(Faunus::json &j) const {
-    double pex = log(duexp.avg()) / dV;
+    double pex = log(duexp.avg()) / dV; // excess pressure
     j = {
         {"dV", dV}, {"Pex/mM", pex/1.0_mM},
         {"Pex/Pa", pex/1.0_Pa}, {"Pex/kT/"+u8::angstrom+u8::cubed, pex}

src/analysis.h

@@ -706,19 +706,18 @@ namespace Faunus {
          */
         class PairFunctionBase : public Analysisbase {
             protected:
-                int dim;
+                int dim=3;
                 int id1=-1, id2=-1; // particle id (mol or atom)
-                double dr;
+                double dr=0;
                 Table2D<double,double> hist;
                 Table2D<double,Average<double>> hist2;
-                std::string name1, name2, file, file2;
-                double Rhypersphere; // Radius of 2D hypersphere
+                std::string name1, name2, file;
+                double Rhypersphere=-1; // Radius of 2D hypersphere
                 Average<double> V;   // average volume (angstrom^3)
                 virtual void normalize();
-            private:
 
+            private:
                 void _from_json(const json &j) override;
-
                 void _to_json(json &j) const override;
 
             public:

src/energy.h

@@ -1150,6 +1150,7 @@ namespace Faunus {
                     std::shared_ptr<POWERSASA::PowerSasa<float,Point>> ps=nullptr;
 
                     void updateSASA(const Tpvec &p) {
+                        assert(ps != nullptr);
                         radii.resize(p.size());
                         std::transform(p.begin(), p.end(), radii.begin(),
                                 [this](auto &a){ return atoms<Tparticle>[a.id].sigma*0.5 + this->probe;});
@@ -1203,7 +1204,7 @@ namespace Faunus {
                         name = "sasa";
                         cite = "doi:10.1002/jcc.21844";
                         probe = j.value("radius", 1.4) * 1.0_angstrom;
-                        conc = j.value("molarity", conc) * 1.0_molar;
+                        conc = j.at("molarity").get<double>() * 1.0_molar;
                         init();
                     }
 
@@ -1220,7 +1221,7 @@ namespace Faunus {
                     double energy(Change &change) override {
                         double u=0, A=0;
                         /*
-                         * ideally we want to call `update` only is `key==NEW` but
+                         * ideally we want to call `update` only if `key==NEW` but
                          * syncronising the PowerSasa object is difficult since it's
                          * non-copyable.
                          */
@@ -1326,7 +1327,7 @@ namespace Faunus {
                                     }
 
                                     if (vec.size()==oldsize)
-                                        std::cerr << "warning: ignoring unknown energy '" << it.key() << "'" << endl;
+                                        throw std::runtime_error("unknown term");
 
                                 } catch (std::exception &e) {
                                     throw std::runtime_error("Error adding energy '" + it.key() + "': " + e.what());

src/move.h

@@ -1193,23 +1193,34 @@ namespace Faunus {
                 Average<double> uavg;
 
                 void init() {
+                    dusum=0;
+                    Change c; c.all=true;
+
                     state1.pot.key = Energy::Energybase::OLD; // this is the old energy (current)
                     state2.pot.key = Energy::Energybase::NEW; // this is the new energy (trial)
                     state1.pot.init();
-                    state2.pot.init();
-                    dusum=0;
-                    Change c; c.all=true;
-                    state2.sync(state1, c);
+
                     uinit = state1.pot.energy(c);
 
+                    state2.sync(state1, c);
+                    state2.pot.init();
+
                     // Hack in reference to state1 in speciation
                     for (auto base : moves.vec) {
                         auto derived = std::dynamic_pointer_cast<Move::SpeciationMove<Tspace>>(base);
                         if (derived)
                             derived->setOther(state1.spc);
                     }
-                    //cout << "s1 = " << uinit << "  s2 = " << state2.pot.energy(c) << endl; 
-                    assert(state1.pot.energy(c) == state2.pot.energy(c));
+#ifndef NDEBUG
+                    double u2 = state2.pot.energy(c);
+                    double error = std::fabs(uinit-u2);
+                    if (uinit!=0)
+                        assert(error/uinit<1e-3);
+                    else
+                        assert(error<1e-6);
+                    //cout << "u1 = " << uinit << "  u2 = " << u2 << endl;
+                    //assert( std::fabs((uinit-u2)/uinit)<1e-3 );
+#endif
                 }
 
             public:

src/tabulate.h

@@ -189,7 +189,7 @@ namespace Faunus
                  * @param d Table data
                  * @param r2 x value
                  */
-                T eval( const typename base::data &d, T r2 ) const
+                inline T eval( const typename base::data &d, T r2 ) const
                 {
                     auto low = std::lower_bound(d.r2.begin(), d.r2.end(), r2);
                     size_t pos = (low - d.r2.begin() - 1);