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

#read_restart pressuretest.rest

read_data graphene1776.conf

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

# ------------- 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

# ------------- 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}*lx*ly/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])/(3*lx*ly*(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