diffusion coefficient for liquid water

[姚懿]

Dear all,

I would like to calculate the quantum self diffusion coefficient for liquid water. My understanding is I need to do cmd and calculate the vvacf or msd from the centroid path.

Is there an example of doing it? I have seen rpmd and cmd in the manual but cannot find examples for them.

Thank you in advance,

Yi

[Michele Ceriotti]

I recommend TRPMD, as it typically gives the same results as (PA)-CMD for zero-frequency properties, and does not require reducing the time step.
It's such a small change with respect to a normal PIMD run that we never bothered add an example but since you ask, I have just pushed an h2o-trpmd example within examples/lammps/ folder
Hope this helps
M

[姚懿]

Hi Prof. Ceriotti,

Thank you for your fast reply.

It is very helpful for me.

Since, it is the first time I heard about TRPMD. I would like to make sure I have the correct reference for it. 

Is this the correct reference for TRPMD? 

====

http://aip.scitation.org/doi/abs/10.1063/1.4883861

How to remove the spurious resonances from ring polymer molecular dynamics

====

Best,

Yi

[Mariana Rossi]

Dear Yi,

That is the correct reference, yes. As you see in that paper, the self-diffusion coefficients of liquid water are essentially the same with CMD, TRPMD and RPMD.

All the best,

Mariana

[Chris]

Hello all,

 

I am also trying to calculate the diffusion coefficient of water, and rather than create a new thread, I decided to continue this one. I’ve run the TRPMD simulations recommended above, but I am having trouble getting realistic results for the diffusion coefficient. I know I need the Kubo-transformed velocity autocorrelation function (as seen in Eq. 43 in https://aip.scitation.org/doi/10.1063/1.2953308 ), so I tried using the getacf.py python script included in the i-PI release, but the results I'm getting do not get anywhere close to the expected results. My results for the radial distribution functions are in good agreement with published results, so I do not think it is an issue with how I'm running my simulation (though that is a possibility). My question is: does i-PI do the necessary Kubo-transform of the velocity autocorrelation function in the getacf python script (I tried looking through the source code but couldn’t determine whether or not this is the case)? If it does not, what steps do I need to take to perform this transform?

 

Thanks,

 

Chris

[Michele Ceriotti]

Venkat can tell you more about the getacf script, that indeed should get you what you are looking for if you run it over the trajectory of the centroids. 

However you should understand that RPMD computes the Kubo transform of the quantum CF, so you do not need to do any transformation to the 

CF you extract from the trajectory.  See here http://dx.doi.org/10.1063/1.4883861 where the explanation is as concise as it gets.

Hope this helps

[Mariana Rossi]

Dear Chris,

I also recommend extracting the mean square displacement instead, as explained for example here: https://pubs.acs.org/doi/abs/10.1021/acs.jpclett.6b01093 . 

It is less prone to noise than the acf for the diffusion coefficient. Especially if the tail of your acf is noisy it is almost impossible to get a good diffusion coefficient.  

All the best,

Mariana

[venkat kapil]

Hello Chris, the RPMD estimator for the velocity velocity autocorrelation is itself an approximation to the Kubo transformed correlation function of the velocity velocity autocorrelation. So there is no need to perform the Kubo transform explicitly. Coming to the script, I think there might be a problem with the units.

I have attached the input.xml that I used. I use a timestep of 0.25 fs and a total of 4 ps for molecular dynamics and a stride of 4 to print out the velocities.  And I used the following to get a quick vvac spectrum. The important thing is how dt is specified. It needs to be the time step at which the velocities are being printed out. Also the final output is in atomic units. Which means that lim_{omega -> 0} vvac is in the units of bohr^2 / atomic_time. To convert to Angstrom^2/picosecond you need to multiply by 11576.76. With that I get  0.238 Angstrom^2/picosecond which is in the right ball park. Of course to get a converged result one needs to average over many short simulations.

python ~/source/i-pi-dev/tools/py/getacf.py -ifile simulation.vc.xyz -mlag 2000 -bsize 4000 -ftpad 4000 -ftwin none -dt "1.0 femtosecond"  -oprefix out

[Chris]

Thank you all for your prompt and informative replies.

 

Michele, thank you for clearing that up.

 

Mariana, from what I can tell, i-PI does not include any functionality for computing the mean square displacement. This isn’t too much of an issue, as many MD codes are capable of doing so on their own. I’ll look into how to compute it using LAMMPS (the client MD code I am using).

 

Venkat, I took the files you provided, ran the i-PI simulation, ran the same python command you provided to get the vvacf, loaded the file into a Python script, performed the integration & divided by 3 (per the formula from the reference). I got a value of 1.49088e-4 [bohr^2 / atomic time], which translates to 1.725 Angstrom^2/ps after multiplying by the unit-conversion prefactor you provided. This is significantly larger than the value you got, so I suspect that there is something that I am neglecting to do. My run did use a different random seed from the file you provided, but I don’t think that would lead to a difference of this magnitude.

I’ve attached a .png file containing the plot of the vvac as well as the Fourier transform of the vvac. As far as I can tell, the Fourier transformed vvac looks correct, so I suspect my issue is not with the i-PI run or the getacf.py command, but rather in how I calculate the diffusion coefficient. The exact Python command I use is:

 

D=numpy.trapz(vvacf, t)/3*11576.76

 

Is this not the correct way to calculate this diffusion coefficient?

 

Thanks,

 

Chris

[Mariana Rossi]

Dear Chris,

If you save the trajectory of your simulation, calculating the MSD is a rather straightforward post-processing. I will try to find some scripts here but it is not dificult to write one's own.

That is how we calculated the water self-diffusion coefficient for TRPMD and CMD that we report in Table I of https://aip.scitation.org/doi/10.1063/1.4883861 . 

There is certainly some problem in the units you are taking for the integration.

All the best,

Mariana

[Chris]

Hi Mariana,

I was able to write up a small script to calculate the MSD from the i-PI trajectories, and after averaging together 10 independent 50 picosecond TRPMD runs of q-TIP4P water, I got D=0.2105 Ang^2/ps, which is in good agreement with the results published in the paper you mentioned. 

I was curious about why I wasn’t able to get the correct answer by integrating the velocity autocorrelation function, so I did some digging. I found there’s a relation for what the vvac should converge to (from “Quantum diffusion in liquid para-hydrogen from ring-polymer molecular dynamics”):

Lim C_vv(t)  = 3/(beta * m)

t->0

doing this for water (where beta = 1/(k_B*T) = 1/(1.380e-23 J/K * 300 K) = 2.415e20 J^-1, and m = 2.9915e-26 kg), I get that

Lim C_vv(t)  = 3/(2.415e20 J^-1* 2.9915e-26 kg) = 4.15225e5 m^2/s^2

t->0

Looking at the vvac file provided by i-PI, I get that C_vv(0) = 7.04074e-06 which, after converting from atomic units, gives C_vv(0) = 3.36969e7 m^2/s^2, which is ~81 times larger than the theoretical value that I calculated. I know that if the tail of the vvacf is particularly noisy it can skew the results for the diffusion coefficient, but maybe this is also related somehow? 

Thanks,

Chris

[venkat kapil]

Dear Chris, thanks for your email. The diffusion coefficient is estimated as a the integral from t=0 infinity of the vacf which means it is the value of the fourier transform (facf) at t=0. The facf has already been divided by 3 so that the value at t=0 gives the diffusion coefficient directly in atomic units which has to be multiplied by  11576.76.

As for calculating it from the vacf there are two things:

1) The pre-factor is not correct since we have already renormalized the facf to give the diffusion coefficient.
2) Even if it were correct, your results would not be correct since you are missing a dt = 1 fs = 0.024188843 a.u.  in 

D = numpy.trapz(vvacf, t)/3*11576.76

[Eslam Ibrahim]

Dear Venkat,

Thank you very much for the useful insights. It is very helpful for what you provided to calculate the quantum diffusion coefficient, but I have a question regarding how to choose the optimal value for the  maximum time lag.

Thank you very much.

Best regards,

Eslam

[Venkat Kapil]

Hey Eslam, you want to keep the maximum time lag as high as possible. However, a longer maxlag can lead to noise in velocity autocorrelation, so you have to check the error bars on the estimated diffusion coefficient.