Tutorial 04 - Optimizing Force Fields Using Partial Charges

[faranaknejat1991-bot]

Hello:

I'm using ubuntu 22.04. I installed OpenFF using "conda install -c conda-forge openff-evaluator".
I have been able to run a successful run from the given example (CCO) and I was able to optimize the Lj parameters.

I am also interested to be able to optimize the partial charges. This is where the issue comes from. I have introduced the "cosmetic_attribute" for the partial charges. The program runs and changes the partial charges, however, since the partial charges do not add up to zero after being changed, the program crashes.

I have looked up the SMIRNOFF format to do this but still the program crashes. I think there has to be a way to constrain the partial charges to be zero after they are being changed. And I don't think this is explained fully in here: "https://openforcefield.github.io/standards/standards/smirnoff/#partial-charge-and-electrostatics-models".

If I use this format for the Force Field which I have specified the parametrize key word inside the Library of Charges:

The Open Force Field Initiative 2025-06-18

Then error message is:
openff.interchange.exceptions.NonintegralMoleculeChargeError: Molecue CCO has a net charge of -0.1315656 compared to a total formal charge of 0.0.

Now, I have looked online and realize the best way to constrain the total charge (Neutrality) is to use the keyword "parameter_eval" under the library of charge tag. The way I a have approached this is to say that the charge of the end hydrogen will be the sum of the negative of the other charges. Then the force field is:

The Open Force Field Initiative 2025-06-18

Then the error message is:
openff.toolkit.utils.exceptions.IncompatibleUnitError: charge_increment=-0.05 dimensionless should have units of elementary_charge.

This error message does not make sense since the unit of the "parameter_eval" is automatically correct.
Please let me know if you need anything else. I have search everywhere and can't find a solution!
I have attached the full env as well.

full-env.yaml

Thanks,
Ab

[faranaknejat1991-bot]

In case you can't see the forcefield, i have them attached here.
The one that say "partial" is for the test that parametrize keyword is defined inside the Library of charges and the other file is foe which I gave used the key word parameter_eval to constrain the charges.

openff-ethanol.txt
openff-ethanol_Partial.txt

[j-wags]

Hi @faranaknejat1991-bot, I'm away from my computer so I can't dig into this, but two recommendations:

• chargeincrementmodel will probably be better than librarycharges here. There's a feature where, if there's one fewer chargeincrement than tagged atoms in the parameter smirks, the last one will be defined as the negative of the sum of the rest. This may be easier than defining all the chargeincrement parameters separately and ensuring they sum to zero using them via parametrize_eval (though you'll still need to use parametrize_eval to ensure symmetric Hs get the same charge)
• using am1bccelf10 as the base charge model is probably unnecessary, since it seems like the goal is to just find optimal charges for your molecule (so starting from zeros would be faster and end up at the same result). You can define partial_charge_method=zeros in the chargeincrementmodel header to get this behavior.

The nonintergralcharge/librarycharges errors were probably due to exactly what the error message said, and I suspect the answer to that is doing something with parameterize_eval. I'm not sure about the unit errors.

[faranaknejat1991-bot]

Hello, Thanks for your advice. I was able to get the program working using the following syntax to define the "Atom" and "ChargeIncrement" tags:

<vdW version="0.4" potential="Lennard-Jones-12-6" combining_rules="Lorentz-Berthelot" scale12="0.0" scale13="0.0" scale14="0.5" scale15="1.0" cutoff="9.0 * angstrom ** 1" switch_width="1.0 * angstrom ** 1" periodic_method="cutoff" nonperiodic_method="no-cutoff">
    <Atom smirks="[#1:1]-[#6X4H3]" epsilon="0.0157 * mole ** -1 * kilocalorie" id="n1" rmin_half="1.546191 * angstrom" charge="0.0447599987188975 * elementary_charge"></Atom>
    <Atom smirks="[#1:1]-[#6X4H2]" epsilon="0.0157 * mole ** -1 * kilocalorie" id="n2" rmin_half="1.546191 * angstrom" charge="0.0173199983934561 * elementary_charge"></Atom>
    <Atom smirks="[#1:1]-[#8X2]" epsilon="0.0157 * mole ** -1 * kilocalorie" id="n3" rmin_half="1.546191 * angstrom" charge="0.39809000367919606 * elementary_charge"></Atom>
    <Atom smirks="[#6X4H3:1]" epsilon="0.1094 * mole ** -1 * kilocalorie" id="n4" rmin_half="1.908 * angstrom" charge="-0.09709999834497769 * elementary_charge"></Atom>
    <Atom smirks="[#6X4H2:1]" epsilon="0.1094 * mole ** -1 * kilocalorie" id="n5" rmin_half="1.908 * angstrom" charge="0.13142999882499376 * elementary_charge"></Atom>
    <Atom smirks="[#8X2:1]" epsilon="0.21 * mole ** -1 * kilocalorie" id="n6" rmin_half="1.6612 * angstrom" charge="-0.6013399971028169 * elementary_charge"></Atom>
</vdW>
<Electrostatics version="0.4" scale12="0.0" scale13="0.0" scale14="0.833333" scale15="1.0" cutoff="9.0 * angstrom ** 1" switch_width="0.0 * angstrom ** 1" periodic_potential="Ewald3D-ConductingBoundary" nonperiodic_potential="Coulomb" exception_potential="Coulomb"></Electrostatics>
<ChargeIncrementModel version="0.4" number_of_conformers="1" partial_charge_method="zero">
    <ChargeIncrement smirks="[#1:1]-[#6X4H2:2]" charge_increment1="0.01 * elementary_charge" parameterize="charge_increment1" />
    <ChargeIncrement smirks="[#1:1]-[#6X4H3:2]" charge_increment1="0.01 * elementary_charge" parameterize="charge_increment1" />
    <ChargeIncrement smirks="[#6X4H3:1]-[#6X4H2:2]" charge_increment1="0.01 * elementary_charge" parameterize="charge_increment1" />
    <ChargeIncrement smirks="[#6X4H2:1]-[#8X2:2]" charge_increment1="0.01 * elementary_charge" parameterize="charge_increment1" />
    <ChargeIncrement smirks="[#8X2:1]#[#1X1:2]" charge_increment1="0.01 * elementary_charge" parameterize="charge_increment1"/>
</ChargeIncrementModel>

As you can see in the above I have specified the partial charges inside the "atom" tags and set partial_charge_method="zero".

This method works however, it is ignoring the initial starting points for the charges for which I get from Gaussian. This is important because I would like to setup a threshold for which partial charges can not change more than that.

So my questions are:

1- Is it possible to enforce the initial starting point for the partial charges and some sort of a constraint?
2- The (charge_increment1="0.01*elementary_charge" ) is the threshold/constraint?

[mattwthompson]

I don't think Evaluator has ever been used to train partial charges, so you're treading into new territory here and there's no guarantee it will work at all.

Also note that most OpenFF infrastructure (and force field science in general) is typically more interesting in covering a range of chemistry than trying to optimize parameters for a single molecule. I think this is part of why Jeff is suggesting you use charge increments; they're much more general than library charges. If all you care about it ethanol and you need to carefully set the charges on the zeroth iteration, I think library charges are better since you can convert a molecule's partial charges into a set of library charge parameters.

I'm not sure how you'd enforce that the library charge parameters can only change by a certain amount (or necessarily that it would be wise to do so). That would be set in the ForceBalance configuration, however, and the ForceBalance docs are a better resource for those options than here.