The Case of The Exploding Water
To run water at 1 atm between carbon sheets, my usual procedure is:
- Simulate a slab of water being “squeezed” and then “expanded” between moving carbon sheets
- Track the middle 20% density during said simulation
- Compare with NPT density of a bulk simulation at specified pressure
It was doing this last step that I saw something very strange: I was (barely) getting a stable simulation – of bulk water!! – with a density one thousandth of what it’s meant to be: 1.5 milligrams per cubic centimeter instead of 1 g/cm^3.
There was nothing wrong with my LAMMPS commands (I’ve attached files below), but there was an issue with my input data file.
The challenge: without visualizing, can you spot the error? (Frankly a visualizer doesn’t help much with this one unless you know what you’re looking for.)
If not, can you at least figure out what physical state is being simulated?
Data file header:
LAMMPS data file via write_data, version 29 Aug 2024, timestep = 500000, units = real
4800 atoms
2 atom types
3200 bonds
1 bond types
1600 angles
1 angle types
0 34.08 xlo xhi
0 39.35219434796489 ylo yhi
-18.8 18.8 zlo zhi
Masses
1 15.994
2 1.008
Pair Coeffs # lj/cut/coul/long
1 0.163406 3.17427
2 0 0
Bond Coeffs # zero
1 0.9789
Angle Coeffs # zero
1 109.47
Atoms # full
1 1 1 -0.447585 1.6014459976669704 3.2757603722079485 14.448500463436522 1 1 -1
2 1 2 0.2237925 1.327766340063397 4.057068563140284 13.926100052673606 1 1 -1
3 1 2 0.2237925 1.8649951332391688 2.553865094584665 13.842155289287193 1 1 -1
4 2 1 -0.447585 4.077378649855302 20.265193327446948 -7.719531104858056 0 1 1
5 2 2 0.2237925 4.243790486345037 21.199474368025154 -7.4793826228196405 0 1 1
6 2 2 0.2237925 4.895302425102546 19.871948474183146 -8.086415833146178 0 1 1
7 3 1 -0.447585 27.005964635076445 35.68058783031934 10.041600207753303 0 0 0
8 3 2 0.2237925 27.41416614821016 36.532291720847645 10.298927346689505 0 0 0
9 3 2 0.2237925 27.689212928151946 34.979593471014944 10.036706865825293 0 0 0
...
Log file snippet (timestep = 1fs, units = real so energies in kcal/mol, temperature in K, pressure in atm)
Step Temp TotEng PotEng KinEng Press Volume Density
0 293.33555 -161.01637 -2958.1502 2797.1338 6901.0076 50426.217 0.94891294
25000 301.25117 2462.0174 -410.59664 2872.6141 80.258098 590675.38 0.081009115
50000 299.8129 2695.8524 -163.04689 2858.8993 9.1856863 4413567.6 0.01084159
75000 299.20019 2689.0934 -163.96338 2853.0567 1.5604506 16676422 0.0028693259
100000 298.00004 2681.7676 -159.84497 2841.6126 1.0117224 24361350 0.0019641805
125000 300.66326 2707.1174 -159.89063 2867.008 0.98702617 26635693 0.001796465
150000 300.37224 2702.788 -161.44498 2864.233 1.0218432 28080397 0.0017040389
175000 297.19658 2673.9455 -160.00557 2833.951 1.0027294 29672859 0.0016125878
200000 303.92651 2735.4474 -162.67773 2898.1251 0.96527529 32107356 0.0014903155
225000 298.71974 2682.788 -165.68738 2848.4754 0.98477773 32042813 0.0014933174
250000 297.79069 2677.8596 -161.75665 2839.6163 1.0346839 32465452 0.0014738772
...
Hint 1: What state is the simulation in?
It’s a gas. You can tell because the volume fluctuations are large but the pressure fluctuations are small, indicating a large compressibility.
Hint 2
Check the charges.
Answer
You can see in the data file that I set hydrogens to a charge of +0.2237925 and oxygens to -0.447585. This is of course exactly half of what OPC3 prescribes, based on me copying wrong from the paper (HTML), but any familiarity with a water model should set alarm bells off – SPC/E is about -0.85 on the oxygens!
Based on the inputs, I was simulating a system of water molecules with half the dipole strength. No wonder it was turning into vapor at 300 K and 1 atm!
Lessons
- In real life, water has an immensely high boiling point (and melting point!) relative to its molecular mass. This fail of mine actually demonstrates the point nicely – simulated with half the usual dipole, water boils at room temperature and pressure! Of course real-life water has hydrogen-bonding as well as dipole-dipole interactions, so the electrostatic dipole of a computational water molecule is probably larger than its actual dipole to compensate.
- When we give advice and say “the simulation does what the force field says”, this is what we mean. A water system with half the usual dipole should vaporize at room temperature. LAMMPS faithfully simulated the physics of my system – I just wasn’t simulating water as it really exists.
- For troubleshooting a failed simulation, always check the log files and understand exactly what the thermodynamic significance of each observable is. For liquid water the pressure of a very small volume should fluctuate significantly. When I was troubleshooting, the very small fluctuations in pressure, and correspondingly small potential energy, was what tipped me off. I double-checked everything that gives water cohesive energy – the LJ parameters? Seemed correct. Charge? That’s when I realized what had gone wrong.
- To verify this for yourself, you can update the charges using two “set type 1 charge …” commands, and run the NPT simulation again to get a much more sensible density of 0.99 g/cm^3.
- Even people who are experienced with LAMMPS make hilarious mistakes. The experience helps us check and recover quickly.
Attachments:
equil-nvt.data (804.6 KB)
npt.lmp (585 Bytes)
(Note: the index variable type for nthermo allows you to run shorter runs from the command line, for example lmp -i npt.lmp -var nthermo 500.
By default, nthermo is set to 25000; the input runs for 20 times nthermo steps and a restart and XTC dump is generated every 4 times nthermo steps.
From memory, dump ... xtc comes from the EXTRA-DUMP package, and you may have to comment that line out if running on a LAMMPS executable without that package.)