Skip to content

Commit

Permalink
Merge #3107
Browse files Browse the repository at this point in the history 
3107: Fix failing grand_canonical sample test and add documentation to samples r=jngrad a=jonaslandsgesell

This PR adds documentation in to the sample files:

* widom_insertion.py
* wang_landau_reaction_ensemble.py
* grand_canonical.py
* reaction_ensemble.py

The PR also fixes partly #3104: The problem with the failing test was that the excess chemical potential did not match the submitted concentration.
I now provide a matching pair of concentration and excess chemical potential.

Co-authored-by: Jonas Landsgesell <[email protected]>
Co-authored-by: Jean-Noël Grad <[email protected]>
  • Loading branch information
@bors @jonaslandsgesell @jngrad
3 people authored Aug 30, 2019
2 parents e2e7c5d + 6e15c45 commit 1098109
Show file tree
Hide file tree
Showing 8 changed files with 284 additions and 270 deletions.
53 changes: 29 additions & 24 deletions samples/grand_canonical.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,19 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
"""
This sample performs a grand canonical simulation of a salt solution.
This example script performs a grand canonical simulation of a system in contact
with a salt reservoir and ensures constant chemical potential.
It takes two command line arguments as input: 1) the reservoir salt concentration
in units of 1/sigma^3 and 2) the excess chemical potential of the reservoir in
units of kT.
The excess chemical potential of the reservoir needs to be determined prior to
running the grand canonical simulation using the script called widom_insertion.py
which simulates a part of the reservoir at the prescribed salt concentration.
Be aware that the reservoir excess chemical potential depends on all interactions
in the reservoir system.
"""
import numpy as np
import sys
import argparse

import espressomd
from espressomd import reaction_ensemble
Expand All @@ -29,18 +38,19 @@
required_features = ["P3M", "EXTERNAL_FORCES", "LENNARD_JONES"]
espressomd.assert_features(required_features)

# print help message if proper command-line arguments are not provided
if len(sys.argv) != 3:
print("\nGot ", str(len(sys.argv) - 1), " arguments, need 2\n\nusage:" +
sys.argv[0] + " [cs_bulk] [excess_chemical_potential/kT]\n")
sys.exit()
parser = argparse.ArgumentParser(epilog=__doc__)
parser.add_argument('cs_bulk', type=float,
help="bulk salt concentration [1/sigma^3]")
parser.add_argument('excess_chemical_potential', type=float,
help="excess chemical potential [kT] "
"(obtained from Widom's insertion method)")
args = parser.parse_args()

# System parameters
#############################################################

cs_bulk = float(sys.argv[1])
excess_chemical_potential_pair = float(
sys.argv[2]) # from widom insertion simulation of pair insertion
cs_bulk = args.cs_bulk
excess_chemical_potential_pair = args.excess_chemical_potential
box_l = 50.0

# Integration parameters
Expand All @@ -54,7 +64,7 @@
system.cell_system.skin = 0.4
temperature = 1.0
system.thermostat.set_langevin(kT=temperature, gamma=.5, seed=42)
system.cell_system.max_num_cells = 2744
system.cell_system.max_num_cells = 14**3


#############################################################
Expand All @@ -78,8 +88,7 @@
for type_1 in types:
for type_2 in types:
system.non_bonded_inter[type_1, type_2].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig,
cutoff=lj_cut, shift="auto")
epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")

RE = reaction_ensemble.ReactionEnsemble(
temperature=temperature, exclusion_radius=2.0, seed=3)
Expand All @@ -96,36 +105,32 @@
p3m = electrostatics.P3M(prefactor=2.0, accuracy=1e-3)
system.actors.add(p3m)
p3m_params = p3m.get_params()
for key in list(p3m_params.keys()):
print("{} = {}".format(key, p3m_params[key]))
for key, value in p3m_params.items():
print("{} = {}".format(key, value))


# Warmup
#############################################################
# warmup integration (with capped LJ potential)
warm_steps = 1000
warm_n_times = 20
# do the warmup until the particles have at least the distance min_dist
# set LJ cap
lj_cap = 20
system.force_cap = lj_cap
system.force_cap = 20

# Warmup Integration Loop
act_min_dist = system.analysis.min_dist()
i = 0
while (i < warm_n_times):
while i < warm_n_times:
print(i, "warmup")
RE.reaction(100)
system.integrator.run(steps=warm_steps)
i += 1
#Increase LJ cap
lj_cap = lj_cap + 10
system.force_cap = lj_cap
# increase LJ cap
system.force_cap += 10

# remove force capping
system.force_cap = 0

#MC warmup
# MC warmup
RE.reaction(1000)

n_int_cycles = 10000
Expand Down
6 changes: 4 additions & 2 deletions samples/reaction_ensemble.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,10 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
"""
This sample simulates the reaction ensemble. It also illustrates how the
constant pH method can be used.
This sample illustrates how to use either the raction ensemble or the constant pH ensemble. You can choose
in which ensemble you want to simulate via either providing --reaction_ensemble or --constant_pH_ensemble as command line argument to the script.
Be aware that in the case of the reaction ensemble the dissociation constant gamma is not the thermodynamic reaction constant K, but rather K*1mol/l and therefore carries a unit!.
In the case of the of the constant pH method gamma is the thermodynamic reaction constant!
"""
import numpy as np
import argparse
Expand Down
10 changes: 9 additions & 1 deletion samples/wang_landau_reaction_ensemble.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,15 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
"""
This sample simulates the Wang-Landau Reaction Ensemble for a harmonic bond.
This example script simulates two reacting monomers which are bonded via a harmonic potential.
The script aborts as soon as the abortion criterion in the Wang-Landau algorithm is met.
The Wang-Landau simulation runs until the Wang-Landau potential is converged and then raises a Warning that it has converged, effectively aborting the simulation.
With the setup of the Wang-Landau algorithm in this script you sample the density of states of a three dimensional reacting harmonic oscillator
as a function of the two collective variables 1) degree of association and 2) potential energy.
The recorded Wang-Landau potential (which is updated during the simulation) is written to the file WL_potential_out.dat
In this simulation setup the Wang-Landau potential is the density of states. You can view the converged Wang-Landau potential e.g. via plotting with gnuplot: splot "WL_potential_out.dat".
As expected the three dimensional harmonic oscilltor has a density of states which goes like sqrt(E_pot).
For a scientific description and different ways to use the algorithm please consult https://pubs.acs.org/doi/full/10.1021/acs.jctc.6b00791
"""
import numpy as np

Expand Down
58 changes: 30 additions & 28 deletions samples/widom_insertion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,11 +17,13 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
"""
This sample measures the excess chemical potential for Widom insertion of
charged particles using the reaction ensemble method.
This example script measures the excess chemical potential of a charged WCA
fluid via Widom's insertion method.
As input this script requires you to provide particle number density in units
of 1/sigma^3.
"""
import numpy as np
import sys
import argparse

import espressomd
from espressomd import code_info
Expand All @@ -33,13 +35,16 @@
required_features = ["LENNARD_JONES", "P3M"]
espressomd.assert_features(required_features)

parser = argparse.ArgumentParser(epilog=__doc__)
parser.add_argument('cs_bulk', type=float,
help="bulk salt concentration [1/sigma^3]")
args = parser.parse_args()

# System parameters
#############################################################
assert len(sys.argv) == 2, "please provide a value for cs_bulk"
cs_bulk = float(sys.argv[1])
box_l = 50.0
N0 = int(cs_bulk * box_l**3)
print("actual cs_bulk", float(N0) / box_l**3)
cs_bulk = args.cs_bulk
N0 = 70
box_l = (N0 / cs_bulk)**(1.0 / 3.0)

# Integration parameters
#############################################################
Expand All @@ -51,7 +56,7 @@
system.cell_system.skin = 0.4
temperature = 1.0
system.thermostat.set_langevin(kT=temperature, gamma=1.0, seed=42)
system.cell_system.max_num_cells = 2744
system.cell_system.max_num_cells = 14**3


#############################################################
Expand All @@ -76,41 +81,39 @@
for type_1 in types:
for type_2 in types:
system.non_bonded_inter[type_1, type_2].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig,
cutoff=lj_cut, shift="auto")
epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")

p3m = electrostatics.P3M(prefactor=0.9, accuracy=1e-3)
p3m = electrostatics.P3M(prefactor=2.0, accuracy=1e-3)
system.actors.add(p3m)
p3m_params = p3m.get_params()
for key in list(p3m_params.keys()):
print("{} = {}".format(key, p3m_params[key]))
for key, value in p3m_params.items():
print("{} = {}".format(key, value))

# Warmup
#############################################################
# warmup integration (with capped LJ potential)
warm_steps = 1000
warm_n_times = 20
# do the warmup until the particles have at least the distance min_dist
# set LJ cap
lj_cap = 20
system.force_cap = lj_cap
system.force_cap = 20

# Warmup Integration Loop
act_min_dist = system.analysis.min_dist()
i = 0
while (i < warm_n_times):
while i < warm_n_times:
print(i, "warmup")
system.integrator.run(steps=warm_steps)
i += 1
# Increase LJ cap
lj_cap = lj_cap + 10
system.force_cap = lj_cap
# increase LJ cap
system.force_cap += 10

# remove force capping
system.force_cap = 0

RE = reaction_ensemble.WidomInsertion(
temperature=temperature, seed=77)

# add insertion reaction
insertion_reaction_id = 0
RE.add_reaction(reactant_types=[],
reactant_coefficients=[], product_types=[1, 2],
product_coefficients=[1, 1], default_charges={1: -1, 2: +1})
Expand All @@ -119,16 +122,15 @@

n_iterations = 100
for i in range(n_iterations):
for j in range(30):
RE.measure_excess_chemical_potential(0) # 0 for insertion reaction
system.integrator.run(steps=2)
for j in range(50):
RE.measure_excess_chemical_potential(insertion_reaction_id)
system.integrator.run(steps=500)
if i % 20 == 0:
print("mu_ex_pair ({:.4f}, +/- {:.4f})".format(
*RE.measure_excess_chemical_potential(0) # 0 for insertion reaction
))
*RE.measure_excess_chemical_potential(insertion_reaction_id)))
print("HA", system.number_of_particles(type=0), "A-",
system.number_of_particles(type=1), "H+",
system.number_of_particles(type=2))

print(RE.measure_excess_chemical_potential(0)) # 0 for insertion reaction
print("excess chemical potential for an ion pair ",
RE.measure_excess_chemical_potential(insertion_reaction_id))
Loading

0 comments on commit 1098109

Please sign in to comment.