Brownian dynamics simulation for planar shear

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 coefficientshear ratey_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


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} []( 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}

The trick in these kind of models is to get the fluid flowing
in the same manner the box is deforming. To do this you
typically need to thermostat with a bias applied to
the thermostat fix that matches the deformation.
So use compute temp/deform as the thermostat
temperature compute. See the howto section
of the manual for details about how to apply
bias to temperature computes and assign
them to thermostats. Then you can use
another compute temp (e.g. profile) to monitor
the temperature with the streaming profile
substratcted. Also you should visualize your
system to verify that the fluid is flowing in sync
with the derfoming box. Sometimes that takes
a while to come to equilibrium. You can also
initialize the fluid with a velocity profile that matches
the box deformation rate. See the velocity command
for how to do that, e.g. the ramp option.