NTerm GLU torsion energies are bad

j-wags commented 4 years ago

GLU

NTerm GLU has a large torsion energy discrepancy (~20 kJ/mol). CTerm ASP is off by 3 kJ/mol. Both discrepancies are in torsions.

Top suspects:

~our C-terminal fixing code being too judicious and modifying carboxylates it shouldn't~
- I don't think it's this, because the SMILES for the OFFMol show a reasonable carboxylate
~The hierarchy is being scrambled and there's some conflict between the protonated/deprotonated SMIRKS~
- I inspect the parameter SMIRKS below, and they do indeed include specific charge/protonation states
Something with SMIRNOFF vs. AMBER improper functional forms
At least for NTerm GLU, a backend toolkit may be taking liberties with how it interprets zwitterions.
Something about the equivalence of the two oxygens in the charged sidechains causing parameter duplication
- Does AMBER apply multiple impropers to that center?
~Is the sidechain somehow chemcially equivalent to an uncapped CTerm?~
- They're chemically distinct. Also, it would have to be ONLY the GLU Cterm, since this problem reproduces when only GLU is used to make the OFFXML

Next step would be to unpack the parameters from each system and see what's being mis-assigned

NTerminal/GLU/GLU
HarmonicBondForce 2.1379218101501465 kJ/mol
HarmonicAngleForce 6.882795810699463 kJ/mol
PeriodicTorsionForce 46.8884162902832 kJ/mol
NonbondedForce -415.76202392578125 kJ/mol
-359.8529052734375 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C([H])([H])C(=O)[O-])[N+]([H])([H])[H]
HarmonicAngleForce 6.882845401763916 kJ/mol
HarmonicBondForce 2.1379215717315674 kJ/mol
NonbondedForce -415.7628173828125 kJ/mol
PeriodicTorsionForce 28.416242599487305 kJ/mol
-378.3258056640625 kJ/mol
67 amber dihedrals (3 impropers) and 73 off dihedrals

Related values for reference:

CTerm isn't that bad:

CTerminal/GLU/GLU
HarmonicBondForce 3.8477184772491455 kJ/mol
HarmonicAngleForce 11.950961112976074 kJ/mol
PeriodicTorsionForce 51.16751480102539 kJ/mol
NonbondedForce 186.67312622070312 kJ/mol
253.63931274414062 kJ/mol
[H][C@@](C(=O)[O-])(C([H])([H])C([H])([H])C(=O)[O-])N([H])C(=O)C([H])([H])[H]
HarmonicAngleForce 11.950912475585938 kJ/mol
HarmonicBondForce 3.8477187156677246 kJ/mol
NonbondedForce 186.671630859375 kJ/mol
PeriodicTorsionForce 50.02157211303711 kJ/mol
252.491943359375 kJ/mol
74 amber dihedrals (4 impropers) and 82 off dihedrals

For comparison, here's mainchain GLU/GLH:

MainChain/GLH/GLH
HarmonicBondForce 2.53347110748291 kJ/mol
HarmonicAngleForce 12.666474342346191 kJ/mol
PeriodicTorsionForce 89.4136734008789 kJ/mol
NonbondedForce -312.77301025390625 kJ/mol
-208.15940856933594 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C([H])([H])C(=O)O[H])N([H])C(=O)C([H])([H])[H]
HarmonicAngleForce 12.666414260864258 kJ/mol
HarmonicBondForce 2.5334713459014893 kJ/mol
NonbondedForce -312.7730712890625 kJ/mol
PeriodicTorsionForce 89.38488006591797 kJ/mol
-208.18833923339844 kJ/mol
105 amber dihedrals (5 impropers) and 115 off dihedrals

MainChain/GLN/GLN
HarmonicBondForce 0.08747541904449463 kJ/mol
HarmonicAngleForce 1.6231834888458252 kJ/mol
PeriodicTorsionForce 78.39984130859375 kJ/mol
NonbondedForce 93.097900390625 kJ/mol
173.2084197998047 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C([H])([H])C(=O)N([H])[H])N([H])C(=O)C([H])([H])[H]
HarmonicAngleForce 1.6231499910354614 kJ/mol
HarmonicBondForce 0.08747542649507523 kJ/mol
NonbondedForce 93.09825134277344 kJ/mol
PeriodicTorsionForce 78.39984893798828 kJ/mol
173.208740234375 kJ/mol
106 amber dihedrals (6 impropers) and 118 off dihedrals

MainChain/GLU/GLU
HarmonicBondForce 2.5307743549346924 kJ/mol
HarmonicAngleForce 5.230216979980469 kJ/mol
PeriodicTorsionForce 62.76970672607422 kJ/mol
NonbondedForce -118.9146728515625 kJ/mol
-48.383941650390625 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C([H])([H])C(=O)[O-])N([H])C(=O)C([H])([H])[H]
HarmonicAngleForce 5.230195999145508 kJ/mol
HarmonicBondForce 2.5307741165161133 kJ/mol
NonbondedForce -118.9154052734375 kJ/mol
PeriodicTorsionForce 62.565921783447266 kJ/mol
-48.5885009765625 kJ/mol
94 amber dihedrals (5 impropers) and 104 off dihedrals

ASP

Mainchain ASP/N/H

MainChain/ASH/ASH
HarmonicBondForce 2.550563097000122 kJ/mol
HarmonicAngleForce 11.37463092803955 kJ/mol
PeriodicTorsionForce 86.83865356445312 kJ/mol
NonbondedForce -358.69512939453125 kJ/mol
-257.9312744140625 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C(=O)O[H])N([H])C(=O)C([H])([H])[H]
HarmonicAngleForce 11.37458324432373 kJ/mol
HarmonicBondForce 2.550563097000122 kJ/mol
NonbondedForce -358.6954040527344 kJ/mol
PeriodicTorsionForce 86.74219512939453 kJ/mol
-258.028076171875 kJ/mol
93 amber dihedrals (5 impropers) and 103 off dihedrals

MainChain/ASN/ASN
HarmonicBondForce 3.6648805141448975 kJ/mol
HarmonicAngleForce 7.699925899505615 kJ/mol
PeriodicTorsionForce 82.42230224609375 kJ/mol
NonbondedForce -492.32086181640625 kJ/mol
-398.5337829589844 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C(=O)N([H])[H])N([H])C(=O)C([H])([H])[H]
HarmonicAngleForce 7.699950218200684 kJ/mol
HarmonicBondForce 3.6648802757263184 kJ/mol
NonbondedForce -492.320556640625 kJ/mol
PeriodicTorsionForce 82.2682876586914 kJ/mol
-398.6873779296875 kJ/mol
94 amber dihedrals (6 impropers) and 106 off dihedrals

MainChain/ASP/ASP
HarmonicBondForce 3.1304473876953125 kJ/mol
HarmonicAngleForce 6.654740333557129 kJ/mol
PeriodicTorsionForce 64.0936279296875 kJ/mol
NonbondedForce -183.40335083007812 kJ/mol
-109.52450561523438 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C(=O)[O-])N([H])C(=O)C([H])([H])[H]
HarmonicAngleForce 6.654727458953857 kJ/mol
HarmonicBondForce 3.1304476261138916 kJ/mol
NonbondedForce -183.40423583984375 kJ/mol
PeriodicTorsionForce 63.85840606689453 kJ/mol
-109.76068115234375 kJ/mol
86 amber dihedrals (5 impropers) and 96 off dihedrals

CTerm ASN/P


CTerminal/ASN/ASN
HarmonicBondForce 4.261889934539795 kJ/mol
HarmonicAngleForce 12.491641998291016 kJ/mol
PeriodicTorsionForce 70.66112518310547 kJ/mol
NonbondedForce -500.82763671875 kJ/mol
-413.4130859375 kJ/mol
[H][C@@](C(=O)[O-])(C([H])([H])C(=O)N([H])[H])N([H])C(=O)C([H])([H])[H]
HarmonicAngleForce 12.491678237915039 kJ/mol
HarmonicBondForce 4.261890411376953 kJ/mol
NonbondedForce -500.82867431640625 kJ/mol
PeriodicTorsionForce 70.4883804321289 kJ/mol
-413.5866394042969 kJ/mol
74 amber dihedrals (5 impropers) and 84 off dihedrals

CTerminal/ASP/ASP
HarmonicBondForce 4.411398410797119 kJ/mol
HarmonicAngleForce 17.078983306884766 kJ/mol
PeriodicTorsionForce 59.746490478515625 kJ/mol
NonbondedForce 229.81549072265625 kJ/mol
311.05230712890625 kJ/mol
[H][C@@](C(=O)[O-])(C([H])([H])C(=O)[O-])N([H])C(=O)C([H])([H])[H]
HarmonicAngleForce 17.078941345214844 kJ/mol
HarmonicBondForce 4.411397457122803 kJ/mol
NonbondedForce 229.813232421875 kJ/mol
PeriodicTorsionForce 56.86469268798828 kJ/mol
308.16827392578125 kJ/mol
66 amber dihedrals (4 impropers) and 74 off dihedrals

Nterm

NTerminal/ASN/ASN
HarmonicBondForce 2.3103301525115967 kJ/mol
HarmonicAngleForce 7.480605125427246 kJ/mol
PeriodicTorsionForce 47.826507568359375 kJ/mol
NonbondedForce -330.39764404296875 kJ/mol
-272.78021240234375 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C(=O)N([H])[H])[N+]([H])([H])[H]
HarmonicAngleForce 7.480668544769287 kJ/mol
HarmonicBondForce 2.3103301525115967 kJ/mol
NonbondedForce -330.3980712890625 kJ/mol
PeriodicTorsionForce 47.77263641357422 kJ/mol
-272.83441162109375 kJ/mol
67 amber dihedrals (4 impropers) and 75 off dihedrals

NTerminal/ASP/ASP
HarmonicBondForce 0.9721102714538574 kJ/mol
HarmonicAngleForce 4.750027656555176 kJ/mol
PeriodicTorsionForce 29.664676666259766 kJ/mol
NonbondedForce -335.2193298339844 kJ/mol
-299.83258056640625 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C(=O)[O-])[N+]([H])([H])[H]
HarmonicAngleForce 4.750013828277588 kJ/mol
HarmonicBondForce 0.9721102714538574 kJ/mol
NonbondedForce -335.2209777832031 kJ/mol
PeriodicTorsionForce 29.531963348388672 kJ/mol
-299.9668884277344 kJ/mol
59 amber dihedrals (3 impropers) and 65 off dihedrals

j-wags commented 4 years ago

The charge/protonation states on the NTerm GLU SMIRKS are indeed explicit, so I don't think there could be confusion with the charged AA being a substructure of the protonated AA

<Proper 
smirks="[H][C@@]([C]=O)(C([H])([H])[C:2]([H])([H:1])[C:3](=[O:4])[O-])[N+]([H])([H])[H]"
periodicity1="2" 
phase1="0.0 * degree" 
id="A14SB-NTerminal-GLU-HC_2C_CO_O2" 
k1="0.0 * mole**-1 * kilocalorie" 
idivf1="1.0"></Proper>

<Improper 
smirks="[H][C@@]([C]=O)(C([H])([H])[C:1]([H])([H])[C:2](=[O:4])[O-:3])[N+]([H])([H])[H]" 
periodicity1="1" 
phase1="180.0 * degree" 
k1="3.5 * mole**-1 * kilocalorie" 
id="A14SB-NTerminal-GLU-2C_O2_CO_O2" 
idivf1="1.0">
</Improper>

j-wags commented 4 years ago

If I rerun with only GLU, the I get the exact same values for NTerm GLU. So this isn't a case of parameters for another AA being incorrectly substituted into the GLU sidechain

j-wags commented 4 years ago

Is it an extreme case of calculating the improper energy differently?

This improper center is very bent, by virtue of its proximity to the positively charged N terminus.

Depending on how I measure around the improper center, I can get angle values ranging from 150.8 and 152.5 degrees.

k*(1+cos(periodicity*theta-phase)) periodicity1="1" phase1="180.0 * degree" k1="3.5 * mole**-1 * kilocalorie"

import math
k=3.5
per=1.
ph=180.
theta = 152.5
k*(1+math.cos(math.radians(per*(theta-ph)))) * 4.184

27.63338664106188
theta=150.8
k*(1+math.cos(math.radians(per*(theta-ph)))) * 4.184
27.427070899539093

This is a substantial difference, but unless I'm missing a factor of 10 in the math, this doesn't explain the observed discrepancy

j-wags commented 3 years ago

It looks like, at some point in the past 40 days, we fixed this. Current output reads:

NTerminal/GLU/GLU
HarmonicBondForce 2.1379218101501465 kJ/mol
HarmonicAngleForce 6.882795810699463 kJ/mol
PeriodicTorsionForce 46.8884162902832 kJ/mol
NonbondedForce -415.76202392578125 kJ/mol
-359.8529052734375 kJ/mol
[H][C@@](C(=O)N([H])C([H])([H])[H])(C([H])([H])C([H])([H])C(=O)[O-])[N+]([H])([H])[H]
HarmonicAngleForce 6.882845401763916 kJ/mol
HarmonicBondForce 2.1379215717315674 kJ/mol
NonbondedForce -415.7628173828125 kJ/mol
PeriodicTorsionForce 46.889320373535156 kJ/mol
-359.85272216796875 kJ/mol
67 amber dihedrals (3 impropers) and 67 off dihedrals

openforcefield / amber-ff-porting