From e8099c63b75f00c849e8d8ef5d8b4da76c356f59 Mon Sep 17 00:00:00 2001 From: Mikael Lund 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 class WidomMolecule : public AnalysisBase @@ -1783,6 +1787,7 @@ namespace Faunus string molecule; Point dir; int molid; + bool absolute_z=false; public: Average 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({1,1,1}) ); // magic! + dir << j.value("dir", vector({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