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.
I stumbled into this thread as I was trying to run a simulation for water using this same force field and also getting a density of 0.993 g/cm^3 for liquid water at 300 K, 1 atm (which is not the one they get in the official publication of the TIP4P/2005f force field: 0.9977 g/cm3 - see http://dx.doi.org/10.1063/1.3663219). For what is worth it, I was not really computing the other properties, such as dipole moment.
Force field information - i.e., bonded and non-bonded potential parameters (charges included), cutoff for the non-bonded interactions, position set up for the virtual site in the tip4p pair_style, pair_modify used - input in my script is literally the same as in @alireza’s simulation, apart from the fact that I only used numbers up to the 3rd decimal case (I doubt any significant difference can come from that).
I did not have the problem @alireza had with the damping constants for the barostat and thermostat (i.e. they are duly defined to the values recommended in the manual for a 0.25 fs timestep). So I did two things:
(i) reran his simulation using the scripts he gave here but fixing the damping constants to the values adviced in the manual, since this was said to be a problem.
(ii) reran his script with my configuration instead as I thought maybe the equilibration I set up was not good (long story short: basically there is just too much space between the molecules to start with due to the way I create the initial configuration and so I initially do some “fix deform” in my own script to shrink the simulation domain. It does not shrink the box up to the exact value of the desired volume as I setup a value slightly smaller - I was relying on a fix npt summoneed later to do reach the equilibrium volume at ambient conditions. Since I thought all of this could be a mess, I tried using my LAMMPS data file with his script).
I got the same value of 0.993 g/cm3 in both cases…. Note that from the official publication, the cutoff for the short range non-bonded interaction seems to be 8 angstroms. I can already tell you that this does not fix the problem either (in fact if you use 8 angstroms you get an ever lower density…..).
I had been using the 29Aug2024 LAMMPS version. Since this post is from the beggining of the year, I guess its not such an outdated version to use in this context and compare. Does anyone have any idea on what the problem might be? Could it be that @alireza was using an older version and somehow some bug related to the tip4p style’ish exists in the new LAMMPS versions?
As a last piece of information, I can tell you that I have the same problem of an underestimated density when using the TIP4P/Ice force field (in which water is rigid), so probably its not a problem related to the molecule being flexible.
I did not try this myself, but I think this kind of small difference is well expected. The original paper is published in 2011, when computational capabilities are way lower than nowadays (so more likely suffer from finite size effects, statistical errors, etc.). Also they are using a different software (which has different default settings in many aspects), with different thermostats and barostats. They are all potential sources of discrepancy.
I think you should make sure that your simulation setup is correct and then just go on. If you want to reproduce the exact numbers then I think the best way is to use GROMACS with all setup same as those in the original paper.
Also,
apart from the fact that I only used numbers up to the 3rd decimal case
If your input only have 3 decimals, then your simulation result is only accurate to that, right? (well, this is not always the case, but if you want to do exact comparison then that may be a problem.)
Some more thoughts from looking at the attached input deck
- when running with fix npt, a Kspace convergence of 1.0e-4 is insufficient in my experience. It results in reasonably well converged forces, but not a well converged virial. I suggest to try 1.0e-6 and see if it makes a difference.
- similarly, the /long styles use tabulation by default for the coulomb forces. You can try
pair_modify table 0to use a slower but more accurate approach and see if it makes a difference - when comparing to other codes, one should check if they do use the tail correction for LJ interactions and adjust
pair_modify tail yes/noaccordingly. - it should not matter for liquids, but since water is strongly hydrogen bonded, there may be a difference between isotropic and anisotropic box adjustments. It should not.
Hi, thanks for reproducing density results for TIP4P/2005f.
Before I dive into what lammps does, let’s get an overview of what TIP4P/2005f is. Unlike rigid bodies, this model is flexible even in the case of dummy atom. From the original paper, you can see that the position of the dummy is dependent on the positions of hydrogen atoms. However, at the equilibrium, it reaches alsmost the same as rigid body in terms of Oxygen - Dummy distance ( as far as I remember). Some say a good approximation is to make it fix but I say wrong answers are mostly free of cost.
Now, let’s go for lammps. LAMMPS does not support flexible TIP4P styles, period, to the best of my knowledge. One way to overcome this challenge is to modify the source code for TIP4P style. Essentially, in this flexible model what makes things a little bit complicated is the calculation of dummy position, because that is very important when it comes to coulomb interactions.
Fo future references, if anyone wants to do flexible models, one has to take care of time steps and masses of hydrogen (or dummy if anyone could re-write the source code for TIP4P flexible models). For flexible models, especially TIP4P/2005f, less than 0.5 fs should be used, while for rigid bodies 1 or 2 fs is used widely (I say 1 fs is better, but that is your choice and your advisor). What makes flexible model sometimes defficult to deal with is the possiblility of hydrogens running away, resulting in a clash. Think of this, there is a dummy atom ( has a mass close to zero like 1e-100), when a force is applied to such atom, what is gonna happen? Next thing you know your simulation is blown up !
One other important note is that you are dealing with a liquid water, if you use aniso or iso coupling pressure approach, the final outcome should not be different at all. Furthermore, even though 993 g/L is not that much different from 997 g/L, it still affects the equation of state for water !
Coming back to the dammping factor, as much as it is not important to some people, but sometimes wrong chosing of damping factor can lead you to incorrect compressibility factor.
Finally, this thread is when I started learning lammps back in the days, so the script that you are seeing is not as accurate as it should.