I am trying to simulate 2 layer graphene sheet(15nmX15nm) inside water box(30nmX30nm) (using lammps version released on 17 Aug 2017). Modelling is done using vmd and merging between water and graphene sheet following the axels tutorial(https://sites.google.com/site/akohlmey/software/topotools/topotools-tutorial—various-tips-tricks).
The interaction potential parameters are taken from references as shown in system.in.settings file. I want to simulate the system at lab conditions of (p = 1 atm and t = 298K), so I use npt after equilibrating the system for 5000 time steps of size 2 femtoseconds. At the end of equilibration pressure is of 5 orders and when npt is used in next phase the size of the box start to increase. I tried following methods:
- Graphene sheet is fixed using rigid single, and langevin thermostat along with nve/limit 0.1 on water for equilibration and then for next phase npt with dilate all option.
- Graphene sheet is fixed using rigid/nve single, and langevin thermostat along with nve/limit 0.1 on water for equilibration and then for next phase npt with dilate all option.
- Graphene sheet is fixed using rigid/nvt single, and langevin thermostat along with nve/limit 0.1 on water for equilibration and then for next phase npt with dilate all option.
- Then in second set, above three cases along with rigid I’ve also used neigh_modify to exclude the interaction between C atoms of graphene sheet and in third set, apart from exclude pair interactions delete_bond multi command for C atoms.
Box expansion was less in cases when only rigid was used and no pair exclusion and delete bond instruction.
I also simulated water box without graphene sheet and found that pressure is too high for the equilibration phase involving langevin and nve/limit, but comes down to 1 atm.
Please go through my input files and suggest where is the glitch.
I came across similar querry on list, where reply was that the reason could be ill geometry or overlapping. if so, how could I figure out this? and solve it.
contents of system.in.init:
boundary p p p
(Hybrid force fields were not necessary but are used for portability.
pair_style hybrid/overlay lj/cut/coul/long 10.0 lj/cut 10.0 # hybrid lj/charmm/coul/long 9.0 10.0 10.0 lj/cut/coul/long 10.0 #hybrid lj/charmm/coul/long 9.0 10.0 10.0
pair_modify tail yes
bond_style hybrid harmonic
angle_style hybrid harmonic
kspace_style pppm 0.00001