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