Hi!
I am trying to simulate confined LJ fluid in a slit pore. I model the slit pore with wall/lj126 interactions. To fill the pore I use GCMC simulation. However, I do something wrong and the simulation is not at the prescribed temperature. In fact the temperature drops suddenly from the prescribed 77 K to less than 1 K and hangs there, such that the fluid freezes.
I tried using the tfac_insert keyword to try to improve the control of temperature, but to no avail. Any advice on what may go wrong would be greatly appreciated. I use the April 8th 2021 version of LAMMPS.
Here is the code that I use:
# Confined LJ fluid in a slit pore with LJ walls
# Define global controls
units real
atom_style full
boundary p f p
# Define variables
variable ext_temp equal 77.0
variable mu equal -1.53
variable epsilon equal 0.7853
variable sigma equal 3.5750
variable cutoff equal 18.0
variable epsilonw equal 1.2250
variable sigmaw equal 3.1700
variable cutoffw equal 18.0
variable Lx equal 20.0
variable Ly equal 100.0
variable Lz equal 20.0
# Create box
region box block 0 ${Lx} 0 ${Ly} 0 ${Lz} units box
lattice sc 4.0
create_box 1 box
create_atoms 1 box subset 2 15021
# Define particle properties
mass 1 28.0134
group fluid type 1
# Define time-integration details
timestep 0.1
variable dt equal dt
variable disp index 0.5
variable Tdamp equal 100.0*${dt}
pair_style lj/cut ${cutoff}
pair_coeff * * ${epsilon} ${sigma} ${cutoff}
velocity fluid create ${ext_temp} 2349852 rot yes dist gaussian
# Define walls
fix WALLS all wall/lj126 ylo EDGE ${epsilonw} ${sigmaw} ${cutoffw} &
yhi EDGE ${epsilonw} ${sigmaw} ${cutoffw} units box
# Define NVT & GCMC commands
fix NVT fluid nvt temp ${ext_temp} ${ext_temp} ${Tdamp} tchain 3
fix_modify NVT dynamic/dof yes
region mcregion block 0.0 20.0 2.0 98.0 0.0 20.0 units box
fix mygcmc fluid gcmc 200 100 0 1 54341 ${ext_temp} ${mu} ${disp} &
region mcregion group fluid
# Define output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable racc equal f_mygcmc[8]/(f_mygcmc[7]+0.1)
compute_modify thermo_temp dynamic/dof yes
compute Tfluid fluid temp
compute_modify Tfluid dynamic/dof yes
variable Tfluid equal c_Tfluid
thermo_style custom step temp v_Tfluid density epair press pe ke etotal enthalpy atoms evdwl ecoul v_iacc v_dacc v_tacc v_racc
thermo 100
restart 5000 mus.restart1 mus.restart2
compute kea fluid ke/atom
compute pea fluid pe/atom
dump 11 fluid custom 1000 xmol id type x y z vx vy vz q c_kea c_pea element
dump_modify 11 sort id
# Run it!
run 1000000
Thank you!