Silent GAFF2 failures can happen easily

[simonaxelrod]

Hi there,

I just wanted to let everyone know about a silent failure mode that I found recently. If you're using GAFF2 to parameterize the ligand, we know that you must provide its leaprc file under systems:

systems:
  # We can call our system anything we want, this example just uses a short name for the receptor hyphenated with the ligand
  PARP1-binder:
    # These names all use the names we defined previously
    receptor: PARP1
    ligand: binder
    solvent: pme
    leap:
      parameters: [leaprc.protein.ff14SB, leaprc.gaff2, leaprc.water.tip4pew] # we knew that this needed to be included

However, what I didn't know was that I must also provide its leaprc file under molecules/binder:

molecules:
  # Define our receptor. We can call it whatever we want so we just use its name here as the directive
  PARP1:
    filepath: protein.pdb
  # Define our ligand molecule
  binder:
    filepath: ligand.mol2
    # Get partial charges
    antechamber:
      charge_method: bcc
      net_charge: 0
    leap:
      parameters: [leaprc.gaff2] # must include here

If you don't, then GAFF atom types will be used, not GAFF2 atom types. This comes from this line in the code. If leaprc.gaff2 is not provided in parameters, then Yank assumes that we're using GAFF.

As far as I've seen, this never leads to an error in actually running the simulation. Instead, the simulation runs without errors, but uses an incorrect ligand force field that's a mix-and-match between GAFF and GAFF2. Therefore it's a silent failure mode and can be dangerous. I suggest that we change this line so that the failure can't happen silently.

Simon