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 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
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}