From what I'm reading, simulation of auto-ionization in liquid substances is hard to do using gas phase DFT or AIMD [1][2][3]. I am using ORCA to do my work and I must confess that I'm shocked to see that when using modern functionals like wB97X-D3 and r2SCAN-D3, and even when going with a triple zeta basis set I wasn't able to see a single auto-ionization at 300 K in a box containing 50 water molecules or 30 DMSO molecules after 2000 1-fs steps or more.

My silly question is: is this calculation possible without using more specialized software like CP2K or am I doing something terribly wrong?

  1. Y. Fu, L. Liu, R.-Q. Li, R. Liu and Q.-X. Guo, J. Am. Chem. Soc., 2004, 126, 814–822.
  2. L. R. Pestana, N. Mardirossian, M. Head-Gordon and T. Head-Gordon, Chem. Sci., 2017, 8, 3554–3565.
  3. P. L. Geissler, C. Dellago, D. Chandler, J. Hutter and M. Parrinello, Science, 2001, 291, 2121–2124.
  • $\begingroup$ You use ORCA for the electronic part, how are you modeling the dynamics? I'm no expert on AIMD, but from a brief glance at the 2nd paper you cite, they do 5 ps (5000 fs) of just equilibration before they proceed to their production runs of 40 ps. So you may just not be running long enough simulations. $\endgroup$
    – Tyberius
    May 3, 2022 at 2:53
  • $\begingroup$ Auto-ionization will be a rare event. You should use an enhanced sampling technique like meta dynamics with one or more collective variables that help to force the proton transfer. You can compute free energy barriers from this approach, if that’s what you’re after. $\endgroup$
    – Stephen
    Mar 9 at 16:19

1 Answer 1


This is because autoionization is very thermodynamically unfavored near room temperature. I don't know much about water in a water-DMSO mixture, but let's take pure water as an example. The equilibrium constant of the autoionization $$ 2\ce{H2O} \to \ce{H3O+} + \ce{OH-} $$ is well-known to be $3.2\times 10^{-18}$ (this differs from the more familiar figure $K_{\rm{w}} = [\ce{H3O+}][\ce{OH-}] = 1.0\times 10^{-14} \rm{M}^2$ by the square of the bulk concentration of liquid water, $55.6 \rm{M}$). The rate constant of the reverse reaction is $1.3\times 10^{11} \rm{M}^{-1} \rm{s}^{-1}$. Thus, from detailed balance, the rate constant of the forward reaction is $4.1\times 10^{-7} \rm{M}^{-1} \rm{s}^{-1}$. Taking into account the bulk concentration of water, this means that one must wait about 12 hours for any given water molecule to be engaged in an autoionization event. Since you have 50 water molecules, the required time is shortened to about 14 minutes, but this is still orders of magnitude longer than the 2 ps simulation that you have done! Obviously no software package, not even CP2K, can do a 14 minute AIMD simulation in the foreseeable future.

The lesson is that, a (perhaps surprisingly) large fraction of reactions cannot be studied by plain AIMD (or by even classical MD) due to the timescale involved, but rather necessitate enhanced sampling techniques. People tend to be interested in reactions that can most conveniently be studied in a wet lab, namely those reactions that take seconds to hours to finish under room temperature, but the timescale of AIMD simulations is over 10 orders of magnitude shorter, such that it does not overlap with the timescale of these reactions at all.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .