[lammps-users] pressure vibration for supercritical CO2 system

LAMMPS LAMMPS Mailing List Mirror
1

I simulated a small system containing only 400 CO2 molecules at 308.15K(NVT).
I get consistent values for coulombic energy and total energies using Ewald
and pppm methods for long potential.
But, the final pressure is jumping all over the place. The box is 32×32×32
with CO2 density 0.89, which is corresponding to the CO2 critical state.
Therefor the expected average pressure should be 300atm finaly,
but my obtained result is 1013.17atm.

The typical results are obtained from last several steps of simulation only.
(1×10^6 steps in all, Ewald)

--------------- Step 998000 ---- CPU = 32206.6205 (sec) ---------------
Total E = 677.5437 Total KE= 1113.6373 Temp = 311.5973
Total PE= -436.0936 E_bond = 194.7543 E_angle = 158.9376
E_dihed = 0.0000 E_impr = 0.0000 E_vdwl = -663.3588
E_coul = 9596.0125 E_long = -9722.4393 Press = 1459.4498

--------------- Step 999000 ---- CPU = 32327.9993 (sec) ---------------
Total E = 627.5608 Total KE= 1074.9995 Temp = 300.7864
Total PE= -447.4387 E_bond = 193.6396 E_angle = 152.9668
E_dihed = 0.0000 E_impr = 0.0000 E_vdwl = -656.2266
E_coul = 9585.1683 E_long = -9722.9868 Press = 39.0421

--------------- Step 1000000 ---- CPU = 32368.4202 (sec) ---------------
Total E = 651.9472 Total KE= 1086.1959 Temp = 303.9192
Total PE= -434.2487 E_bond = 193.8684 E_angle = 166.6106
E_dihed = 0.0000 E_impr = 0.0000 E_vdwl = -654.6532
E_coul = 9583.0347 E_long = -9723.1093 Press = 1469.1984

The input script I use is below.(lammps 2001)
My input date file is attached to this email.

I appreciate any input.

#real system with NVT ensemble

units real
extra memory 1.5 1.5 4.0 1.5
maximum cutoff 18
neighbor 1 1 1 5 1
dimension 3
periodicity 0 0 0

read data data.co2

coulomb style ewald 10.0 1.0E-4
dielectric 1.0
nonbond style lj/cutoff 15.0 1
nonbond coeff 1 1 0.055639162 2.757 10.0
nonbond coeff 1 2 0.094128308 2.895 10.0
nonbond coeff 2 2 0.159242846 3.033 10.0

bond style harmonic
bond coeff 1 450 1.149

angle style harmonic
angle coeff 1 147.7 180

temp control nose/hoover 308.15 308.15 0.001

thermo flag 1000
thermo style 2

dump atoms 1000 dump.lj.nvt

restart 10000 2 file1 file2

timestep 1

run 1000000

data-isotropic.co2 (72.3 KB)

2

dear xzj(??),

I simulated a small system containing only 400 CO2 molecules at 308.15K(NVT).
I get consistent values for coulombic energy and total energies using Ewald
and pppm methods for long potential.
But, the final pressure is jumping all over the place. The box is 32��32��32
with CO2 density 0.89, which is corresponding to the CO2 critical state.
Therefor the expected average pressure should be 300atm finaly,
but my obtained result is 1013.17atm.

there are two concerns: 1) is the 300atm the experimental pressure or
for the potential parameters that you are using (or in other words, are
the parameters that you are using optimized for near-critical
conditions). 2) a typical property of materials at the critical point
are large density fluctuations, which make 'small' classical MD
simulations usually victims of serious finite size effects, so that
critical properties are usually approximated by aproaching the critical
conditions from above and below (following coexistance curves) and
then extrapolate to the critical point.

cheers,
  axel.

3

A 400 molecule simulation is so small that it is normal to see
large pressure fluctuations. If you are cutting off the LJ interactions
then you are missing a tail correction to the pressure. See the
pair_modify tail command for how to include the correction. There
is also a thermo_style custom pave option which will time-average
the pressure values printed by thermo (also see thermo_modify window),
which might help you get smoother results.

Steve