GCMC with CO2

Dear Aidan,

I am Tuan Ho from geochemistry department.

I performed gcmc simulation using lammps for supercritical CO2 into an empty box to test the model and methodology. The problem is that for a given chemical potential the density of CO2 and pressure are very very high. I tested the same CO2 model (TRaPPE) and chemical potential using TOWHEE and it gave me very good results (i.e., mu=-8Kcal/mol, T=338K will yield density of 0.6384g/ml and pressure 198atm). I also tested with MD NPT simulation using lammps. At 338K and 198atm it gave me the same density as TOWHEE did.

There is one warning because of using fix gcmc: WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (…/neighbor.cpp:461).Because the Coulombic interactions will affect the energy calculation I think it might affect the virial term in the pressure calculation. This might be the reason for the problem that I have (I guess).

equi1.lammps (1.21 KB)

out.data (18.4 KB)

CO2.txt (442 Bytes)

The fact that you are able to get the same density from TOWHEE GCMC and LAMMPS NPT is a good sign. It means that both codes are running the same model. You did not mention the density that you obtained using LAMMPS GCMC. It might just be an issue with how chemical potential is defined. It could also be an issue with fix rigid not working with fix gcmc. fix gcmc has special code dealing with fix shake. You might have better look if you run your model with SHAKE, or with a fully flexible CO2 model.

Aidan

Also, try doing your testing on the gas phase, it will go a lot quicker.

Dear Aidan,

Also, try doing your testing on the gas phase, it will go a lot quicker.

I try for the gas phase and it worked ok with the rigid or flexible CO2 GCMC LAMMPS. For the same chemical potential both LAMMPS GCMC and TOWHEE produced consistent results. Please see below for more on the supercritical CO2

The fact that you are able to get the same density from TOWHEE GCMC and LAMMPS NPT is a good sign. It means that both codes are running the same model. You did not mention the density that you obtained using LAMMPS GCMC.

The density I obtained for the same chemical potential in LAMMPS GCMC is ~0.9g/ml with pressure of ~2000 atm.

It might just be an issue with how chemical potential is defined.

I run for the gas phase for the same chemical potential LAMMPS GCMC and TOWHEE produced the same results. I don’t think it is the case.

It could also be an issue with fix rigid not working with fix gcmc. fix gcmc has special code dealing with fix shake.

I performed one simulation only used fix gcmc (without fix nvt, or fix rigid). This didn’t work for dense phase. Fix shake does not work for linear molecule.

You might have better look if you run your model with SHAKE, or with a fully flexible CO2 model.

I did tested with the fully flexible CO2 using only fix gcmc (without fix nvt or fix rigid). Again it worked for gas phase but not for dense phase. There are several reasons that make me believe that the insertion of molecule to the dense phase is the problem here. First, in some cases I received this warning in the middle of the run: WARNING: System is not charge neutral, net charge = 0.35 (…/kspace.cpp:297)… The warning indicates that it inserted C-O (C charge =0.70, O charge = -0.35), not O-C-O. Second, visualization of the snapshot, I found some CO2 molecules have C-O bond and angle reduce from 1.16 to 0.96A. angle =145. I understood that because it is flexible but in lammps description: Note that fix GCMC does not use configurational bias MC or any other kind of sampling of intramolecular degrees of freedom. Inserted molecules can have different orientations, but they will all have the same intramolecular configuration, which was specified in the molecule command input. Why is this the case?

I tried another approach: run the CO2 model with large bond and angle force constant. This did not work as well.

Thank you very much for your time.

Tuan

Dear Aidan,

Also, try doing your testing on the gas phase, it will go a lot quicker.

I try for the gas phase and it worked ok with the rigid or flexible CO2 GCMC LAMMPS. For the same chemical potential both LAMMPS GCMC and TOWHEE produced consistent results. Please see below for more on the supercritical CO2

The fact that you are able to get the same density from TOWHEE GCMC and LAMMPS NPT is a good sign. It means that both codes are running the same model. You did not mention the density that you obtained using LAMMPS GCMC.

The density I obtained for the same chemical potential in LAMMPS GCMC is ~0.9g/ml with pressure of ~2000 atm.

It might just be an issue with how chemical potential is defined.

I run for the gas phase for the same chemical potential LAMMPS GCMC and TOWHEE produced the same results. I don’t think it is the case.

It could also be an issue with fix rigid not working with fix gcmc. fix gcmc has special code dealing with fix shake.

I performed one simulation only used fix gcmc (without fix nvt, or fix rigid). This didn’t work for dense phase. Fix shake does not work for linear molecule.

You might have better look if you run your model with SHAKE, or with a fully flexible CO2 model.

I did tested with the fully flexible CO2 using only fix gcmc (without fix nvt or fix rigid). Again it worked for gas phase but not for dense phase. There are several reasons that make me believe that the insertion of molecule to the dense phase is the problem here. First, in some cases I received this warning in the middle of the run: WARNING: System is not charge neutral, net charge = 0.35 (…/kspace.cpp:297)… The warning indicates that it inserted C-O (C charge =0.70, O charge = -0.35), not O-C-O. Second, visualization of the snapshot, I found some CO2 molecules have C-O bond and angle reduce from 1.16 to 0.96A. angle =145. I understood that because it is flexible but in lammps description: Note that fix GCMC does not use configurational bias MC or any other kind of sampling of intramolecular degrees of freedom. Inserted molecules can have different orientations, but they will all have the same intramolecular configuration, which was specified in the molecule command input. Why is this the case?

I tried another approach: run the CO2 model with large bond and angle force constant. This did not work as well.

Thank you very much for your time.

Tuan

Hi Aidan,

Thank you for quick response.

  1. When using fix rigid and fix gcmc: bond angle and bond length does change (some time it looks like water molecule than CO2).

  2. Because in the previous email you suggested not to used fix rigid so I just used fix gcmc alone (because fix shake does not work for CO2 so there is no way I can check if fix gcmc alone (with shake keywords) can work with rigid model). BTW, there is one thing I don’t understand: fix gcmc doesn’t sample the intraconfiguration. Why do we need ‘shake’ keyword in the fix gcmc? Does the description imply what the code actually is?

  3. Now I strongly believe that the intramolecular configuration cause all the problems with dense phase. I don’t think TOWHEE takes angle and bond terms into account. That is why it yielded very good result.

So the problem now is how to keep molecule configuration as it is in the molecule file in fix gcmc. No bond and angle terms involes in the gcmc simulation.

Please let me know if you have any suggestion to go futher.

Tuan

Hi Aidan,

Thank you for quick response.

  1. When using fix rigid and fix gcmc: bond angle and bond length does change (some time it looks like water molecule than CO2).