Charge Parameter Discrepancies between LigParGen Server and LigParGen v2.1

[shehan807]

Hi,

I'm creating a new issue similar to #9--using LigParGen (LigParGen v2.1, BOSS v5.0) and comparing it to the online server, I get a case where the charges are different in the resulting OpenMM XML files. For example, if I run...

ligpargen -n TMP -p TMP -r TMP -c +1 -o 3 -cgen CM1A -s 'C[P+](C)(C)C'

the <NonBondedForce> (and <AtomTypes> section, for reference) from the TMP.openmm.xml looks like

<!--       Generated with LigParGen   --> 
<!--        William L. Jorgensen Lab  -->
<!--     Author: israel.cabezadevaca@yale.edu  -->
<!--    OPLS Force Field with CM1A derived Atomic Charges  -->
<ForceField>
<AtomTypes>
<Type name="opls_800" class="P800" element="P" mass="30.974"/>
<Type name="opls_801" class="C801" element="C" mass="12.011"/>
<Type name="opls_802" class="C802" element="C" mass="12.011"/>
<Type name="opls_803" class="C803" element="C" mass="12.011"/>
<Type name="opls_804" class="C804" element="C" mass="12.011"/>
<Type name="opls_805" class="H805" element="H" mass="1.008"/>
<Type name="opls_806" class="H806" element="H" mass="1.008"/>
<Type name="opls_807" class="H807" element="H" mass="1.008"/>
<Type name="opls_808" class="H808" element="H" mass="1.008"/>
<Type name="opls_809" class="H809" element="H" mass="1.008"/>
<Type name="opls_810" class="H810" element="H" mass="1.008"/>
<Type name="opls_811" class="H811" element="H" mass="1.008"/>
<Type name="opls_812" class="H812" element="H" mass="1.008"/>
<Type name="opls_813" class="H813" element="H" mass="1.008"/>
<Type name="opls_814" class="H814" element="H" mass="1.008"/>
<Type name="opls_815" class="H815" element="H" mass="1.008"/>
<Type name="opls_816" class="H816" element="H" mass="1.008"/>
</AtomTypes>

...

<NonbondedForce coulomb14scale="0.5" lj14scale="0.5">
<Atom type="opls_800" charge="2.128122" sigma="0.374000" epsilon="0.836800"/>
<Atom type="opls_801" charge="-0.726714" sigma="0.350000" epsilon="0.276144"/>
<Atom type="opls_802" charge="-0.758699" sigma="0.350000" epsilon="0.276144"/>
<Atom type="opls_803" charge="-0.752600" sigma="0.350000" epsilon="0.276144"/>
<Atom type="opls_804" charge="-0.775626" sigma="0.350000" epsilon="0.276144"/>
<Atom type="opls_805" charge="0.152006" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_806" charge="0.152006" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_807" charge="0.152006" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_808" charge="0.157288" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_809" charge="0.157288" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_810" charge="0.157288" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_811" charge="0.156856" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_812" charge="0.156856" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_813" charge="0.156856" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_814" charge="0.162355" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_815" charge="0.162355" sigma="0.250000" epsilon="0.125520"/>
<Atom type="opls_816" charge="0.162357" sigma="0.250000" epsilon="0.125520"/>
</NonbondedForce>

whereas if I input the exact same smiles string (and optimization and charge model settings) into the server, the resulting xml file looks like

<ForceField>
<AtomTypes>
<Type name="opls_811" class="H811" element="H" mass="1.008000" />
<Type name="opls_810" class="H810" element="H" mass="1.008000" />
<Type name="opls_812" class="H812" element="H" mass="1.008000" />
<Type name="opls_816" class="H816" element="H" mass="1.008000" />
<Type name="opls_813" class="H813" element="H" mass="1.008000" />
<Type name="opls_804" class="C804" element="C" mass="12.011000" />
<Type name="opls_809" class="H809" element="H" mass="1.008000" />
<Type name="opls_803" class="C803" element="C" mass="12.011000" />
<Type name="opls_801" class="P801" element="P" mass="30.973762" />
<Type name="opls_814" class="H814" element="H" mass="1.008000" />
<Type name="opls_808" class="H808" element="H" mass="1.008000" />
<Type name="opls_815" class="H815" element="H" mass="1.008000" />
<Type name="opls_807" class="H807" element="H" mass="1.008000" />
<Type name="opls_800" class="C800" element="C" mass="12.011000" />
<Type name="opls_806" class="H806" element="H" mass="1.008000" />
<Type name="opls_802" class="C802" element="C" mass="12.011000" />
<Type name="opls_805" class="H805" element="H" mass="1.008000" />
</AtomTypes>

...

<NonbondedForce coulomb14scale="0.5" lj14scale="0.5">
<Atom type="opls_806" charge="0.162200" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_804" charge="-0.870400" sigma="0.350000" epsilon="0.276144" />
<Atom type="opls_810" charge="0.162200" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_815" charge="0.162300" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_808" charge="0.162200" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_803" charge="-0.870400" sigma="0.350000" epsilon="0.276144" />
<Atom type="opls_816" charge="0.162300" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_802" charge="-0.870500" sigma="0.350000" epsilon="0.276144" />
<Atom type="opls_800" charge="-0.870500" sigma="0.350000" epsilon="0.276144" />
<Atom type="opls_812" charge="0.162300" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_811" charge="0.162300" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_813" charge="0.162300" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_814" charge="0.162300" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_809" charge="0.162200" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_801" charge="2.535000" sigma="0.374000" epsilon="0.836800" />
<Atom type="opls_805" charge="0.162200" sigma="0.250000" epsilon="0.125520" />
<Atom type="opls_807" charge="0.162200" sigma="0.250000" epsilon="0.125520" />
</NonbondedForce>

You can see that for the phosphorus atom, the charge differs by ~0.4 e, 2.535500 from the server vs 2.128122 from the code. For the case of C[P+](C)(C)C, I did additionally get a runtime warning UFFTYPER: Unrecognized charge state for atom: 1, which seems to be an rdkit issue. This isn't necessarily the source of the inconsistency because I also tried testing C[N+](C)(C)C which suppresses the warning but results in inconsistent xml files. It would be very helpful to know if (1) this is reproducible for others and (2) where the discrepancy is coming from, thanks!

[iGulitch]

Have you tried locally with BOSS v.5.1?

[shehan807]

For BOSS v.5.1 (April 2024 Release), the output appears to be the same as BOSS v.5.0 for the case tested above, i.e., my .XML file still looks different than the server output.

Is the server using an older version of BOSS? If, say I use a error-prone "DUMMY" SMILES string as an input for the server, I get the following error:

 export BOSSdir=/var/www/html/ligpargen/apps/boss-4.9;/var/www/html/ligpargen/apps/anaconda2/bin/python2.7 /var/www/html/ligpargen/apps/ligpargenCode/Converter.py -s 'DUMMY' -r UNK -o 0 -c 0 > /tmp/errorServer.log 

This seems to suggest the LigParGen server is currently using BOSS v.4.9. This might explain the discrepancy since v.4.9 (Released Jan. 2013) would not include the parameter changes in OPLS/2020.

[Isra3l]

Hey,

Yes, the server was initially released in 2017 and uses BOSS4.9. The GitHub version of LigParGen has since been updated and rewritten from scratch to address compatibility issues with newer Linux distributions and to align with the latest BOSS version.

[shehan807]

Is there any comment on which parameters are "correct" and/or suggested use given the discrepancy, i.e., should users revert to the GitHub version when possible? Thanks so much!!

[Isra3l]

Hey Shehan,

I recommend using this version, as it includes the latest OPLS updates. There is not a “correct” or “incorrect” version -both have been reliable over time, though the latest version should ideally offer improved parameters. To my knowledge, there are no plans to update the server. In any case, for many molecules, the parameters will remain the same.

By the way, have you added the BOSS executable to the Docker container? I am unsure if that is allowed under the BOSS license. I won’t be able to accept any commit unless I am completely certain it complies with the software license.

Best

[shehan807]

I see, in that case, I won't use the BOSS executable through Docker to err on the side of caution. I'm not so sure if there is a work around at the moment, however.