I setup a simulation to calculate CO2 properties using the known EPM2-flexible model. My system has 1000 CO2 molecules.
The forcefield parameters are as below:
In system.in.init:
units real
atom_style full
bond_style harmonic
angle_style harmonic
pair_style lj/cut/coul/long 10.0
kspace_style pppm 1.0e-5
pair_modify mix arithmetic
In system.in.settings:
bond_coeff 1 1283.381 1.149
angle_coeff 1 147.5 180
pair_coeff 1 1 0.055898 2.757
pair_coeff 2 2 0.159984 3.033
group co2 type 1 2
special_bonds lj/coul 0 0 0
Yet, my calculated self-diffusivity is around 2.3 times larger than the reported diffusivity for EPM2-flexible model with the same forcefield parameters. It also deviates almost the same from the experimental data. The viscosity deviates only around 12%, though. I calculate the fluid properties at the following conditions: 50 Mpa and 424 K. Do you have any idea why my simulation cannot reproduce the data reported in the literature despite using the same forcefield? I ran my simulation using two other systems with 2000 and 8000 number of molecules but got the same results.
The paper specifies a cutoff of 14 angstroms but you have 10 angstroms instead. Please note I have not exhaustively looked for any other differences; this is simply the first one I saw and most obviously likely to impact your results.
Thanks for the reply. I also changed the cutoff to 14 A but the results didn’t change. They also used velocity rescaling for thermostat while I used nose-hoover.
I just looked through your input files. You did not note in your initial post that you performed an NPT simulation. The paper describes obtaining their results with NVT simulations.
This matters tremendously, because firstly the barostat by itself changes the observed molecular dynamics. Secondly, a fluctuating periodic box makes the “unwrapped” displacement ambiguous to define or calculate. A recent paper gives a good overview of the difficulties involved and proposes a method to overcome them.
Long story short, try running with NVT instead. The paper is reasonably descriptive about how they did things and while there are some points that are sadly missing, you should be able to reproduce their results with some care.
Thanks for the helpful comment. An idea popped into my mind. Instead of defining the volume based on the experimental density at the given T and P, can I do NPT first and then do NVT and calculate Ds and viscosity in the second NVT?
Dear @srtee, I would like to appreciate your feedback kindly. Using NVT and predetermined density solved my problem. However, as expected, the pressure oscillates very much!
Well yes, you are simulating a very small block of a supercritical fluid. There are enough discussions on “pressure fluctuations”, if you search this forum, for someone to write a very long (and very repetitive) book.
Congratulations for successfully identifying the cause of the discrepancy though!
…and supercritical fluids are well known for their large pressure (and local density!) fluctuations, especially near the critical point. I recall doing experiments in physical chemistry lab in college on determining the critical point where the clear liquid turning opaque was an indication of reaching the critical point.
After reading this paper, I now wonder what coordinates LAMMPS uses to calculate the diffusivity—wrapped or unwrapped? I guess it should be unwrapped coordinates, as they are mostly affected by the barostat. Also, are other properties affected by these unbound displacements? properties like enthalpy of vaporization, density, RDF, solubility, heat capacity, etc.?
The LAMMPS source code is quite readable if you would like to work out the answer to these questions yourself (certainly more than other comparable molecular dynamics source codes). In addition, it is easy to continuously print out coordinates and velocities and then import those into Python (using MDAnalysis, for example) and then write your own property-calculation scripts.
LAMMPS uses unwrapped coordinates to calculate MSD, clearly. What the paper is telling you is that how to unwrap coordinates is not straightforward in a simulation with changing box size.
Now consider how you would calculate, say, a density from a molecular dynamics simulation. Do you need to unwrap particle coordinates to calculate density? Concretely speaking – let’s say you have a snapshot of your simulation and one of the CO2 molecules is very near the box edge. If you were calculating the MSD, would it make a difference whether the CO2 had gotten there by crossing a box boundary? What about if you were calculating the density?
Please note that people have been discussing various artifacts of PBCs for a very long time – such as size dependence of diffusivities (like the Yeh and Hummer correction) and even NPT diffusivities (Hummer had a paper about it ten years ago – I happened to link this one because it’s more complete and recent). Molecular dynamics is sometimes half art, but it’s definitely not black magic. There are sensible procedures to work out whether aspects of your simulation setup cause artificial changes to observed properties.
A final point to consider is that if you are simulating a supercritical fluid it is not surprising that NPT and NVT ensemble averages may show very different fluctuations and convergence times. The density is liquid-like so particles have very short mean free paths and collide often – but the compressibility is gas-like, so with a barostat the box size will fluctuate significantly. You should survey the literature on simulating supercritical fluids – or more generally gas-liquid-supercritical phase transitions – if you want to be confident in your results.