LAMMPS System Instability Switching from Berendsen to Nose-Hoover Thermostat

Hello all,

I am trying to simulate a graphene-like sheet with Napthalene - like hydrophobic molecules all in a periodic box of TIP3P water. I am treating the graphene sheet as a rigid body which cannot translate or rotate in the box, and the hydrophobic molecules as independent rigid bodies which have translational and rotational degrees of freedom.

After minimization, I run an initial equilibration simulation (NPT using Berendsen thermostat and Nose-Hoover barostat) for 500 ps. After the equilibration run, I switch to Nose-Hoover thermostat for the production run (using restart files). However, upon switching thermostats, the system destabilizes and expands. I am testing two similar systems (see details below) and have attached graphs of the system volume for both systems, as well as the input script.

Does anyone have ideas as to what could be causing the instability in the system upon switching thermostats in a (seemingly) equilibrated system? Are the restart files an issue or is something else the matter?

I am testing two different sets of LJ parameters on the graphene sheet. For one set of parameters (call this system A), the system expands ~10% and then returns to the equilibrated volume within ~200ps. For the more hydrophobic set of parameters (call this system B), the system expands by > 100% and does not appear to return to the equilibrated volume. To me, the significant difference in the system stability for two different sets of LJ parameters is surprising.

Thanks for all the help,

Ryan DeFever
Clemson University
Chemical Engineering '14

Further details, input scripts, and figures:

I have tested the restart files by just restarting the simulation under the same conditions (NPT, Berendsen thermostat, Nose-Hoover barostat) as the equilibration run. This appears to work fine, as the volume remains stable. The two simulations use the same starting configurations. All I change is the LJ pair coefficients for the carbon in the sheet.

Input Script for production run. Restarting from equilibration run. (7Jun13 Build of LAMMPS)

log log.####sim
units real
atom_style full
pair_style lj/cut/coul/long 10.0 10.0
pair_modify mix geometric
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0.0 0.0 0.5

read_restart ######.rst

angle_coeff 1 55 104.52
kspace_style pppm 0.0001
neigh_modify delay 0 every 1 check yes

group water type 23 24
group npth type 19 20 21 22
group sheet type 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

compute temp_sheet sheet temp
compute temp_npth npth temp

thermo 10
velocity sheet set 0 0 0 units box # This is to keep the sheet from having a slight rotation left over from the restart files.
timestep 2

dump 1 all xyz 10000
dump 2 all xtc 100 dump.######.xtc
dump_modify 1 sort id element C1 C11 C2 C12 C3 C13 C4 C14 C5 C15 O1 O11 O2 O12 H1 H11 H2 H12 CN1 CN2 HN1 HN2 Ow Hw

fix 1 all shake 0.0001 10 0 b 1 a 1
fix 2 water npt temp 300.0 300.0 500 iso 1.0 1.0 1000.0
fix 3 sheet rigid/nve single torque * off off off force * off off off
fix 4 npth rigid/nve molecule
fix 5 all momentum 100 linear 1 1 1 angular

thermo_style custom step temp press vol etotal ke pe evdwl ecoul elong

restart 10000 #####.rst #####.rst1
run 5250000

Links to figures (attached, but in case those do not work):

System A:
System B:

System A:
Thermostat issues 2 (A).png
System B:

Thermostat issues (B).png