From c1e1f0fba57f32b9447c6d368aedb6ce544b6d31 Mon Sep 17 00:00:00 2001 From: Jonas Landsgesell Date: Wed, 28 Aug 2019 15:54:22 +0200 Subject: [PATCH 01/17] add documentation --- samples/grand_canonical.py | 7 +++++-- samples/reaction_ensemble.py | 6 ++++-- samples/wang_landau_reaction_ensemble.py | 4 +++- samples/widom_insertion.py | 6 +++--- 4 files changed, 15 insertions(+), 8 deletions(-) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index 98f089f438f..2b47b8571df 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -17,7 +17,10 @@ # along with this program. If not, see . # """ -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. +From the Widom insertion script you obtain the excess chemical potential you need to provide for this script. """ import numpy as np import sys @@ -32,7 +35,7 @@ # 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.argv[0] + " [cs_bulk [1/sigma^3]] [excess_chemical_potential [kT]]\n") sys.exit() # System parameters diff --git a/samples/reaction_ensemble.py b/samples/reaction_ensemble.py index 9d3a615d19b..6f8852e4ef4 100644 --- a/samples/reaction_ensemble.py +++ b/samples/reaction_ensemble.py @@ -17,8 +17,10 @@ # along with this program. If not, see . # """ -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 diff --git a/samples/wang_landau_reaction_ensemble.py b/samples/wang_landau_reaction_ensemble.py index ab2792a61c3..91ddb7f0e78 100644 --- a/samples/wang_landau_reaction_ensemble.py +++ b/samples/wang_landau_reaction_ensemble.py @@ -17,7 +17,9 @@ # along with this program. If not, see . # """ -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 potential which is updated during the simulation is written to the file WL_potential_out.dat. """ import numpy as np diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index df832981c80..3bd6b5c614f 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -17,8 +17,8 @@ # along with this program. If not, see . # """ -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 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 @@ -131,4 +131,4 @@ 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(0)) # 0 for insertion reaction From abd3e9e8a70baa5c66057b40879de3a2a8ec3ae0 Mon Sep 17 00:00:00 2001 From: Jonas Landsgesell Date: Wed, 28 Aug 2019 16:02:29 +0200 Subject: [PATCH 02/17] fix test --- testsuite/scripts/samples/test_grand_canonical.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/scripts/samples/test_grand_canonical.py b/testsuite/scripts/samples/test_grand_canonical.py index af46713722e..48673ac0f32 100644 --- a/testsuite/scripts/samples/test_grand_canonical.py +++ b/testsuite/scripts/samples/test_grand_canonical.py @@ -20,7 +20,7 @@ sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@SAMPLES_DIR@/grand_canonical.py", n_int_cycles=51, n_int_steps=5, - warm_n_times=10, warm_steps=50, cmd_arguments=[0.01, 0.01]) + warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659461370000200e-04, -3.365753246398797693e-01]) @skipIfMissingFeatures From 9a2eb6e81cc03166e3a4f1b6f87371ca30174e4f Mon Sep 17 00:00:00 2001 From: Jonas Landsgesell Date: Wed, 28 Aug 2019 16:12:57 +0200 Subject: [PATCH 03/17] up --- samples/widom_insertion.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index 3bd6b5c614f..e579c44645c 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -131,4 +131,7 @@ system.number_of_particles(type=1), "H+", system.number_of_particles(type=2)) -print("excess chemical potential for an ion pair ",RE.measure_excess_chemical_potential(0)) # 0 for insertion reaction +print( + "excess chemical potential for an ion pair ", + RE.measure_excess_chemical_potential( + 0)) # 0 for insertion reaction From e4e416c9a658eb41eaa9998aa181998a3699d4e6 Mon Sep 17 00:00:00 2001 From: Jonas Landsgesell Date: Thu, 29 Aug 2019 12:40:41 +0200 Subject: [PATCH 04/17] add doc --- samples/wang_landau_reaction_ensemble.py | 8 +++++++- samples/widom_insertion.py | 12 +++++------- testsuite/scripts/samples/test_grand_canonical.py | 2 +- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/samples/wang_landau_reaction_ensemble.py b/samples/wang_landau_reaction_ensemble.py index 91ddb7f0e78..44b049d496a 100644 --- a/samples/wang_landau_reaction_ensemble.py +++ b/samples/wang_landau_reaction_ensemble.py @@ -19,7 +19,13 @@ """ 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 potential which is updated during the simulation is written to the file WL_potential_out.dat. +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 which are introduced are the degree of association and the potential energy. +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. +The recorded Wang-Landau potential (which is updated during the simulation is written to the file WL_potential_out.dat) +contains the density of states. You can view the converged Wang-Landau potential e.g. via plotting with gnuplot: sp "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 of the algorithm please consult https://pubs.acs.org/doi/full/10.1021/acs.jctc.6b00791 """ import numpy as np diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index e579c44645c..6249593cc87 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -17,7 +17,7 @@ # along with this program. If not, see . # """ -This example script measures the excess chemical potential via Widom's insertion 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 @@ -37,9 +37,8 @@ ############################################################# 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) +N0 = 70 +box_l = (float(N0)/cs_bulk)**(1.0/3.0) # Integration parameters ############################################################# @@ -79,7 +78,7 @@ 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()): @@ -119,9 +118,8 @@ n_iterations = 100 for i in range(n_iterations): - for j in range(30): + for j in range(50): RE.measure_excess_chemical_potential(0) # 0 for insertion reaction - system.integrator.run(steps=2) system.integrator.run(steps=500) if i % 20 == 0: print("mu_ex_pair ({:.4f}, +/- {:.4f})".format( diff --git a/testsuite/scripts/samples/test_grand_canonical.py b/testsuite/scripts/samples/test_grand_canonical.py index 48673ac0f32..2fa1565fb30 100644 --- a/testsuite/scripts/samples/test_grand_canonical.py +++ b/testsuite/scripts/samples/test_grand_canonical.py @@ -20,7 +20,7 @@ sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@SAMPLES_DIR@/grand_canonical.py", n_int_cycles=51, n_int_steps=5, - warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659461370000200e-04, -3.365753246398797693e-01]) + warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659461370000200e-04, -0.336]) @skipIfMissingFeatures From 9874d6e601a6126b061536c5c8bef4826780ae1b Mon Sep 17 00:00:00 2001 From: Jonas Landsgesell Date: Thu, 29 Aug 2019 12:44:50 +0200 Subject: [PATCH 05/17] wording --- samples/wang_landau_reaction_ensemble.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/samples/wang_landau_reaction_ensemble.py b/samples/wang_landau_reaction_ensemble.py index 44b049d496a..a0f329bfadc 100644 --- a/samples/wang_landau_reaction_ensemble.py +++ b/samples/wang_landau_reaction_ensemble.py @@ -19,13 +19,13 @@ """ 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. -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 which are introduced are the degree of association and the potential energy. 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. -The recorded Wang-Landau potential (which is updated during the simulation is written to the file WL_potential_out.dat) -contains the density of states. You can view the converged Wang-Landau potential e.g. via plotting with gnuplot: sp "WL_potential_out.dat". +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 of the algorithm please consult https://pubs.acs.org/doi/full/10.1021/acs.jctc.6b00791 +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 From f5650ffaddc591d9e08caa48ce5b00054690d994 Mon Sep 17 00:00:00 2001 From: Jonas Landsgesell Date: Thu, 29 Aug 2019 12:46:10 +0200 Subject: [PATCH 06/17] fmt --- samples/widom_insertion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index 6249593cc87..f772a4d28ad 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -38,7 +38,7 @@ assert len(sys.argv) == 2, "please provide a value for cs_bulk" cs_bulk = float(sys.argv[1]) N0 = 70 -box_l = (float(N0)/cs_bulk)**(1.0/3.0) +box_l = (N0 / cs_bulk)**(1.0 / 3.0) # Integration parameters ############################################################# From e73a2cbd12fc924aa6f1828f6c79285e2d8f1fe2 Mon Sep 17 00:00:00 2001 From: Jonas Landsgesell Date: Thu, 29 Aug 2019 18:14:59 +0200 Subject: [PATCH 07/17] truncate and add comment to how one obtains the magic parameters --- testsuite/scripts/samples/test_grand_canonical.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/testsuite/scripts/samples/test_grand_canonical.py b/testsuite/scripts/samples/test_grand_canonical.py index 2fa1565fb30..fbc07ade012 100644 --- a/testsuite/scripts/samples/test_grand_canonical.py +++ b/testsuite/scripts/samples/test_grand_canonical.py @@ -20,7 +20,9 @@ sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@SAMPLES_DIR@/grand_canonical.py", n_int_cycles=51, n_int_steps=5, - warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659461370000200e-04, -0.336]) + warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659e-04, -0.336]) #note the command line arguments which are submitted are 1) the salt concentration in 1/sigma^3 and the excess chemical potential in kT. + #the excess chemical potential needs to be determined using the Widom insertion method. E.g. with the example script widom_insertion.py + @skipIfMissingFeatures From f4d673d0b528d6f89ce55d0cb43deeadc4b2d9a9 Mon Sep 17 00:00:00 2001 From: Jonas Landsgesell Date: Thu, 29 Aug 2019 18:28:01 +0200 Subject: [PATCH 08/17] fmt --- testsuite/scripts/samples/test_grand_canonical.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/testsuite/scripts/samples/test_grand_canonical.py b/testsuite/scripts/samples/test_grand_canonical.py index fbc07ade012..0d21a07a690 100644 --- a/testsuite/scripts/samples/test_grand_canonical.py +++ b/testsuite/scripts/samples/test_grand_canonical.py @@ -20,11 +20,11 @@ sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@SAMPLES_DIR@/grand_canonical.py", n_int_cycles=51, n_int_steps=5, - warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659e-04, -0.336]) #note the command line arguments which are submitted are 1) the salt concentration in 1/sigma^3 and the excess chemical potential in kT. - #the excess chemical potential needs to be determined using the Widom insertion method. E.g. with the example script widom_insertion.py + warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659e-04, -0.336]) # note the command line arguments which are submitted are 1) the salt concentration in 1/sigma^3 and the excess chemical potential in kT. + # the excess chemical potential needs to be determined using the Widom + # insertion method. E.g. with the example script widom_insertion.py - @skipIfMissingFeatures class Sample(ut.TestCase): system = sample.system From 876ed64492bc1cb60a53d140273e8ce9417f4ab3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 30 Aug 2019 16:17:02 +0200 Subject: [PATCH 09/17] Cleanup minimization code --- samples/grand_canonical.py | 14 +++++--------- samples/widom_insertion.py | 12 ++++-------- 2 files changed, 9 insertions(+), 17 deletions(-) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index 2b47b8571df..f2dc02b60cd 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -108,27 +108,23 @@ # 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 diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index f772a4d28ad..44008f33375 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -89,21 +89,17 @@ # 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 From 5a8ebf3e47dae0d9d4c8bb41a9d5a93b4e7556f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 30 Aug 2019 16:17:58 +0200 Subject: [PATCH 10/17] Increase test tolerance Running this test 5000 times resulted in 65 occurences of a deviation larger than 10%, with a maximal deviation of 13.6%, a mean of 3.3% and a median of 2.9%. Setting the test threshold to 15% should significantly reduce the frequency of random failures. --- testsuite/scripts/samples/test_grand_canonical.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/testsuite/scripts/samples/test_grand_canonical.py b/testsuite/scripts/samples/test_grand_canonical.py index 0d21a07a690..9748f03b20c 100644 --- a/testsuite/scripts/samples/test_grand_canonical.py +++ b/testsuite/scripts/samples/test_grand_canonical.py @@ -30,8 +30,8 @@ class Sample(ut.TestCase): system = sample.system def test_deviation_to_target_concentration(self): - # deviation < 10% - self.assertLess(abs(sample.deviation), 10.0) + # deviation < 15% + self.assertLess(abs(sample.deviation), 15.0) if __name__ == "__main__": From 21a7d5da5f0c6fd44e171aabf652df957244606a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 30 Aug 2019 17:00:41 +0200 Subject: [PATCH 11/17] Document origin of hardcoded values --- samples/grand_canonical.py | 12 +++++------ samples/widom_insertion.py | 20 +++++++++---------- .../scripts/samples/test_grand_canonical.py | 10 ++++++---- 3 files changed, 22 insertions(+), 20 deletions(-) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index f2dc02b60cd..9e474aad8a2 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -35,15 +35,15 @@ # 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 [1/sigma^3]] [excess_chemical_potential [kT]]\n") + sys.argv[0] + " [cs_bulk [1/sigma^3]] [excess_chemical_potential [kT]]\n" + "The excess chemical potential needs to be determined using Widom's insertion method\n") sys.exit() # System parameters ############################################################# cs_bulk = float(sys.argv[1]) -excess_chemical_potential_pair = float( - sys.argv[2]) # from widom insertion simulation of pair insertion +excess_chemical_potential_pair = float(sys.argv[2]) box_l = 50.0 # Integration parameters @@ -57,7 +57,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 ############################################################# @@ -99,8 +99,8 @@ 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 diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index 44008f33375..5a688201e87 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -50,7 +50,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 ############################################################# @@ -81,8 +81,8 @@ 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 ############################################################# @@ -106,6 +106,9 @@ 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}) @@ -115,17 +118,14 @@ n_iterations = 100 for i in range(n_iterations): for j in range(50): - RE.measure_excess_chemical_potential(0) # 0 for insertion reaction + 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( - "excess chemical potential for an ion pair ", - 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)) diff --git a/testsuite/scripts/samples/test_grand_canonical.py b/testsuite/scripts/samples/test_grand_canonical.py index 9748f03b20c..e42a03a0ff4 100644 --- a/testsuite/scripts/samples/test_grand_canonical.py +++ b/testsuite/scripts/samples/test_grand_canonical.py @@ -20,10 +20,12 @@ sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@SAMPLES_DIR@/grand_canonical.py", n_int_cycles=51, n_int_steps=5, - warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659e-04, -0.336]) # note the command line arguments which are submitted are 1) the salt concentration in 1/sigma^3 and the excess chemical potential in kT. - # the excess chemical potential needs to be determined using the Widom - # insertion method. E.g. with the example script widom_insertion.py - + warm_n_times=10, warm_steps=50, cmd_arguments=[8.523659e-04, -0.336]) +# Note: the command line arguments which are submitted are 1) the salt +# concentration in 1/sigma^3 and 2) the excess chemical potential in kT. +# The excess chemical potential needs to be determined using the Widom +# insertion method, e.g. with the example script widom_insertion.py + @skipIfMissingFeatures class Sample(ut.TestCase): From e68fb6792eb723c4a8245145624a4d76d08056d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 30 Aug 2019 17:15:47 +0200 Subject: [PATCH 12/17] Use argparse module --- samples/grand_canonical.py | 19 ++++++++++--------- samples/widom_insertion.py | 10 +++++++--- 2 files changed, 17 insertions(+), 12 deletions(-) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index 9e474aad8a2..46e428b5fb3 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -23,7 +23,7 @@ From the Widom insertion script you obtain the excess chemical potential you need to provide for this script. """ import numpy as np -import sys +import argparse import espressomd from espressomd import reaction_ensemble @@ -32,18 +32,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 [1/sigma^3]] [excess_chemical_potential [kT]]\n" - "The excess chemical potential needs to be determined using Widom's insertion method\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]) +cs_bulk = args.cs_bulk +excess_chemical_potential_pair = args.excess_chemical_potential box_l = 50.0 # Integration parameters diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index 5a688201e87..96b39c22076 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -21,7 +21,7 @@ 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 @@ -33,10 +33,14 @@ 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]) +cs_bulk = args.cs_bulk N0 = 70 box_l = (N0 / cs_bulk)**(1.0 / 3.0) From 428b51abfd8c1ee07329724c8260ce7eca38c923 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 30 Aug 2019 17:25:20 +0200 Subject: [PATCH 13/17] Remove repetitions in the documentation --- samples/grand_canonical.py | 1 - src/core/reaction_ensemble.cpp | 18 ++++++------------ 2 files changed, 6 insertions(+), 13 deletions(-) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index 46e428b5fb3..ca0f8c56b16 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -20,7 +20,6 @@ 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. -From the Widom insertion script you obtain the excess chemical potential you need to provide for this script. """ import numpy as np import argparse diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 84d7c0210c9..b7a94a70f31 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -164,7 +164,6 @@ void ReactionAlgorithm::add_reaction( * set. */ void ReactionAlgorithm::check_reaction_ensemble() { - /**checks the reaction_ensemble struct for valid parameters */ if (reactions.empty()) { throw std::runtime_error("Reaction system not initialized"); } @@ -369,7 +368,6 @@ double ReactionEnsemble::calculate_acceptance_probability( std::map &old_particle_numbers, int dummy_old_state_index, int dummy_new_state_index, bool dummy_only_make_configuration_changing_move) { - /**calculate the acceptance probability in the reaction ensemble */ const double factorial_expr = calculate_factorial_expression(current_reaction, old_particle_numbers); @@ -607,7 +605,6 @@ int ReactionAlgorithm::delete_particle(int p_id) { * range becoming excessively huge. */ - /**deletes the particle with the provided id */ int old_max_seen_id = max_seen_particle; if (p_id == old_max_seen_id) { // last particle, just delete @@ -1197,16 +1194,13 @@ int WangLandauReactionEnsemble::initialize_wang_landau() { } /** - * Calculates the expression which occurs in the Wang-Landau acceptance - * probability. + * Calculates the expression in the acceptance probability of the Wang-Landau + * reaction ensemble */ double WangLandauReactionEnsemble::calculate_acceptance_probability( SingleReaction ¤t_reaction, double E_pot_old, double E_pot_new, std::map &old_particle_numbers, int old_state_index, int new_state_index, bool only_make_configuration_changing_move) { - /**determine the acceptance probabilities of the reaction move - * in Wang-Landau reaction ensemble - */ double beta = 1.0 / temperature; double bf; if (do_not_sample_reaction_partition_function || @@ -1756,15 +1750,15 @@ int ConstantpHEnsemble::do_reaction(int reaction_steps) { return 0; } +/** + * Calculates the expression in the acceptance probability of the constant pH + * method. + */ double ConstantpHEnsemble::calculate_acceptance_probability( SingleReaction ¤t_reaction, double E_pot_old, double E_pot_new, std::map &dummy_old_particle_numbers, int dummy_old_state_index, int dummy_new_state_index, bool dummy_only_make_configuration_changing_move) { - /** - *Calculates the expression in the acceptance probability of the constant pH - *method. - */ double ln_bf; double pKa; const double beta = 1.0 / temperature; From 28f0df94285569e916abe4fff9f648aad380053a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 30 Aug 2019 17:56:31 +0200 Subject: [PATCH 14/17] Re-shuffle Doxygen blocks --- src/core/reaction_ensemble.cpp | 39 ++++++---------------------------- src/core/reaction_ensemble.hpp | 21 ++++++++++++++++++ 2 files changed, 28 insertions(+), 32 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index b7a94a70f31..5af8d9e5ed9 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -16,15 +16,6 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ -/** reaction ensemble method according to smith94x for the reaction ensemble at - *constant volume and temperature, for the reaction ensemble at constant - *pressure additionally employ a barostat! NOTE: a chemical reaction consists of - *a forward and backward reaction. Here both reactions have to be defined - *separately. The extent of the reaction is here chosen to be +1. If the - *reaction trial move for a dissociation of HA is accepted then there is one - *more dissociated ion pair H+ and A- - */ - /** @file */ #include "reaction_ensemble.hpp" @@ -76,10 +67,10 @@ template bool is_in_list(T value, const std::vector &list) { return false; } -void EnergyCollectiveVariable::load_CV_boundaries( - WangLandauReactionEnsemble &m_current_wang_landau_system) { /**save minimum and maximum energies as a function of the other collective * variables under min_boundaries_energies, max_boundaries_energies **/ +void EnergyCollectiveVariable::load_CV_boundaries( + WangLandauReactionEnsemble &m_current_wang_landau_system) { m_current_wang_landau_system.do_energy_reweighting = true; // load energy boundaries from file @@ -416,7 +407,6 @@ void WangLandauReactionEnsemble::on_end_reaction(int &accepted_state) { update_wang_landau_potential_and_histogram(accepted_state); } -bool ReactionAlgorithm::generic_oneway_reaction(int reaction_id) { /** *generic one way reaction *A+B+...+G +... --> K+...X + Z +... @@ -430,6 +420,7 @@ bool ReactionAlgorithm::generic_oneway_reaction(int reaction_id) { *the box. Matching particles simply change the types. If there are more *reactants than products, old reactant particles are deleted. */ +bool ReactionAlgorithm::generic_oneway_reaction(int reaction_id) { SingleReaction ¤t_reaction = reactions[reaction_id]; current_reaction.tried_moves += 1; @@ -568,9 +559,7 @@ void ReactionAlgorithm::replace_particle(int p_id, int desired_type) { * Hides a particle from short ranged interactions and from the electrostatic * interaction. Additional hiding from interactions would need to be implemented * here. - */ -void ReactionAlgorithm::hide_particle(int p_id, int previous_type) { - /** + * *remove_charge and put type to a non existing one --> no interactions anymore *it is as if the particle was non existing (currently only type-based *interactions are switched off, as well as the electrostatic interaction) @@ -581,6 +570,7 @@ void ReactionAlgorithm::hide_particle(int p_id, int previous_type) { *there would be a need for a rule for such "collision" reactions (a reaction *like the one above). */ +void ReactionAlgorithm::hide_particle(int p_id, int previous_type) { auto part = get_particle_data(p_id); double d_min = distto(partCfg(), part.r.p, p_id); @@ -595,7 +585,6 @@ void ReactionAlgorithm::hide_particle(int p_id, int previous_type) { set_particle_type(p_id, non_interacting_type); } -int ReactionAlgorithm::delete_particle(int p_id) { /** * Deletes the particle with the given p_id and stores if it created a hole * at that position in the particle id range. This method is intended to @@ -604,7 +593,7 @@ int ReactionAlgorithm::delete_particle(int p_id) { * avoid the id * range becoming excessively huge. */ - +int ReactionAlgorithm::delete_particle(int p_id) { int old_max_seen_id = max_seen_particle; if (p_id == old_max_seen_id) { // last particle, just delete @@ -1319,9 +1308,9 @@ int WangLandauReactionEnsemble::do_reaction(int reaction_steps) { // boring helper functions + /**increase the Wang-Landau potential and histogram at the current nbar */ void WangLandauReactionEnsemble::update_wang_landau_potential_and_histogram( int index_of_state_after_acceptance_or_rejection) { - /**increase the Wang-Landau potential and histogram at the current nbar */ if (index_of_state_after_acceptance_or_rejection >= 0) { if (histogram[index_of_state_after_acceptance_or_rejection] >= 0) { histogram[index_of_state_after_acceptance_or_rejection] += 1; @@ -1689,20 +1678,6 @@ int ConstantpHEnsemble::get_random_valid_p_id() { return random_p_id; } -/** - * Constant-pH Ensemble, for derivation see Reed and Reed 1992 - * For the constant pH reactions you need to provide the deprotonation and - * afterwards the corresponding protonation reaction (in this order). If you - * want to deal with multiple reactions do it multiple times. Note that there is - * a difference in the usecase of the constant pH reactions and the above - * reaction ensemble. For the constant pH simulation directily the - * **apparent equilibrium constant which carries a unit** needs to be provided - * -- this is equivalent to the gamma of the reaction ensemble above, where the - * dimensionless reaction constant needs to be provided. Again: For the - * constant-pH algorithm not the dimensionless reaction constant needs to be - * provided here, but the apparent reaction constant. - */ - /** *Performs a reaction in the constant pH ensemble */ diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 99bd5f46957..6fe1d0928c1 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -226,6 +226,14 @@ class ReactionAlgorithm { ////////////////////////////////////////////////////////////////actual /// declaration of specific reaction algorithms +/** reaction ensemble method according to smith94x for the reaction ensemble at + *constant volume and temperature, for the reaction ensemble at constant + *pressure additionally employ a barostat! NOTE: a chemical reaction consists of + *a forward and backward reaction. Here both reactions have to be defined + *separately. The extent of the reaction is here chosen to be +1. If the + *reaction trial move for a dissociation of HA is accepted then there is one + *more dissociated ion pair H+ and A- + */ class ReactionEnsemble : public ReactionAlgorithm { public: ReactionEnsemble(int seed) : ReactionAlgorithm(seed) {} @@ -341,6 +349,19 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { double delta_CV); }; +/** + * Constant-pH Ensemble, for derivation see Reed and Reed 1992 + * For the constant pH reactions you need to provide the deprotonation and + * afterwards the corresponding protonation reaction (in this order). If you + * want to deal with multiple reactions do it multiple times. Note that there is + * a difference in the usecase of the constant pH reactions and the above + * reaction ensemble. For the constant pH simulation directily the + * **apparent equilibrium constant which carries a unit** needs to be provided + * -- this is equivalent to the gamma of the reaction ensemble above, where the + * dimensionless reaction constant needs to be provided. Again: For the + * constant-pH algorithm not the dimensionless reaction constant needs to be + * provided here, but the apparent reaction constant. + */ class ConstantpHEnsemble : public ReactionAlgorithm { public: ConstantpHEnsemble(int seed) : ReactionAlgorithm(seed) {} From 8d0081918380f9438fa5720ab7035f9c1580b48e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 30 Aug 2019 18:10:34 +0200 Subject: [PATCH 15/17] Update documentation --- src/core/reaction_ensemble.cpp | 125 ++++++++++++++++----------------- src/core/reaction_ensemble.hpp | 20 +++--- 2 files changed, 72 insertions(+), 73 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 5af8d9e5ed9..f36232674ab 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -16,6 +16,7 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . */ + /** @file */ #include "reaction_ensemble.hpp" @@ -67,8 +68,9 @@ template bool is_in_list(T value, const std::vector &list) { return false; } - /**save minimum and maximum energies as a function of the other collective - * variables under min_boundaries_energies, max_boundaries_energies **/ +/** Save minimum and maximum energies as a function of the other collective + * variables under min_boundaries_energies, max_boundaries_energies + */ void EnergyCollectiveVariable::load_CV_boundaries( WangLandauReactionEnsemble &m_current_wang_landau_system) { @@ -407,19 +409,18 @@ void WangLandauReactionEnsemble::on_end_reaction(int &accepted_state) { update_wang_landau_potential_and_histogram(accepted_state); } - /** - *generic one way reaction - *A+B+...+G +... --> K+...X + Z +... - *you need to use 2A --> B instead of A+A --> B since in the last case you - *assume - *distinctness of the particles, however both ways to describe the reaction - *are equivalent in the thermodynamic limit (large particle numbers) further - *it is crucial for the function in which order you provide the reactant and - *product types since particles will be replaced correspondingly! If there are - *less reactants than products, new product particles are created randomly in - *the box. Matching particles simply change the types. If there are more - *reactants than products, old reactant particles are deleted. - */ +/** + * Generic one way reaction + * A+B+...+G +... --> K+...X + Z +... + * You need to use 2A --> B instead of A+A --> B since in the last case you + * assume distinctness of the particles, however both ways to describe the + * reaction are equivalent in the thermodynamic limit (large particle numbers). + * Furthermore, it is crucial for the function in which order you provide the + * reactant and product types since particles will be replaced correspondingly! + * If there are less reactants than products, new product particles are created + * randomly in the box. Matching particles simply change the types. If there + * are more reactants than products, old reactant particles are deleted. + */ bool ReactionAlgorithm::generic_oneway_reaction(int reaction_id) { SingleReaction ¤t_reaction = reactions[reaction_id]; @@ -516,9 +517,9 @@ bool ReactionAlgorithm::generic_oneway_reaction(int reaction_id) { for (int p_ids_created_particle : p_ids_created_particles) { delete_particle(p_ids_created_particle); } - // 2)restore previously hidden reactant particles + // 2) restore previously hidden reactant particles restore_properties(hidden_particles_properties, number_of_saved_properties); - // 2)restore previously changed reactant particles + // 3) restore previously changed reactant particles restore_properties(changed_particles_properties, number_of_saved_properties); reaction_is_accepted = false; @@ -560,16 +561,17 @@ void ReactionAlgorithm::replace_particle(int p_id, int desired_type) { * interaction. Additional hiding from interactions would need to be implemented * here. * - *remove_charge and put type to a non existing one --> no interactions anymore - *it is as if the particle was non existing (currently only type-based - *interactions are switched off, as well as the electrostatic interaction) - *hide_particle() does not break bonds for simple reactions. as long as there - *are no reactions like 2A -->B where one of the reacting A particles occurs - *in the polymer (think of bond breakages if the monomer in the polymer gets - *deleted in the reaction). This constraint is not of fundamental reason, but - *there would be a need for a rule for such "collision" reactions (a reaction - *like the one above). - */ + * Removing the particle charge and changing its type to a non existing one + * deactivates all interactions with other particles, as if the particle was + * inexistent (currently only type-based interactions are switched off, as well + * as the electrostatic interaction). + * This function does not break bonds for simple reactions, as long as there + * are no reactions like 2A -->B where one of the reacting A particles occurs + * in the polymer (think of bond breakages if the monomer in the polymer gets + * deleted in the reaction). This constraint is not of fundamental reason, but + * there would be a need for a rule for such "collision" reactions (a reaction + * like the one above). + */ void ReactionAlgorithm::hide_particle(int p_id, int previous_type) { auto part = get_particle_data(p_id); @@ -585,14 +587,12 @@ void ReactionAlgorithm::hide_particle(int p_id, int previous_type) { set_particle_type(p_id, non_interacting_type); } - /** - * Deletes the particle with the given p_id and stores if it created a hole - * at that position in the particle id range. This method is intended to - * only - * delete unbonded particles since bonds are coupled to ids. This is used to - * avoid the id - * range becoming excessively huge. - */ +/** + * Deletes the particle with the given p_id and stores the id if the deletion + * created a hole in the particle id range. This method is intended to only + * delete unbonded particles since bonds are coupled to ids. This is used to + * avoid the id range becoming excessively huge. + */ int ReactionAlgorithm::delete_particle(int p_id) { int old_max_seen_id = max_seen_particle; if (p_id == old_max_seen_id) { @@ -728,19 +728,13 @@ int ReactionAlgorithm::create_particle(int desired_type) { // set velocities set_particle_v(p_id, vel); double d_min = distto(partCfg(), pos_vec, p_id); - if (d_min < exclusion_radius) - particle_inside_exclusion_radius_touched = - true; // setting of a minimal - // distance is allowed to - // avoid overlapping - // configurations if there is - // a repulsive potential. - // States with very high - // energies have a probability - // of almost zero and - // therefore do not contribute - // to ensemble averages. - + if (d_min < exclusion_radius) { + // setting of a minimal distance is allowed to avoid overlapping + // configurations if there is a repulsive potential. States with + // very high energies have a probability of almost zero and + // therefore do not contribute to ensemble averages. + particle_inside_exclusion_radius_touched = true; + } return p_id; } @@ -862,16 +856,17 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type( // symmetric } - // //correct for enhanced proposal of small radii by using the Metropolis - // hastings algorithm for asymmetric proposal densities. - // double - // old_radius=std::sqrt(std::pow(particle_positions[0]-cyl_x,2)+std::pow(particle_positions[1]-cyl_y,2)); - // double - // new_radius=std::sqrt(std::pow(new_pos[0]-cyl_x,2)+std::pow(new_pos[1]-cyl_y,2)); - // bf=std::min(1.0, - // bf*exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius); - ////Metropolis-Hastings Algorithm for asymmetric proposal density + // // correct for enhanced proposal of small radii by using the + // // Metropolis-Hastings algorithm for asymmetric proposal densities + // double old_radius = + // std::sqrt(std::pow(particle_positions[0]-cyl_x,2) + + // std::pow(particle_positions[1]-cyl_y,2)); + // double new_radius = + // std::sqrt(std::pow(new_pos[0]-cyl_x,2)+std::pow(new_pos[1]-cyl_y,2)); + // bf = std::min(1.0, + // bf*exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius); + // Metropolis-Hastings Algorithm for asymmetric proposal density if (m_uniform_real_distribution(m_generator) < bf) { // accept m_accepted_configurational_MC_moves += 1; @@ -943,7 +938,7 @@ int WangLandauReactionEnsemble::get_flattened_index_wang_landau( int index = -10; // negative number is not allowed as index and therefore // indicates error std::vector individual_indices(nr_collective_variables); // pre result - // individual_indices.resize(nr_collective_variables,-1); //initialize + // individual_indices.resize(nr_collective_variables,-1); //initialize // individual_indices to -1 // check for the current state to be an allowed state in the [range @@ -1308,7 +1303,7 @@ int WangLandauReactionEnsemble::do_reaction(int reaction_steps) { // boring helper functions - /**increase the Wang-Landau potential and histogram at the current nbar */ +/** Increase the Wang-Landau potential and histogram at the current nbar */ void WangLandauReactionEnsemble::update_wang_landau_potential_and_histogram( int index_of_state_after_acceptance_or_rejection) { if (index_of_state_after_acceptance_or_rejection >= 0) { @@ -1705,10 +1700,10 @@ int ConstantpHEnsemble::do_reaction(int reaction_steps) { for (int reaction_i = 0; reaction_i < reactions.size(); reaction_i++) { SingleReaction ¤t_reaction = reactions[reaction_i]; for (int reactant_i = 0; reactant_i < 1; - reactant_i++) { // reactant_i<1 since it is assumed in this place - // that the types A, and HA occur in the first place - // only. These are the types that should be switched, - // H+ should not be switched + reactant_i++) { // reactant_i < 1 since it is assumed in this place + // that the types A and HA occur in the first place + // only. These are the types that should be + // switched, H+ should not be switched if (current_reaction.reactant_types[reactant_i] == type_of_random_p_id) { list_of_reaction_ids_with_given_reactant_type.push_back(reaction_i); @@ -1765,9 +1760,9 @@ WidomInsertion::measure_excess_chemical_potential(int reaction_id) { for (int p_ids_created_particle : p_ids_created_particles) { delete_particle(p_ids_created_particle); } - // 2)restore previously hidden reactant particles + // 2) restore previously hidden reactant particles restore_properties(hidden_particles_properties, number_of_saved_properties); - // 2)restore previously changed reactant particles + // 3) restore previously changed reactant particles restore_properties(changed_particles_properties, number_of_saved_properties); std::vector exponential = { exp(-1.0 / temperature * (E_pot_new - E_pot_old))}; diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 6fe1d0928c1..3b9095a8a8a 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -106,6 +106,7 @@ struct DegreeOfAssociationCollectiveVariable : public CollectiveVariable { } }; +/** Base class for reaction ensemble methods */ class ReactionAlgorithm { public: @@ -226,13 +227,14 @@ class ReactionAlgorithm { ////////////////////////////////////////////////////////////////actual /// declaration of specific reaction algorithms -/** reaction ensemble method according to smith94x for the reaction ensemble at - *constant volume and temperature, for the reaction ensemble at constant - *pressure additionally employ a barostat! NOTE: a chemical reaction consists of - *a forward and backward reaction. Here both reactions have to be defined - *separately. The extent of the reaction is here chosen to be +1. If the - *reaction trial move for a dissociation of HA is accepted then there is one - *more dissociated ion pair H+ and A- +/** Reaction ensemble method according to smith94x. + * Works for the reaction ensemble at constant volume and temperature. For the + * reaction ensemble at constant pressure additionally employ a barostat! + * NOTE: a chemical reaction consists of a forward and backward reaction. + * Here both reactions have to be defined separately. The extent of the + * reaction is here chosen to be +1. If the reaction trial move for a + * dissociation of HA is accepted then there is one more dissociated ion + * pair H+ and A- */ class ReactionEnsemble : public ReactionAlgorithm { public: @@ -246,6 +248,7 @@ class ReactionEnsemble : public ReactionAlgorithm { bool dummy_only_make_configuration_changing_move) override; }; +/** Wang-Landau reaction ensemble method */ class WangLandauReactionEnsemble : public ReactionAlgorithm { public: WangLandauReactionEnsemble(int seed) : ReactionAlgorithm(seed) {} @@ -350,7 +353,7 @@ class WangLandauReactionEnsemble : public ReactionAlgorithm { }; /** - * Constant-pH Ensemble, for derivation see Reed and Reed 1992 + * Constant-pH Ensemble, for derivation see Reed and Reed 1992. * For the constant pH reactions you need to provide the deprotonation and * afterwards the corresponding protonation reaction (in this order). If you * want to deal with multiple reactions do it multiple times. Note that there is @@ -377,6 +380,7 @@ class ConstantpHEnsemble : public ReactionAlgorithm { int get_random_valid_p_id(); }; +/** Widom insertion method */ class WidomInsertion : public ReactionAlgorithm { public: WidomInsertion(int seed) : ReactionAlgorithm(seed) {} From 8fb80978fa951c6fb62e74731cb0b28e4919dc73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 30 Aug 2019 18:21:41 +0200 Subject: [PATCH 16/17] Formatting and numpydoc style --- samples/grand_canonical.py | 16 +- samples/widom_insertion.py | 9 +- src/python/espressomd/reaction_ensemble.pyx | 197 ++++++++++---------- 3 files changed, 117 insertions(+), 105 deletions(-) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index ca0f8c56b16..923bde2f5dd 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -17,9 +17,16 @@ # along with this program. If not, see . # """ -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. +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 argparse @@ -81,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) diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index 96b39c22076..da01deecfe2 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -17,8 +17,10 @@ # along with this program. If not, see . # """ -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. +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 argparse @@ -79,8 +81,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") p3m = electrostatics.P3M(prefactor=2.0, accuracy=1e-3) system.actors.add(p3m) diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index be2034a8edb..01219a75329 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -13,10 +13,10 @@ class WangLandauHasConverged(Exception): cdef class ReactionAlgorithm: """ - This class provides the base class for Reaction Algorithms like the Reaction Ensemble algorithm, the Wang-Landau - Reaction Ensemble algorithm and the constant pH method. Initialize the - reaction algorithm by setting the standard pressure, temperature, and the - exclusion radius. + This class provides the base class for Reaction Algorithms like the Reaction + Ensemble algorithm, the Wang-Landau Reaction Ensemble algorithm and the + constant pH method. Initialize the reaction algorithm by setting the + standard pressure, temperature, and the exclusion radius. Note: When creating particles the velocities of the new particles are set according the Maxwell-Boltzmann distribution. In this step the mass of the @@ -26,16 +26,15 @@ cdef class ReactionAlgorithm: Parameters ---------- temperature : :obj:`float` - The temperature at which the reaction is performed. + The temperature at which the reaction is performed. exclusion_radius : :obj:`float` - Minimal distance from any particle, within which new - particle will not be inserted. This is useful to avoid - integrator failures if particles are too close and there - is a diverging repulsive interaction, or to prevent two - oppositely charged particles from being placed on top of - each other. The Boltzmann factor :math:`\exp(-\\beta - E)` gives these configurations a small contribution to - the partition function, therefore they can be neglected. + Minimal distance from any particle, within which new particle will not + be inserted. This is useful to avoid integrator failures if particles + are too close and there is a diverging repulsive interaction, or to + prevent two oppositely charged particles from being placed on top of + each other. The Boltzmann factor :math:`\exp(-\\beta E)` gives these + configurations a small contribution to the partition function, + therefore they can be neglected. """ cdef object _params cdef CReactionAlgorithm * RE @@ -47,8 +46,7 @@ cdef class ReactionAlgorithm: return "temperature", "exclusion_radius", "seed" def _set_params_in_es_core(self): - deref(self.RE).temperature = self._params[ - "temperature"] + deref(self.RE).temperature = self._params["temperature"] # setting a volume is a side effect, sets the default volume of the # reaction ensemble as the volume of the cuboid simulation box. this # volume can be altered by the command "reaction ensemble volume @@ -57,8 +55,7 @@ cdef class ReactionAlgorithm: # only contained in the volume of the cylinder) if(deref(self.RE).volume < 0): deref(self.RE).set_cuboid_reaction_ensemble_volume() - deref(self.RE).exclusion_radius = self._params[ - "exclusion_radius"] + deref(self.RE).exclusion_radius = self._params["exclusion_radius"] def set_cylindrical_constraint_in_z_direction(self, center_x, center_y, radius_of_cylinder): @@ -70,11 +67,11 @@ cdef class ReactionAlgorithm: Parameters ---------- center_x : :obj:`float` - x coordinate of center of the cylinder. + x coordinate of center of the cylinder. center_y : :obj:`float` - y coordinate of center of the cylinder. + y coordinate of center of the cylinder. radius_of_cylinder : :obj:`float` - radius of the cylinder + radius of the cylinder """ deref(self.RE).cyl_x = center_x @@ -135,14 +132,14 @@ cdef class ReactionAlgorithm: """ Sets the particle type for non-interacting particles. Default value: 100. - This is used to temporarily hide - particles during a reaction trial move, if they are to be deleted after - the move is accepted. - Please change this value if you intend to use the type 100 for some other - particle types with interactions. Please also note that particles in the current - implementation of the Reaction Ensemble are only hidden with respect to - Lennard-Jones and Coulomb interactions. Hiding of other interactions, - for example a magnetic, needs to be implemented in the code. + This is used to temporarily hide particles during a reaction trial + move, if they are to be deleted after the move is accepted. Please + change this value if you intend to use the type 100 for some other + particle types with interactions. Please also note that particles + in the current implementation of the Reaction Ensemble are only + hidden with respect to Lennard-Jones and Coulomb interactions. Hiding + of other interactions, for example a magnetic, needs to be implemented + in the code. """ deref(self.RE).non_interacting_type = non_interacting_type @@ -160,24 +157,22 @@ cdef class ReactionAlgorithm: Parameters ---------- gamma : :obj:`float` - Equilibrium constant of the reaction, :math:`\gamma` - (see the User guide, section 6.6 for the definition and further details). + Equilibrium constant of the reaction, :math:`\gamma` (see the User + guide, section 6.6 for the definition and further details). reactant_types : list of :obj:`int` - List of particle types of reactants in the reaction. + List of particle types of reactants in the reaction. reactant_coefficients : list - List of stoichiometric coefficients of the - reactants in the same order as the list of - their types. + List of stoichiometric coefficients of the reactants in the same + order as the list of their types. product_types : list - List of particle types of products in the reaction. + List of particle types of products in the reaction. product_coefficients : list - List of stoichiometric coefficients of - products of the reaction in the same order as - the list of their types + List of stoichiometric coefficients of products of the reaction in + the same order as the list of their types default_charges : dictionary - A dictionary of default charges for types that occur in the provided reaction. + A dictionary of default charges for types that occur in the provided reaction. check_for_electroneutrality : Boolean - Check for electroneutrality of the given reaction if True. + Check for electroneutrality of the given reaction if True. """ self._params["check_for_electroneutrality"] = True @@ -267,7 +262,7 @@ cdef class ReactionAlgorithm: Parameters ---------- reaction_steps : :obj:`int`, optional - The number of reactions to be performed at once, defaults to 1. + The number of reactions to be performed at once, defaults to 1. """ deref(self.RE).do_reaction(int(reaction_steps)) @@ -275,10 +270,11 @@ cdef class ReactionAlgorithm: def displacement_mc_move_for_particles_of_type(self, type_mc, particle_number_to_be_changed=1): """ - Performs a displacement Monte Carlo move for particles of given type. New positions - of the displaced particles are chosen from the whole box with a uniform probability distribution. - If there are multiple types, that are being moved in a simulation, they should be moved in a - random order to avoid artefacts. + Performs a displacement Monte Carlo move for particles of given type. + New positions of the displaced particles are chosen from the whole box + with a uniform probability distribution. If there are multiple types, + that are being moved in a simulation, they should be moved in a random + order to avoid artefacts. Parameters ---------- @@ -317,8 +313,12 @@ cdef class ReactionAlgorithm: for i in range(deref(self.RE).reactions[single_reaction_i].product_types.size()): product_coefficients.append( deref(self.RE).reactions[single_reaction_i].product_coefficients[i]) - reaction = {"reactant_coefficients": reactant_coefficients, "reactant_types": reactant_types, "product_types": product_types, "product_coefficients": - product_coefficients, "reactant_types": reactant_types, "gamma": deref(self.RE).reactions[single_reaction_i].gamma} + reaction = {"reactant_coefficients": reactant_coefficients, + "reactant_types": reactant_types, + "product_types": product_types, + "product_coefficients": product_coefficients, + "reactant_types": reactant_types, + "gamma": deref(self.RE).reactions[single_reaction_i].gamma} reactions.append(reaction) return {"reactions": reactions, "temperature": deref(self.RE).temperature, "exclusion_radius": deref(self.RE).exclusion_radius} @@ -434,11 +434,12 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): def reaction(self, reaction_steps=1): """ - Performs reaction_steps reactions. Sets the number of reaction steps which are - performed at once. Do not use too many reaction steps - steps consecutively without having conformation - changing steps in between (especially important for the Wang-Landau reaction ensemble). Providing a number for the parameter reaction steps reduces the need for the interpreter to be - called between consecutive reactions. + Performs reaction_steps reactions. Sets the number of reaction steps + which are performed at once. Do not use too many reaction steps steps + consecutively without having conformation changing steps in between + (especially important for the Wang-Landau reaction ensemble). Providing + a number for the parameter reaction steps reduces the need for the + interpreter to be called between consecutive reactions. """ status_wang_landau = deref( @@ -455,14 +456,13 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): Parameters ---------- associated_type : :obj:`int` - Particle type of the associated state of the reacting species. + Particle type of the associated state of the reacting species. min : :obj:`float` - Minimum value of the collective variable. + Minimum value of the collective variable. max : :obj:`float` - Maximum value of the collective variable. + Maximum value of the collective variable. corresponding_acid_types : list - List of the types of the version of the - species. + List of the types of the version of the species. """ for k in kwargs: @@ -498,23 +498,23 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): Parameters ---------- filename : :obj:`str` - Filename of the energy boundary file which provides the - potential energy boundaries (min E_pot, max E_pot) tabulated - for all degrees of association. Make sure to only list the - degrees of association which are used by the degree of - association collective variable within this file. The energy - boundary file can be created in a preliminary energy run. By - the help of the functions - :meth:`update_maximum_and_minimum_energies_at_current_state` - and :meth:`write_out_preliminary_energy_run_results`. This - file has to be obtained before being able to run a - simulation with the energy as collective variable. + Filename of the energy boundary file which provides the + potential energy boundaries (min E_pot, max E_pot) tabulated + for all degrees of association. Make sure to only list the + degrees of association which are used by the degree of + association collective variable within this file. The energy + boundary file can be created in a preliminary energy run. By + the help of the functions + :meth:`update_maximum_and_minimum_energies_at_current_state` + and :meth:`write_out_preliminary_energy_run_results`. This + file has to be obtained before being able to run a + simulation with the energy as collective variable. delta : :obj:`float` - Provides the discretization of the potential energy range. Only - for small enough delta the results of the energy reweighted - averages are correct. If delta is chosen too big there are - discretization errors in the numerical integration which occurs - during the energy reweighting process. + Provides the discretization of the potential energy range. Only + for small enough delta the results of the energy reweighted + averages are correct. If delta is chosen too big there are + discretization errors in the numerical integration which occurs + during the energy reweighting process. """ for k in kwargs: @@ -546,25 +546,17 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): Parameters ---------- final_wang_landau_parameter : :obj:`float` - Sets the final Wang-Landau parameter, which is the Wang-Landau parameter after which the simulation should stop.). + Sets the final Wang-Landau parameter, which is the Wang-Landau + parameter after which the simulation should stop.). full_path_to_output_filename : :obj:`str` - Sets the path to the output file of the - Wang-Landau algorithm which contains the - Wang-Landau potential + Sets the path to the output file of the Wang-Landau algorithm which + contains the Wang-Landau potential do_not_sample_reaction_partition_function : :obj:`bool` - Avoids sampling the - Reaction ensemble partition - function in the Wang-Landau - algorithm. Therefore this - option makes all degrees of - association equally - probable. This option may - be used in the sweeping - mode of the reaction - ensemble, since the - reaction ensemble partition - function can be later added - analytically. + Avoids sampling the Reaction ensemble partition function in the + Wang-Landau algorithm. Therefore this option makes all degrees of + association equally probable. This option may be used in the + sweeping mode of the reaction ensemble, since the reaction ensemble + partition function can be later added analytically. """ for k in kwargs: @@ -638,10 +630,18 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): def displacement_mc_move_for_particles_of_type(self, type_mc, particle_number_to_be_changed=1): """ - Performs an MC (Monte Carlo) move for particle_number_to_be_changed particle of type type_mc. Positions for the particles are drawn uniformly random within the box. The command takes into account the Wang-Landau terms in the acceptance probability. - If there are multiple types, that - need to be moved, make sure to move them in a random order to avoid - artefacts. For the Wang-Landau algorithm in the case of energy reweighting you would also need to move the monomers of the polymer with special moves for the MC part. Those polymer configuration changing moves need to be implemented in the case of using Wang-Landau with energy reweighting and a polymer in the system. Polymer configuration changing moves had been implemented before but were removed from espresso. + Performs an MC (Monte Carlo) move for particle_number_to_be_changed + particle of type type_mc. Positions for the particles are drawn + uniformly random within the box. The command takes into account the + Wang-Landau terms in the acceptance probability. + If there are multiple types, that need to be moved, make sure to move + them in a random order to avoid artefacts. For the Wang-Landau algorithm + in the case of energy reweighting you would also need to move the + monomers of the polymer with special moves for the MC part. Those + polymer configuration changing moves need to be implemented in the + case of using Wang-Landau with energy reweighting and a polymer in the + system. Polymer configuration changing moves had been implemented + before but were removed from espresso. """ use_wang_landau = True @@ -651,7 +651,9 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): cdef class WidomInsertion(ReactionAlgorithm): """ - This class implements the Widom insertion method in the canonical ensemble for homogeneous systems, where the excess chemical potential is not depending on the location. + This class implements the Widom insertion method in the canonical ensemble + for homogeneous systems, where the excess chemical potential is not + depending on the location. """ @@ -693,7 +695,10 @@ cdef class WidomInsertion(ReactionAlgorithm): def measure_excess_chemical_potential(self, reaction_id=0): """ - Measures the excess chemical potential in a homogeneous system. Returns the excess chemical potential and the standard error for the excess chemical potential. It assumes that your samples are uncorrelated in estimating the standard error. + Measures the excess chemical potential in a homogeneous system. + Returns the excess chemical potential and the standard error for the + excess chemical potential. It assumes that your samples are + uncorrelated in estimating the standard error. """ return self.WidomInsertionPtr.get().measure_excess_chemical_potential(int(reaction_id)) From 6e15c453649e99cdae2af68c76b99c9bfb962b69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 30 Aug 2019 18:34:27 +0200 Subject: [PATCH 17/17] Fix typos and Sphinx/Doxygen issues --- src/core/reaction_ensemble.cpp | 2 +- src/core/reaction_ensemble.hpp | 3 +- src/python/espressomd/reaction_ensemble.pyx | 32 ++++++++++----------- 3 files changed, 18 insertions(+), 19 deletions(-) diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index f36232674ab..ec63e5ac772 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -58,7 +58,7 @@ double average_list_of_allowed_entries(const std::vector &vector) { } /** - * Checks whether a number is in a std:vector of numbers. + * Checks whether a number is in a std::vector of numbers. */ template bool is_in_list(T value, const std::vector &list) { for (int i = 0; i < list.size(); i++) { diff --git a/src/core/reaction_ensemble.hpp b/src/core/reaction_ensemble.hpp index 3b9095a8a8a..c9f50dae8bd 100644 --- a/src/core/reaction_ensemble.hpp +++ b/src/core/reaction_ensemble.hpp @@ -224,8 +224,7 @@ class ReactionAlgorithm { get_random_position_in_box_enhanced_proposal_of_small_radii(); }; -////////////////////////////////////////////////////////////////actual -/// declaration of specific reaction algorithms +///////////////////////////// actual declaration of specific reaction algorithms /** Reaction ensemble method according to smith94x. * Works for the reaction ensemble at constant volume and temperature. For the diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index 01219a75329..6eb1bafbc95 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -32,7 +32,7 @@ cdef class ReactionAlgorithm: be inserted. This is useful to avoid integrator failures if particles are too close and there is a diverging repulsive interaction, or to prevent two oppositely charged particles from being placed on top of - each other. The Boltzmann factor :math:`\exp(-\\beta E)` gives these + each other. The Boltzmann factor :math:`\\exp(-\\beta E)` gives these configurations a small contribution to the partition function, therefore they can be neglected. """ @@ -157,21 +157,21 @@ cdef class ReactionAlgorithm: Parameters ---------- gamma : :obj:`float` - Equilibrium constant of the reaction, :math:`\gamma` (see the User + Equilibrium constant of the reaction, :math:`\\gamma` (see the User guide, section 6.6 for the definition and further details). reactant_types : list of :obj:`int` List of particle types of reactants in the reaction. - reactant_coefficients : list + reactant_coefficients : list of :obj:`int` List of stoichiometric coefficients of the reactants in the same order as the list of their types. - product_types : list + product_types : list of :obj:`int` List of particle types of products in the reaction. - product_coefficients : list + product_coefficients : list of :obj:`int` List of stoichiometric coefficients of products of the reaction in the same order as the list of their types default_charges : dictionary A dictionary of default charges for types that occur in the provided reaction. - check_for_electroneutrality : Boolean + check_for_electroneutrality : :obj:`bool` Check for electroneutrality of the given reaction if True. """ @@ -356,7 +356,7 @@ cdef class ReactionEnsemble(ReactionAlgorithm): if k in self._valid_keys(): self._params[k] = kwargs[k] else: - raise KeyError("%s is not a vaild key" % k) + raise KeyError("%s is not a valid key" % k) self._set_params_in_es_core() @@ -379,7 +379,7 @@ cdef class ConstantpHEnsemble(ReactionAlgorithm): if k in self._valid_keys(): self._params[k] = kwargs[k] else: - raise KeyError("%s is not a vaild key" % k) + raise KeyError("%s is not a valid key" % k) self._set_params_in_es_core() @@ -425,7 +425,7 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): if k in self._valid_keys(): self._params[k] = kwargs[k] else: - raise KeyError("%s is not a vaild key" % k) + raise KeyError("%s is not a valid key" % k) self.WLRptr.reset(new CWangLandauReactionEnsemble(int(self._params["seed"]))) self.RE = self.WLRptr.get() @@ -435,8 +435,8 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): def reaction(self, reaction_steps=1): """ Performs reaction_steps reactions. Sets the number of reaction steps - which are performed at once. Do not use too many reaction steps steps - consecutively without having conformation changing steps in between + which are performed at once. Do not use too many reaction steps + consecutively without having conformation-changing steps in-between (especially important for the Wang-Landau reaction ensemble). Providing a number for the parameter reaction steps reduces the need for the interpreter to be called between consecutive reactions. @@ -461,7 +461,7 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): Minimum value of the collective variable. max : :obj:`float` Maximum value of the collective variable. - corresponding_acid_types : list + corresponding_acid_types : list of :obj:`int` List of the types of the version of the species. """ @@ -547,7 +547,7 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): ---------- final_wang_landau_parameter : :obj:`float` Sets the final Wang-Landau parameter, which is the Wang-Landau - parameter after which the simulation should stop.). + parameter after which the simulation should stop. full_path_to_output_filename : :obj:`str` Sets the path to the output file of the Wang-Landau algorithm which contains the Wang-Landau potential @@ -598,9 +598,9 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm): Records the minimum and maximum potential energy as a function of the degree of association in a preliminary Wang-Landau reaction ensemble simulation where the acceptance probability includes the factor - :math:`\exp(-\\beta \\Delta E_{pot})`. The minimal and maximal + :math:`\\exp(-\\beta \\Delta E_{pot})`. The minimal and maximal potential energies which occur in the system are needed for the energy - reweighting simulations where the factor :math:`\exp(-\\beta \\Delta E_{pot})` + reweighting simulations where the factor :math:`\\exp(-\\beta \\Delta E_{pot})` is not included in the acceptance probability in order to avoid choosing the wrong potential energy boundaries. @@ -689,7 +689,7 @@ cdef class WidomInsertion(ReactionAlgorithm): if k in self._valid_keys(): self._params[k] = kwargs[k] else: - raise KeyError("%s is not a vaild key" % k) + raise KeyError("%s is not a valid key" % k) self._set_params_in_es_core()