Big Water Hole in Bilayer Graphene Oxide Simulations by a NPT run

LAMMPS LAMMPS Mailing List Mirror
#1

Dear all,

I used LAMMPS to perform NPT simulations of bilayer graphene oxide in water.

The water model is TIP4P/2005 (Abascal, J. L. F. & Vega, C. *The Journal of
Chemical Physics* *123*, 234505-234512 (2005).)
Graphene edges is functionalized by hydroxyl groups. The initial
coordinates of chemical groups are taken from previous DFT study (Yan,
J.-A. & Chou, M. Y. *Physical Review B* *82*, 125403 (2010).). The force
field for graphene oxide I use in MD simulations is OPLS-AA.

The time step is 0.1 fs. The total simulation time is 100 ps. The
boundaries in three directions is all periodic.
The ensemble is NPT of 298K thermostat and 1atm isotropic barostat.
(fix 4 all npt temp 298.0 298.0 100.0 iso 1.0 1.0 1000.0)

The initial state is a non-equilibrium configuration. The two layers of
graphene oxide tend to be stuck together. The remaining space should be
filled by surrounding water molecules. However, a large vacuum hole emerges
and the simulation box won't contract its size to decrease the vacuum hole.
The attached PNG images are the initial and final frames of my simulation.
The formation of the vacuum hole seems to be unrealistic. Is there anyone
have the same problems in non-equilibrium simulations within water
molecules? I am wondering how to resolve this issue.

I found a paper (Patel, R. Y. & Balaji, P. V. *The Journal of Physical
Chemistry B* *109*, 14667-14674 (2005)) raise a question for different MD
results of isotropic and anisotropic pressure coupling methods. I can
change my pressure coupling method from "iso" to "aniso", which the xyz
pressures are "couple none". (fix 4 all npt temp 298.0 298.0
100.0 aniso 1.0 1.0 1000.0). Is it the key to solve this problem?

Thank you in advance.

Best,
-Wenpeng

The input script and log file of my simulation are attached as below:
########### input script ##########
# System and potential definition
units real
boundary p p p
atom_style full
# processors 2 2 2

read_data ./data.GO-WM

# Pair Coeffs

int.png

final.png

#2

Dear all,

I used LAMMPS to perform NPT simulations of bilayer graphene oxide in water.

The water model is TIP4P/2005 (Abascal, J. L. F. & Vega, C. *The Journal of
Chemical Physics* *123*, 234505-234512 (2005).)
Graphene edges is functionalized by hydroxyl groups. The initial
coordinates of chemical groups are taken from previous DFT study (Yan,
J.-A. & Chou, M. Y. *Physical Review B* *82*, 125403 (2010).). The force
field for graphene oxide I use in MD simulations is OPLS-AA.

The time step is 0.1 fs. The total simulation time is 100 ps. The

0.1fs seems very short. if you run shake not only on the water, but also
on the hydroxyl groups, you should be able to go for 1fs.

a 100ps simulation is _nothingness_.
have you equilibrated the water system
before inserting the the object. different
water models have different equilibrium
densities, which can be away from the
experimental density of water.

boundaries in three directions is all periodic.
The ensemble is NPT of 298K thermostat and 1atm isotropic barostat.
(fix 4 all npt temp 298.0 298.0 100.0 iso 1.0 1.0 1000.0)

The initial state is a non-equilibrium configuration. The two layers of
graphene oxide tend to be stuck together. The remaining space should be
filled by surrounding water molecules. However, a large vacuum hole emerges
and the simulation box won't contract its size to decrease the vacuum hole.
The attached PNG images are the initial and final frames of my simulation.
The formation of the vacuum hole seems to be unrealistic. Is there anyone
have the same problems in non-equilibrium simulations within water
molecules? I am wondering how to resolve this issue.

this kind of hole can happen often when setting up a composite system.

one way to "solve" it would be to first _only_ equilibrate the solvent
and not time integrate the inserted object at all (so it remains where
it is). what you have happening as a consequence of your setup is
akin to a massive nanoscale implosion. you probably should use
aniosotropic cell relaxation and run at a higher target pressure than
1 atmosphere (which is almost nothing on such a small scale), so
that you "squeeze" your system a bit. once you are sure there is
no vacuum bubble left, you can reset the pressure target and equilibrate
some more. it may also be helpful to reinitialize the velocities a
few times, to try dissipate any shockwaves, or even use a langevin
thermostat in combination with a berendsen barostat at the beginning
of your equilibration.

once you have your water box properly equilibrated around your
inserted object, you can start time integrating all atoms and
continue equilibration.

please note, that similar problems are very common in simulations
of solvated proteins, so you may want to look at the typical simulation
set up instructions and tutorials for such simulations to get an idea
how to go about these things.

cheers,
    axel.

I found a paper (Patel, R. Y. & Balaji, P. V. *The Journal of Physical
Chemistry B* *109*, 14667-14674 (2005)) raise a question for different MD
results of isotropic and anisotropic pressure coupling methods. I can
change my pressure coupling method from "iso" to "aniso", which the xyz
pressures are "couple none". (fix 4 all npt temp 298.0 298.0
100.0 aniso 1.0 1.0 1000.0). Is it the key to solve this problem?

you also will need to simulate longer and

#3

Another path that you could use could be to create a box with periodic
boundaries in x and y but not in z. Then fix the first graphene layer
at the bottom of the box, place all the H2O molecules on top of it and
enclose it at the top with the second graphene layer that will be
rigid but mobile (note that the simulation box dimensions in x and y
would have to be tuned to the graphene lattice structure). Next, apply
a force to the top rigid layer that corresponds to your desired
pressure and wait until equilibration is reached by measuring the
center of mass z-coordinate of the top graphene wall. When
equilibration is attained, dump the data and cut appropriately the
extra vacuum space above the mobile wall by redefining the box
dimension in z. Also, delete (if needed) graphene edge atoms within
the skin region about the x, y boundaries. Finally, use this new
config as the input for new NPT simulations where the z dimension is
now periodic as well. That should work. Yet, I don't know if you would
find it to be an easier approach.
Carlos

#4

carlos,

your suggestion has one flaw.

since you need to cut and stitch the non-periodic
boundary back together, you create the same
problem as before: you need to avoid overlap,
so you need to leave extra space and you get
an implosion resulting in a shockwave again.
and in this case it is probably worse, since it
is directional and not as easily diffused as the
one from embedding into a well equilibrated
periodic water box. also, you'll need to run
for quite a while to "heal" the volume, where
to combine the water back together.

axel.

#5

Hi Axel,
I knew that my solution only works perfectly for the infinite layer. That’s why I wrote “if needed” between brackets within my response. It wasn’t clear to me from the post that the finite size of the object mattered. And in many cases the people posting here don’t fully disclose what they are after. Anyways, yes you are correct that my proposal is best for the infinite layer.
Carlos