negative vapor pressure

Hi all,

I want to use LAMMPS to simulate hydrocarbon liquid-vapor interfaces, to estimate critical points and calculate surface tension as a function of temperature.

My simulation setup is a liquid slab in the center of a 50x50x250 A simulation box, with the plane of the liquid-vapor interface perpendicular to the z axis. I’m calculating the surface tension as 0.5Lz(-0.5*(+)), where Lz is the box length and Pii are pressure tensor components. Pzz is the vapor pressure.

When I use the OPLS-AA force field with pppm/disp to calculate long-range interactions, I get surface tensions and liquid densities in reasonable agreement with experiment, but the vapor pressures are very negative, on the order of -50 atm depending on the molecule and temperature. I’ve tried different system sizes and permutations of my input file, and the only thing that seems to give positive vapor pressures is turning off SHAKE and using a 0.25 fs timestep. I’m not sure why this is happening; bulk liquid NPT simulations with SHAKE and a 2 fs timestep give correct densities, and energy is conserved in NVE interface simulations with SHAKE and a 2 fs timestep.

I feel like I’m using a command incorrectly, or am overlooking something else obvious, but I haven’t been able to figure out what. A sample input file is below, and I can send a data file if needed. I’m using the 11 Aug 2017 version of LAMMPS. Any suggestions would be much appreciated.

Thanks,

Brian

atom_style full

units real

boundary p p p

pair_style lj/long/coul/long long long 10.0

pair_modify mix geometric

kspace_style pppm/disp 1.0e-6

kspace_modify force/disp/real 0.0001

kspace_modify force/disp/kspace 0.002

bond_style harmonic

angle_style harmonic

dihedral_style opls

special_bonds lj/coul 0.0 0.0 0.5

read_data …/…/data.decalin1000.360.5K.lammps

variable dt equal 2.0

timestep ${dt}

variable tdamp equal ${dt}*100

thermo_style custom step temp pe etotal press pxx pyy pzz

thermos 500

fix press_tensor all ave/time 1 100 100 c_thermo_press[1] c_thermo_press[2] c_thermo_press[3] file press_tensor.dat

fix bal all balance 100 1.0 shift xyz 20 1.0

fix thermostat all nvt temp 360.5 360.5 ${tdamp}

fix const all shake 0.0001 20 10 b 2 a 2 3

run 1000000

Saw that this message was unanswered. Possibly Aidan has

some ideas (CCd).

Steve