Dear lammps users
I want to simulate the Couette flow. I have explicit two layer wall atoms at Z direction (one at the top, another at the bottom) and they are defined as a group. They are fixed at the initial position.There is another one group of atoms defined as flow and the angle and bonds of water molecules are fixed by “fix constrain”. I minimize my system first, and then create the velocity to the flow and then I equilibrate the flow . Finally I make the top wall to move to generate the Couette flow. When the Couette flow is generated, the thermostat is only applied to the flow direction perpendicular to the wall by new “fix nvt” command. Although the simulation can run at the equilibrium stage, the output the temperature,press and total energy output are “-nan”. I suppose there is something wrong with the thermostat, but I can not identify the problem. The following is the input file and output.
#nanoscale water flow squeezed between two substrate planes
units real
atom_style full
dimension 3
boundary p p s #not sure about the boundary
neighbor 2.0 bin
neigh_modify delay 5
read_data water_graphiteatom.data
#define each layer and water molecular as different groups, the sequence of wall follows the Z direction
group bot_wall id 1:864
group top_wall id 865:1728
group water molecule 2:1501
group wall subtract all water
#L-J potentials
pair_style hybrid lj/cut/coul/cut 15 lj/cut 15
pair_coeff 1 1 lj/cut 0.05570 3.4
pair_coeff 2 2 lj/cut/coul/cut 0.00000 0.000
pair_coeff 3 3 lj/cut/coul/cut 0.15535 3.166
pair_coeff 1 3 lj/cut 0.09302 3.283
pair_coeff 2 3 lj/cut/coul/cut 0.00000 0.000
pair_coeff 1 2 lj/cut 0.00000 0.000
bond_style harmonic
bond_coeff 1 554.1349 1.0
angle_style harmonic
angle_coeff 1 45.7696 109.47
#minimize the system
min_style sd
min_modify dmax 0.005
minimize 1.0e-4 1.0e-6 10000 1000
#fix the bond and angle of water moleculars
fix constrain water shake 1.0e-4 100 0 b 1 a 1
#initial velocities
velocity water create 300.0 104112 dist gaussian
fix NHER1 water nvt temp 300 300 100
#equlibrating flow
velocity bot_wall set 0.0 0.0 0.0 units box
velocity top_wall set 0.0 0.0 0.0 units box
fix fix_wall wall setforce 0.0 0.0 0.0
dump 1 all custom 1000 flow_solid_NHER_equi.lammpstrj mol type x y z ix iy iz mass element
dump_modify 1 element C H O
timestep 1
thermo 100
thermo_modify flush yes
run 300
#Couette flow
velocity bot_wall set 0.0 0.0 0.0 units box
velocity top_wall set 0.7 0.0 0.0 units box
fix still wall setforce 0.0 0.0 0.0
#thermostating
unfix NHER1
compute flow_temperature water temp/partial 0 0 1
fix NHER2 water nvt temp 300 300 100
fix_modify NHER2 temp flow_temperature
thermo_modify temp flow_temperature
#profile of the flow
compute chunky water chunk/atom bin/1d z lower 0.2
fix profile1 water ave/chunk 1 1000 1000 chunky vx density/number density/mass file Couette.profile Couette.profile
dump 2 all custom 10000 flow_solid_lng.lammpstrj mol type x y z ix iy iz mass element
dump_modify 2 element C H O
run_style verlet
run 300
output
Step Temp E_pair E_mol TotEng Press Volume
1000 190.33213 132942.61 0 135624.45 -nan 62591.886
1100 -nan 126709.3 0 -nan -nan 62591.886
1200 -nan 126709.3 0 -nan -nan 62591.886
1300 -nan 126709.3 0 -nan -nan 62591.886
Any help will be appreciated.
Thanks
Fan Li