Hello,

We are currently trying to learn how to use LAMMPS and we are having a begginer problem with the viscosity case in.wall.2d.

When we try to calculate the pressure manually (by reducing the result of compute atom/stress), we do not obtain the same results as the compute pressure command when the imposed shear rate on the top plate is not 0. However, when this shear rate is 0, the computed pressure is identical to the one we manually obtain. Furthermore, Pyy is always equal to the one we calculate manually.

What is the reason for this difference and how can we manully recompute the pressure tensor correctly when the shear rate is not 0?

Best regards,

Bruno

Joined LAMMPS input script :

# sample LAMMPS input script for viscosity of 2d LJ liquid

# use shearing wall, thermostat via fix langevin

# settings

variable x equal 20

variable y equal 20

variable ylo equal -2.5

variable yhi equal 23

variable rho equal 0.6

variable t equal 1.0

variable rc equal 2.5

variable srate equal 2.7

# problem setup

units lj

dimension 2

atom_style atomic

neigh_modify delay 0 every 1

lattice sq ${rho}

region simbox block 0 x {ylo} ${yhi} -0.1 0.1

create_box 3 simbox

create_atoms 1 box

mass * 1.0

pair_style lj/cut ${rc}

pair_coeff * * 1 1

region lower block INF INF INF 0.0 INF INF

region upper block INF INF $y INF INF INF

group lower region lower

group upper region upper

set group lower type 2

set group upper type 3

group wall union lower upper

group flow subtract all wall

mass 2 2.0

mass 3 2.0

velocity flow create t 97287 velocity upper set {srate} 0.0 0.0 units box

compute thermal flow temp/partial 0 1 0

compute flow flow temp

fix 1 all nve

fix 2 flow langevin $t $t 0.1 498094

fix_modify 2 temp thermal

fix 3 wall setforce 0.0 0.0 0.0

fix 4 flow ave/spatial 20 500 10000 y center 0.05 vx &

units reduced file profile.wall.2d

fix 5 all enforce2d

# equilibration run

variable ybox equal y*ylat compute tilt flow temp/ramp vx 0 {srate} y 0 ${ybox} units box

thermo 1000

thermo_style custom step temp c_tilt epair etotal press pxy

dump 1 all image 1000 image.*.jpg mass type zoom 1.6 adiam 1.2

run 3000

# data gathering run

compute atomstress all stress/atom

compute newpxx all reduce sum c_atomstress[1]

compute newpyy all reduce sum c_atomstress[2]

compute newpxy all reduce sum c_atomstress[4]

variable visc equal -pxy/(v_srate/ly)

variable newpxy equal -c_newpxy/(vol)

variable visc2 equal -v_newpxy/(v_srate/ly)

variable newpyy equal -c_newpyy/vol

variable newpxx equal -c_newpxx/vol

variable normal equal v_newpyy-v_newpxx

fix vave all ave/time 1000 1 1000 v_visc v_visc2 v_normal file visc.out

thermo_style custom step temp pxy v_newpxy pyy v_newpyy pxx v_newpxx v_visc v_visc2 v_normal

thermo_modify temp tilt

thermo 1

# only need to run for 5000 steps to make a good 100-frame movie

#dump 1 all custom 50 dump.wall.2d id type x y z vx

dump 2 all image 1000 image.*.jpg mass type zoom 1.6 adiam 1.2

#dump_modify 2 pad 5 amap 0.0 ${srate} ca 0.0 2 min blue max red

run 50