GCMC of confined LJ fluid


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!

Your script contains a subtle error in how the thermostat is defined, the correct version is:

#fix_modify NVT dynamic/dof yes
compute_modify NVT_temp dynamic/dof yes

Note that this use case is demonstrated in …/examples/gcmc/in.gcmc.h2o


1 Like

Thank you, Aidan! Now the simulation works perfectly as intended.

~ Filip