GCMC+NVT Density of flexible CO2 at dense phase

Dear all,

I am Liang Huang from University of California, Berkeley.

Our group is trying to validate the density of flexible CO2 model in bulk phase at high pressure with the hybrid GCMC +NVT method. The problem is that with flexile CO2 model at 300 atm and 338 K, the density we got from GCMC-NVT (~0.33 g/cm3) is far smaller than the NIST data (~0.81 g/cm3). The obtained temperature is around 338 K, but the averaged pressure is around 434 atm, much larger than 300 atm.

We have tested the flexible CO2 model with NPT MD simulation, which can give us consistent density with the NIST data. So, we believe the problem is in the GCMC part.

We are performing MPI parallel run with Kspace and tail yes (full energy option on). The LAMMPS version is 23 Oct 2017. We used pressure and fugacity coefficient, the fugacity coefficient is from the Peng -Robison equation.

We have looked through the email list and tried some possible solutions.

(1)Intra_energy option.

As suggested, we have added the intra_energy option. But there were no big differences using the intra_energy or not (it may because that we get a very small intra_energy for CO2). We got this value (-0.00556 Kcal/mol) for one CO2 molecule at zero step by performing NVE simulation with consistent parameters as GCMC+NVT. Hopefully, the way we got it is correct.

(2) Initial equilibrium

We have performed minimization and initial NVT equilibrium before the hybrid GCMC+NVT method. But we got no improvement on the CO2 density.

(3) Maniputing N, X, M paraments

We have set M to be 0, and tested both 100, 10, 0 and 10000 2500 0. The density results for both cases are still much smaller than targeted value.

(4)Full energy option no

We have made another attempt to close the full energy option by deleting Kspace and tail correction. Again, it didn’t work.

(5)Delete-bonds command

We have used the delete-bonds command (delete_bonds co2 multi remove) to close the intramolecular interaction in the template molecule. It didn’t work too.

(6)Pure GCMC

We have also tested the results of pure GCMC. Although it was said that pure GCMC in LAMMPS is not efficient. The density is increased to be 0.65 g/cm3, but still smaller than the target one (0.81 g/cm3).

(7)Force filed

We have tried several force fields. Even the commercial PCFF+ of Medea and COMPASS of Material Studio cannot give a good density with the hybrid GCMC+NVT method.

Now, we are totally lost. I have attached our input file, CO2 molecule file and Log file. Could you give us some advices on what we miss in our simulation? Any of your comment and suggestion will be highly appreciated.



co2exam.in (3.7 KB)

co2mol_gcmcnvt.dat (429 Bytes)

log.lammps (12.3 KB)

1 Like

please note, that there are known issues with the kind of simulation you are running when running in parallel. in more recent LAMMPS versions, this way of running is not allowed anymore for more than one MPI rank. thus you should update to the latest version to avoid computing incorrect trajectories even is the differences are likely small.

i don’t think, however, this issue of running in parallel is connected to the mismatch of density you are observing, and i am not enough of an expert in using fix gcmc to comment on that. i am copying Aidan, who is the resident LAMMPS expert for this.


Dear Axel,

Thanks so much for your comment.

We have noted the latest discussion on the issues with this kind of combination. Thus, we have set the translation/rotation parameter M in GCMC to be 0.

As you comment, we also don’t think this kind of combination can not even give a good CO2 density at bulk phase. Aidan has addressed the problem associated with GCMC+rigid/nvt/small for CO2 density at dense phase. The flexible model should be more mature than the rigid CO2 model. So, we believe that we must have missed something in our simulation. Hopefully, Aidan can give us some comments and help us out.



Dear all,

We have also tried to run our hybrid GCMC+NVT simulation with only one MPI rank. However, the obtained CO2 density at 300 atm and 338 K (~0.3 g/cm3) is still much smaller than the targeted one (~0.81 g/cm3).

Moreover, we have taken apart the hybrid simulation into iterative GCMC +MD+GCMC+MD….The density we got with parallel run is around 0.4 g/cm3.

We really need your help. What can we do to increase the CO2 density ? We have been stuck in the problem for such a long time.



CCing this to Aidan …