Hello Axel,
Well, as you have told my script is just syntactically correct! Here is the condensed form of my script
I’m trying to siumulate water flow between graphene sheets.
My setup consists of 2 graphene sheets along each of the zboundaries (totally 4), and water molecules in between the channels. I’ve been using the p p f boundary conditions.
My input script works in the following way :

Initially set the regions and group them into walls and the liquid.

Assign the force fields :
pair_style hybrid lj/cut/coul/long 10.0 lj/cut 8.0
pair_coeff 1 1 lj/cut/coul/long 0.15535 3.166 # where 1 is O atom, 2 is H atom, 3 is C atom
pair_coeff 1 2 lj/cut/coul/long 0.03411 2.0845
pair_coeff 2 2 lj/cut/coul/long 0.0077 0.98
pair_coeff 1 3 lj/cut 0.3135 3.19
pair_coeff 2 3 lj/cut 0.253 3.28
pair_coeff 3 3 lj/cut 0.068443 3.407
2a) setting the bond and angle interactions, SHAKE algorithm
bond_style hybrid harmonic
bond_coeff 1 harmonic 1000 1.0
special_bonds lj/coul 0.0 0.0 0.5
angle_style hybrid harmonic
angle_coeff 1 harmonic 1000 109.47
group spce type 1 2
fix fshake spce shake 0.0001 10 100 b 1 a 1
kspace_style pppm 1e04
kspace_modify slab 3.0
pair_modify mix arithmetic
neighbor 2.0 bin
neigh_modify delay 1 every 1 check yes
neigh_modify exclude group upper upper
neigh_modify exclude group lower lower

I then minimize the energy of the system :
unfix fshake
minimize 1.0e4 1.0e6 100 1000
fix fshake spce shake 0.0001 10 100 b 1 a 1 
I use the following commands for the equilibration :
velocity all create $T 102486 mom yes rot yes dist gaussian # $T = 298K
fix NVT liquid nvt temp $T $T 10 drag 0.2
After equilibration, I use the following commands (temp/rescale) to set the temperatures on the graphene walls :
fix templwal1a lowera temp/rescale 10 298.0 298.0 0.02 1.0
fix templwal1b lowerb temp/rescale 10 298.0 298.0 0.02 1.0
fix templwal2ba uppera temp/rescale 10 298.0 298.0 0.02 1.0
fix templwal2bb upperb temp/rescale 10 298.0 298.0 0.02 1.0
fix liquid liquid temp/rescale 10 298.0 298.0 0.02 1.0
5)After this I try to calculate the thermal conductivity of water using the GK formula and the compute heat flux :
reset_timestep 0
compute myKE liquid ke/atom
compute myPE liquid pe/atom
compute myStress liquid stress/atom virial
compute flux liquid heat/flux myKE myPE myStress
variable Jx equal c_flux[1]/v_liquidvol
variable Jy equal c_flux[2]/v_liquidvol
variable Jz equal c_flux[3]/v_liquidvol
fix JJ liquid ave/correlate $s $p d &
c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running
variable scale1 equal {convert}/${kB}/$T/$T/v_liquidvol*s*{dt}
variable k11 equal trap(f_JJ[3])*v_scale1
variable k22 equal trap(f_JJ[4])*v_scale1
variable k33 equal trap(f_JJ[5])*v_scale1
variable tempk equal (v_k11+v_k22+v_k33)/3.0
thermo_style custom step temp v_liquidvol v_k11 v_k22 v_k33
 I then run the script file, but the results are unphysical. I get the thermal conductivity of water to be 10.35W/mK which is not correct.
Could you please have a glance at the above script, and suggest me how to go about?
time step for equilibration : 0.5fs., no.of steps for equilibration : 20000
for thermal conductivity : 1000000
Thanks and Regards,
Deepak