Hello all,
I hope you all are doing great.
Recently, I have focused on the bulk properties of water using a flexible water model (TIP4P/2005f). To do so, I read the original article for this water model and noticed some differences and similarities such as O-M distance is the same as TIP4P/2005, full morse function is used for TIP4P/2005f, the angle coefficient and the teta is also different.
I also figured out I should use a smaller time step for flexible models, so I used 0.25 fs. I performed simulations in NPT (T=298.15K and P=1 atm) for 8000000 time steps (2 ns). Besides, the number of water molecules was 521. I have attached the input files to this topic.
To check if my model is correct, I computed a couple of properties, including dipole moment, density, temperature, pressure, average r-OH, and average teta.
However, what I am observing is not close to what I expect. Although temperature was maintained at 298.15 K, pressure was around -17 atm, density was 0.993 g/cm^3. Furthermore, the dipole moment was 17.9 (C.A) which was off the one reported in the main article. I also calculated and from dumped trajectories, although I got the average r-OH as 0.96 A which was close to the reported value, the average teta was 102.603, slightly off.
I have cross checked my parameters and commands I used, but I did not notice any mistake. I am thinking I am missing something very serious.
I was wondering if you could give me a hint.
Cheers,
Al.
data.sim (133.3 KB)
in.sim (1.3 KB)
Hi a quick update on this.
I made a silly mistake, which was damping.
the one that I chose was too small, that is why I got that pressure.
Also, for dipole, it should be done using dipole/tip4p/chunk
We all make straightforward errors from time to time.
Thank you for documenting what you did – this will help future users very much.
Could you put up a reply with what the old damping constant was and what pressure you got, and what the new damping constant was and the new pressure? This will not just help other users – in the (very unlikely) event that you have not kept your own notes, this will help you too!
Hi,
Yes, sometimes I get lost in my scripts that is why I may make an obvious mistake.
In the old setup, I set timestep 0.25. Since I am using real units, then it is 0.25 fs. Well, damping is technically determining at what time the coupling should be done. So for example, if someone chooses 1000 (step) in the fix npt command or any fix command that requires damping, it is indeed t = timestep * step, in my case it became 0.025 ps (0.25 fs * 100 * 0.001) for temperature coupling and 0.25 ps (0.25 fs * 1000 * 0.001) for pressure coupling.
MD is about fluctuations and averages and other fun stuff in statistical thermodynamics. Such small time coupling leads to huge fluctuations. With this setting, as I mentioned in my first post, I got around -17 atm. This was basically telling me the box wanted to shrink more which was unreal.
And also, equilibration was completely done because the system was quite small and it was run for 2ns.
Now, if we choose a 4000 step in damping, we would have 1 ps. Ideally, 1 ps is good enough for both temperature and pressure. My experience with GROMAC says temperature can be coupled faster. So even if someone chooses 2000 step (0.5 ps) for temperature, and 4000 step (1 ps) for pressure could be the best choice.
Note, this damping is working for water in the case of my system in the liquid state. Depending on systems, one needs to use the right damping. For example, for n-alkanes, I would use somewhere between 5 ps and 10 ps.
Finally, when I used 4000 step (1 ps) with a timestep of 0.25, after getting block averages and discarding 1000 ps from outpute, I got the pressure around -0.8 atm. If did more sampling, I would be around 1 atm at the end.