From e8099c63b75f00c849e8d8ef5d8b4da76c356f59 Mon Sep 17 00:00:00 2001
From: Mikael Lund <mikael@Mikaels-MacBook-Pro.local>
Date: Mon, 9 Oct 2017 11:42:14 +0200
Subject: [PATCH] WidomMolecules: added `absz` keyword for rods on surface and
 added external energy on group upon insertion

---
 include/faunus/analysis.h | 20 +++++++++++++++++---
 scripts/yason.py          |  2 +-
 2 files changed, 18 insertions(+), 4 deletions(-)

diff --git a/include/faunus/analysis.h b/include/faunus/analysis.h
index 965a3bd21..9f53658a5 100755
--- a/include/faunus/analysis.h
+++ b/include/faunus/analysis.h
@@ -1763,14 +1763,18 @@ namespace Faunus
          * the system and calculate a Widom average to measure
          * the free energy of the insertion process. The position
          * and molecular rotation is random.
+         * For use with rod-like particles on surfaces, the `absz`
+         * keyword may be used to ensure orientations on only one
+         * half-sphere.
          *
          * JSON input:
          *
          * Keyword       | Description
          * :------------ | :-----------------------------------------
-         * `dir`         | Inserting direction array. Default [1,1,1]
+         * `dir`         | Inserting direction array. Default `[1,1,1]`
          * `molecule`    | Name of molecule to insert
          * `ninsert`     | Number of insertions per sample event
+         * `absz`        | Apply `std::fabs` on all z-coordinates of inserted molecule (default: `false`)
          */
         template<typename Tspace>
             class WidomMolecule : public AnalysisBase
@@ -1783,6 +1787,7 @@ namespace Faunus
                 string molecule;
                 Point dir;
                 int molid;
+                bool absolute_z=false;
 
             public:
                 Average<double> expu;
@@ -1798,7 +1803,15 @@ namespace Faunus
                     for ( int i = 0; i < ninsert; ++i )
                     {
                         auto pin = rins(spc->geo, spc->p, spc->molecule[molid]); // ('spc->molecule' is a vector of molecules
+                        assert( !pin.empty() );
+
+                        if (absolute_z)
+                            for (auto &p : pin)
+                                p.z() = std::fabs(p.z());
+
                         double u = pot->v2v(pin, spc->p); // energy between "ghost molecule" and system in kT
+                        Group g = Group(0, pin.size()-1 );
+                        u += pot->g_external(pin, g);
                         expu += exp(-u); // widom average
                     }
                 }
@@ -1816,7 +1829,6 @@ namespace Faunus
                             << pad(SUB, w, "Particle density") << WidomMolecule::rho.avg() << angstrom + superminus + cubed
                             << "\n"
                             << pad(SUB, w, "Excess chemical potential") << excess << kT << "\n";
-                        //todo : print rho in mol/L
                     }
                     return o.str();
                 }
@@ -1834,6 +1846,7 @@ namespace Faunus
                                 { "insertions", expu.cnt },
                                 { "rho", rho.avg() },
                                 { "rho_unit", angstrom + superminus + cubed },
+                                { "absz", absolute_z },
                                 { "mu",
                                     {
                                         {"excess", excess },
@@ -1851,8 +1864,9 @@ namespace Faunus
             {
                 name = "Widom Molecule";
                 ninsert = j.at("ninsert");
-                dir << j.value("dir", vector<double>({1,1,1}) );  // magic!
+                dir << j.value("dir", vector<double>({1,1,1}) );
                 molecule = j.at("molecule");
+                absolute_z = j.value("absz", false);
                 // look up the id of the molecule that we want to insert
                 molid = -1;
                 for ( unsigned long i = 0; i < spc.molecule.size(); ++i )
diff --git a/scripts/yason.py b/scripts/yason.py
index 972075d7b..d9f196843 100755
--- a/scripts/yason.py
+++ b/scripts/yason.py
@@ -3,7 +3,7 @@
 import sys
 import json, sys, argparse
 try:
-    from ruamel import yaml
+    import ruamel_yaml as yaml
 except:
     import yaml