TinkerTools / tinker

Tinker: Software Tools for Molecular Design
https://dasher.wustl.edu/tinker/
Other
130 stars 61 forks source link

Potential typo in MM3 silicon equilibrium angle #142

Closed philipturner closed 6 months ago

philipturner commented 9 months ago

I am working on a forcefield that combines MM4 and MM3 parameters, and one of the many cross-references I look at is your repository. Allinger made a typo in the MM3 silicon paper (1997) specifying effectively:

torsion5      1    1   19    1      0.000 0.0 1   0.800 180.0 2   0.000 0.0 3
torsion5      1    1   19    1      0.000 0.0 1   0.000 180.0 2   0.250 0.0 3

He corrected the 1-1-19-1 duplicate parameter in the MM3(2000) parameters he handed to you. The second "1-1-19-1" was split into two similar ones with an alkane carbon replaced with an alkene. "0.250" for V3 materialized in the latter. I did some sketches of mixed C and Si cyclopentane rings, and the removal of the V3 term would coincide with the ideal torsion angle shifting from 120° to 90° around the silicon.

Replacement of V3 with V2 only happened before when Csp2 was involved (1-2-2-19). The absolute value 0.800 for 1-1-19-1 does match the 0.850 from the 1-1-1-19 parameter, as both changed from 0.167 to ~0.800. However, 2-1-1-19 had 0.167 in both torsion6 and torsion5 (no difference), and 2-1-1-19 torsion5 wasn't present in the original paper. There seems something fishy going on with X-1-1-19 and Csp2. Potentially, all 5-ring alkane-only parameters in that section could actually involve a Csp2 somewhere. Allinger only corrected some of the parameters in MM3(2000), or at least they're questionable now. I'll stick with the ones for 6-membered rings instead.

torsion5      1    1   19    1      0.000 0.0 1   0.800 180.0 2   0.000 0.0 3
torsion5      1    1   19    2     -0.500 0.0 1   0.000 180.0 2   0.000 0.0 3
torsion5      2    1   19    1      0.000 0.0 1   0.000 180.0 2   0.250 0.0 3

I have also corrected what I strongly think is a typo with his MM4 formula (from MM4-alkene as I abbreviate)[^1] for torsion-stretch. Tinker does not run the full MM4 forcefield, only MM3, but the point is, I've seen more than one typo in his papers previously. An especially severe one for my line of work (crystalline carbon and silicon) is the angle-bending parameter for quaternary Sisp3 (19-19-19). I think it is actually 109.5°, and the paper's 118° is a typo. Your repo contains the potentially faulty 118°.

https://github.com/philipturner/molecular-renderer/commit/6585df047349c698a2b9c70967f6f46946399745

[^1]: Molecular Mechanics (MM4) Calculations on Alkenes* NEYSA NEVINS, KUOHSIANG CHEN, and NORMAN L. ALLINGER Computational Center for Molecular Structure and Design, Department of Chemistry, University of Georgia, Athens, Georgia 30602-2556 Received 25 March 1995; accepted 29 September 1995

jayponder commented 9 months ago

Thanks for your email! After some quick looks at structures and QM calculations on silanes and cyclosilanes, it seems clear that the Si-Si-Si angle parameter is a "typo" in the original Allinger MM3 parameters. (As you note, Lou sent me a file at some point in the early 2000's that I translated into Tinker format.)

I'm sure enough of this typo that I'll replace:

angle 19 19 19 0.250 118.00 110.80 111.20

with the following:

angle 19 19 19 0.250 109.50 110.80 111.20

As to the torsions, it seems much less clear to me whether that's a typo or what Allinger intended. I do know that he "cared" about the silicon parameters, and my suspicion is that what was in his parameter file (and now in Tinker) is what he intended. So, I'm inclined to leave those "torsion5" parameters as is.

Also, do you have a full implementation of MM4? At one time I was very interested in trying to implement MM4 fully in Tinker. But it appeared to me that Allinger never really had enough of a parameterization completed to allow a full force field. And I didn't want to end up guessing at some large percentage of values. At some point, I heard rumors that there was a code around that ran MM4- but it was not freely available, and eventually I lost interest. If you have a full MM4 parameter file, and would be willing to share, I'd like to take a look.

Best wishes, Jay P.

philipturner commented 9 months ago

What I have, is agglomerating all the parameters from his MM4 papers. With MM4, the forcefield is no longer about parameters, but about all the heuristics and the algorithms embedded into those research papers' text. For example, the method for generating Electronegativity Effect and Bohlmann Band Effect corrections to fluorines, with multiple connected carbons, requires reading the paper very carefully to understand. There's also some cross-terms that only apply to specific edge cases, and one must be careful about understanding to which molecules they apply. For example, I dropped the torsion-torsion term meant to fix up fluoro-cyclohexane conformational shapes, as it's likely not transferrable to bulk diamond.

Right now, I'm writing an MM4 forcefield for sp3 H, C, F, and Si (with a restriction of no bonds, angles, or torsions containing F and Si simultaneously). The silicon parameters come from MM3, but using some updated algorithms or combination heuristics from the MM4 papers. I'm also using some Morse parameters that someone got for MM3 using CCSD(T), although the paper explained that the zero-point energy is combinable with the original MM3's bond lengths and stiffnesses. So I'm inserting the MM4 bond lengths and stiffnesses there.

I suspect that the MM4 program was actually MM3, but with a fraction of the parameters updated using the MM4 papers. Allinger wrote the MM3 silicon paper with MM4 in mind, but didn't see the need to make MM4 cross-terms for it. The data set was so small, using MM4 cross-terms might be "overfitting" (in the sense used in Machine Learning).

https://github.com/philipturner/molecular-renderer/tree/main/Sources/MM4/Parameters

As to the torsions, it seems much less clear to me whether that's a typo or what Allinger intended. I do know that he "cared" about the silicon parameters, and my suspicion is that what was in his parameter file (and now in Tinker) is what he intended. So, I'm inclined to leave those "torsion5" parameters as is.

After some quick looks at structures and QM calculations on silanes and cyclosilanes

Do you have any results from quantum chemistry calculations showing the shape of cyclosilanes? Specifically the -C-C-C-C-Si- ring and the -C-C-Si-C-Si- ring. Perhaps if the atom nuclei could be exported as raw 3D coordinates into a text file, so I can measure some angles.

jayponder commented 9 months ago

Below are the atomic coordinates for minimum energy structures at the MP2/6-311G(1d,1p) level for the -C-C-C-C-Si- ring and the -C-C-Si-C-Si- ring. Of course, it's dangerous to fit force field parameters to a single structure :) I'm going to close this issue here on GitHub. If you have other questions or comments, email me at ponder@dasher.wustl.edu.

-C-C-C-C-Si- Ring

     C            0.374173574134    -0.038215551663     0.029145704021
     C            0.042208845189     0.048173538152     1.534234181840
     C            1.361498934513    -0.018980953310     2.315643962314
     C            2.312087681507     1.052449628780     1.739994005084
     Si           1.963073497412     0.983260933736    -0.121194257135
     H           -0.450970884385     0.305930079906    -0.599725255994
     H            0.578396582086    -1.080456385801    -0.241414542186
     H           -0.650101611915    -0.740842053665     1.846013871606
     H           -0.442582735493     1.009311463738     1.748675237095
     H            1.801521737141    -1.014502527347     2.173832182598
     H            1.199157166904     0.113634248537     3.390401627570
     H            3.354868454340     0.875450216430     2.015780146801
     H            2.033349111174     2.037454902073     2.131459294476
     H            3.023113895097     0.298365207660    -0.894804602536
     H            1.750548752298     2.319091252776    -0.722552555554

-C-C-Si-C-Si- Ring

     C            0.436559989542    -0.021509834557    -0.070008725009
     C           -0.175030171743     0.040664578954     1.352031202458
     Si           1.270072281487    -0.070457393573     2.567966751376
     C            2.592470798994     0.941963695435     1.674120078942
     Si           1.923771464832     1.149968957044    -0.082753592247
     H           -0.301082787522     0.215712235450    -0.843118691980
     H            0.790877052560    -1.040136410532    -0.273958643718
     H           -0.928542417136    -0.737929543456     1.507145436186
     H           -0.681533994528     1.003948511432     1.494778642656
     H            1.705503094093    -1.480021075740     2.699020954727
     H            0.939801326509     0.446486514825     3.914300446718
     H            3.559494280506     0.431547810184     1.681437207276
     H            2.737610891177     1.914388357462     2.153497345019
     H            2.912045637587     0.812839717055    -1.130868413736
     H            1.468325553642     2.542657880016    -0.298100998669
philipturner commented 6 months ago

Norman was right. The correct parameter is 118°.

https://github.com/philipturner/molecular-renderer/blob/main/Sources/HardwareCatalog/Materials/NonCarbonElements/README.md#elemental-silicon (permalink)

jayponder commented 6 months ago

Yes, Norman "Lou" Allinger was a really smart and very wise fellow :)