Add initial cyclodextrin files

[andrrizzi]

Here are the initial host-guest files.

Note that this works only with the dev YANK currently in branch solvent-parameters.

Could you take a look at the cyclodextrin.yaml script and see if there's something that we should change in the protocol? In particular:

1. Which atom should we use to compute the centroids to restrain? For now, I went with the carbon atoms of the ligand and the receptor.
2. Should we use HMR? Any particular time step?
3. I've parametrized everything with AM1-BCC/GAFF (the host as well) rather than using the original RESP charges in the mol2 files because in Niel Henriksen and Michael Gilson's paper it seems to have given the best statistics overall for binding free energies.
4. The script is configured to assign OpenEye AM1-BCC charges to guests and antechamber bcc charges to hosts. This is because OpenEye on the host gave me this omega error:

   Warning: force field setup failed due to missing parameters for molecule [dodecaholmiooxy-pentakis(holmiooxymethyl)-BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2}-dodecaosmaheptacyclo[BLAH]dotetracontanyl]methoxyholmium

@davidlmobley this script should give you a good starting point to setup your calculations. Tweak the script as needed. It's still untested (I'll run 5-10 iterations of all the combinations in the next couple of days), but the setup stage should work on the YANK version I mentioned at the top. Note that the online analysis (options.online_analysis_target_error) and the automatic selection of the alchemical path (protocols.binding-auto.complex.alchemical_path=auto) are very recent features, and we haven't run simulations with these activated yet, so you may want to just specify a finite number of iterations (options.number_of_iterations) and your own alchemical path for now.

If you just want to create and inspect the system files, run this

from yank.experiment import ExperimentBuilder

exp_builder = ExperimentBuilder('cyclodextrin.yaml')
exp_builder.setup_experiments()

[jchodera]

Could you take a look at the cyclodextring.yaml script and see if there's something that we should change in the protocol? In particular: Which atom should we use to compute the centroids to restrain? For now, I went with the carbon atoms of the ligand and the receptor.

I'd just use all heavy atoms: (mass > 1.5). I'm surprised mdtraj doesn't have selections for heavy or not hydrogen.

Should we use HMR? Any particular time step?

I think we're still awaiting the results of @maxentile's benchmarking here to provide guidance. We might add a host-guest system to openmmtools.testsystems to help him with system-specific benchmarking to provide specific guidance for future calculations.

For now, we should use BAOAB with 2 fs timestep and standard hydrogen masses.

I've parametrized everything with AM1-BCC/GAFF (the host as well) rather than using the original RESP charges in the mol2 files because in Niel Henriksen and Michael Gilson's paper it seems to have given the best statistics overall for binding free energies.

We have to be careful of the goals here. Tagging @davidlmobley to be aware of this.

• For purposes of comparing the performance of different codes, we need to produce reference calculations with the exact same parameters and charges. That would mean we either have to use the RESP charges Niel's mol2 files or replace Niel's mol2 files with yours. I'd prefer to replace his files, but we should make sure that @davidlmobley is OK with this.
• For purposes of comparing with experiment, we should use a consistent approach to parameterizing and charging our molecules, ideally with what we believe to be our current "best practices". That means:
  • The latest GAFF 1.x from AmberTools 17 and/or GAFF 2.0. (Do we intend to do both?)
  • The canonical AM1-BCC charging scheme
  • For proteins, AMBER99SB-ILDN

The script is configured to assign OpenEye AM1-BCC charges to guests and antechamber bcc charges to hosts.

That's almost assuredly not what we want to do here.

This is because OpenEye on the host gave me this omega error: Warning: force field setup failed due to missing parameters for molecule [dodecaholmiooxy-pentakis(holmiooxymethyl)-BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2},BLAH$l^{2}-dodecaosmaheptacyclo[BLAH]dotetracontanyl]methoxyholmium

Have you reported this error to OpenEye? What did they say?

@davidlmobley this script should give you a good starting point to setup your calculations.

This should be fine to start @davidlmobley's comparison of restraint convergence statistics, but ideally we would standardize the parameterized systems in the benchmark set paper for comparing convergence rates of different approaches with different restraints.

[jchodera]

ionic_strength: 50*millimolar

Is 50 mM the experimental ionic strength?

@andrrizzi : Can you be sure to toss a README.md in there with information on the provenance of these files, the protocol you used to set them up, and the journal reference for experimental data and conditions?

[davidlmobley]

I'm going to be tied up for the next few hours, so I can just comment really quickly that, since the purpose here is to benchmark agreement (and efficiency) versus existing gold standard data from Niel, and NOT so much to compare against experiment, then you would presumably want to run with exactly the same charges he is using rather than reassigning charges.

@nhenriksen -- can you comment on why you're depositing the charges you are in the files for the benchmark sets? There is some discussion here about what charges we should be using. One answer is that we should be using whatever you deposited since that's going to be the "standard" benchmark, but then the question is, "But why didn't you deposit what worked best for you?"

[davidlmobley]

Also, @andrrizzi -- what are you generating a conformer for which fails? The host? I would NOT generate conformers for the host; it's too big. I'd just use the provided structure.

[davidlmobley]

Finally, @andrrizzi - it looks to me like your error has "holmium" in the name, which is suspicious, since nothing here contains the element holmium. I suspect you are processing a mol2 file containing the atom type "Ho" or "ho" (the latter being a GAFF atom type) without setting the input flavor to "forcefield" so it is processing the atom types as SYBYL types and reading "ho" (hyrogen bonded to oxygen) as "holmium), which is not what you want.

[nhenriksen]

@davidlmobley The charges (and other parameters) I used in the benchmarks are what I concluded would be a "default" approach for running these calculations. It uses the Q4MD-CD (Link) force field for the CDs. Q4MD-CD uses RESP charges and FF99/Glycam LJs and bonded terms. I used RESP charges on the guest to be consistent, and then GAFF v1.7 for LJs and bonded terms.

I was presuming the benchmark files are mostly to compare code and not to represent the best FF for the job. In fact, of the force fields I've tested, it's a little unclear what the best FF is, because it depends on what data you're comparing to. The best FF for binding free energy (that I tested), which used AM1-BCC/GAFF v1.7 for CD/guest and TIP3P water, seemed to sample incorrect conformations. For example, there was essentially no water in the cavity of the unbound CD, which doesn't seem to be correct. Also the binding enthalpy results were terrible.

[davidlmobley]

Thanks, @nhenriksen .

@andrrizzi - is there an easy way for you to just use Niel's parameters? If not, then I'd suggest you use his charges but your FF if you have to.

[jchodera]

I'm going to be tied up for the next few hours, so I can just comment really quickly that, since the purpose here is to benchmark agreement (and efficiency) versus existing gold standard data from Niel, and NOT so much to compare against experiment, then you would presumably want to run with exactly the same charges he is using rather than reassigning charges.

@davidlmobley: I think we are once again running into a problem with terminology.

I want to be very clear so that there is no potential for confusing the different goals of (1) validation (reproducing the results of existing codes for the exact same parameters), (2) performance (comparing the time required to reach a given error threshold for a given set of parameters), and (3) accuracy (assessing how accurately experimental data is reproduced for a given protocol).

The goals of our YANK paper include (1) validation and (3) accuracy, but not (2) benchmarking.

On the other hand, your work on comparison of restraints methods focuses on (1) validation and (2) benchmarking, but not (3) accuracy.

Our primary concern for this PR, in the yank-benchmark repository (which we should maybe call yank-accuracy so as to eliminate any potential for confusion), is to evaluate (3) accuracy.

[davidlmobley]

Ah, OK. If you are after accuracy of computed free energies then I agree you do not want to run with the exact files Niel provided and you instead want something closer to what Andrea has set up (assuming free energy is your primary metric for "accuracy", as Niel noted).

[andrrizzi]

Sorry for the delay, I didn't have much time to answer during my program retreat.

@jchodera

I'd just use all heavy atoms: (mass > 1.5) ... we should use BAOAB with 2 fs timestep and standard hydrogen masses

Ok! If we'll switch the default integrator to BAOAB instead of Langevin we won't be able to run the validation purely in YAML though. We'll have to add a Python script modifying the integrator move to be Langevin.

and/or GAFF 2.0. (Do we intend to do both?)

No, I think we decided to go with GAFF 1.0. I remember we had the discussion with @davidlmobley too by email.

That's almost assuredly not what we want to do here.

Is it because you don't want to use antechamber or because you want them to be consistent? I can send an email to OpenEye, the host is pretty big for conformer generations though, as @davidlmobley says. In Henriksen & Gilson's paper antechamber was used for am1-bcc charges for both host and guests.

Is 50 mM the experimental ionic strength?

Yep!

Can you be sure to toss a README.md

Sure! I was planning to do it.

@davidlmobley

what are you generating a conformer for which fails? The host?

Correct.

processing a mol2 file containing the atom type "Ho" or "ho" (the latter being a GAFF atom type) without setting the input flavor to "forcefield"

Thanks! I'll look into it tonight.

is there an easy way for you to just use Niel's parameters?

Yes, it should be fairly easy if you have the correct leap parameter files.

• To use RESP charges you'd just have to remove the openeye.quacpac key, and set antechamber.charge_method: null to make antechamber use the charges in the file.
• For the solvent phase, I believe you just need to make sure that your leaprc.gaff installed is the same version as the one used in benchmarksets (1.7 I think). I would double check the parameters for water and ions though.
• If you want to regenerate the complex phase too, you'll have to substitute the ff leap parameters files in the script with the ones used in benchmarksets. Just keep in mind that, differently from APR calculations, YANK sometimes adds extra ions that have to be alchemically annihilated with the ligand to keep the solvent box neutral. If I remember correctly, this host is neutrally charged, so this should not be the case, but it mat be good to double-check.

[jchodera]

Ok! If we'll switch the default integrator to BAOAB instead of Langevin we won't be able to run the validation purely in YAML though. We'll have to add a Python script modifying the integrator move to be Langevin.

I think we either want to make BAOAB the integrator default or use LangevinIntegrator as default and expose the splitting string as an optional argument that defaults to BAOAB.

[jchodera]

No, I think we decided to go with GAFF 1.0. I remember we had the discussion with @davidlmobley too by email.

We should figure out exactly which GAFF 1.x we want to use, then. Probably the one that comes with AmberTools17?

[jchodera]

Is it because you don't want to use antechamber or because you want them to be consistent? I can send an email to OpenEye, the host is pretty big for conformer generations though, as @davidlmobley says.

We should always report thing that don't work if we think they should, otherwise there's really no way to get them fixed and we'll always be cobbling together hacks instead of actually evaluating a protocol for accuracy.

In Henriksen & Gilson's paper antechamber was used for am1-bcc charges for both host and guests.

The aim of our exercise is to evaluate a specified protocol to see how accurate it is. We can't do that if the protocol is a bunch of exceptions or special cases that apply to only one system. I think that we should stick to a single charging scheme/code throughout our entire benchmark. We had wanted to use the canonical AM1-BCC charging scheme because the developer of the method (@cbayly13) feels this is the best practice for using these charges, but if there's no way we can make it work universally (even with @cbayly13's help), we can shift everything to Antechamber-generated charges and use that scheme throughout.

Oddly enough, I had exactly the opposite experience when charging CB7 hosts by hand: I failed to get Antechamber's MOPAC-based scheme to work, while the OpenEye scheme seemed to work fine. As @davidlmobley says, however, he suspects you were using the mol2 files from Niel that contained GAFF atom types rather than Tripos atom types since holonium is not something we would expect to find in CB7 but instead could come from accidentally having an ho GAFF type in your molecule.

[jchodera]

is there an easy way for you to just use Niel's parameters?

For the validation set, where our goal is to show we can reproduce Niel's calculations to demonstrate the code is working as intended, that would be ideal. For the accuracy benchmark comparison to experiment (which this repository is concerned with), we must use a uniform protocol for all systems in order to assess the error in applying that protocol in reality.

[andrrizzi]

I've added the readme files and switched the restrained atoms to all heavy atoms as you suggested.

expose the splitting string as an optional argument that defaults to BAOAB

Good idea! I think we should be able to add this quickly.

he suspects you were using the mol2 files from Niel that contained GAFF atom types rather than Tripos atom types

The script that downloads the files also standardize the atom names, but I've modified YANK to always load mol2 with the forcefield flavor as @davidlmobley suggested and it works now! Thanks! It seems to take an awful lot of memory though. I wasn't able to find the charges on my laptop, but maybe it'll work on the cluster.

[jchodera]

I'm not sure how using the forcefield flavor would help. OpenEye doesn't know anything about GAFF atom or bond types, so I'm not certain how that would let it perceive the correct bond orders and Tripos types. I'm pretty sure we simply never want to be feeding antechamber mol2 files into yank or OpenEye tools since they are non-compliant with the mol2 format. Niel also provided Tripos mol2 files and SDF files, so why not use those?

[andrrizzi]

Sorry if I wasn't clear! By standardizing the atom names, I meant that I've converted them to tripos names, so that has never been the problem that was causing the OpenEye parametrization issue. Setting the flavor to forcefield did solve it.

[andrrizzi]

I've set another couple of flavors though (generic default and mol2 default), so that may be the change responsible for getting past the problem instead of the forcefield flavor.

[davidlmobley]

Using the forcefield flavor normally does allow OpenEye to correctly process molecules with GAFF types, @jchodera. We use this routinely. Chris Bayly added that for precisely this use case.

[jchodera]

Can you explain how this works, @davidlmobley? Does it actually know GAFF atom and bond types? If so, why is this not documented?

[davidlmobley]

I believe it infers the elements from the atom names. And the bind orders in GAFF files are still useful, so this results in obtaining a full correct representation of the molecular graph even without understanding the GAFF types.

I think there may be even more under the hood, though. For example, it can also write back or via this flavor, and the atom types are retained in that case.

I presume this is documented (at least, I think I've seen documentation on the FF flavor reader/writer) but if not then there should be a support request submitted. As I mentioned we use this routinely and have for probably two years ,since it was added. Most of the examples in SMIRNOFF use this so that mol2 files process properly whether or not they are "standard".

[jchodera]

Let's just say I'm wary of relying explicitly on undocumented behavior.

The documentation for the OEIFlavor_MOL2.Forcefield flavor suggests only that the element type is inferred from the atom typename:

This specifies a commonly used variant of MOL2, wherein the chemical element for an atom is inferred from the atom name, but not using the second character if it is capitalized. For example:

CA specifies carbon (not calcium)
Ca specifies calcium
NA specifies nitrogen (not sodium)
Na specifies sodium
HA specifies hydrogen (not hahnium)
Exceptions were made for atom names beginning with FE, ZN, MG, MN, CL, and BR: these specify chemical elements Fe, Zn, Mg, Mn, Cl, and Br, respectively.

No mention of bond orders is made. I don't see how one could infer bond orders in a foolproof manner from element names alone without knowing the formal charge as well, but without documentation of what is done, it is impossible to know how reliable this might be.

We should not use Antechamber format mol2 files---we should only use compliant Tripos mol2 files.

[nhenriksen]

Can you clarify what you mean by "compliant Tripos mol2 files"? Is it just using SYBYL atom types? or are there additional criteria which need to be met?

just to clarify, I'm assuming the format should be as defined here. But just checking if there are additional criteria beyond those specifications.

[davidlmobley]

@jchodera

No mention of bond orders is made. I don't see how one could infer bond orders in a foolproof manner from element names alone without knowing the formal charge as well, but without documentation of what is done, it is impossible to know how reliable this might be.

The bond orders are in the mol2 files even in Antechamber format. Antechamber uses a different kekulization scheme, but they are still there, and are still useful. There is no need to infer them. All that has to be done is to read them and correctly perceive the elements and one has a complete representation of the molecule. This is FAR FAR less dangerous than what, say, Antechamber does when trying to get TO a mol2 file from a PDB file (with or without CONECT entries), or what OpenEye tools do in a similar situation, because (a) you DO have a complete representation of the bonding pattern, and (b) you DO have a complete and correct representation of the bond orders, even if the kekulization scheme is different.

(I'm excluding the case of mol2 files which have incorrect bond orders because of user and/or input error -- e.g. if you provide antechamber with a PDB file and it incorrectly infers bond orders and produces a mol2 file with incorrect bond orders, they will be incorrect. But this would be true whether you use Tripos mol2 files with SYBYL types or whether you use GAFF types.)

Let's just say I'm wary of relying explicitly on undocumented behavior.

What specifically is the "undocumented behavior" you are worried about here? The forcefield flavor infers the element types, and uses the bond information to create an OpenEye molecule. As you point out, the document explains how it infers the element types (and my explanation was incorrect). What's undocumented about this? I guess you could argue that it's "undocumented" that it uses the bond orders to figure out what molecule it is -- but this is just standard OpenEye toolkit behavior; whatever information you provide about the molecule is used to represent that molecule.

(Regarding failure of conformer generation) We should always report thing that don't work if we think they should, otherwise there's really no way to get them fixed and we'll always be cobbling together hacks instead of actually evaluating a protocol for accuracy.

I suspect the real reason @andrrizzi 's conformer generation fails is just because he provided a forcefield flavor mol2 file without using the forcefield flavor reader (conformer generation via Omega is not supposed to work for 'methoxyholmium...' molecules). It seems likely to me that if the host is processed correctly, conformer generation will work or at least attempt to work. It's still not obvious to me that trying to generate a conformation for the host is a good idea, but it may work once this is fixed.

@nhenriksen

Can you clarify what you mean by "compliant Tripos mol2 files"? Is it just using SYBYL atom types? or are there additional criteria which need to be met?

Yes, in OpenEye's view a mol2 file was invented by Tripos and so their specs determine what can be in it and how. This means SYBYL atom types, so a mol2 file with AMBER or GAFF atom types is technically not a Tripos mol2 file. It used to be that their readers could not handle these at all, but eventually they had to add a forcefield "flavor" reader because these "non-compliant" mol2 files are so common in the molecular modeling field. The main thing the forcefield flavor reader is there to cover is cases where the atom types are non-SYBYL.

[jchodera]

The bond orders are in the mol2 files even in Antechamber format. Antechamber uses a different kekulization scheme, but they are still there, and are still useful. There is no need to infer them. All that has to be done is to read them and correctly perceive the elements and one has a complete representation of the molecule.

Correctly perceiving bond orders from an Antechamber mol2 file requires that OEChem translate Antechamber bond orders (which are not well documented) to Tripos mol2 or some other type of bond orders. The OEChem docs do not specify that this translation is performed or how it is done. If you use this feature, you are relying on undocumented behavior of OEChem that is required for your translation, end of story. This is dangerous at best.

If you think we should ever allow Antechamber mol2 files to be processed this way, we need to press OpenEye to fully document this feature.

[jchodera]

What specifically is the "undocumented behavior" you are worried about here? The forcefield flavor infers the element types, and uses the bond information to create an OpenEye molecule.

It's the bond order interpretation that I'm worried about. The GAFF bond order types produced by Antechamber are not well documented. I'm concerned about this because I ran into the exact same problem when writing this function and am not at all convinced my version of the code works properly under all circumstances. It's not a theoretical problem, but a very concrete one.

[jchodera]

I suspect the real reason @andrrizzi 's conformer generation fails is just because he provided a forcefield flavor mol2 file without using the forcefield flavor reader (conformer generation via Omega is not supposed to work for 'methoxyholmium...' molecules). It seems likely to me that if the host is processed correctly, conformer generation will work or at least attempt to work.

Agreed, but since Niel provided Tripos format mol2 files and SDF files (in addition to the Antechamber ones), we should just use the mol2 files that are compliant to the Tripos format instead of relying on undocumented behavior of the forcefield flavor!

It's still not obvious to me that trying to generate a conformation for the host is a good idea, but it may work once this is fixed.

The canonical AM1-BCC scheme from @cbayly13 ensures that the original conformation is included, so even if Omega fails to generate new conformations as well, the protocol should still produce useful charges.

[davidlmobley]

OK, yes, if you think this is tricky and needs to be documented then I think it's worth submitting a support request asking them to document it. It's extremely useful to have OpenEye be able to process typical force field-style mol2 files without intermediate conversion (since these are SO OFTEN provided by people in our field without alternatives) and they added this functionality for precisely this use case, so if you're uncomfortable with using it as is you should definitely ask them to document it.

[jchodera]

To make this work robustly, we'd need to have a different flavor for each forcefield, wouldn't we?

My understanding was that the original reason for the forcefield flavor was NOT to allow OpenEye to correctly perceive molecular structure, but simply to allow it to RETAIN the entries in the atom and bond type categories instead of overwriting them. It would be awesome if it could deal with GAFF natively, but I really think that would require a GAFF-specific flavor, which is what my support request would be for.

I honestly just don't think it's possible to automagically detect forcefield type and the meaning of integral bond types in a robust way. It seems like wishful thinking.

[andrrizzi]

@Lnaden just a heads up that before starting the calculations, you may need to regenerate the CD input files by calling scripts/prepare_input_files.py:prepare_cyclodextrin_files() since the files have been updated in MobleyLab/benchmarksets#47 with SYBYL atom types and different positions.

[jchodera]

Should we merge this PR, @andrrizzi ?