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