Wolfram Language Paclet Repository
Community-contributed installable additions to the Wolfram Language
Extending Coverage of the CFF
$CFFPublishedPatternsParameters | association of SMARTS patterns with symbolic parameters and parameter rules with numeric values and units for the published force field |
$CFFAuxiliaryPatternsParameters | association of SMARTS patterns with symbolic parameters and parameter rules with numeric values and units for extensions to the published force field |
CFFParameterize | generate energy lists and assign parameters |
Data and functions 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, non-bond, and cross 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).
There are two collections of both parts of the parameterization. The first collection is for those elements that have appeared in peer reviewed publications (), and the second collection is for extensions that have not undergone peer review (). In general, these two collections should be disjoint (orthogonal). Users can provide their own patterns and parameters via the options to , and more permanently by editing the paclet file CFFAuxilliaryPatternsParameters.wl which contains the definition of .
$CFFAuxiliaryPatternsParameters
$CFFAuxiliaryPatternsParameters
SMARTS Patterns
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}
Only one instantiation of an unsymmetrical SMARTS pattern needs to be included. For example, for the bending of a CC—C bond angle the pattern suffices and inclusion of the reversed pattern is not necessary.
"[CX3]=[CX3]-[CX4]"
"[CX4]-[CX3]=[CX3]"
Data Structure
<| "Patterns" -> <| "Stretch" -> <| ("[*]~[*]" -> {fc, r0})... |>, "Bend" -> <| ("[*]~[*]~[*]" -> {fc, θ0})... |>, "Torsion" -> <| ("[*]~[*]-[*]~[*]" -> {v, p, 1 | -1})... |>, "Twist" -> <| ("[*](~[*])(~[*])=[*](~[*])~[*]" -> {v})... |>, "OutOfPlane" -> <| ("[*](=[*])(~[*])~[*]" -> {fc})... |>, "Nonbond" -> <| ("[*]" -> {lj, r})... |>, "StretchStretch" -> <| ({, } -> fc)... |>, "StretchBend" -> <| ({r0, θ0} -> fc)... |>, "BendBend" -> <| ({atm, , } -> fc)... |>, "BendBendTwist" -> <| ({, } -> fc)... |>, "OopOop" -> <| ({, } -> fc)... |> |>, "Parameters" -> ((par -> [val, units]) | ( :> ))... |>
r0
i
r0
j
θ0
i
θ0
j
θ0
i
θ0
j
fχ
1
fχ
2
par
i
par
j
Examples
The current implementation of auxiliary patterns and parameters allows coverage of methane, cyclobutanes, cyclopropanes, and alkyl-substituted benzene rings. None of these parameters have been optimized against data.
Methane
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. By eliding the auxiliary patterns and parameters and using we can see what diagonal and non-bond patterns and parameters are needed.
CFFParamterize
In[60]:=
Block
=True,
="Patterns","Parameters"{},
[Molecule["C"]]//KeyTake[{"Stretch","Bend","Torsion","Twist","Nonbond","OutOfPlane"}]
 | $AllowSymbolicCFFEnergy |
| $CFFAuxiliaryPatternsParameters |
 | CFFParameterize |
Out[60]=
Stretch{{{1,2},Missing[],Missing[]},{{1,3},Missing[],Missing[]},{{1,4},Missing[],Missing[]},{{1,5},Missing[],Missing[]}},Bend{2,1,3},Missing,Missing[],{2,1,4},Missing,Missing[],{2,1,5},Missing,Missing[],{3,1,4},Missing,Missing[],{3,1,5},Missing,Missing[],{4,1,5},Missing,Missing[],Torsion{},Twist{},Nonbond{},OutOfPlane{}
0,_
0,_
0,_
0,_
ℋ
Blank[θ]
Blank[θ]
0
ℋ
Blank[θ]
Blank[θ]
0
ℋ
Blank[θ]
Blank[θ]
0
ℋ
Blank[θ]
Blank[θ]
0
ℋ
Blank[θ]
Blank[θ]
0
ℋ
Blank[θ]
Blank[θ]
0
So, one stretch parameter and one bend parameter are needed. Assuming that the stretching parameters for methane are the same as for methyl seems reasonable. However, we cannot do that for the bending parameters because the bending reference angle is 109.47°. The SMARTS pattern for the carbon atom in methane is "[CX4H4]", and can be used to extract the patterns and symbolic parameters we chose.
SMARTS pattern and symbolic parameters for methane C—H stretch :
In[126]:=
KeySelect
["Patterns","Stretch"],StringContainsQ["[CX4H4]"]
 | $CFFAuxiliaryPatternsParameters |
Out[126]=
[CX4H4]-[H]{,}
r
0,r
SMARTS pattern and symbolic parameters for methane H—C—H bend:
In[15]:=
KeySelect
["Patterns","Bend"],StringContainsQ["[CX4H4]"]
| $CFFAuxiliaryPatternsParameters |
Out[15]=
[H]-[CX4H4]-[H]{,}
ℋ
λ
λ
0
The numeric stretch parameters can be extracted from the published parameters, and the numeric bend parameters can be extracted from the auxiliary parameters.
Numeric parameter values for methane C—H stretch and H—C—H bend:
In[67]:=
FilterRules
["Parameters"],,//Association//Echo//ReplaceRepeated
["Parameters"]
 | $CFFPublishedPatternsParameters |
""
"r"
""
0,"r"
| $CFFPublishedPatternsParameters |
»
,
r
681.5/(mol)
kcal
th
2
Å
0,r
0,d
Out[67]=
,
r
681.5/(mol)
kcal
th
2
Å
0,r
1.105
Å
Numeric parameter values for methane H—C—H bend:
In[66]:=
FilterRules
["Parameters"],{,}//Association//Echo//ReplaceRepeated
["Parameters"]
| $CFFAuxiliaryPatternsParameters |
"ℋ"
"λ"
"λ"
0
| $CFFPublishedPatternsParameters |
»
,
ℋ
λ
ℋ
α
λ
0
109.471
°
Out[66]=
,
ℋ
λ
79./mol
kcal
th
λ
0
109.471
°
Using the same technique we can see what cross term parameters were needed.
In[70]:=
Block
=True,
="Patterns","Parameters"{},
Molecule["C"],"AdditionalPatterns""Stretch""[CX4H4]-[H]",,"Bend""[H]-[CX4H4]-[H]"{,}//KeyTake[{"StretchStretch","StretchBend","BendBend","BendBendTwist","OopOop"}]
| $AllowSymbolicCFFEnergy |
| $CFFAuxiliaryPatternsParameters |
| CFFParameterize |
""
"r"
""
0,"r"
"ℋ"
"λ"
"λ"
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 are 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:
So, here are the energy lists for methane:
We can check the normal modes to see if we made a good choice for the parameter values.
Herzberg gives the following assignments for the fundamental frequencies (first two columns), and the CFF frequencies are in the last column.
Cyclobutane
Cyclobutane
The published CFF predicts a planar structure for cyclobutane, which is counter to observation. By eliding the auxiliary patterns and parameters we can obtain results with just the published parameters.
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 SMARTS pattern for the carbon atom in cyclobutane is "[C;r4]", and can be used to extract the patterns and symbolic parameters we chose.
SMARTS pattern and symbolic parameters for cyclobutane C—C—C bend:
Numeric parameter values for cyclobutane C—C—C bend:
The energy report shows that bend energy is greatly reduced, as is the torsion energy due to the now somewhat staggered C—H bonds.
Alkyl benzene
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 [Ermer81Aspekte].
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 that have been implemented:
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.
Water
Water
Create a force field that works only for water.
Next, define the substructure patterns for the stretch and bend terms along with their symbolic force constants and reference values, and define the symbolic parameters for the stretch-stretch and stretch-bend cross-terms:
Next, define the nominal numeric values for the parameters:
Now, parameterize a water molecule to check that everything has been taken care of:
Compute the energy: