I am running an alchemical transformation simulation and calculation the absolute solvation free energy a ion pair. The ion pair consists of a metal cation and a nitrate anion. As I learnt from the previous post, there is no way to decouple two molecules at the same time. So I make an ion-pair itp file which contains two ions. But when I run the simulation, I see the nitrate is kind of glued to metal ion — the nitrate plane surface is patched to the metal ion, which does not make sense physically. As I turned off the free energy decouple and run the regular simulation, such unphysical behavior disappeared.
My current speculation is that, when the decoupling proceeds, the intramolecular interaction is not handled properly as these two ions are defined as one molecule. One work-around is to decouple them separately, but that will obviously increase the simulation cost. I am wondering if anyone has come to similar issue and have a better resolution for this.
One solution is to do manually write the topology for decoupling the ions. Then you need to define a dummy atom type with LJ parameters set to zero and use this in either the A or B state for the atom parameters in your molecule types, along with setting the charge to zero in that same state.
But putting both ions in the same moleculetype should also work. I don’t understand why that doesn’t work for you. Did you maybe define a bond between the ions?
I do not have any bond between these two ions. I am attaching my .top file for the ion pair.
Another thing came into my notion is that the ion pair carrying a +1 net charge. Although the problem seems to be short range interactions between two ions, can the long range PPPM calculation cause this trouble?
Theoretically, the alchemical transformation should not be affected by how I phase in and out ions. (Now I am phasing out +1 ion-pair and appearing in +1 cation separately.) But practically, will this cause trouble? I guess a more thorough way to do this to define them as A and B states to do the decoupling together.
I recall now that there actually can be issues depending on how you phase out ions. For example, ions which charge +0.5 and -0.5 like each other a lot, but don’t really like water, so they will tend to cluster together.
But you probably want to keep the net charge of the system zero and therefore simultaneously change the two ions.
Can you please elaborate more on this issue? Why are charge +0.5 and -0.5 not dispersed by water molecules? I also carries out a similar alchemical transformation in other software package. I did not see such clustering behavior there. There I used slow-growth, but I think that should not be the case because I did not see any incident like this in the trajectory.
Hydrogen bonds between water molecules are stronger than the interaction with an ion with 0.5 charge. Thus the two ions will stick togther, if they manage to find each other, when the charge is lower than something like 0.8.
What happens precisely depends on how you interpolate the Hamiltonians and on if you use soft-core interactions and how they look like.