Dear all
I want to calculate diffusion coefficient VS time in Gromacs (like figure below)to know at what time the system reaches to a stable diffusion state. What is the best way to do such analysis?
also I have another question regarding this figure. How does the diffusion coefficient calculated in different times? Is it like this: (0 ns to 1 ns), (1 ns to 2ns), (2 ns to 3n),... or like (0 ns to 1ns), (0 ns to 2ns), (0 ns to 3 ns),...?
Thank you very much
Dear Mohammad Rezayani
It depends on how you calculate the coefficient itself. One of the most common methods is using the mean square displacement (MSD) of the center of mass of the molecule/atom of your interest (see for example https://en.wikipedia.org/wiki/Mean_squared_displacement), however you should pay special attention to the condition of applicability of the formula you want to use to your particular system. For example the description of the diffusion of a water molecule in pure water is well described by Einstein's equation (i.e. the one in wikipedia page), but the diffusion of a water molecule across a lipid membrane is not.
If you want to use Einstein's formula, the best thing to do in my opinion is to extract from you data the coordinates of interest for the given molecule and then calculate the MSD from there. By using the relation in wiki's link, you then plot the MSD divided by (2*n*t) as a function of time, where n is the dimension of your MSD and t is the time itself. That is the value for the diffusion coefficient in that moment. You will see that at the beginning it fluctuates a lot and then it should reach a (reasonably) constant value. You can do all of this via GROMACS' tool gmx msd (http://manual.gromacs.org/archive/5.0.2/programs/gmx-msd.html) that calculates the mean square displacement of the chosen group and also does a linear fit. Read the documentation, it's a useful tool.
However, it depends on the formula you are using and the diffusion condition you are in. Einstein's formula works well if MSD is linear in time, that is if your molecule behaves as a n-D random walker. Take also into account that these formulas usually work as mean values, so it is better to have very long runs to collect a lot of data or (better) to average on a large ensemble of molecule. Probably starting with calculating the MSD is always a good option: if it is linear in time, then Einstein's formula is likely to be a good choice to calculate the diffusion.
About the plot, you should read the description in the paper. I think it depends on what they are studying.
Hope this helps,
Nicola
You need the coordinates (trajectory file) of the atomic groups of the molecules in your system at each time frame ( .xtc/ .trr/ .cpt/ .trj/ .gro/ .g96/ .pdb/ tng file), then use the GROMACS tool gmx msd. You have to choose a specific atomic group (http://manual.gromacs.org/archive/5.0.2/programs/gmx-msd.html). This will calculate the msd as a function of time, and also gives you the diffusion co-efficient by fitting the long-time tail with a linear equation (Einstein formula). MSD has three different dynamical regimes: ballistic (B), trapped (T), and diffusive (D). For further information see equation 1, figure 2 of Article Connection of large amplitude angular jump motions with temp...
Calculate mean square displacement from the simulation as a function of time. Then take the derivative of the data. D=1/6 lim(t->infinity) d MSD/dt
So plot 1/6 d MSD/dt , which is plotted in the above figure
Dear Debasish Das Mahanta
Thank you very much for your full description. I got it depply.
Thanks agin
Dear Kolattukudy P. Santo
Thank you very much for your reply. It got more clear for me.
Dear Nicola Piasentin
Thank you very much for your full description. It made me more clear to understand diffusion mechanism.
Thanks again
Dear Debasish Das Mahanta
I read your article but I got a little confused. As I understand diffusion coefficient is sum of rotational and transnational diffusion. Then diffusion coefficient obtained from NMR is sum of them? If it is true the D obtained from MSD could not be the real D obtained from NMR. Then to get a exact D I have to calculate both of them. is it right?
And the second question is that how you obtain rotational diffusion coefficient? Did you use a specific script or we can obtain it directly from MD software like Gromacs?
Finally, What is the exact reason for ballistic diffusion? change in volume during npt level?
Thank you in advance
Dear Mohammad Rezayani
1. Yes, the contributions of the real diffusion coefficient come from all types of motion involved. However, it is not a very easy task to deconvolute the translational and rotational motions experimentally. In the simulation study, you can do it with some careful analysis.
2. To calculate RMSD one has to map every single particle trajectory (XYZ coordinates) in the angular space as a function of time and then calculate the angular path spanned by any particle at a particular time interval. You can find the definition there in the paper (explained in details). It's very easy. I can write your own code easily.
- Badges
- Science topic