Dear LAMMPS team,
I am investigating water structures on metal surfaces using LAMMPS and wanted to compute the msd of only one oxygen atom (not the whole oxygen group) over the simulation.
I first tried to compute msd for a group of oxygen atoms using :
compute msd ox msd com yes , where ox is the group of oxygen atoms
LAMMPS will output me the msd. However, when I only considered one oxygen atom the simulation using:
compute msd oxy_1 msd com yes , where oxy_1 is a specific atom , it will give me a MSD of 0 the whole simulation duration.
It really didn’t make sense to me. So I tried to use msd without “com yes”:
compute msd oxy_1 msd com no
and I get non-zero results.
Why can’t I use com yes to calculated the msd of one oxygen atom, whereas with the whole group of oxygen atoms I get some results.
Sorry if this problem is too trivial for you, I would be grateful for any tips and answers 
Kind regards,
Oskar
You should not be needing the “com yes” option in either case. If your system has a drift, you need to address what is causing it. The purpose of this option is mainly to determine self-diffusion in flow simulations.
Of course, the displacement of a single atom is equal to the displacement of its center of mass and thus your zero msd is expected output.
Axel
1 Like
Thank you so much for the fast response.
One last question I still have to that issue: I used the SPC water model to run my simulation.
Is the fix shake command necessary when using water models on metal surfaces?
I understand that the bonds and angles will be constrained every timestep, however I get quite different results for msd for both fix shake and without fix shake (see attachements)
Again, thank you in advance.
Kind regards,
Oskar
with_fix_shake_group_oxygen_msd_time_10ns_15_H2O_with_fix_shake.eps (38.7 KB)
no_fix_shake_oxygen_group_msd_time_10ns_15_H2O.eps (37.8 KB)
Thank you so much for the fast response.
One last question I still have to that issue: I used the SPC water model to run my simulation.
Is the fix shake command necessary when using water models on metal surfaces?
the (standard) SPC water model is a rigid model and thus fix shake is required. However, there are many water models around, rigid and flexible. of course each of those will produce different results. similarly, the way how you model the metal and particularly how you describe the water/metal interactions have a crucially important impact on the results. modeling water/metal interfaces is not a new topic. There were already plenty of water models and different ways to model these interfaces (and a wealth of literature on the subject) over 25 years ago, when I learned MD simulations as an undergrad. There must be a plethora now. I suggest you do some research on what is the most suitable and accurate for your needs.
I understand that the bonds and angles will be constrained every timestep, however I get quite different results for msd for both fix shake and without fix shake (see attachements)
you are supposed to get different results as mentioned above, however, the differences you see in the graphs are obviously not really because of the model but because of (very) bad statistics. to get converged MSD results for a single atom you need to do many (non correlated!) runs and average over them. compute msd is typically used for bulk systems where the MSD is computed for each atom in the group and then averaged over those atoms. and even there you may need to average over multiple runs depending on the sampling of states you get.
Axel.
Thank you sir.
It helped me a lot.
I reran some simulation of water layers on metal surfaces with the tip3p water model and try to understand the “compute msd com yes” issue about the drifting of the group of oxygen atoms (from water) on metal surfaces. As you mentioned earlier, in the file tip3p_group_oxygen_msd_time_10ns_15_H2O.eps (msd for group of oxygen without com yes) we see a very random “trend”, which according to my understanding makes sense, because we only ran one simulation and is bad statistics.
However if you look at drift_tip3p_group_oxygen_msd_time_10ns_15_H2O.eps (msd for group of oxygen with com yes) we see a “normal” msd vs time curve. Does that mean (with regard to com yes setting) that after a certain point of time the water layer moves as a whole unit? (because we are subtracting the drift from the centre of mass). Can I somehow use this result even though we only ran one simulation and only changed the setting to com yes?
Thank you again.
Kind regards,
Oskar Cheong
drift_tip3p_group_oxygen_msd_time_10ns_15_H2O.eps (40.5 KB)
tip3p_group_oxygen_msd_time_10ns_15_H2O.eps (40.2 KB)
Thank you sir.
It helped me a lot.
I reran some simulation of water layers on metal surfaces with the tip3p water model and try to understand the “compute msd com yes” issue about the drifting of the group of oxygen atoms (from water) on metal surfaces. As you mentioned earlier, in the file tip3p_group_oxygen_msd_time_10ns_15_H2O.eps (msd for group of oxygen without com yes) we see a very random “trend”, which according to my understanding makes sense, because we only ran one simulation and is bad statistics.
However if you look at drift_tip3p_group_oxygen_msd_time_10ns_15_H2O.eps (msd for group of oxygen with com yes) we see a “normal” msd vs time curve. Does that mean (with regard to com yes setting) that after a certain point of time the water layer moves as a whole unit? (because we are subtracting the drift from the centre of mass). Can I somehow use this result even though we only ran one simulation and only changed the setting to com yes?
Thank you again.
Kind regards,
Oskar Cheong
drift_tip3p_group_oxygen_msd_time_10ns_15_H2O.eps (40.5 KB)
tip3p_group_oxygen_msd_time_10ns_15_H2O.eps (40.2 KB)
What you are asking about is not really a LAMMPS issue by a question of what it is that you want to measure. I have already given an explanation how I see that option and how and where it should be used. That is ultimately something you need to discuss with your adviser not me, since this concerns your research and not LAMMPS.
A few more general remarks.
If a graph looks like expected, it does not automatically mean it is correct. It can still be very wrong. Neither of your attached graphs looks good and well converged for a MSD graph. They are far too noisy and inconsistent. And there are very well established protocols to determine the statistical relevance and convergence of data.
If you have a single and very long trajectory, it usually possible to analyze it in chunks and then compare the result for the total with two results for halves, four results for quarters and so on until the individual sections become too noisy. for “short” runs you usually need to average over multiple decorrelated runs and then the convergence just follows the usually sqrt()-like behavior for noisy data.
Axel.