Chemical potential and fugacity problems in fix_gcmc

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I want to conduct GCMC simulations to investigate methane adsorption on kerogen by using fix gcmc.

Before that I want to try the simulation of bulk methane using TraPPE force field.

I encountered some problems when calculating the density profile:

  1. When I specified Chemical potential (data from NIST database), the density after equilibrium is much higher than expected. Temperature continues increasing, higher than the specified value. Pressure is decreasing but far higher than pressure specified (Several orders of magnitude).

  2. When I specified Pressure and Fugacity coefficient (also from NIST database),

For pressure less than 10 bar, the density profile agrees well with NIST database;

When pressure higher than 10 bar, the density calculated is much lower than expected and the error gets higher with the increase of pressure. Temperature and pressure are both lower than specified. e.g. T specified is 300 K, p specified is 200 bar, fugacity coefficient is 0.75842, density expected is 0.155 g/cm3, while T calculated is around 160 K, p calculated is around 140 bar, density calculated is 0.084 g/cm3. Attached please find the results using Pressue and Fugacity coefficient.

  1. I guess the problem might be that the temperature was not thermostated at the specified value. So I added the fix_nvt command to the input file. The fix nvt was combined with fix_modify dynamic/dof yes and fix_modify temp as the LAMMPS Manuals says. However, an error was produced as follows: ERROR: Lost atoms: original 10036 current 9730 (…/thermo.cpp:442).

I am not familiar with gcmc at LAMMPS. I have looked through the LAMMPS Mailing List and solved many problems but I haven’t found the answers to the above three problems.

The version of LAMMPS I am using is lammps-17Nov16. I have also tried the version of lammps-11Aug17. But the problem remains unsolved and the speed of calculation slows down and I haven’t found the reason.

Any help would be helpful.

Thank you very much!

Best regards,

Juan

#####################INPUT SCRIPT############################

units real

boundary p p p

atom_style atomic

region mybox block -60.0 60.0 -60.0 60.0 -60.0 60.0

create_box 1 mybox

create_atoms 1 random 10000 3576738 mybox

group mygas type 1

mass 1 16.0425

pair_style lj/cut 10.0

pair_coeff 1 1 0.148 3.730

pair_modify tail yes

velocity mygas create 300.0 4928459 mom yes rot yes dist gaussian

variable chempot equal -1.201599293

variable fugacoeff equal 0.75842

variable pressatm equal 200*0.9869233

#compute mdtemp mygas temp

#compute_modify mdtemp dynamic/dof yes

#fix mdnvt mygas nvt temp 300.0 300.0 100.0

#fix_modify mdnvt temp mdtemp

#fix 2 mygas gcmc 10 1000 1000 1 3456543 300.0 {chempot} 0.1 pressure {pressatm} fugacity_coeff ${fugacoeff}

fix 2 mygas gcmc 10 1000 1000 1 3456543 300.0 ${chempot} 0.1

timestep 1

thermo 500

thermo_style custom step temp press density pe ke etotal

run 200000

results using Pressue and Fugacity coefficient.jpg

2

  1. Chemical potential always includes an arbitrary constant. The LAMMPS choice is explained on the doc page. There is no reason to expect LAMMPS and NIST to use the same constant.

  2. If the temperature is different, of course the pressure will be different

  3. I suggest following the example examples/gcmc/in.gcmc.lj. No dynamics is used and the temperature comes out correctly. If you do dynamics, you have to be more careful. See the H2O and CO2 examples, which both use the comand:

compute_modify thermo_temp dynamic/dof yes

Aidan