Wolfram Language Paclet Repository
Community-contributed installable additions to the Wolfram Language
Extending Coverage of the CFF
$ConsistentForceFieldPatterns | association of SMARTS patterns and symbolic parameters for the published force field |
$ConsistentForceFieldParameters | list of rules of symbolic parameters and values for the published force field |
$CFFAuxiliaryPatterns | association of SMARTS patterns and symbolic parameters |
$CFFAuxiliaryParameters | list of rules of symbolic parameters and values |
Data for extending the Ermer-Lifson Consistent Force Field.
The parameterization of the Ermer-Lifson Consistent Force Field (CFF) in Wolfram Language is handled in two parts. The first is comprised of SMARTS patterns and corresponding symbolic parameters for the molecular substructures for the various stretch, bend, torsion, etc. diagonal force field terms, and pairs of symbolic parameter reference values for the off-diagonal terms (cross-terms). The second part is comprised of the symbolic parameters and their numeric values and units ( objects).
The SMARTS patterns are written keeping in mind Einstein's aphorism "Everything should be kept as simple as possible, but no simpler." To that end, the patterns are written with some generality for economy, and made specific enough to exclude molecules that were not intended to be modeled by the force field. Additionally, the rules are ordered so that more specific patterns come before more general ones. For example, the simplest SMARTS pattern for twist about a CC double bond is . However, that pattern will match the double bonds in hexa-2,4-diene and conjugated double bonds are not part of the force field.
"[#6]-[#6]=[#6]-[#6]"
Out[84]=
In[31]:=
Needs["RobertNachbar`ConsistentForceField`"]
Find the substructures using the simplest double bond SMARTS pattern:
In[32]:=
FindMoleculeSubstructure["CC=CC=CC",MoleculePattern["[#6]-[#6]=[#6]-[#6]"],All]
Out[32]=
{11,22,33,44,13,24,35,46}
A better SMARTS pattern is .
"[CX4]-[CX3]=[CX3]-[CX4]"
Find the substructures using the a more specific double bond SMARTS pattern:
In[33]:=
FindMoleculeSubstructure["CC=CC=CC",MoleculePattern["[CX4]-[CX3]=[CX3]-[CX4]"],All]
Out[33]=
{}
And it will still match unconjugated double bonds, as in hepta-2,5-diene.
Out[88]=
Test the SMARTS pattern on unconjugated double bonds:
In[34]:=
FindMoleculeSubstructure["CC=CCC=CC",MoleculePattern["[CX4]-[CX3]=[CX3]-[CX4]"],All]
Out[34]=
{11,22,33,44,14,25,36,47}
The current implementation of auxiliary parameters allows coverage of methane, cyclobutanes, and alkyl-substituted benzene rings. None of these parameters have been optimized against data.
Methane
Methane was not among the compounds used for calibration, and depiction of the C—H stretch and H—C—H bend terms were limited to methyl and methylene groups. We can use to see what parameters are needed.
CFFParamterize
In[35]:=
Block
=True,
[Molecule["C"],"PrimaryParameters"{}]["MetaInformation"]["CFFEnergyLists"]
 | $AllowSymbolicCFFEnergy |
 | CFFParameterize |
Out[35]=
For the stretch term, we chose the same parameters as for methyl, so only the SMARTS pattern had to be included.
SMARTS pattern for methane C—H stretch :
In[36]:=
KeySelect
["Stretch"],StringContainsQ["[CX4H4]"]
 | $CFFAuxiliaryPatterns |
Out[36]=
[CX4H4]-[H]{,}
r
0,r
Numeric parameter values for methane C—H stretch:
In[37]:=
FilterRules
,//Association//Echo//ReplaceRepeated
| $ConsistentForceFieldParameters |
""
"r"
""
0,"r"
 | $ConsistentForceFieldParameters |
»
,
r
681.5/(mol)
kcal
th
2
Å
0,r
0,d
Out[37]=
,
r
681.5/(mol)
kcal
th
2
Å
0,r
1.105
Å
For the bend term we chose the same force constant as for methyl, but the reference angle needed to 109.47°. Therefore we defined new parameters for both.
SMARTS pattern for methane H—C—H bend:
In[38]:=
KeySelect
["Bend"],StringContainsQ["[CX4H4]"]
| $CFFAuxiliaryPatterns |
Out[38]=
[H]-[CX4H4]-[H]{,}
ℋ
λ
λ
0
Numeric parameter values for methane H—C—H bend:
In[39]:=
FilterRules
,//Echo//ReplaceRepeated
| $CFFAuxiliaryParameters |
"ℋ"
"λ"
"λ"
0
| $ConsistentForceFieldParameters |
»
,
ℋ
λ
ℋ
α
λ
0
109.471
°
Out[39]=
,
ℋ
λ
79./mol
kcal
th
λ
0
109.471
°
Using the option in will allow us to see which cross-term parameters were needed.
"AdditionalPatterns"
In[40]:=
Block
=True,
Molecule["C"],"PrimaryParameters"{},"AdditionalPatterns""Stretch""[CX4H4]-[H]",,"Bend""[H]-[CX4H4]-[H]"{,}["MetaInformation"]["CFFEnergyLists"]
| $AllowSymbolicCFFEnergy |
| CFFParameterize |
""
"r"
""
0,"r"
"ℋ"
"λ"
"λ"
0
Out[40]=
Atoms{C,H,H,H,H},Stretch{{{1,2},,},{{1,3},,},{{1,4},,},{{1,5},,}},Bend{{{2,1,3},,},{{2,1,4},,},{{2,1,5},,},{{3,1,4},,},{{3,1,5},,},{{4,1,5},,}},Torsion{},Twist{},OutOfPlane{},Nonbond{},StretchStretch{{{2,1,3},,,},{{2,1,4},,,},{{2,1,5},,,},{{3,1,4},,,},{{3,1,5},,,},{{4,1,5},,,}},StretchBend{2,1,3},,,,{2,1,4},,,,{2,1,5},,,,{3,1,2},,,,{3,1,4},,,,{3,1,5},,,,{4,1,2},,,,{4,1,3},,,,{4,1,5},,,,{5,1,2},,,,{5,1,3},,,,{5,1,4},,,,BendBend{2,1,3,4},,,,{2,1,3,5},,,,{2,1,4,5},,,,{3,1,2,4},,,,{3,1,2,5},,,,{3,1,4,5},,,,{4,1,2,3},,,,{4,1,2,5},,,,{4,1,3,5},,,,{5,1,2,3},,,,{5,1,2,4},,,,{5,1,3,4},,,,BendBendTwist{},OopOop{},Constraint{}
r
0,r
r
0,r
r
0,r
r
0,r
ℋ
λ
λ
0
ℋ
λ
λ
0
ℋ
λ
λ
0
ℋ
λ
λ
0
ℋ
λ
λ
0
ℋ
λ
λ
0
ℱ
r,r
0,r
0,r
ℱ
r,r
0,r
0,r
ℱ
r,r
0,r
0,r
ℱ
r,r
0,r
0,r
ℱ
r,r
0,r
0,r
ℱ
r,r
0,r
0,r
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
ℱ
r,λ
0,r
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
†
ℱ
λ,λ
λ
0
λ
0
Since the reference value parameter for the methane C—H stretch is the same as for methyl, a new stretch-stretch parameter was not needed. Stretch-bend and bend-bend parameters were needed, and we used a value of 0, like the other cross-terms involving H in the published CFF.
Numeric parameter values for methane cross-terms:
In[41]:=
FilterRules
,
| $CFFAuxiliaryParameters |
"ℱ"
"r","λ"
†
"ℱ"
"λ","λ"
Out[41]=
,
ℱ
r,λ
0/(Åmol)
kcal
th
†
ℱ
λ,λ
0/mol
kcal
th
In[42]:=
methane=
Molecule["C"],"AdditionalPatterns"
,"AdditionalParameters"
;methane["MetaInformation"]["CFFEnergyLists"]
| CFFParameterize |
| $CFFAuxiliaryPatterns |
| $CFFAuxiliaryParameters |
Out[43]=
Atoms{C,H,H,H,H},Stretch{{{1,2},681.5,1.105},{{1,3},681.5,1.105},{{1,4},681.5,1.105},{{1,5},681.5,1.105}},Bend{{{2,1,3},79.,1.91063},{{2,1,4},79.,1.91063},{{2,1,5},79.,1.91063},{{3,1,4},79.,1.91063},{{3,1,5},79.,1.91063},{{4,1,5},79.,1.91063}},Torsion{},Twist{},OutOfPlane{},Nonbond{},StretchStretch{},StretchBend{},BendBend{},BendBendTwist{},OopOop{},Constraint{}
In[44]:=
 | CFFEnergy |
Out[44]=
0.223336/mol
kcal
th
In[45]:=
methane=
[methane]
| CFFOptimizeGeometry |
Out[45]=
Molecule
 |  |
|
In[46]:=
| CFFEnergyReport |
Out[46]=
Type | Energy (kcal/mol) |
Stretch | 0 |
Bend | 0 |
Torsion | 0 |
Twist | 0 |
Out-of-Plane | 0 |
Nonbond | 0 |
Stretch-stretch | 0 |
Stretch-bend | 0 |
Bend-bend | 0 |
Bend-bend-twist | 0 |
Oop-oop | 0 |
Constraint | 0 |
Total | 0 |
RMS Gradient (kcal/mol/Å) | 0 |
We can check the normal modes to see if we made a good choice for the parameter values.
In[262]:=
modes=
[methane];
| CFFNormalModes |
In[319]:=
Keys@modes
Out[319]=
{Frequencies,NormalModes,NormalModeLabels,NormalModeSymmetries,InfraredActive,CoordinateLabels}
In[321]:=
| NormalModeTable |
Out[321]=
,
,
S[1] | S[2] | S[3] | |
λ | 2984.94 wavenumbers | 2984.94 wavenumbers | 2984.94 wavenumbers |
Symmetry | C 1 | C 1 | C 1 |
IR active | True | True | True |
r[1,2] | -0.488839 Å | 0.311525 Å | -0.700341 Å |
r[1,3] | 0.85855 Å | 0.288761 Å | -0.0774476 Å |
r[1,4] | -0.354518 Å | 0.308701 Å | 0.778145 Å |
r[1,5] | -0.0151925 Å | -0.908987 Å | -0.000356198 Å |
θ[2,1,3] | -3.43048 ° | -5.56996 ° | 7.21697 ° |
θ[2,1,4] | 7.82537 ° | -5.75497 ° | -0.721928 ° |
θ[2,1,5] | 4.67682 ° | 5.54375 ° | 6.50165 ° |
θ[3,1,4] | -4.67682 ° | -5.54375 ° | -6.50165 ° |
θ[3,1,5] | -7.82537 ° | 5.75497 ° | 0.721928 ° |
θ[4,1,5] | 3.43048 ° | 5.56996 ° | -7.21697 ° |
S[4] | S[5] | S[6] | |
λ | 2823.83 wavenumbers | 1507.02 wavenumbers | 1507.02 wavenumbers |
Symmetry | T d | D 2 | D 2 |
IR active | False | False | False |
r[1,2] | -0.498055 Å | 0 Å | 0 Å |
r[1,3] | -0.498055 Å | 0 Å | 0 Å |
r[1,4] | -0.498055 Å | 0 Å | 0 Å |
r[1,5] | -0.498055 Å | 0 Å | 0 Å |
θ[2,1,3] | 0 ° | -32.0621 ° | 40.4934 ° |
θ[2,1,4] | 0 ° | -19.0373 ° | -48.0133 ° |
θ[2,1,5] | 0 ° | 51.0993 ° | 7.51988 ° |
θ[3,1,4] | 0 ° | 51.0993 ° | 7.51988 ° |
θ[3,1,5] | 0 ° | -19.0373 ° | -48.0133 ° |
θ[4,1,5] | 0 ° | -32.0621 ° | 40.4934 ° |
S[7] | S[8] | S[9] | |
λ | 1345.46 wavenumbers | 1345.46 wavenumbers | 1345.46 wavenumbers |
Symmetry | C 1 | C 1 | C 1 |
IR active | True | True | True |
r[1,2] | -0.00747661 Å | -0.0148425 Å | -0.0272915 Å |
r[1,3] | 0.00745227 Å | -0.0217059 Å | 0.0222338 Å |
r[1,4] | -0.0249871 Å | 0.0172673 Å | 0.00992506 Å |
r[1,5] | 0.0250115 Å | 0.0192811 Å | -0.00486742 Å |
θ[2,1,3] | -0.037151 ° | -55.7737 ° | -7.71807 ° |
θ[2,1,4] | -49.5404 ° | 3.70033 ° | -26.5015 ° |
θ[2,1,5] | 26.7586 ° | 6.7733 ° | -49.0752 ° |
θ[3,1,4] | -26.7586 ° | -6.7733 ° | 49.0752 ° |
θ[3,1,5] | 49.5404 ° | -3.70033 ° | 26.5015 ° |
θ[4,1,5] | 0.037151 ° | 55.7737 ° | 7.71807 ° |
Herzberg gives the following assignments for the fundamental frequencies (first two columns), and the CFF frequencies are in the last column.
Cyclobutane
The CFF predicts a planar structure for cyclobutane, which is counter to observation.
Compute the CFF-optimized structure of cyclobutane:
Compute the symmetry of the result:
Compute the ring bond lengths, bond angles, torsion angles, and out-of-plane pucker angle:
A survey of the literature reveals these structural and energy details for cyclobutane:
The bend energy is too high because the internuclear angles in 4-membered rings are not a good measure of the angle between the bonding orbitals (the so-called banana bonds).
Compute the energy report for the contributions of the various terms:
To remedy this problem, we chose to introduce new symbolic bend parameters for 4-membered rings with their own numeric values. The process used is summarized below.
The energy report shows that bend energy is greatly reduced, as is the torsion energy due to the now somewhat staggered C—H bond.
Alkyl benzene
Ermer developed a set of parameters for alkyl-substituted benzene molecules so that he could study 1,5-dibenzocycloocta-1,5-diene and other similar molecules. These parameters are described in his Habilitationsschrift.
Many of the parameters for alkene double bonds were applied also to the benzene ring, so only a few new parameters were needed. These are the SMARTS patterns needed:
These are the symbolic parameters for the new values introduced by Ermer:
The numeric values are:
Another molecule studied was 1,4,7-tribenzocyclonona-1,4,7-triene:
Ermer did not calibrate these parameters, so they must be used with caution. Geometries are probably acceptable for undistorted benzene rings, as are relative energies. The normal modes and vibrational frequencies are less certain. The table below shows the observed infrared and Raman frequencies for the fundamental modes of benzene and perdeuterio benzene, as reported by Hertzberg, and the corresponding frequencies computed with CFF. The agreement is only qualitatively acceptable.