Hi

I am trying to set up a simulation to study rheology of a system particles. I would like to use *implicit solvent and planar shear*. I don’t want to use SLLOD as it would require solvent explicitly.

To this end, I use the following fixes:

Fix nve

Fix langevin

Fix deform

To add a linear streaming velocity profile v_x(y), I add a force (x-component only) equal to

f_x = drag coefficient*shear rate*y_coordinate, to all the particles.

I use *compute temp/profile* so that the bias added due to the imposed streaming velocity can be subtracted while calculating thermal properties.

However, the temperature output from thermo_style is off from the fixed temperature. It is about twice!

Then I tried to use *fix_modify* *langevin* to remove bias using *temp/profile*, now the difference was even more!

I also find that the *damp*_parameter used in the Langevin thermostat affect the viscosity measured as

-pxy/shear_rate, and the steady state temperature!

I would appreciate if anyone could tell me what is going wrong here. I attach below the script.

Best regards

Ethaya

# settings

variable x equal 8

variable y equal 8

variable z equal 8

variable rho equal 0.8442

variable t equal 0.725

variable rc equal 2.5

variable rco equal 2.51

variable xyrate equal 1.0 # strain rate

variable sdrag equal 1.0

#variable damp equal 1./{sdrag}
variable damp equal 0.1
variable eqstep equal 4200
variable sstep equal 4200
variable accstep equal {eqstep}+{sstep}
variable nestep equal 42000
variable vxstep equal 4200
variable vistep equal 4200
variable dumrun equal 2/{xyrate}/0.001

variable dumstep equal 0.02/${xyrate}/0.001

# problem setup

units lj

dimension 3

atom_style atomic

neigh_modify delay 0 every 1

# problem setup

lattice fcc ${rho}

region simbox prism 0 $x 0 $y 0 $z 0.0 0.0 0.0

create_box 1 simbox

create_atoms 1 box

pair_style lj/smooth {rc} {rco}

pair_coeff * * 1 1

mass * 1.0

velocity all create $t 97287

#velocity all ramp vx -6.7184 6.7184 y 0.0 13.4368

timestep 0.001

fix 1 all nve

fix 2 all langevin $t t {damp} 498094

# equilibration run

thermo 1000

run ${eqstep}

unfix 1

unfix 2

# turn on NEMD shear and equilibrate some more

#velocity all scale $t

# shear rate defined relative to perpendicular dimension

timestep 0.001

variable dforce atom {sdrag}*{xyrate}*y

compute def all temp/profile 1 1 1 xyz 10 10 10

fix 1 all nve

fix 2 all langevin $t t {damp} 4749574

#fix_modify 2 temp def

fix 3 all addforce v_dforce 0. 0. every 1

fix 4 all deform 1 xy erate ${xyrate} remap v

compute layers all chunk/atom bin/1d y center 0.05 units reduced

fix 5 all ave/chunk 100 10 {vxstep} layers vx file profile.nemd.2d
#compute prof all temp/deform
thermo_style custom step temp epair etotal press pxy
thermo_modify temp def
thermo {vxstep}

run ${sstep}

# data gathering run

thermo {vistep}
variable visc equal -pxy/{xyrate}

fix vave all ave/time 100 10 {vistep} v_visc ave running start {accstep}

thermo_style custom step temp press pxy v_visc f_vave

thermo_modify temp def

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

# set 54K steps above in equil and 5400 here simply to make good movie

# 54K and 5400 are multiples of box-swap periodicity = 2700 steps

run ${nestep}

#run for visualization

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

dump 1 all custom {dumstep} [dump.xyz](http://dump.xyz) id type x y z
dump_modify 1 sort id
#dump 2 all image 50 image.*.jpg vx type zoom 1.2 adiam 1.2
#dump_modify 2 pad 5 amap 0.0 {srate} ca 0.0 2 min blue max red

run ${dumrun}