Error in MOF + CO2 system with gcmc and fix rigid/nvt/small

LAMMPS LAMMPS Beginners
,
1

Hello,

I’m trying to simulate a MOF + CO2 system using gcmc. For the co2 molecules im using the fix rigid/nvt/small. The problem is that when a molecule is inserted in the system (sometimes in the thermo output the number of accepted mc insertions is still zero), I get errors such as “Out of range atoms - cannot compute PPPM”. I understand this could be because the molecules are inserted too close to the MOF structure and the forces/energies are too high.

I tried using overlap_cutoff, but I get the warning:

Energy of old configuration in fix gcmc is > MAXENERGYTEST

if the cutoff is too big (more than 2 angstroms), and it just makes the co2 amount drop (I guess it just rejects every insertion attempt). Otherwise, if smaller than that, it doesnt solve the problem.

Im using the 17 Apr 2024 version (I also tried with 2 Aug 2023 - Update 3). I tried with different MOF structures (with bigger pore size) and force fields (Dreiding and UFF, with parameters given by the lammps_interface software), using rigid/nve/small with Langevin thermostat, smaller timesteps and tfac values… neither worked.

Here is a small version of the input script, not using a fix on the mof (rigid structure). In this particular case, two molecules are added to the system (the first one is a good example of why the system can be breaking, it “pops out” when inserted, good to see in ovito) before it breaks. This is the best case I could get to, otherwise not a single molecule would be inserted.

######################   SYSTEM CREATION   ######################
units           real
atom_style      full
boundary        p p p

pair_style      lj/cut/coul/long 12
bond_style      harmonic
angle_style     hybrid cosine/squared harmonic
dihedral_style  charmm
improper_style  umbrella

special_bonds   dreiding
dielectric      1.0
box tilt        large
read_data       thermo.data

neigh_modify    every 1 delay 0 check yes
kspace_style    pppm 1e-6
pair_modify     tail yes mix arithmetic

molecule co2_mol CO2.txt

######################   SIMULATION   ######################

group           co2     type 6 7
compute         co2T co2 temp
compute_modify  co2T dynamic/dof yes

dump dp1 co2 atom 100 dump.lammpstrj
dump dp2 all atom 100 all.lammpstrj

fix     co2nvt co2 rigid/nvt/small molecule temp 293 293 100 mol co2_mol
fix_modify co2nvt temp co2T

fix     co2gcmc co2 gcmc 200 5 0 0 5342 293 -5 0.5 mol co2_mol full_energy tfac_insert 1.6 group co2 rigid co2nvt overlap_cutoff 1.5

variable	iacc equal f_co2gcmc[4]/(f_co2gcmc[3]+0.1)
variable	dacc equal f_co2gcmc[6]/(f_co2gcmc[5]+0.1)

thermo_style    custom step atoms c_co2T v_iacc v_dacc
thermo 200
timestep 0.1

run 100000

The files for CO2.txt and thermo.data are in my drive (sorry I cant upload files here yet).

Any suggestion on how to improve the input? I’m still learning how to use LAMMPS, so maybe I’m missing something. Thanks a lot.

PS: This problem disappears when I use the regular nvt fix (flexible co2) and works perfectly fine, even without the overlap keyword, I don’t completely understand why. But I would like to use the rigid model if possible, that’s why im asking for advise here.
Also, sorry for opening the topic and closing, I thought I fixed it for a moment and I wanted to check.

2

Update:
First, something I didnt mention before, I had already used fix rigid/nvt/small with co2 in bulk (no mof structure) and gcmc, it worked perfectly. It was adding the framework that makes the problem arise. Either the structure being flexible or rigid.

Now, I could make a “bad” solution to the problem. Instead of using both fixes at the same time, which looks like its what is producing the problem with the mof (if its only one or the other there is no problem), now there is one run for the rigid/nvt/small and another for the gcmc. This way it works, for some reason.

Here is a version of the input, the thermo.data and CO2.txt files are the same than before:

######################   SYSTEM CREATION   #################################
units           real
atom_style      full
boundary        p p p

pair_style      lj/cut/coul/long 12
bond_style      harmonic
angle_style     hybrid cosine/squared harmonic
dihedral_style  charmm
improper_style  umbrella

special_bonds   dreiding
dielectric      1.0
box tilt        large
read_data       thermo.data
neigh_modify    every 1 delay 0 check yes

kspace_style    pppm 1e-6
pair_modify     tail yes mix arithmetic


#******* molecules ******
molecule co2_mol CO2.txt

############################   SIMULATION   #################################
group           mof     type 1 2 3 4 5
group           co2     type 6 7

compute         co2T co2 temp
compute_modify  co2T dynamic/dof yes
compute         mofT mof temp

variable N_co2  equal count(co2)/3

dump dco2 co2 atom 100 dump_co2.lammpstrj
dump dall all atom 100 dump_all.lammpstrj

timestep 0.1

label loop
    variable a loop 200
    ## NVT run
    fix     co2nvt co2 rigid/nvt/small molecule temp 298 298 100 mol co2_mol
    fix_modify co2nvt temp co2T

    thermo_style    custom step v_N_co2 c_co2T c_mofT
    thermo 500
    run 9900
    unfix co2nvt

    ## GCMC run
    fix     co2gcmc co2 gcmc 1 5 0 0 1263 298 -5 0.5 mol co2_mol full_energy tfac_insert 1.6 group co2 overlap_cutoff 1.5

    variable	iacc equal f_co2gcmc[4]/(f_co2gcmc[3]+0.1)
    variable	dacc equal f_co2gcmc[6]/(f_co2gcmc[5]+0.1)

    thermo_style    custom step v_N_co2 c_co2T c_mofT v_iacc v_dacc
    thermo 10

    run 100
    unfix co2gcmc

    if "${a} == 200" then "jump SELF mainrun"
    next a
    jump SELF loop


label mainrun

I have barely used the jump command so it might not be applied the best way. Neither have I found the best values for how long each of the nvt/gcmc runs should be for a good convergence (I only tried a few different values for now), so there are still things to improve, but it seems to work well. I tried with both rigid and flexible mof structures as well.

For now the only really bad thing is that, due to some warnings that repeat over time, after a while they get “omitted” for being too much of them… that could be dangerous.

I don’t really understand why this is making the script work, given that in theory is doing something very similar than using both fixes at the same time, some gcmc attempts every “n” nvt steps…

3

I may have a suggestion, but I am not sure it solves the problem (I have never done an gcmc for rigid molecules).

So, in general (outside of gcmc context), I have learnt that if you use fix rigid without turning off the interactions inside the rigid body, there can be some bogus calculation of pressure (and maybe energy?), which I thought may be causing the issue, since you dont observe it when you consider flexible CO2 molecules. I suppose that in this case, since the intra-molecular non-bonded interactions should be excluded thanks to the special bonds command, the issue may be the bonded interactions (I see that you have set very large values of K, which I understand should not really matter in principle since the molecule is rigid anyways).

Maybe you can try setting these values of K to 0 in bond and angle interactions to see if it solves the problem? Maybe you can even explicitly use commands to turn off the non-bonded and bonded interactions within CO2 anyways. Something like:

delete_bonds	co2 multi
neigh_modify	exclude molecule/intra co2
4

The large values for K are for simulating almost “rigid” molecules when using the flexible model (I thought as well that it shouldn’t affect the simulation when using then fix/rigid).

So, I tried the things you proposed and it didnt work neither. But I think the same as you, that the problem (and its solution) may be related to the energy of the system…

Thanks a lot for the help!

1 Like
5

Hey

Can you provide full your working version of this script please?

I am a student and my task is to obtain the isotherms of carbon dioxide adsorption in the MOF-5 structure, and your code would help me a lot!