Large remain pressure in graphene/oxygen after equilibration by NPT

LAMMPS LAMMPS Mailing List Mirror
1

Dear All,

I use lammps-9Dec14 version.

I have a system of a graphene sheet, plus oxygen atoms on both sides of

the sheet. The oxygen atoms are placed not closer than 1 nm from the sheet.

See attached figure, graphene-oxygen.jpg for the initial structure.

The reactive force field potential is used. NPT is used for graphene equilibration at 300K and zero pressure in all directions. Oxygen atoms are kept stationary during equilibration.

I find that the output pressure in the log file remains to be large after equilibration (the log file is attached as equi.lammps). I past below the thermpdynamic info near the end of

simulation.

Step tempMobi pressMob Temp Press Volume

9200 298.69932 -312.07759 199.08861 -312.07759 413496.43
9300 295.30714 -905.86599 196.82766 -905.86599 413556.73
9400 301.45809 -87.312429 200.92738 -87.312429 413566.18
9500 295.21937 -1037.9511 196.76916 -1037.9511 413605.69
9600 291.60771 892.00661 194.36192 892.00661 413297.37
9700 305.64637 103.38666 203.71895 103.38666 413364.74
9800 308.01659 -454.75043 205.29874 -454.75043 413508.28
9900 300.22849 -1020.2742 200.10783 -1020.2742 413478.3
10000 285.29167 774.90448 190.15216 774.90448 413418.08

Could anyone help me in this problem?
I paste below my input script. The structure data is attached (CHO.dat).

Thanks,

Jie

units real
atom_style charge
boundary p p p

read_data CHO.data

pair_style reax/c lmp_control
pair_coeff * * ffield.reax.cho C O

neighbor 2 bin
neigh_modify every 10 delay 0 check no

timestep 0.2

group mobile type 1
group immobile type 2

fix center all recenter INIT INIT INIT units box

dump position all atom 1000 qui.lammpstrj

thermo 100

neigh_modify exclude group immobile immobile

fix 1 all qeq/reax 1 0.0 10.0 1e-6 param.qeq

compute tempMobile mobile temp
compute pressMobile all pressure tempMobile

thermo_style custom step c_tempMobile c_pressMobile temp press vol

fix Fmovestuff mobile npt temp 300.0 300.0 10.0 iso 0.0 0.0 200.0 dilate mobile

fix_modify Fmovestuff temp tempMobile

run 10000

write_restart system_after_fix_rigid.rst

CHO.data (130 KB)

graphene-oxygen.jpg

equi.lammps (9.42 KB)

2

Dear All,

I use lammps-9Dec14 version.

I have a system of a graphene sheet, plus oxygen atoms on both sides of
the sheet. The oxygen atoms are placed not closer than 1 nm from the sheet.
See attached figure, graphene-oxygen.jpg for the initial structure.

The reactive force field potential is used. NPT is used for graphene
equilibration at 300K and zero pressure in all directions. Oxygen atoms are
kept stationary during equilibration.

I find that the output pressure in the log file remains to be large after
equilibration (the log file is attached as equi.lammps). I past below the
thermpdynamic info near the end of
simulation.

Step tempMobi pressMob Temp Press Volume

    9200 298.69932 -312.07759 199.08861 -312.07759 413496.43
    9300 295.30714 -905.86599 196.82766 -905.86599 413556.73
    9400 301.45809 -87.312429 200.92738 -87.312429 413566.18
    9500 295.21937 -1037.9511 196.76916 -1037.9511 413605.69
    9600 291.60771 892.00661 194.36192 892.00661 413297.37
    9700 305.64637 103.38666 203.71895 103.38666 413364.74
    9800 308.01659 -454.75043 205.29874 -454.75043 413508.28
    9900 300.22849 -1020.2742 200.10783 -1020.2742 413478.3
   10000 285.29167 774.90448 190.15216 774.90448 413418.08

Could anyone help me in this problem?

you request isotropic volume change, but your system is not isotropic.
also, how can you relax pressure contributions from your oxygens, if
they are kept immobile?

also, it is a *very* bad idea to use pairwise interaction exclusions
like you do with a manybody potential. you create inconsistent
handling of manybody interactions. if you want to keep the oxygens
immobile do not time integrate them or set their resulting forces to
zero with fix setforce and initialize their velocities to zero, too.

axel.

3

Dear Axel,

Thanks.

4

Dear Axel,

Thanks.

you request isotropic volume change, but your system is not isotropic.
also, how can you relax pressure contributions from your oxygens, if
they are kept immobile?

I replace iso by x, y, z component,

fix Fmovestuff mobile npt temp 300.0 300.0 10.0 x 0.0 0.0 200.0 y 0.0 0.0
200.0 z 0.0 0.0 200.0 dilate mobile

Yes. oxygens are kept immobile. What I want is to relax graphene so that
pressure in graphene plane (in x and y) is zero. The pressure contributions
from your oxygens can be finite. Is there a way to do that?

the only way i see that this can be done cleanly would be in a
simulation that *only* contains the graphene.
after that system is relaxed, you add the oxygen atoms.

axel

5

Dear Axel,

Thanks.

For curious. After I add
"
velocity immobile zero angular
velocity immobile zero linear
fix 2 immobile setforce 0.0 0.0 0.0

" for oxygen atoms,

and delete iso volume change by setting

"

fix Fmovestuff mobile npt temp 300.0 300.0 10.0 x 0.0 0.0 200.0 y 0.0 0.0 200.0 z 0.0 0.0 200.0 couple none dilate mobile" for graphene,

and delete

“neigh_modify exclude group immobile immobile”,

I run the new input script, the pressure in the out put is still around 1000,

"Step tempMobi Press Volume

9200 273.52866 -503.97012 329780.24
9300 312.31974 -429.30177 329761.54
9400 310.36913 -656.97742 329864.66
9500 302.72363 205.82644 329690.58
9600 285.41184 685.77908 329769.21
9700 299.75422 -1350.3062 329779.21
9800 305.65174 -784.1244 329817.13
9900 292.8356 -1019.6085 329918.59
10000 314.19726 802.96647 329702.25
"

So the large pressure is from freeze oxygen atoms?

Thanks,

Jie