Openff-toolkit: Antechamber/sqm fails on about 1% of input

[j-wags]

In the first 300 molecules of our test_molecule database, antechamber/sqm has failed three times. These are used in calculating partial charges using the bcc and mul charge methods.

ZINC00335972
ZINC35326424
ZINC00160396

Both failures had C#N motifs. However, C#N motifs are not unique to the failure, as ZINC28286664 (which has one) was processed successfully. This bug may have something to do with the C#N motifs being connected to an aromatic system.

Larger and smaller molecules were processed successfully, as were singly and doubly charged molecules.

OE QuacPac did not have a problem calculating partial charges for the same molecules.

The error message in sqm.out was similar in each case:

QMMM: ERROR! QMMM: Unable to achieve self consistency to the tolerances specified QMMM: No convergence in SCF after 1000 steps. QMMM: E = -0.4629E+06 DeltaE = 0.5704E-08 DeltaP = 0.9306E-05 QMMM: Smallest DeltaE = 0.3434E-08 DeltaP = 0.1597E-04 Step = 986

Possible solutions:

• Check the protonation/geometry by eye to see if the problem is in file conversion
• Loosen the convergence criteria or increase the number of iterations allowed in sqm

[j-wags]

For now, I've added a check for this error in the compute_partial_charges function and moved on. However, this might be a good issue to tackle sooner rather than later.

[jchodera]

Another useful thing would be to report these failures to the Amber mailing list to make sure this issue is on their radar.

[jchodera]

Tagging @dgasmith who might have a recommended replacement for sqm for AM1.

[davidlmobley]

That's interesting. I've never run into this and I've used antechamber for charging a lot. BUT, I've also never done molecules with Antechamber which had C#N motifs connected to an aromatic system that I can think of.

[dgasmith]

Looks like MOPAC has AM1: http://openmopac.net/manual/semiempirical_theory.html

[jchodera]

Is mopac conda-installable and available without a license agreement @dgasmith?

[dgasmith]

Unfortunately no, only program I could think of however.

[j-wags]

Did some searching through the mailing list, and found some options that we could vary to see if they help us get convergence.

http://archive.ambermd.org/201704/0379.html

One keyword I'll try is ndiis_attempts=700

[j-wags]

I relaxed (?) many settings in the SQM run for ZINC00160396, but still didn't get a non-error run. (I'm not a QM person, so maybe this was nonsensical)
Original

Run semi-empirical minimization
 &qmmm
    qm_theory='AM1', grms_tol=0.0005,
  scfconv=1.d-10,   qmcharge=0,
 /

Relaxed (?)

Run semi-empirical minimization
 &qmmm
    qm_theory='AM1', grms_tol=0.002, maxcyc=2000,
  scfconv=1.d-8,   qmcharge=0,   ndiis_attempts=1000
 /

A few more failures in the first 1300 molecules from our test set. Most have nitrile moieties.

Nitriles:
ZINC13130634
ZINC28114756
ZINC10022957
ZINC14143999
ZINC24528707
ZINC01615096 (not aromatic!)

No nitrile, but has a net charge of -3 (I think it's ok to fail on these extreme cases)
ZINC22013758

[j-wags]

Per John's suggestion, I'll submit this to the AMBER mailing list.

@dgasmith (and others?), do we have plans to use Psi4 as a partial charge / Wiberg bond order calculation method in the future? The interface we currently have with antechamber goes through the CLI and relies on temporary files, which seems unstable to me. It might be good to start thinking of an alternative.

[dgasmith]

Psi4 can do partial charges / Wiberg bond orders now, but will be much slower than SQM/MOPAC.

[davidlmobley]

I think we should be OK with slow, @dgasmith and @j-wags . And I really don't think it should be Jeff's priority to deal with problems with Antechamber/sqm. My perspective is that we should:
a) allow people to provide their own charges
b) support OE's good charging package as a fast and good solution
c) allow them to use antechamber
d) support any other reliable and open options which are recommended to us ( psi4 we obviously should support)

but given that folks can do a and b, it's not our job to ensure antechamber is as robust/reliable. Our funding isn't to create open source semiempirical QM codes, nor to optimize existing ones. There are good fast options (e.g. OE Quacpac) and if people don't want to pay for them then they can help us by bringing forward alternatives.

[jchodera]

I disagree on this, @davidlmobley. Anything that affects force field predictions must be part of the force field definition, otherwise accuracy will be highly variable for users. We need an open, reliable way to assign partial charges. We will ultimately need to pair with a quantum chemist who can deliver this.

[davidlmobley]

I agree with you on this, @jchodera :

We will ultimately need to pair with a quantum chemist who can deliver this.

Key words there being "ultimately" and "quantum chemist", i.e., not necessarily right now, and not us. :)

It's not Jeff's job to deal with this. And, whether we invest Consortium resources in trying to create a good open quantum solution would be something we'd have to sort out with the board.

That said, after discussions yesterday I think it seems likely Michael Schauperl will end up creating a good open solution for this, probably in collaboration with Lee-Ping's group and with support from Jeff.

[jchodera]

It's not Jeff's job to deal with this. And, whether we invest Consortium resources in trying to create a good open quantum solution would be something we'd have to sort out with the board.

Totally agree---right now, we just offer an open source alternative (even if it fails some fraction of the time) and move on until we get the resources to provide a more robust solution. But I think retreating from the principle that the force field uniquely encodes everything required to compute the energies would be a major mistake.

[davidlmobley]

Yes, agreed.

[j-wags]

I've followed this conversation and won't consider this to be an immediate priority for the near future. Thanks for clearing this up.

For anyone who takes over this task later, my email to the AMBER users list is here: http://archive.ambermd.org/201811/0010.html

One response suggests that the wavefunction might be divergent. I'm not an expert on quantum chemistry, but I interpreted this to mean that the nuclei are in a place where the SDF's defined bond orders/electron sharing don't make sense. Judging by a quick Pymol view of the situation, there isn't anything horribly wrong with the atom geometry or bond orders.

David Case also responded that increasing ndiis_attempts to 700 fixed the problem on his machine. However, this didn't fix it for me, and sqm still throws the same non-convergence error when I set it all the way up to 1000.

My next guesses would be

• Try it in AmberTools 18. Right now, I'm using AmberTools 16 from the ambermini conda repo.
• See if something messed up the molecule when it was read into/out of the openforcefield molecule format.

[j-wags]

As a follow-up: This is resolved in AmberTools 18 and newer. So if people run into this again, the solution is conda create -n ambertools -c ambermd ambertools (which currently pulls down AmberTools 19)

[jchodera]

From Ian Craig:

Just as reported here (https://github.com/openforcefield/openforcefield/issues/130#issuecomment-531276506) we're also finding that sqm fails to converge in a small-single-digit-% of cases which are often nitriles. Since moving to the sqm in ambertools 19 fixes the issue (as described in issue #130), would it be possible to either: (i) update ambermini to match ambertools19, or (ii) have ambertools 19 installed along with openforcefield instead of ambermini?

[jchodera]

We do have an AmberTools 19 (ambertools 19.9) package built via conda-forge here, but it's currently only available for linux:
https://anaconda.org/conda-forge/ambertools/files
Here's the feedstock:
https://github.com/conda-forge/ambertools-feedstock

The current issue is that it doesn't yet build on osx for unknown reasons. There is an open PR to add this, so it may be valuable to see if we can get this to work, since it would drastically simplify our ability to stay up to date with AmberTools: When AmberTools releases a new patch update, all that's required is to bump this line and the build logic will automatically apply the appropriate patch level, where the versioning scheme here is {major_release}.{patch_level}.

[jchodera]

As a stopgap fix, we could also copy the official AmberTools 19.0 release from the AmberMD channel using anaconda copy, but this makes it harder to keep up to date with new version-level or patch-level releases.

[davidlmobley]

Note that we also have an issue now -- #473 -- where charging with antechamber fails with NEWER versions of antechamber for carboxylates. It seems we will have to suffer one problem or another.

I will check in with MolSSI as to whether they have any updates on alternative open charging solutions.

[j-wags]

Can Omnia handle channel designations in recipes (like ambermd::ambertools)? Or is that a special conda-forge syntax?

If we can't force a channel designation, I agree with @jchodera that our best current option is to copy the latest AT package onto Omnia. I can coordinate updating that with changing the CLI behavior in the OFFTK AmberToolsToolkitWrapper too.

Once AmberTools is up on conda-forge, we can remove the stopgap.

Taking the full AmberTools package will also close #153, which is the decision we were trending towards anyway. So, two birds with one stone!

[jchodera]

Can Omnia handle channel designations in recipes (like ambermd::ambertools)? Or is that a special conda-forge syntax?

Is this supported by conda-forge in recipes? Can you point to some docs?

[dgasmith]

The :: is not supported by cf by choice. It is supported by env files, and I think build files too.

[jchodera]

Can someone point me to the docs for what you're suggesting? Is this a new anaconda feature that was recently introduced to allow channels to be specified for build dependencies in recipes?

[jaimergp]

I am afraid the only official docs is the class docstring. The new syntax was introduced in conda 4.6, I think.

[jaimergp]

I don't think you can use it in recipes yet, though. See here: https://github.com/conda/conda-build/issues/532