viscosity calculation thermostat problem

Dear Lammps users,
Sorry to ask an old question. I have searched the email list but still can not get a good explanation.
I am currently trying to use the shearing the wall method to calculate the shear viscosity of water which is confined between two SAM-functionalized gold substrate. If I only use fix nve for the water, I can get a linear velocity profile from the simulation, but the temperature of my water increases a lot, which I think is not right.
The problem is when I try to use (fix nve + fix langevin) for the water, I could not get the linear velocity profile. The velocity points for my water is just randomly dispersed around 0.
I would be really grateful if you can give me some hints or explanation.

Below is my lammps script (fix nve + fix langevin for my water)

read_restart npt.restart
restart 1000 left.restart right.restart
pair_style hybrid lj/class2/coul/long 10.0 10.0 lj/cut/coul/long 10.0 10.0 morse 8
pair_modify mix arithmetic
kspace_style pppm 1.0e-5
bond_style hybrid class2 harmonic
angle_style hybrid class2 harmonic
dihedral_style class2
improper_style class2

include ‘pair_coeff’
include ‘bond_CH3’

group SAM1 id <> 1 1280
group SAM2 id <> 6081 7360
group Gold_1 id <> 1281 6080
group Gold_2 id <> 7361 12160
group Water id <> 12161 18160
group AuS type 3 5
group Au union Gold_1 Gold_2
group bulk1 union Gold_1 SAM1
group bulk2 union Gold_2 SAM2
group wall union bulk1 bulk2
group moving union Water bulk1

#--------------- eighborlist------------------------------------
neighbor 2.5 bin
neigh_modify delay 0 every 1 check yes
#--------------variable definition and initilization------------
variable vtimpstep equal 0.25
timestep {vtimpstep} variable Pdamp equal {vtimpstep}*1000
variable Tdamp equal ${vtimpstep}*100

variable shearrate equal 0.0015

variable atom2Pa equal 101325
variable fs2s equal 1.0e-15
variable ly equal 34
variable npt_run_time_1 equal 3000
variable npt_run_time_2 equal 3000
variable runtime equal 2000000
variable PaS2mPaS equal 1000

variable Text equal 300
variable Pext equal 1.0
variable length equal 38
fix constrain Water shake 1.0e-4 100 0 b 6 a 10

fix 2 Water nve

compute Water_thermo Water temp/profile 0 1 1 x 20

fix 1 Water langevin 300 300 ${Tdamp} 123456 tally yes
fix_modify 1 temp Water_thermo

compute 1 Water stress/atom Water_thermo
compute 11 Water reduce sum c_1[1]
compute 12 Water reduce sum c_1[2]
compute 13 Water reduce sum c_1[3]
compute 14 Water reduce sum c_1[4]
compute 15 Water reduce sum c_1[5]
compute 16 Water reduce sum c_1[6]

variable waterpress equal -(c_11+c_12+c_13)/(3lxly*38)

fix move bulk1 move linear ${shearrate} 0.0 0.0 units box

compute layers Water chunk/atom bin/1d z lower 0.02 units reduced
fix chunkavg Water ave/chunk 10 500 5000 layers vx c_1[4] c_1[5] c_1[6] file profile.wall

variable viscxz equal ((c_15*{atom2Pa}/(lx*ly*{length}))/({shearrate}/({length}*{fs2s})))*{PaS2mPaS}

fix visave all ave/time 10 100 1000 v_viscxz ave running file avg_vis.profile

thermo 1000
thermo_style custom step c_Water_thermo temp press v_waterpress pxy c_14 c_15 c_16 f_1 v_viscxz f_visave
thermo_modify temp Water_thermo

dump 2 all custom 4000 shearCH3.lammpstrj id type x y z vx vy vz fx fy fz

run 4000000

write_restart shear.restart


Hi Dezhao,

The issue with your script seems to be how you are applying shear to your system. By using fix move on ‘bulk1’ you effectively move this entire group as a single body ‘without regard to forces on the atoms’ and, from what I can tell, ‘bulk2’ is not time integrated at all (fix Langevin doesn’t do this). As a result, thermostatting these ‘immobile’ groups won’t work and hence you get a large temperature rise in the sheared region.

Its not surprising that you don’t get a linear velocity profile when using a stochastic Langevin themostat on the water layer including the streaming direction. In any case, it is usually not preferable to directly thermostat confined fluids NEMD simulations for several reasons, see e.g. J. Chem. Phys. 132 (24), 244706 .

You need to re-think how you apply shear to your system so that the wall regions you are thermostatting can actually dissipate energy - there are plenty of examples to follow both in the LAMMPS/examples directory and in the mailing list.

Kind regards,


Thank you very much for your kind explanation!
That helps a lot!