Skip to content
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 cyclobutane parameters which produce puckered ring #65

Closed
bannanc opened this issue Nov 18, 2017 · 17 comments
Closed

Add cyclobutane parameters which produce puckered ring #65

bannanc opened this issue Nov 18, 2017 · 17 comments

Comments

@bannanc
Copy link
Collaborator

bannanc commented Nov 18, 2017

@davidlmobley Noticed that when running simulations with smirnoff99Frosst on molecules with cyclobutyl groups the ring is flat, which we don't think is correct and isn't true in GAFF (or parm@frosst).

When we talked about this in person, David mentioned there are cyclobutyl groups in the AlkEthOH molecule set which means the carbons should be assigned CT. I'm going to start by double checking which parameters the molecules in AlkEthOH get with smirnoff99Frosst and which atom types those parameters are assigned. Then I'll look to see if the parameters agree.

I'm guessing that there are some kind of angle or torsion parameters with @ for ring bonds that might be messing the cyclobutane when they shouldn't.

Mostly for my own reference, David has code for minimizing and running gas phase simulations with SMIRNOFF in the MobleyLab/SMIRNOFF_paper_code or we also have minimizations in the MobleyLab/off-ffcompare repo which is public.

@bannanc
Copy link
Collaborator Author

bannanc commented Nov 21, 2017

Here's an update:

I did a minimization with the off-ffcompare pipeline for cyclobutane.
Here are the minimized structures:

877e4f86-064c-45aa-95f8-d289e7e40e61

  • cyan/white - smirnoff99Frosst
  • yellow - MMFF94S
  • red - GAFF

I also tried running compare_molecule_energies.py with cyclobutane. When I used the smirnoff tailored to exactly replicate parm@frosst for AlkEthOH (Frosst_AlkEthOH_parmAtFrosst.ffxml) all of the energies agreed. When I used smirnoff99Frosst from this repository, I got the following error message:

AssertionError: Error: Harmonic Angle (1, 3, 2) has force constant values of 0.101962 and 0.127452 kJ/(mol degree**2), respectively.

@davidlmobley
Copy link
Collaborator

So, in other words, you've confirmed that there is a problem, @bannanc , and particularly that the problem might be with angles?

(Also, you may only be seeing the FIRST different term because of when the error is raised; it's possible you want to try it with a flag set so it doesn't raise assertion errors so you can go ahead and see the energies/which other terms disagree. My recollection is that it will let you do that...)

@jchodera
Copy link
Member

That cyan structure looks exactly planar. Could it be that there hasn't been anything to break the symmetry, and that it's just "stuck" in the planar energy maximum? If you feed in the MMFF minimized geometry, does it go back to planar?

@jchodera
Copy link
Member

Ah, I see this was drawn from a simulation, so it probably already had its symmetry broken and went back to being planar during minimization. Sorry!

@bannanc
Copy link
Collaborator Author

bannanc commented Nov 21, 2017

@jchodera the off-ffcompare pipeline minimizes with a series of force fields all starting from the same input tripos mol2 file so they all started with the same input coordinates. David has run into this problem before so we're fairly certain there is an issue with parameters in smirnoff.

@davidlmobley yes, basically I've confirmed the problem is what we thought it was.
The only difference in parameters is that in parm@frosst there are three different angles around the carbons:

HC-CT-HC    35.0      109.50
CT-CT-HC    50.0      109.50
CT-CT-CT    40.0      109.50

and in smirnoff there are only two.

[*:1]~[#6X4:2]-[*:3]      50.  109.5   parm99 generic Csp3
[#1:1]-[#6X4:2]-[#1:3]    35.  109.5   parm99 generic Csp3

In GAFF there have a special carbon type for 4-membered rings (cy) and that has a different set of angles:

cy-cy-hc      44.83      114.79
cy-cy-cy      69.75       88.40
hc-cy-hc      39.22      108.98

It was easier for me to add a parameter to smirnoff and rerun in the pipeline I'm using then to find parm@frosst input files. If I add an angle parameter corresponding to CT-CT-CT, making the smirnoff parameters:

[*:1]~[#6X4:2]-[*:3]      50.  109.5   parm99 generic Csp3
[#6X4:1]-[#6X4:2]-[#6X4:3]  40.  109.5 CCB adding to match CT-CT-CT
[#1:1]-[#6X4:2]-[#1:3]    35.  109.5   parm99 generic Csp3

I still get a flat center ring:
08ef3f15-6a3a-4078-ace4-184b2d91202f

  • cyan/white - smirnoff here
  • pink - smirnoff with added CT-CT-CT angle

I'm working on a parm@frosst minimization right now to confirm, but I think this is an inherited issue.

@davidlmobley
Copy link
Collaborator

OK, so if this is an issue also present in parm@Frosst we would want to look at what GAFF/GAFF2 have done to get this right.

@bannanc
Copy link
Collaborator Author

bannanc commented Nov 21, 2017

Confirmed, parm@frosst also minimizes to a flat ring, I can try it with a substituent to see if that changes anything.

44d1052d-76ef-4a82-94b8-88deae00c540
cyan - smirnoff (here)
purple - parm@frosst

I'm pulling GAFF2 parameters since that is where we have gotten estimates for other additions to the force field. Since this is a completely new atomtype, there is a lot to parse, but as I see it these are the most general types we would want to add:

  • ring carbon - ring carbon - non ring atom (different for hydrogen or heavy atom?)
    • "[#6X4r4:1]-;@[#6X4r4:2]-;!@[*:3]"
    • "[#6X4r4:1]-;@[#6X4r4:2]-;!@[#1:3]"
  • ring carbon - ring carbon - ring carbon
    • "[#6X4r4:1]-;@[#6X4r4:2]-;@[#6X4r4:3]"
  • non ring atom - ring carbon - non ring atom (differences for heavy atoms and hydrogen?)
    • "[*:1]-;!@[#6X4r4:2]-;!@[*:3]"
    • "[#1:1]-;!@[#6X4r4:2]-;!@[*:3]"
    • "[#1:1]-;!@[#6X4r4:2]-;!@[#1:3]"

The full set of GAFF2 angles around cyclobutane rings is long so I'm not going to copy it in here right now. I'll try to get it sorted into the categories above and propose what we should add later this week.

@bannanc
Copy link
Collaborator Author

bannanc commented Nov 21, 2017

Oh, just as a note, there aren't many cy parameters for torsions, but they all seem to be the same as c3 (GAFF version of CT, tetrahedral/sp3 carbon).

@bannanc
Copy link
Collaborator Author

bannanc commented Nov 28, 2017

Here are my recommendations, the only changes are in angles, but I'll outline everything I looked at for posterity.
I will be comparing parameters for cy (GAFF2 for sp3 carbon in 4 membered ring) to c3 (GAFF2 for general sp3 carbons).

Non-Bond

Unsurprisingly the LJ parameters for cy are identical to c3. We do not need to make changes here.

Torsions

The only dihedral parameters GAFF2 has for cy are copied directly from c3 so we have no evidence that new parameters are necessary.
There are no improper torsions for that atom type.

Bond

cy-cy  218.51   1.556  
c -cy  222.38   1.551
c3-cy  236.93   1.532
ce-cy  250.06   1.516
cf-cy  250.06   1.516
ca-cy  250.49   1.515
cv-cy  250.66   1.515
cx-cy  250.83   1.515
c2-cy  255.91   1.509
cd-cy  259.34   1.505
cc-cy  259.34   1.505

Almost all parameters here are the same as the corresponding values for c3 in GAFF2. However, I notices that in general, SMIRNOFF99Frosst has larger force constants than GAFF. For example:
SMIRNOFF: [#6X4:1]-[#6X4:2] 310.0 1.526 parm99 CT-CT
GAFF2: c3-c3 232.52 1.538 SOURCE1_SOURCE5 88072 0.0058
The smirnoff values come from the frcmod file so this isn't related to the factor of 2 that amber doesn't show.

One thing to note here is that the bond between two ring carbons has a slightly lower force constant. I am not convinced we need another parameter for this, at least a this point in time, but might be worth keeping in mind.

Angle

I started by organizing parameters from GAFF2 based on the categories below. For each I list the category, my proposed SMIRNOFF parameter and then show the relevant lines from gaff2.dat.

internal ring angle
SMIRNOFF: [#6r4]-;@[#6r4]-;@[#6r4] 72. 88.
GAFF2:

cy-cy-cy   72.002      88.400      
cv-cy-cy   73.645      86.710     cyclobutene (internal ring angle)

external angle with only heavy atoms
SMIRNOFF: [!#1:1]-[#6r4:2]-;!@[!#1:3] 63. 119.
GAFF2 (2 ring atoms):

cy-cy-ss   60.567     118.400    
cy-cy-s6   60.829     117.220    
br-cy-cy   61.847     119.270    
cc-cy-cy   62.477     121.220    cc = Sp2 carbons in non-pure aromatic systems
cd-cy-cy   62.477     121.220    cd = Sp2 carbons in non-pure aromatic systems
c3-cy-cy   62.618     118.700
cx-cy-cy   68.163     101.230     cx = sp3 carbon in cyclopropyl group 
c2-cy-cy   68.568     100.400     c2 = sp2 carbon
cl-cy-cy   69.535     117.010
c -cy-cy   73.514      85.090     c = Sp2 C carbonyl group
cy-cy-no   79.836     115.430
cy-cy-na   80.160     119.560
cy-cy-oh   82.786     114.190
cy-cy-os   83.069     111.770
cy-cy-f    85.702     112.870
cy-cy-n4   90.361      89.940    n4 = sp3 N with 4 connected atoms
cy-cy-n    91.445      90.630    n  = sp2 nitrogen in amide group
cy-cy-n3   92.778      87.580   n3 = sp3 N with 3 connected atoms

GAFF2 (1 ring atom):

ce-cy-ss   60.454     120.540
cf-cy-ss   60.454     120.540
br-cy-cy   61.847     119.270
c -cy-c3   63.261     116.710  
c3-cy-c3   65.146     111.440
ss-cy-ss   68.296      95.040 
cl-cy-cy   69.535     117.010 
c -cy-cl   71.094     112.180
c -cy-n    80.522     117.340
n -cy-ss   81.655     105.130
c -cy-os   82.061     115.000
c3-cy-n3   82.158     113.650
n -cy-s6   82.361     103.180
c3-cy-os   83.640     112.310
c3-cy-n    86.069     104.130
cf-cy-n    94.195      87.940
ce-cy-n    94.195      87.940
f -cy-py   99.595     113.190
n -cy-os  109.575     107.420
f -cy-f   120.355     108.560

@cbayly13 this is a section I could use some help. In general I expect the external angles to be slightly larger than 109.5º. I defined our type based on that assumption, but there are some notable exceptions here, I've tried to define the atom type in those cases. Do you think it is necessary at this point to track down these exceptions and make another angle parameter for them?

external angle with 1 heavy atom and 1 hydrogen
SMIRNOFF: [!#1:1]-[#6r4:2]-;!@[#1:3] 45. 115.
GAFF2 (ring-ring-hydrogen):

cy-cy-h2   44.908     116.780
cy-cy-hx   45.080     115.920
cy-cy-hc   45.290     114.790
cy-cy-h1   45.598     113.240

GAFF2 (no ring - ring - hydrogen):

h1-cy-s6   41.611     111.620
h1-cy-ss   41.667     111.560
h2-cy-s6   41.906     110.010
h2-cy-ss   41.997     109.770
cx-cy-hc   45.715     118.300
c -cy-h1   45.784     113.030
ce-cy-h2   45.874     117.400
cf-cy-h2   45.874     117.400
c -cy-hc   46.135     111.320
ce-cy-h1   46.242     115.490
cf-cy-h1   46.242     115.490
cf-cy-hc   46.374     114.840
ca-cy-hc   46.450     114.530
c3-cy-h1   46.561     111.780
c3-cy-hc   46.907     110.140
c2-cy-hc   46.979     112.800
cl-cy-hc   47.632     114.000
cd-cy-hc   48.302     107.200
cc-cy-hc   48.302     107.200
cl-cy-h1   49.261     106.590
hx-cy-n4   58.977     110.620
h1-cy-nh   59.704     113.860
h1-cy-n3   59.591     113.930
h2-cy-n    59.772     114.500
h1-cy-n    60.615     111.280
h1-cy-os   62.029     110.130
h2-cy-os   62.315     109.190
h1-cy-oh   62.604     111.490
h1-cy-na   62.857     106.380

external angle with 2 hydrogens
SMIRNOFF: don't add anything as this agrees with the existing parameter
[#1:1]-[#6X4:2]-[#1:3] 35. 109.5 parm99 generic Csp3
GAFF2:

h1-cy-h1   38.688     109.460
hc-cy-hc   38.780     108.980
hx-cy-hx   38.598     110.800

@bannanc
Copy link
Collaborator Author

bannanc commented Nov 28, 2017

TL;DR
I want to add these parameters

[#6r4:1]-;@[#6r4:2]-;@[#6r4:3] 72. 88.
[!#1:1]-[#6r4:2]-;!@[!#1:3] 63. 119.
[!#1:1]-[#6r4:2]-;!@[#1:3] 45. 115.

below this line in the frcmodish file for SMIRNOFF

@davidlmobley
Copy link
Collaborator

@bannanc - do you want to see what happens if you add those and run some dynamics, e.g. comparing geometries for these types of molecules resulting from new and old cases?

@bannanc
Copy link
Collaborator Author

bannanc commented Nov 28, 2017

I added the parameters and did a minimization (no dynamics).

image
cyan/white - current smirnoff
pink - smirnoff + cyclobutan parameters proposed above.

@bannanc
Copy link
Collaborator Author

bannanc commented Nov 28, 2017

@davidlmobley and @cbayly13

@davidlmobley
Copy link
Collaborator

🎉 🎉
Excellent. Seems like a great start then!

@cbayly13
Copy link
Collaborator

Caught up on the issue; I think @bannanc 's reasoning and actions are sound. I am surprised our parm@Frosst cyclobutane is planar; we knew it should be puckered and I thought we put in the effort to make it so but I guess not. I am surprised it is modulated all by bond angles and not torsions but there it is. The GAFF2 parameters look overfit but divide into broad categories as @bannanc has done; with heteroatom substituents like N and S it looks like a couple of fairly general child parameters might be warranted. parm@Frosst was definitely wrong here; GAFF2 has led the way. Nice catch @davidlmobley and @bannanc .

@bannanc
Copy link
Collaborator Author

bannanc commented Nov 29, 2017

Thanks @cbayly13
I think further exploration is necessary, but it seemed worth getting at least a generic fix in there.

@davidlmobley davidlmobley changed the title Cyclobutane parameters Add cyclobutane parameters which produce puckered ring Nov 29, 2017
@davidlmobley
Copy link
Collaborator

davidlmobley commented Nov 29, 2017

Closed via openforcefield/openff-toolkit#83 and #66

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants