Issue tracker for Tutorial 12

[jngrad]

Regarding the version with WCA active and P3M active:

• it takes 1 hour to run the tutorial
• the tutorial only says "Observe how the ionization equilibrium is affected by various interactions" without interpreting the plot
  • the sigmoid is shifted to higher pH values, i.e. the degree of dissociation is smaller at pH > pKa, probably because once several monomers are basic, protons lost by neighboring monomers have a higher chance of being recaptured by basic monomers due to the Coulomb interaction
  • it's difficult to assess if the resulting plots are correct; are analytical formula available?
• the test for the WCA+P3M case doesn't check numerical results

[jonaslandsgesell]

If there is low electrostatic interactions present, then the ionization of one group suppresses the ionization of other groups (because of the like-charge repulsion of charged monomers along the chain). Therefore at same pH (for a given acid which has a certain pK), the ionization is lower (if you turn on electrostatics). One can compute this when using the thermodynamic framework, starting from dF=sum_i nu_i mu_i=0 (for constant volume and temperature) and inserting the chemical potential. For the excess part of the chemical potential you need a model, e.g. the Debye-Hückel model. For more about the chemical equilibrium condition or the partition functuons of the reaction ensemble or the constant pH ensemble, see https://pubs.rsc.org/no/content/articlehtml/2019/sm/c8sm02085j as well as https://link.springer.com/article/10.1140/epjst/e2016-60324-3 for understanding the difference between the constant pH method and the reaction ensemble method (in short: the constant pH method does not take care of ionic strength changing with pH, i.e. the constant pH method uses implicit protons which results in wrong electrostatic interactions at extreme pH).

Note: if there is strong electrostatic interactions (i.e. low epsilon_r) which induce strong ionic correlations then the dissociation can also be enhanced at same pH and pK, so the dissociation increases (if you turn on electrostatics).

You really should separate between the pH (a measure for the chemical potential of H+) and pK (a property of the acidic/basic monomers of the polymer). Plotting stuff as a function of pH-pK suggests that this was some unique order parameter and it would not matter whether you change pH or pK. But his is only true for ideal systems! For systems with interactions changing pK or changing pH is very different. Changing pH changes e.g. the electrostatic screening in the system very different to changing pK. This should be stated clearly. Only in intermediate pH ranges changing pH or pK is approximately behaving similar. See https://link.springer.com/article/10.1140/epjst/e2016-60324-3

If you have a model for how the chemical potentials change with pH you can predict (from dF=0=sum_i nu_i mu_i) how alpha changes with pH at least for a monomeric system. But already here there is (to my knowledge) mostly empirical formulas available for this because computing the chemical potential for electrostatically interacting species would require to compute mu_i=partial F/partial N which is not possible because F=-kTln(Z) cannot be computed analytically because the partition function Z cannot be computed analytically. All this gets even more intractable for polymeric systems where bonds need to be taken into account in the excess part of the partition function.

I hope this answers the first part of your question.

[kosovan]

I think that an earlier version of the tutorial was running faster and I agree that one-hour runtime is too long. My suggestion is to decrease the number of pH values for which we calculate the result with electrostatics and/or decrease the accuracy of electrostatics, or replace it with DH. However, it should be explained in the tutorial that we sacrificed accuracy for speed, and that this should not be done in production runs.

Interpreting the plot is indeed desired. I plan to work on fixing some other issues in this tutorial, which showed up when I was teaching it. Hopefully in November 2020, when I will be going through the tutorials again with my new PhD students.

Testing can be done by generating a fixed set of numerical results as a reference. Jonas is right that there is no theory which could predict the correct result.

[kosovan]

I propose to address the following issues in the constant-pH tutorial ( I plan to do it during the coding day):

• Title of the notebook "The Reacton-Ensemble Method" does not correspond to the filename and to the content
• Specify the expected prior knowledge:
  • Chemistry, chemical reaction equilibria
  • Espresso - Electrostatics, LJ tutorials
• Refernce cocnentration with ominus or with zero
• Unit of energy should not be denoted at Delta E or Delta x
• Reduced unit of energy of epsilon, unit of length is sigma, refer to ug for reduced units
• BOX_L_UNITLESS and C_ACID_UNITLESS is very confusing, should be BOX_L.magnitude
• Similarly, "to("sim_length * sim_energy / sim_charge^2")" and "kT=KT.to("sim_energy").magnitude" might be confusing but it is better
• Test how to speed-up the calculations with electrostatics:
  • Lower accuracy of P3M
  • Use DH
  • Explain that poor-quality electrostatics is not suitable for production simulations
  • Provide a pre-computed titration curve with electrostatics, let the user compute just 2-3 additional data points and then plot them together with the pre-computed data
• Explain PROB_REACTION - it is not "probability of accepting the reaction"
• The nested cycles setting up WCA interaction might be difficult to understand: "for type_1, type_2 in ((x, y) for x in TYPES.values() for y in TYPES.values())"
• User experience: Setting up the electrostatics was easier to follow in the earlier version of the tutorial (?)
• Explain why we use nsquare cell system
• RE.temperature shoudl be renamed to KT, same as in the thermostat
• Avoid the term 'simulation units' - some statements are ambiguous: "make sure to provide the temperature and exclusion_radius in simulation units!"
• Explain that no other values of coefficients are possible; in future, remove this arguments from the constant-pH
• Rename equilibrate_pH to equilibrate_reaction
• Functions equlibrate_pH and perform_sampling should not depend on global variables
• Modify perform_sampling - if there is no WCA, then nothing happens in 50% of the moves
  • the function should not depend on global variables PROB_REACTION, N_ACID, TYPES, ...
  • PROB_REACTION should be 1.0 if no MD integration is performed
  • consider inner loop, so that we do several reaction moves and several configuration moves within one "sample"
• Desired_block_size should be the desired number of independent samples per block, not the number of data points
• N_BLOCKS is not used in the function
• Setting of parameters that we want to vary should be done within one single cell.
• All relevant simulation parameters should be set within one block of the code, outside the function calls that pass the parameter values to espressomd.system()

[kosovan]

I will try to finish the changes to the constant-pH tutorial during the coding day today.