Corrected reaction ensemble documentation, small changes in RE syntax… #1741

kosovan · 2017-12-28T21:52:16Z

…, and corrections to teststuite and samples to comply with the current syntax.

Partly Fixes #1650 Some issues still need to be clarified.

Description of changes:

fixed the description of the standard_pressure
rephrased documentation of several other functions and parameters
changed RE.add() to RE.add_reaction() to make the semantics consistent with RE.add_collective_variable..()
adopted the samples to the current state of Reaction Ensemble code

PR Checklist

Tests?
- Interface
- Core
Docs?

…, and corrections to teststuite and samples to comply with the current syntax.

codecov · 2017-12-28T22:32:03Z

Codecov Report

Merging #1741 into python will increase coverage by <1%.
The diff coverage is 75%.

@@           Coverage Diff           @@
##           python   #1741    +/-   ##
=======================================
+ Coverage      59%     60%   +<1%     
=======================================
  Files         416     416            
  Lines       28815   28965   +150     
=======================================
+ Hits        17288   17386    +98     
- Misses      11527   11579    +52

Impacted Files	Coverage Δ
src/core/reaction_ensemble.hpp	`76% <ø> (ø)`	⬆️
src/core/reaction_ensemble.cpp	`74% <75%> (ø)`	⬆️
src/core/rattle.cpp	`91% <0%> (-2%)`	⬇️
src/core/p3m-dipolar.cpp	`79% <0%> (-2%)`	⬇️
src/core/integrate.cpp	`73% <0%> (-1%)`	⬇️
src/core/forces_inline.hpp	`64% <0%> (-1%)`	⬇️
src/core/metadynamics.cpp	`6% <0%> (-1%)`	⬇️
src/core/energy_inline.hpp	`51% <0%> (-1%)`	⬇️
src/core/magnetic_non_p3m_methods.cpp	`7% <0%> (-1%)`	⬇️
src/core/topology.cpp	`3% <0%> (-1%)`	⬇️
... and 20 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update f74064d...fca320a. Read the comment docs.

jonaslandsgesell · 2017-12-30T17:04:42Z

src/python/espressomd/reaction_ensemble.pyx

+                        Note that if you provide the reaction constant :math:`K_c`
+                        as described in the user guide, this parameter should be kept at its default value.
+                        This paramater is retained for backward compatibility with earlier
+                        implementation, where the pressure-based constaht :math:`K_p` was


The user guide only deals with a dimensionless reaction constant K (K_p is not mentioned there). Please also name it this way here.

Sorry i overlooked parts. But how is K=Exp(...) related to K_p?

I don't see that there is a statement K=K_p although it is used in that way. After all this only holds for ideal gases but all the reaction constants which are provided in the reaction ensemble are ideal Gas reaction constants. Adding Interactions to the system effectively alters reaction constant K_c which is observed in the system then. Can we introduce a statement like this or so you see Problems with it?

You are right, K_p is not defined anywhere. At the same time, though K_p = prod_i (p_i / p_standard)^{nu_i} is well defined only if the system is a gas, and it is equal to K only if the system is an ideal gas. Regarding the message of this note, saying that standard_pressure is needed only in gas-phase reactions, I would leave it as is. I expect that most users need not bother, and those who need to study gas-phase reactions should be familiar with K_p.

jonaslandsgesell · 2017-12-30T17:08:09Z

src/python/espressomd/reaction_ensemble.pyx

@@ -400,12 +397,12 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):

    def add_collective_variable_degree_of_association(self, *args, **kwargs):
        """
-        Adds a reaction coordinate of the type degree of association.
+        Adds degree of association as the collective variable (reaction coordinate) for the Wang-Landau Reaction Ensemble.


You can have multiple reaction coordinates of this type. Therefore it was previsously "adds a reaction coordinate of the type..." implying that you can have more than one reaction coordinates of this type.

Therefore I would explicitly state that you can have more than one reaction coordinate of this type.

Agreed. I will add a note that several collective variables can be set simultaneously.

…docs_reaction_ensemble

…ow it is implemented

…ous conventional definitions of equilibrium constant.

…docs_reaction_ensemble

* Replaced equilibrium_constant with Gamma in testsuite. * Modified the scripts such that they do not set the standard pressure but still produce the same numbers as before

… everywhere in the code.

kosovan · 2018-01-17T16:11:00Z

@jonaslandsgesell and @helvrud please check my changes and comment if any further changes are needed.

I think that together with #1767 this is the final fix for #1650 After these two PRs are merged, this issue can be closed.

jonaslandsgesell · 2018-01-17T19:27:18Z

Looks good to me

jonaslandsgesell · 2018-01-17T19:07:00Z

doc/sphinx/run.rst

+
+where :math:`q_i` and :math:`\Lambda_i` are the single-particle partition function and
+thermal wavelength of species :math:`i`. In the absence of intermolecular interactions
+(ideal gas or ideal solution), :math:`\Gamma` proportional to the concentration-based reaction constant


Gamma is the concentration based reaction constant in units of 1/(sigma^3)^nu_bar, if no interactions are present. Proportional is too weak for my taste.

I agree but I will not specify the units because these are defined by the simulation script. The units have to be consistent with other units in the simulation script. Can be mol/L or 1/sigma^3 or anything else.

It all depends on the unit of the volume that is provided that is correct. However for the case that the user does not provide a volume (which is standard and V=box_l_xbox_l_ybox_l_z, V is in units of sigma^3) and then as a standard the unit of gamma is 1/(sigma^3)^nubabr

The sentence before says: "Note that the dimension of :math:\Gamma is :math:V^{\bar\nu}, therefore its units must be consistent with the units in which Espresso measures the box volume, i.e. :math:\sigma^3." Isn't this precisely the same thing as you wrote, just spelled differently?

jonaslandsgesell · 2018-01-17T19:16:04Z

src/core/reaction_ensemble.hpp

@@ -89,7 +89,8 @@ class ReactionAlgorithm {

  std::vector<SingleReaction> reactions;
  std::map<int, double> charges_of_types;
-  double standard_pressure_in_simulation_units = -10.0;
+  //double standard_pressure_in_simulation_units = -10.0;
+  //double standard_pressure_in_simulation_units = 1.0;


Please remove the comment

jonaslandsgesell · 2018-01-17T19:20:43Z

testsuite/constant_pH.py

    # avoid extreme regions in the titration curve e.g. via the choice
    # (np.random.random()-0.5)
+#    target_alpha=0.5; # choose target alpha close to 0.5 to get good statistics in a small number of steps
+#    Ka=target_alpha*target_alpha/(1.0-target_alpha)*c0; # determine Ka which corresponds to the desired c0 and alpha
+#    pKa=-np.log10(Ka);


I did not choose 0.5 because this value was quite robust regarding errors from my experience for typical errors which happened during coding. But of course in principle any point on the curve should be suitable for testing, so I don't have a preference.

jonaslandsgesell · 2018-01-17T19:21:35Z

testsuite/constant_pH.py

    pKa_minus_pH = 1
    pH = 2  # or randomly via: 4*np.random.random()
    pKa = pKa_minus_pH + pH
    # could be in this test for example anywhere in the range 0.000001 ... 9
    K_HA_diss_apparent = 10**(-pKa)
+    #K_HA_diss_apparent = 10**(-pKa)*0.00108; # first, ensure that by removing standard pressure we obtain the same hard-coded numbers as before
+    #print("pKa:", pKa, "pk_pKa:", -np.log10(K_HA_diss_pk));


Please remove

jonaslandsgesell · 2018-01-17T19:23:00Z

testsuite/reaction_ensemble.py

    exclusion_radius = 1.0
    # could be in this test for example anywhere in the range 0.000001 ... 9,
    # chosen for example like 8.8*np.random.random()+0.2
-    K_HA_diss = 8.8 * 0.5 + 0.2
+    # FIXME choose desired alpha and then calculate Ka based on the desired alpha
+    # use alpha close to 0.5 to get good statistics in few steps, 


Please don't leave comments.

jonaslandsgesell · 2018-01-22T10:03:13Z

testsuite/reaction_ensemble.py

-            pK_a, pH)
+        pK_a = -np.log10(K_HA_diss)
+        #ideal=ReactionEnsembleTest.ideal_degree_of_association( pK_a, pH);
+        ideal=0.05080885682460401; # hard-coded value from the original setup by Jonas


I would not hard code these values but use the function that determines these values

jonaslandsgesell · 2018-01-22T10:04:20Z

testsuite/reaction_ensemble.py

+        pK_a = -np.log10(K_HA_diss)
+        #ideal=ReactionEnsembleTest.ideal_degree_of_association( pK_a, pH);
+        ideal=0.05080885682460401; # hard-coded value from the original setup by Jonas
+        target=0.05024999999999979; # hard-coded value from the original setup by Jonas


what is target used for? It could be removed. That s where the Mac Os test is currently failing.

…cs_reaction_ensemble

jonaslandsgesell · 2018-02-01T15:33:35Z

testsuite/constant_pH.py

                cls.type_HA], reactant_coefficients=[1], product_types=[
                cls.type_A, cls.type_H], product_coefficients=[
                1, 1], default_charges={cls.type_HA: 0, cls.type_A: -1, cls.type_H: +1})
        cls.RE.constant_pH = cls.pH

    @classmethod
-    def ideal_degree_of_association(cls, pH):
+    def ideal_alpha(cls, pH):


this is the formula for the degree of association= 1-degree of dissociation=1-alpha, please dont use alpha but a more descriptive name!

Looks to me like exchanging the degree of association and the degree of dissociation works here because both are 0.5 here. Therefore I would choose a point in the parameter space where both are distinct.

jonaslandsgesell · 2018-02-01T15:33:56Z

testsuite/constant_pH.py

-    pH = pKa_minus_pH  # or randomly via: 4*np.random.random()
-    pH = 2  # or randomly via: 4*np.random.random()
+    pKa_minus_pH = 0
+    pH = pKa_minus_pH  


this line is overwritten by the next line

jonaslandsgesell · 2018-02-01T15:37:41Z

testsuite/reaction_ensemble.py

+        average_NA /= num_samples
+        average_NHA /= num_samples
+        average_alpha = average_NA / float(N0)
+        print("average_NH:", average_NH,


please dont print in tests

jonaslandsgesell · 2018-02-01T15:39:49Z

testsuite/reaction_ensemble.py

@@ -81,7 +81,8 @@ def setUpClass(cls):
            reactant_types=cls.reactant_types,
            reactant_coefficients=cls.reactant_coefficients,
            product_types=cls.product_types,
-            product_coefficients=cls.product_coefficients, default_charges={cls.type_HA: 0, cls.type_A: -1, cls.type_H: +1}, check_for_electroneutrality=True)
+            product_coefficients=cls.product_coefficients, 
+            default_charges={cls.type_HA: 0, cls.type_A: -1, cls.type_H: +1}, check_for_electroneutrality=True)

    @classmethod
    def ideal_alpha(cls, Gamma,N0,V,nubar):


jonaslandsgesell · 2018-02-02T14:38:53Z

The reaction ensemble method was moved by @KonradBreitsprecher to the appendix since it was not obvious that the code is documented. If you add a minimal example to the documentation on how to setup a reaction in the rst you can move it to the advanced methods section.

…opefully no more failures will follow.

jonaslandsgesell

Looks good to me

fweik · 2018-02-06T15:07:42Z

https://gitlab.icp.uni-stuttgart.de/espressomd/espresso/pipelines/1843

Corrected reaction ensemble documentation, small changes in RE syntax…

Corrected reaction ensemble documentation, small changes in RE syntax…

62386f9

…, and corrections to teststuite and samples to comply with the current syntax.

hmenke approved these changes Dec 29, 2017

View reviewed changes

jonaslandsgesell suggested changes Dec 30, 2017

View reviewed changes

fweik assigned jonaslandsgesell Jan 4, 2018

kosovan and others added 5 commits January 8, 2018 21:34

Merge branch 'python' of https://github.com/espressomd/espresso into …

614ec07

…docs_reaction_ensemble

Merge branch 'python' of https://github.com/espressomd/espresso into …

d5abe8e

…docs_reaction_ensemble

Some RE docs fixes addressing comments by Jonas

e8269dc

Merge branch 'python' into docs_reaction_ensemble

444a418

Modified the order of reactants so that the reaction corresponds to h…

ed39290

…ow it is implemented

fweik added Python ApiChange labels Jan 12, 2018

kosovan added 5 commits January 12, 2018 12:18

Fixed syntax error

53f0b0e

Renamed RE equilibrium_constant to Gamma to avoid confusion with vari…

e92d6fb

…ous conventional definitions of equilibrium constant.

Merge branch 'python' of https://github.com/espressomd/espresso into …

3c72255

…docs_reaction_ensemble

Modified testsuite

8d71531

* Replaced equilibrium_constant with Gamma in testsuite. * Modified the scripts such that they do not set the standard pressure but still produce the same numbers as before

Removed standard_pressure and renamed equilibrium_constant with Gamma…

8c41ea6

… everywhere in the code.

jonaslandsgesell reviewed Jan 17, 2018

View reviewed changes

Small change in definition of Gamma, removed comments.

61b32f4

jonaslandsgesell approved these changes Jan 21, 2018

View reviewed changes

jonaslandsgesell reviewed Jan 22, 2018

View reviewed changes

kosovan added 5 commits February 1, 2018 13:34

new reaction ensemble test

4250140

Merge branch 'python' of git://github.com/espressomd/espresso into do…

24fedd2

…cs_reaction_ensemble

Resolved some problems from previous merge

03d1f43

Adopted constant_pH test to the same form as reaction_ensemble test

131b311

fixed a typo

6645e25

jonaslandsgesell reviewed Feb 1, 2018

View reviewed changes

Addressing comments from Jonas

0fb0287

kosovan and others added 6 commits February 3, 2018 21:09

Merge branch 'python' into docs_reaction_ensemble

eeae734

Removed one line forgotten during the merge

8ddceff

Dropped rel_error_alpha from msg and increased the number of steps. H…

ce290ae

…opefully no more failures will follow.

Changed number of samples

6c8325f

Cast all values in msg to str()

55d1e54

Corrected the accuracy

fca320a

jonaslandsgesell approved these changes Feb 6, 2018

View reviewed changes

fweik merged commit b5cf1f1 into espressomd:python Feb 6, 2018

ashreyaj pushed a commit to ashreyaj/espresso that referenced this pull request Aug 9, 2018

Merge pull request espressomd#1741 from kosovan/docs_reaction_ensemble

9f0e105

Corrected reaction ensemble documentation, small changes in RE syntax…

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Corrected reaction ensemble documentation, small changes in RE syntax… #1741

Corrected reaction ensemble documentation, small changes in RE syntax… #1741

kosovan commented Dec 28, 2017

codecov bot commented Dec 28, 2017 •

edited

Loading

jonaslandsgesell Dec 30, 2017

jonaslandsgesell Jan 2, 2018

jonaslandsgesell Jan 2, 2018 •

edited

Loading

kosovan Jan 11, 2018

jonaslandsgesell Dec 30, 2017

kosovan Jan 11, 2018

kosovan commented Jan 17, 2018

jonaslandsgesell commented Jan 17, 2018

jonaslandsgesell Jan 17, 2018

kosovan Jan 20, 2018

jonaslandsgesell Jan 20, 2018

kosovan Jan 20, 2018

jonaslandsgesell Jan 17, 2018

jonaslandsgesell Jan 17, 2018

jonaslandsgesell Jan 17, 2018

jonaslandsgesell Jan 17, 2018

jonaslandsgesell Jan 22, 2018

jonaslandsgesell Jan 22, 2018 •

edited

Loading

jonaslandsgesell Feb 1, 2018 •

edited

Loading

jonaslandsgesell Feb 1, 2018 •

edited

Loading

jonaslandsgesell Feb 1, 2018

jonaslandsgesell Feb 1, 2018

jonaslandsgesell Feb 1, 2018

jonaslandsgesell commented Feb 2, 2018

jonaslandsgesell left a comment

fweik commented Feb 6, 2018

Corrected reaction ensemble documentation, small changes in RE syntax… #1741

Corrected reaction ensemble documentation, small changes in RE syntax… #1741

Conversation

kosovan commented Dec 28, 2017

PR Checklist

codecov bot commented Dec 28, 2017 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jonaslandsgesell Jan 2, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kosovan commented Jan 17, 2018

jonaslandsgesell commented Jan 17, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jonaslandsgesell Jan 22, 2018 • edited Loading

Choose a reason for hiding this comment

jonaslandsgesell Feb 1, 2018 • edited Loading

Choose a reason for hiding this comment

jonaslandsgesell Feb 1, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jonaslandsgesell commented Feb 2, 2018

jonaslandsgesell left a comment

Choose a reason for hiding this comment

fweik commented Feb 6, 2018

codecov bot commented Dec 28, 2017 •

edited

Loading

jonaslandsgesell Jan 2, 2018 •

edited

Loading

jonaslandsgesell Jan 22, 2018 •

edited

Loading

jonaslandsgesell Feb 1, 2018 •

edited

Loading

jonaslandsgesell Feb 1, 2018 •

edited

Loading