-
Notifications
You must be signed in to change notification settings - Fork 93
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add ChargeIncrementModelHandler #471
Conversation
The-SMIRNOFF-force-field-format.md
Outdated
</ChargeIncrementModel> | ||
``` | ||
The sum of formal charges for the molecule or fragment will be used to determine the total charge the molecule or fragment will possess. | ||
|
||
`<ChargeIncrementModel>` provides several optional attributes to control its behavior: | ||
* The `number_of_conformers` attribute (default: `"10"`) is used to specify how many conformers will be generated for the molecule (or capped fragment) prior to charging. | ||
* The `quantum_chemical_method` attribute (default: `"AM1"`) is used to specify the quantum chemical method applied to the molecule or capped fragment. | ||
* The `partial_charge_method` attribute (default: `"CM2"`) is used to specify how uncorrected partial charges are to be generated from the quantum chemical wavefunction. Later additions will add restrained electrostatic potential fitting (RESP) capabilities. | ||
* The `partial_charge_method` attribute (default: `"Mulliken"`) is used to specify how uncorrected partial charges are to be generated from the quantum chemical wavefunction. Later additions will add restrained electrostatic potential fitting (RESP) capabilities. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In our last discussion, I think we wanted to combine quantum_chemical_method
and partial_charge_method
to instead have a general two-stage procedure:
- Stage 1: Produce initial charges via some QM or ML method: AM1-Mulliken, graph convolutional network, VCharge, etc.
- Stage 2: Correct the initial charges with BCCs
That will require more significant revision of the proposed ChargeIncrementModel spec.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok, I thought we had been trending toward this decision, but I couldn't find it in writing anywhere. I'm going to have partial_charge_method
be the winner here (so, for example, partial_charge_method="AM1-Mulliken"
). For the ToolkitRegistry, each toolkit will have a supports_charge_method(method)
function which takes a method
string and returns a bool indicating whether it can calculate charges using that method.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The partial charge method may itself have other parameters, so it might be useful to split this into two separate tags for the different stages.
When we had discussed this, I think we discussed generalizing it into perhaps two tags: an initial charging tag (that might take optional parameters) and a BondChargeIncrement tag (that describes how to modify charges).
Happy to chat further if you want to discuss ideas for making this flexible and simple to implement.
This pull request introduces 1 alert when merging 643a750 into e3a8856 - view on LGTM.com new alerts:
|
This pull request introduces 1 alert and fixes 1 when merging 3d3c9ef into 3de67ca - view on LGTM.com new alerts:
fixed alerts:
|
Reviewing! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Great work on this @j-wags! It's a very well-put-together implementation, with a lot of well-thought-out details.
Most of my inline comments are minor. Some additional notes below:
-
Should we add an explicit test for
ToolkitWrapper._check_n_conformers
? In fact, that method could be made into astaticmethod
, since it really doesn't depend onself
at all. A test of this should check that all the logic flow control works as expected. -
When are we planning to deprecate
ToolkitAM1BCC
? -
What do we expect to happen if the
ChargeIncrementModel
section in a SMIRNOFF spec hasnumber_of_conformers=0
? Does our machinery currently handle this case, as we'd expect for thegasteiger
method? -
Remove the warning in the SMIRNOFF spec that
ChargeIncrementModel
is not yet implemented.
if partial_charge_method == "am1bccnosymspt": | ||
optimize = False | ||
symmetrize = False | ||
quacpac_status = oequacpac.OEAssignCharges(oemol, charge_method['oe_charge_method'](optimize, symmetrize)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can you add a small comment on why this partial charge method requires additional parameters?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Added
# The OpenFF toolkit has always supported a version of AM1BCC with no geometry optimization
# or symmetry correction. So we include this keyword to provide a special configuration of quacpac
# if requested.
result = toolkit_registry.call('from_object', other, allow_undefined_stereo=allow_undefined_stereo) | ||
except NotImplementedError: | ||
pass | ||
result = toolkit_registry.call('from_object', |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What is the intention behind this block? It's not clear to me.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As in, what is the behavior we want to support here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh, I thought about this so much, I didn't write down what the successful logic means. D'oh!
Added
# Each ToolkitWrapper may provide a from_object method, which turns some particular type(s)
# of object into OFFMols. For example, RDKitToolkitWrapper's from_object method will
# return an OFFMol if provided with an RDMol, or raise a ValueError if it is provided
# an OEMol (or anything else). This makes the assumption that any non-ValueError errors raised
# by the toolkit _really are_ bad and should be raised immediately, which may be a bad assumption.
top = Topology.from_molecules([create_ethanol(), create_reversed_ethanol()]) | ||
sys = ff.create_openmm_system(top) | ||
nonbonded_force = [force for force in sys.getForces() if isinstance(force, openmm.NonbondedForce)][0] | ||
expected_charges = [0.2, -0.15, -0.05, 0., 0., 0., 0., 0., 0., |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd like to walk through this one with you, just to make sure I'm thinking about it correctly @j-wags.
assert -1.e-5 < charge_sum.value_in_unit(unit.elementary_charge) < 1.e-5 | ||
|
||
@pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') | ||
@pytest.mark.parametrize("partial_charge_method", ['am1bcc', 'am1-mulliken']) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is there a reason we aren't including gasteiger
here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes -- Gasteiger isn't conformer-dependent, so this test (ensuring that charge methods generate different charges for different conformations) doesn't apply to it. Added
Skip Gasteiger because it isn't conformer dependent.
to the test docstring
|
||
@pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), | ||
reason='RDKitToolkit and AmberToolsToolkit not available') | ||
@pytest.mark.parametrize("partial_charge_method", ['am1bcc', 'am1-mulliken']) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Same question as above: why no gasteiger
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is it because gasteiger
doesn't require conformers at all?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yup. Added the same text to the docstring as above.
This pull request introduces 1 alert and fixes 1 when merging 982fa49 into d0318bf - view on LGTM.com new alerts:
fixed alerts:
|
Co-authored-by: David Dotson <[email protected]>
This pull request introduces 1 alert and fixes 1 when merging d7e092a into d0318bf - view on LGTM.com new alerts:
fixed alerts:
|
This pull request introduces 2 alerts and fixes 1 when merging fff6051 into d0318bf - view on LGTM.com new alerts:
fixed alerts:
|
Worked on the inline comments today. I'll pick back up on your "overall" comments tomorrow! |
Regarding the tests, good idea. I've added them Regarding converting to staticmethod, I don't think it's a good long-term idea. We're going to convert ToolkitWrappers to have a cache of previously-performed operations, which would forbid us from using staticmethods and classmethods (since the only place a static or classmethod could store a cache is on the class itself, which is undesirable). For that reason, we're moving toward making all ToolkitWrapper methods into instance methods.
Since it's in the already-published Parsley and smirnoff99Frosst FFs, never.
I've updated the
Good catch. Removed. This should address everything you pointed out. I'll be leading the 0.7.0 release, so let's work through #572 first and I'll resolve any conflicts with this once it's merged. |
This pull request introduces 1 alert and fixes 1 when merging e363a70 into d0318bf - view on LGTM.com new alerts:
fixed alerts:
|
This pull request introduces 1 alert and fixes 1 when merging e910488 into a2f1ff3 - view on LGTM.com new alerts:
fixed alerts:
|
This pull request introduces 1 alert and fixes 1 when merging efe98bb into 9553e9c - view on LGTM.com new alerts:
fixed alerts:
|
raise_exception_types
kwarg toToolkitRegistry.call
, enabling it to optionally intercept or pass through certain exceptions. This is needed in this PR, as many charge-calculation methods will be accessed through theassign_partial_charges(partial_charge_method)
method in different toolkits, and each toolkit needs to be able signal to theToolkitRegistry
that a charge method is unrecognized but that other toolkits might support it.AmberToolsToolkitWrapper.assign_partial_charges(mol, partial_charge_method)
for supported methodsOpenEyeToolkitWrapper.assign_partial_charges(mol, partial_charge_method)
for supported methodsBuiltInToolkitWrapper
class, withassign_partial_charges(mol, partial_charge_method)
supportingzeroes
andformal_to_partial_charge
compute_partial_charges_am1bcc
seems Unable to assign charges for some ligands #492temporary_cd
back because it turns out we actually do use it (Drop temporary_directory/temporary_cd for tempfile #595)