I am currently running MD simulations using the pGFN-FF model and have observed some unusual bending of CO2 molecules confined within a framework. Although not all molecules are affected, there are cases where the O-C-O angle drops to about 165°. The temperatures sampled are relatively low (up to 250 K), and this behavior is observed in both NVT and NPT simulations.
I have explored a few possible solutions without success. First, I attempted to treat the CO2 molecules as rigid, but in the version I am using (GULP 6.1.2), rigid molecules do not appear to be available for MD simulations.
I also specified the C=O bonds as double bonds using “bondtype C9 O9 double”, and tried imposing linearity with a lin3 potential. However, neither approach seems to affect the MD behavior. My impression is that these definitions may be ignored once the pGFN-FF keywords are activated.
After reviewing the available pGFN-FF options, I could not identify any other obvious way to constrain the O-C-O angle. Would you have any suggestions on how to prevent or reduce these deviations from linearity when using your pGFN-FF implementation?
Any advice would be greatly appreciated.
Best regards,
One of the possible reasons for what you observe is due to the nature of GFN-FF method and therefore pGFN-FF (which is largely just a periodic implementation of this). It’s important to remember that this approach defines a bonding topology based on the initial structure and this is never updated. This means although bonds can break, in principle, they can never form. If your initial structure happened to have CO2 located in an environment where it formed bonds other than the obvious C=O ones within the molecule, due to close proximity to the framework, then these molecules will remain parameterised as though bonded to the framework for the rest of the simulation. This means it’s best to ensure that the CO2 molecules are well separated from each other and the framework in the initial configuration & then it’s possible you won’t see significant non-linearity. Imposing additional potentials should be allowed in GULP, though I wouldn’t recommend this as it invalidates the parameterisation of the rest of the force field. Rigid molecule MD isn’t currently supported as GULP is primarily a lattice dynamics code rather than MD one. However, pGFN-FF is not intended for rigid molecule simulations - it’s a waste to use an advanced model that allows for flexibility and does all that extra calculation when you’re fixing all the intramolecular degrees of freedom.
Thank you for your quick response, very helpful as always. Your comment is probably pointing to the root cause, and I really appreciate your help in diagnosing what was going wrong.
I also completely agree that imposing rigidity in pGFN-FF would be a waste.