LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] viscosity calculation thermostat problem
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] viscosity calculation thermostat problem

From: James Ewen <j.ewen14@...68...>
Date: Fri, 20 Apr 2018 11:57:07 +0100

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,


On Thu, Apr 19, 2018 at 10:17 PM, Dezhao Huang <dhuang2@...46...> wrote:
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)/(3*lx*ly*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


Dezhao Huang
2nd year Ph.D. student
University of Notre Dame
Department of Aerospace and Mechanical engineering
150D CoMSEL Fitzpatrick Hall
Notre Dame, IN 46556-5710

Check out the vibrant tech community on one of the world's most
engaging tech sites,!
lammps-users mailing list