Background
There are three types of partial charge interpolations in a hybrid topology transformation:
-
lambda_electrostatics_core - this value is linearly interpolated across the global lambda schedule. These core atom partial charges are varied using an offset of (new_charge - old_charge) * lambda_electrostatic_core.
-
lambda_electrostatics_insert - this value is 0.0 until 0.5, and then 2.0 * (global_lambda - 0.5). These new particle partial charges are varied using an offset of (charge) * lambda_electrostatic_insert.
-
lambda_electrostatics_delete - the value is 2.0 * global_lambda until 0.5 and then 1.0. These old particle partial charges are varied using an offset of (charge) * lambda_electrostatic_delete.
The issue
Because of the way our lambda schedule is set up, we effectively end up with unit charges at the end states, but non-zero charge amounts in the middle states.
Example: take the example of a Benzene (q_core = -0.13, q_delete = 0.13) to Toluene (q_core = -0.078, q_insert = 0.078).
The partial charge will vary in the following manner across our typical 11 lambda windows:
Lambda Total charge
0 0
0.1 -0.0208
0.2 -0.0416
0.3 -0.0624
0.4 -0.0832
0.5 -0.104
0.6 -0.0832
0.7 -0.0624
0.8 -0.0416
0.9 -0.0208
1.0 0
This also happens (and is even more pronounced!) when you have an alchemical ion for a charge transformation. For example, take the case of Benzene (q_core = -0.13, q_delete = 0.13) to Benzoic acid (q_core = -0.95, q_insert = -0.05). Here you will have a core transformation from oxygen to Chloride + two dissappearing hydrogens.
Lambda Total charge
0 0
0.1 -0.0914
0.2 -0.1828
0.3 -0.2742
0.4 -0.3656
0.5 -0.457
0.6 -0.3656
0.7 -0.2742
0.8 -0.1818
0.9 -0.0914
1.0 0
Solution -- when dealing with net charges
One thing we can do is add multiple particle parameter offsets on the ion to water transformation and make it so that the water->ion transformation varies by the same amount as the ligand charge transformation.
I have an incoming PR that should demonstrate how this could happen.
Solution -- when dealing with non-net charge transformations
Because the excess charge is smaller, the size of the effect isn't as impactful, but it is still something that should eventually get dealt with.
This one is a little bit harder. One option is that we stuff the change in charge on a water Oxygen (somewhere in bulk solvent) - the water would essentially go from a net charge of 0 to a non-zero net charge and then back to 0, which should yield a free energy change of zero.
We would do this both in solvent & in complex to and it should cancel out (not worked through the whole thermodynamic cycle yet).
That being said, there would be no way to deal with this in relative hydration free energies.
How critical is this bug?
In pratice, when dealing with net charges in PME, a neutralizing background correction is applied. Given that the thermodynamic cycle essentially cancels out errors, we're unlikely to see major changes in accuracy - especially when dealing with non-net charge transformations.
The impact on net charge transformations is likely to be higher, because of the magnitude of the charge excess, but the exact size of the error is unknown.
Could we just change the lambda schedule?
As far as I could come up with, there's no way to get this to work with our current formalism without causing a situation where there are naked charges.
Background
There are three types of partial charge interpolations in a hybrid topology transformation:
lambda_electrostatics_core - this value is linearly interpolated across the global lambda schedule. These core atom partial charges are varied using an offset of (new_charge - old_charge) * lambda_electrostatic_core.
lambda_electrostatics_insert - this value is 0.0 until 0.5, and then 2.0 * (global_lambda - 0.5). These new particle partial charges are varied using an offset of (charge) * lambda_electrostatic_insert.
lambda_electrostatics_delete - the value is 2.0 * global_lambda until 0.5 and then 1.0. These old particle partial charges are varied using an offset of (charge) * lambda_electrostatic_delete.
The issue
Because of the way our lambda schedule is set up, we effectively end up with unit charges at the end states, but non-zero charge amounts in the middle states.
Example: take the example of a Benzene (q_core = -0.13, q_delete = 0.13) to Toluene (q_core = -0.078, q_insert = 0.078).
The partial charge will vary in the following manner across our typical 11 lambda windows:
Lambda Total charge
0 0
0.1 -0.0208
0.2 -0.0416
0.3 -0.0624
0.4 -0.0832
0.5 -0.104
0.6 -0.0832
0.7 -0.0624
0.8 -0.0416
0.9 -0.0208
1.0 0
This also happens (and is even more pronounced!) when you have an alchemical ion for a charge transformation. For example, take the case of Benzene (q_core = -0.13, q_delete = 0.13) to Benzoic acid (q_core = -0.95, q_insert = -0.05). Here you will have a core transformation from oxygen to Chloride + two dissappearing hydrogens.
Lambda Total charge
0 0
0.1 -0.0914
0.2 -0.1828
0.3 -0.2742
0.4 -0.3656
0.5 -0.457
0.6 -0.3656
0.7 -0.2742
0.8 -0.1818
0.9 -0.0914
1.0 0
Solution -- when dealing with net charges
One thing we can do is add multiple particle parameter offsets on the ion to water transformation and make it so that the water->ion transformation varies by the same amount as the ligand charge transformation.
I have an incoming PR that should demonstrate how this could happen.
Solution -- when dealing with non-net charge transformations
Because the excess charge is smaller, the size of the effect isn't as impactful, but it is still something that should eventually get dealt with.
This one is a little bit harder. One option is that we stuff the change in charge on a water Oxygen (somewhere in bulk solvent) - the water would essentially go from a net charge of 0 to a non-zero net charge and then back to 0, which should yield a free energy change of zero.
We would do this both in solvent & in complex to and it should cancel out (not worked through the whole thermodynamic cycle yet).
That being said, there would be no way to deal with this in relative hydration free energies.
How critical is this bug?
In pratice, when dealing with net charges in PME, a neutralizing background correction is applied. Given that the thermodynamic cycle essentially cancels out errors, we're unlikely to see major changes in accuracy - especially when dealing with non-net charge transformations.
The impact on net charge transformations is likely to be higher, because of the magnitude of the charge excess, but the exact size of the error is unknown.
Could we just change the lambda schedule?
As far as I could come up with, there's no way to get this to work with our current formalism without causing a situation where there are naked charges.