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