I’ve dumped the velocities with a fix in a text file. The particles velocities do not change due to the commands trying to move the upper block, and when I apply the shake fix command, the velocities turn to nan. Can it be because temp/rescale doesn’t integrate in time? I’ve tried different velocities for the velocity set command, and all I get is that the total energy changes, and only this. I’m trying to apply Nose- Hoover but then I become an error message about out of range atoms. Am I using any command non-applicable to non-equilibrium simulations?
Here is the code:
"
dimension 3
units real
boundary p p f
atom_style full
kspace_style pppm 1.0e-5
kspace_modify slab 3
bond_style harmonic
angle_style harmonic
#Create atoms and boxes
lattice fcc 4.08
read_data water100Auv03.dat
region box block -1.765539 1.765539 -1.765539 1.765539 -3.000000 3.000000
#Pair_style hybrid definition
pair_style hybrid lj/cut/coul/long 10.0 eam lj/cut 2.5
#Definition of pair coefficients
#Water
pair_coeff 1 2 lj/cut/coul/long 0.0 0.0 10.0 #pair coefficients for OH
pair_coeff 1 1 lj/cut/coul/long 0.102 3.188 10.0 #pair coefficients between OO
pair_coeff 2 2 lj/cut/coul/long 0.0 0.0 10.0 #pair coefficients between HH
bond_coeff 1 450.000000 0.9572
angle_coeff 1 55.000000 104.52
#Gold
pair_coeff 3 3 eam Au_u3.eam #pair coefficient for planes
#Water+ planes
pair_coeff 1 3 lj/cut 1.0 1.0 2.5
pair_coeff 2 3 lj/cut 1.0 1.0 2.5
#running conditions
neighbor 0.3 bin
neigh_modify delay 5
#neigh_modify check yes
#delete overlaped atoms
delete_atoms overlap 0.1 all all
#create region
region upper block INF INF INF INF 1.8 INF
region lower block INF INF INF INF INF -1.8
region imob union 2 upper lower
group upper region upper
group lower region lower
group imob region imob
group box region box
group flow subtract box imob
set group imob type 3
#Dynamics
compute tv all temp/partial 0 1 1
velocity all create 350 482748 temp tv
variable K atom “mass[]*(vy[]^2)+(vz[]^2)/2”
fix 1 all ave/spatial 1 1 100 x 0.0 0.6 v_K x y z vx vy vz file water100Auv08.txt
#stabilizing T=300K
fix 2 all temp/rescale 10 350 300 100 5 0.5
fix_modify 2 temp tv
fix 3 flow shake 0.00001 200 10 b 1 a 1
thermo 100
run 1000
#moving it
velocity upper set 0.3 0 0 sum yes
#unfix 2
#fix 4 all nvt 300.0 300.0 1.0
write_restart restart.water100Auv08
dump 1 all atom 1000 dump.water100Auv08
dump 2 all xyz 1000 dump.water100Auv08.*.xyz
#timestep 50
run 5000
"
2009/11/30 Steve Plimpton <sjplimp@…24…>