# The internal pressure isn't equal to the external pressure

In my system, there are some methane molecules inside two graphite walls, so it’s a nano-channel simulation with the p p f boundary condition.

Fixing an atome on the lower wall in x y z direction and an atome on the upper wall in x y direction (free in z direction like a piston) and applying 10bar pressure both on the upper and lower wall, I expect that the pressure of methane molecules is equal to 10bar, but after 4e6 steps, the pressure converge to about 60 bar. I have examined my competing pressure method in a square box which just contain methane molecules with the p p p boundary condition, the result is identity with the default calculated value.

Meng

# Methane&Graphene equilibrium MD simulations

atom_style atomic

units metal

dimension 3

boundary p p f

neighbor 0.2 bin

neigh_modify delay 10 every 1 check yes

mass 1 12.0 # C

mass 2 16.0 # CH4

variable Ctype equal 1

variable CH4type equal 2

variable upperfixID equal 1

variable lowerfixID equal 449

variable lowerlength equal 4.0

variable upperlength equal 81.0

variable extpress equal 6.2415E-6 #atm to 6.2415E-7eV/A^3 = 1bar

variable kB equal 8.6173324E-5 # eV/K

variable p equal 1000 # correlation length

variable s equal 10 # sample interval

variable d equal \$p*\$s # dump interval

variable T equal 148.1

variable dt equal 0.001

# Define interaction parameters

pair_style hybrid airebo 3.0 1 1 lj/cut 10.0

# AIREBO:

pair_coeff * * airebo CH.airebo C C

pair_coeff {Ctype} {CH4type} lj/cut 0.00555 3.6 10

pair_coeff {CH4type} {CH4type} lj/cut 0.01276 3.8 10

# ------------- Regions and Groups ----------------

group graphite type \${Ctype} # all graphite atoms

group methane type \${CH4type} # all methane atoms

group fixpoints id {lowerfixID} {upperfixID}

group lowerfix id \${lowerfixID}

group upperfix id \${upperfixID}

region 1 block INF INF INF INF INF \${lowerlength} units box

group lowerwall region 1

region 2 block INF INF INF INF \${upperlength} INF units box

group upperwall region 2

variable envforce equal \${extpress}lxly/count(upperwall)

compute peratom methane stress/atom NULL

compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3]

variable prmethane equal -(c_p[1]+c_p[2]+c_p[3])/(3lxly*(bound(methane,zmax)-bound(methane,zmin)))

variable lengthz equal (bound(methane,zmax)-bound(methane,zmin))

thermo \$d

thermo_style custom step temp lx ly lz v_lengthz v_prmethane

fix pressure methane ave/time \$s \$p \$d v_prmethane file pressure.time

# ------------- Initial velocity ----------------

velocity graphite create \$T 48748 units box

velocity graphite zero linear

velocity methane create \$T 488 units box

velocity methane zero linear

# ------------- Boundary conditions ----------------

fix prlow lowerwall addforce 0.0 0.0 \${envforce}

fix prup upperwall addforce 0.0 0.0 -\${envforce}

velocity lowerfix set 0.0 0.0 0.0

velocity upperfix set 0.0 0.0 NULL

fix 13 lowerfix setforce 0.0 0.0 0.0

fix 12 upperfix setforce 0.0 0.0 NULL

fix wnve methane nve

fix cnvt graphite nvt temp \$T \$T 0.1

run 4000000

write_restart pressuretest4e6_1776.rest

In my system, there are some methane molecules inside two graphite walls, so
it’s a nano-channel simulation with the p p f boundary condition.

Fixing an atome on the lower wall in x y z direction and an atome on the
upper wall in x y direction (free in z direction like a piston) and applying
10bar pressure both on the upper and lower wall, I expect that the pressure
of methane molecules is equal to 10bar, but after 4e6 steps, the pressure
converge to about 60 bar. I have examined my competing pressure method in a
square box which just contain methane molecules with the p p p boundary
condition, the result is identity with the default calculated value.

i don't understand what that is supposed to prove. you are comparing
the proverbial apples and oranges.

i have some general comments about the input below. however, it is to
convoluted and missing the data file to specifically try to understand
what you are actually doing with your walls.

Meng

[...]

neighbor 0.2 bin

neigh_modify delay 10 every 1 check yes

hmm... don't these settings give you "dangerous builds"?
the skin distance seems very small and the delay in contrast very large.

[...]

# Define interaction parameters

pair_style hybrid airebo 3.0 1 1 lj/cut 10.0

# ------------- Coefficients ----------------

# AIREBO:

pair_coeff * * airebo CH.airebo C C

pair_coeff \{Ctype\} {CH4type} lj/cut 0.00555 3.6 10

pair_coeff \{CH4type\} {CH4type} lj/cut 0.01276 3.8 10

this is an inconsistent definition of the interactions. why don't you use:

pair_coeff * * airebo CH.airebo C NULL

[...]

# ------------- Boundary conditions ----------------

fix prlow lowerwall addforce 0.0 0.0 \${envforce}

fix prup upperwall addforce 0.0 0.0 -\${envforce}

wouldn't "fix aveforce" be the fix to be used here?

[...]

velocity lowerfix set 0.0 0.0 0.0

velocity upperfix set 0.0 0.0 NULL

fix 13 lowerfix setforce 0.0 0.0 0.0

fix 12 upperfix setforce 0.0 0.0 NULL

i don't understand this part at all. how can you have immobilized
atoms and mobile walls in the same input and they are all interacting?
that doesn't make any sense. furthermore, you are computing the
pressure you monitor for the entire system, even immobile atoms and
those that are supposed to generate the pressure. this is very
inconsistent, again.

axel.