Hello everyone,
I am parameterizing a protonated molecule for molecular dynamics simulations using the CHARMM36/CGenFF force field.
My workflow was as follows:
- I used the MP2/cc-pVDZ optimized geometry of the molecule reported a literature.
- I performed a single-point B3LYP/6-311++G(d,p) calculation with Pop=CHELPG in Gaussian-09 to obtain ESP-derived atomic charges.
- I generated a .mol2 file from the optimized .xyz geometry using Open Babel.
- I replaced the Open Babel charges with my Gaussian CHELPG charges.
- I submitted the molecule to CGenFF (ParamChem).
However, I noticed that the partial charges assigned by CGenFF are quite different from my CHELPG charges, particularly for several aromatic carbons and the protonated ammonium group.The param pentaly is coming 0.9 and charge pentalty is 0.65.
My questions are:
- Is it expected that CGenFF will assign different charges even if the input MOL2 already contains QM-derived CHELPG charges?
- Does ParamChem ignore the charges present in the input MOL2 and assign its own CHARMM-compatible charges?
- For simulations with the CHARMM36 force field, should I use the charges assigned by CGenFF, or is it acceptable to manually replace them with my CHELPG charges?
- If I want to use my own CHELPG charges, would I also need to re-optimize other force-field parameters (bonded terms, Lennard-Jones parameters, etc.), or is replacing only the charges considered acceptable?
I would appreciate any advice from researchers who have experience with CGenFF parameterization or QM-derived charge fitting.
Hi Somdatta,
I’ll briefly answer all of your questions below:
- Yes. It is absolutely expected that CGenFF will assign different partial charges to your ligand’s atoms than input .mol2 partial charges. The assigned partial charges by CGenFF correspond with each atom’s “assumed” analogous assignment in the CHARMM general force field; once the parameter generator parses your throughput molecule, you’re explicitly telling the program to refine your throughput to hypothetically best-match CHARMM’s parameter set, NOT QM data.
- As described above, yes and no. It respects the charges in that they are leveraged when parsing your throughput molecule while initially matching component ligand atoms with analogous atoms within the CHARMM general force field, however, these charges are NOT necessarily compatible with the CHARMM force fields, and are rightfully discarded.
- Absolutely do NOT ever under any circumstances replace CGenFF-fitted partial charges with your own CHELPG charges. Please read the CHARMM36m and CGenFF methodology publications to understand why.
- You cannot use your own CHELPG charges. Raw QM-derived CHELPG charges are not fine-tuned/calibrated to accurately recapitulate a molecule’s behaviour with the CHARMM36m force field ensemble. Please read the CHARMM36m and CGenFF methodology publications to understand why.
At the moment, you need to pursue force field parameter refinement. You cannot blindly and manually override the CGenFF assignments because you will essentially be replacing charge data that is optimized to at least work with the other components of your MD simulation (critically water), with charge data that is almost guaranteed to lead to unphysical ligand behaviour in your simulation. The most rigorous and well-described tools for validating and optimizing poor CGenFF parameters, including partial charges and electrostatics, are “FFParam” (from the MacKerell group) and “ffTK” (from a third-party group). I recommend reading into both of these options, then deciding which is most applicable for your molecular system and your technical expertise. I’ve used both, and they are quite comparable with respect to workflow, rigour, and being user-friendly. If you’re just looking to optimize parameters for use with CHARMM36m; it’s just a matter whether you prefer a Python (FFParam) or TCL-VMD (ffTK) GUI/API
Here is also the DOI corresponding to the FFParam v2.0 article, which makes for a great read to enrich your understanding of parameter optimization, while also outlining why rigorous parameter optimization (not replacement) is required to maintain the validity of your simulations: doi/10.1021/acs.jpcb.4c01314.
But perhaps most importantly though, are you sure you need to optimize your parameters at all??? Both of the penalty scores you’ve explicitly mentioned here are <1.0, which are actually very great and STRONGLY suggest that CGenFF was able to find a very reliable analogous parameter set to describe the atoms in your molecule. Provided the CHELPG charges do not illustrate an entirely different charge distribution than your CGenFF-assigned charges, based on what you’ve said I don’t immediately see any cause for concern or need for parameter optimization whatsoever… Again, you should NOT see a magnitude match between QM-derived CHELPG charges and CGenFF partial charges; these are NOT 1:1 quantitatively comparable. Deviations between the two charge sets are both expected and even good signs that your molecule has been parameterized to behave accurately within the CHARMM force field ensemble.
Best of luck,
Kyle
Thank you, Kyle, for your detailed response.
Your explanation has really helped clarify my doubts. Since I’m quite new to force-field parameterization, I was hesitant about directly using the CGenFF-assigned charges for my simulations and was worried that I might need to replace them with my QM-derived CHELPG charges.
Your explanation of why CGenFF assigns its own optimized partial charges, and why manually replacing them would be inappropriate, has given me much more confidence in the parameterization process. Based on your advice, I’ll proceed with the CGenFF-generated files for my MD simulations.
I’ll also take a closer look at FFParam and ffTK, as you suggested, to better understand the parameter optimization workflow and the underlying methodology. I think that will help me build a stronger understanding of force-field parameterization.
Thank you again for taking the time to explain everything so clearly. I really appreciate your help.
Regards,
Somdatta