LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] negative vapor pressure
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] negative vapor pressure

From: Brian Morrow <brianhmorrow@...24...>
Date: Wed, 23 Aug 2017 17:54:47 -0400

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.5*Lz*(<Pzz>-0.5*(<Pxx>+<Pyy>)), 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. 


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