Problem calculating stress manually in simple viscosity case


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,

Joined LAMMPS input script :

sample LAMMPS input script for viscosity of 2d LJ liquid

use shearing wall, thermostat via fix langevin


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

thermo_modify temp tilt

I imagine the difference is due to this line,

and which affects the kinetic energy term in the stress/pressure.

With that line, you are telling pxx in the thermo command
to subtract out the tilt velocity before computing the KE term
in the stress.

But there is no way currently to subtract that out from the per-atom
stress. So they become different for non-zero stress.

So either don’t use the tilt temperature in the global pressure, or don’t include
the kinetic energy term in either the global or per-atom stress/pressure
and I’m guessing they will again agree.


The 12Feb14 patch adds a temperature compute to the compute stress/atom
command, same as compute pressure. So it should now be possible
to adjust the KE contribution to the stress/atom to get it to
match the global pressure, i.e. using the same temperature compute
with a bias in both cases.