Skip to content

Commit

Permalink
WidomMolecules: added absz keyword for rods on surface and added
Browse files Browse the repository at this point in the history
external energy on group upon insertion
  • Loading branch information
Mikael Lund authored and Mikael Lund committed Oct 9, 2017
1 parent 1871ab9 commit e8099c6
Show file tree
Hide file tree
Showing 2 changed files with 18 additions and 4 deletions.
20 changes: 17 additions & 3 deletions include/faunus/analysis.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -1783,6 +1787,7 @@ namespace Faunus
string molecule;
Point dir;
int molid;
bool absolute_z=false;

public:
Average<double> expu;
Expand All @@ -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
}
}
Expand All @@ -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();
}
Expand All @@ -1834,6 +1846,7 @@ namespace Faunus
{ "insertions", expu.cnt },
{ "rho", rho.avg() },
{ "rho_unit", angstrom + superminus + cubed },
{ "absz", absolute_z },
{ "mu",
{
{"excess", excess },
Expand All @@ -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 )
Expand Down
2 changes: 1 addition & 1 deletion scripts/yason.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
import sys
import json, sys, argparse
try:
from ruamel import yaml
import ruamel_yaml as yaml
except:
import yaml

Expand Down

0 comments on commit e8099c6

Please sign in to comment.