Pressure keyword is not working correctly in GCMC with reaxFF force field

LAMMPS LAMMPS General Discussion
1

Hi all,
I am currently running a GCMC +MD simulation to do co2 adsorption using ReaxFF force field (file attached below) in lammps version: LAMMPS (29 Aug 2024 - Update 1). I want to run the isotherms under different pressure for co2. I test with pure co2 in empty box firstly. I insert co2 into box at pressure 5 atm and T =300K with GCMC and do MD at the same time with mainly command lines:

units real

variable tfac equal 5.0/3.0
fix 3 gas gcmc 2000 100 0 0 1219 300.0 0.0 0.5 mol co2 region simbox pressure 5.0 overlap_cutoff 1.0 group gas tfac_insert ${tfac}
fix 4 gas langevin 300.0 300.0 50 562489
fix 5 gas nve

There is a problem that when I set up P values directly to replace chemical potential in gcmc, like when I set up like 5 atm, the inserted co2 numbers will be about 10000, which is much larger than the theoretical co2 number ( theoretical is about 100 co2 molecules based on my box size and T = 300K), and the pressure in my box will also be pretty large, which is about 2000 atm. While when I set up chemical potential −108.42 Kcal/mol directly (corresponding to 5atm), it will generate the co2 number and density pretty closer to the ideal gas theoretical values. The results align with the ideal gas and GCMC works well.

And I also have run the example gcmc co2 simulation script that is given in lammps source example, which is using LJ potential. It works pretty well under pressure value and chemicial potential keywords, both got CO2 molecules number and density close to the theoretical values of ideal gas.

So I was wondering if there is problem for my force field apply or there is a problem of P to chemical potential conversion with reaxFF in gcmc. If you have any idea about what is happening or any suggestions, please leave a message,
great appreciation ahead! :grinning:
ffield.reax.new (17.8 KB)

best,
Li

2

Perhaps @athomps has some suggestion(s) for you.

3

The ‘pressure’ keyword provides an alternative way to specify the chemical potential of the ideal gas reservoir, as described on the doc page. For sufficiently small specified pressure (and correspondingly large box size), the ratio of the pressure in the real gas to the specified pressure will approach a constant, which is supposed to be unity. We do not have any examples for the pressure keyword, but I believe it was only designed and tested for the LJ potential, so it makes sense that LJ worked for you as expected. It is possible that an extra factor was omitted in the code that matters when you switch from lj units to real units, or you switch from lj/cut to reaxff, or you add mol, region, or tfac_insert keyword. Can you do some more testing and find which things ‘break’ the pressure option?

1 Like
4

Thanks for your reply. For the pressure keyword is not working correctly for CO2 model under reaxFF force field. Under P from 1 to 5atm, 300k, using ReaxFF force field, CO2 is not behaving like the ideal gas, I need to consider it as real gas, so if just set up P directly in my gcmc simulation, the chemical potential it yields to the gas reservoir is actually not matching real co2 chemical potentials values needed at the same pressure values. The fugacity keywords should taken into account, which results smaller chemical potential values than the ideal gas.

Also, pressure keyword is working well with N2 gas at the same force field, T and P conditions. since N2 is behaving as ideal gas.

Thanks for you hints and help!