This is an old revision of the document!
Diffusion constant, viscosity and size effects
When simulating liquids or solids under periodic boundary conditions, we are making two fundamental approximations:
- We simulate an infinite system, thus neglecting the fact that any real-world system has surfaces. This approximation becomes problematic, when the real-world system to be studied consists only of a few simulation cells.
- We impose the condition that the properties of the system under study repeat exactly from one simulation cell to the next. The quality of this approximation depends on the system under study and the quantity of interest.
Here, we want to calculate the diffusion constant of water at room temperature ($T=300\,\text{K}$). Since we are interested in a property of bulk water, we don't need to worry about the first approximation. But we need to pay attention to the second one. The theory derived in 10.1021/jp0477147 allows us to estimate the finite size effects in the diffusion constant $D_{pbc}(L)$, calculated under periodic boundary conditions with cell size $L$. With this information at hand, we will be able to extrapolate the results for finite cell sizes to the diffusion constant $D=\lim\limits_{L\rightarrow \infty}D_{pbc}(L)$, effectively getting rid of the second approximation.
Calculating transport properties typically requires lots of sampling. Start the MD simulation for 32 water molecules and see how far you can get (aim at least for 200ps).
- While the job is running, check output of CP2K to verify that all is fine. What is the average temperature?
- We want to simulate diffusion at room temperature. Why aren't we using the NVT ensemble? Hint: Think about how thermostats work.
- Use the provided script
./get_t_sigma abc.enerto calculate the standard deviation of the temperature for your simulation as well as for the provided simulations of larger cells containing 64, 128 and 256 water molecules. - How are temperature fluctuations expected to depend on system size? Use gnuplot's fitting functionality to check whether they follow the corresponding law.
The mean squared displacement (msd) $$\text{msd}(t) = \langle |r(t)-r(0)|^2 \rangle$$ is related to the diffusion constant (see eq. 17 in the article).
We have provided a Fortran program that calculates the msd from a trajectory in a .xyz file in units of $\unicode{x212B}^2$ as a function of time in fs.
gfortran msd.f90 -o msd.x # compile msd.x executable ./msd.x < msd.in > msd-64.out # run it with msd.in input file
- While your simulation is running, calculate the msd for the provided simulations of 64, 128 and 256 water molecules, modifying
msd.inas needed. Note:msd.xmay run up to 30min for the largest cell.
do we have a progress indicator? - Plot the msd as a function of time using a double logarithmic scale. Can you identify different regimes? Note: With our method, the number of samples collected for $\text{msd}(t_0)$ decreases like $1/t_0$.
- Obtain the diffusion constant $D_{pbc}$ by fitting a line through the mean square displacement data in the range $2-10$ ps.
- Compare against the values in Table I. Note: We are using a slightly different force field, but the values should be similar. If not, check your units!
On the next day, when your MD simulation has finished.
- Calculate also $D_{PBC}(L)$ for the 32 water molecules.
- Plot $D_{PBC}$ as a function of $1/L$, where $L$ is the length of the edge of the simulation box.
- Perform a linear fit of this curve to obtain the diffusion constant $D=D_{pbc}(L=\infty)$
- Use Eq. 12 in the article to calculate the viscosity.
- Compare the result to data in the paper.