Find the First and third parts here:​
​1st: Consistent Force Field: alkylbenzenes, alkanes, alkenes. Energy, geometry optimization, vibrations.​
​3rd: Consistent Force Field: alkylbenzenes, alkanes, alkenes. Conformational analysis of molecular gears

Table of Contents


Introduction
In the previous installment we introduced the ConsistenForceField paclet and demonstrated the basic uses of many of the functions. In this post we’ll show how they can be used for conformational analysis applications, in particular mapping conformational space and conformational reaction coordinates of molecular rotors and gears.
Conformational analysis is a broad topic with a long history (almost 150 years [2]), and has many applications in organic, inorganic, and physical chemistry. The earliest applications of conformational analysis studied the distortion from the “ideal” geometry of molecules. For example, cyclohexane can adopt a perfect geometry with all bond angles at the tetrahedral angle of
ArcCos-
1
3
≈109.471°
and the ring C—C—C—C torsion angles of ±60°. The shape of the ring is called a “chair” with the carbon atoms alternately above and below the mean plane of the ring. However, because of van der Waals repulsion between the axial hydrogen atoms the ring is slightly flattened, resulting in somewhat larger C—C—C bond angles and somewhat smaller C—C—C—C torsion angles.
∡
6,1,2
111.034
°
∡
1,2,3
111.034
°
∡
2,3,4
111.034
°
∡
3,4,5
111.034
°
∡
4,5,6
111.034
°
∡
5,6,1
111.034
°
∡
6,1,2,3
-55.9541
°
∡
1,2,3,4
55.9541
°
∡
2,3,4,5
-55.9541
°
∡
3,4,5,6
55.9541
°
∡
4,5,6,1
-55.9541
°
∡
5,6,1,2
55.9541
°
The relative stability of different conformations was also of interest. Methylcyclohexane exists as a mixture of two conformations, one with the methyl group equatorial and the other with the methyl group axial. There is a 0.81 kcal/mol difference in the steric energy of the two conformations.
equatorial methyl
axial methyl
5.42227
kcal
th
/mol
6.23199
kcal
th
/mol
The earliest studies, of course, relied on experimental observations (infrared and Raman vibrational spectroscopy, microwave rotational spectroscopy, electron diffraction, and x-ray diffraction). Later, nuclear magnetic resonance (NMR) spectroscopy became available, and variable temperature NMR experiments allowed the accurate determination of energy barriers. Computational methods began to be applied almost 80 years years ago (see the Historical background section of the earlier post). In addition to giving optimized geometries and computed energies of ground and transition states, computational methods allow us to study the detailed mechanisms (within the limits of the theory) of these dynamic conformational changes. That will be one of the aims of this post.

Symmetry distortion coordinates

In the discussion below we will make use of “symmetry distortion coordinates” that are linear combinations of internal coordinates (torsion angles in our case) that transform appropriately under the operations of the symmetry group of the molecule under study [7, 9, 8, 13]. The molecules we’ll be covering are flexible, so a higher symmetry group that encompasses all the conformations accessible is often needed. The nonrigid symmetry groups of Longuet-Higgins [1] and permutation-inversion groups [4, 5] are examples. In the formulas for the symmetry distortion coordinates, torsion angles defined by four atoms are denoted by
ω
i,j,k,l
where
i
,
j
,
k
, and
l
are atom indices. For symmetrical rotors, such as methyl and phenyl groups, we’ll use a rotational coordinate ϕ for convenience. For 2-fold rotors there are two torsional angles
ω
i,j,k,
l
1
and
ω
i,j,k,
l
2
and we can define
ϕ=
ω
i,j,k,
l
1
+
ω
i,j,k,
l
2
-180°2
were
ω
i,j,k,
l
1
<
ω
i,j,k,
l
2
. For 3-fold rotors there are three torsion angles
ω
i,j,k,
l
1
,
ω
i,j,k,
l
2
, and
ω
i,j,k,
l
3
and we can define
ψ=
ω
i,j,k,
l
1
+
ω
i,j,k,
l
2
+
ω
i,j,k,
l
3
-360°2
were
ω
i,j,k,
l
1
<
ω
i,j,k,
l
2
<
ω
i,j,k,
l
3
. Both ϕ and ψ rotational angles can be computed by
CFFTwoFoldTorsion
and
CFFThreeFoldTorsion
, respectively.
The general practice in structural chemistry is to express bond angles, torsion angles and dihedral angles in degrees instead of the usual radians employed by mathematicians. The functions
MoleculeValue
,
CFFAngle
,
CFFTorsion
,
CFFDihedral
,
CFFTwoFoldTorsion
, and
CFFThreeFoldTorsion
all return results as
Quantity[value, "AngularDegrees"]
. So for consistency, we’ll use degrees in this post.
The symmetry distortion coordinates are denoted as
S
i
where
i
is an integer if it transforms as a 1-dimensional irreducible representation of the symmetry group or as a pair
S
i,a
and
S
i,b
if it transforms as a 2-dimensional irreducible representation. We’ll just be presenting the symmetry distortion coordinates as given, so you won’t have to work with character tables!

Getting started

First, install the ConsistenForceField paclet, and then load it.
In[]:=
PacletInstall["RobertNachbar/ConsistentForceField",UpdatePacletSitesTrue]
In[]:=
Needs["RobertNachbar`ConsistentForceField`"]
One 1-Fold Rotor: n-Butane
n-Butane is one of the molecules first studied by organic chemistry students in the section on conformational analysis. The molecule is not rigid, and relatively free rotation about the C—C bonds occurs at room temperature. The potential energy of the molecule varies as the bonds rotate due to changing interatomic distances (non-bonded or van der Waals strain) and changing overlap of vicinal molecular orbitals (Pitzer strain).
In[]:=
MoleculePlot3D[Molecule["CCCC"],ImageSize->Tiny]
Out[]=

Traditional energy profile method

The rotation of the methyl groups at each of the molecule is less interesting than that of the central C—C bond. As is typically done, we can use
MoleculeModify
to set the C1—C2—C3—C4 torsion to a series of values to sample the potential energy of the molecule. At each setting, the geometry of all the atoms is optimized to minimize the energy, and a penalty constraint is applied to keep the torsion angle close to the set value. The computation is done in parallel to reduce the waiting time. This operation is sometimes call a “potential energy scan.”
In[]:=
LaunchKernels[8];
In[]:=
ParallelEvaluate[Needs["RobertNachbar`ConsistentForceField`"]];
In[]:=
profile=ParallelMap[{#->CFFOptimizeGeometry[MoleculeModify[Molecule["CCCC"],{"SetTorsionAngle",{1,2,3,4}->#}],"Constraints"->{{"Torsion",{1,2,3,4},#,100}}]}&,Range[-180,180,5]]//Association;//EchoTiming
⌚
41.4661
In[]:=
CloseKernels[];
The
CFFEnergy
function can be mapped over the list of molecules and then plotted to see the energy profile. Over the full rotation there are three minima and three maxima (the barriers to rotation).
In[]:=
ListPlot[CFFEnergy/@profile,FrameLabel->{"torsion (°)","CFF energy (kcal/mol)"},PlotLabel->"Constrained energy profile for n-butane torsion"]
Out[]=
The minima occur when the six bonds radiating from the two central carbon atoms are staggered, and the maxima occur when they are eclipsed.
Note that there are two mirror image pairs: center and right in the top row, and left and right in the bottom row. Taking only one member from those pairs and the two unique conformations gives a total of four. Using the spatial orientation of the C2—C1 and C3—C4 bonds, we can call them “anti” (
ω=180°
), “gauche” (
ω≈65°
), “eclipsed” (
ω≈120°
), and “syn” (
ω=0°
).
The energies shown in the profile figure are not entirely accurate because a penalty constraint was used in the geometry optimization which slightly inflates the computed energy and the constrained torsion angle is not exactly at the set value. To address these issues, we’ll use a newer method.

Conformation space mapping

This method is not really new (we used it more than 40 years ago), but it’s use does not seem to have found wide use in the area of EFF calculations. In the area of electronic structure calculations (molecular orbital or quantum mechanics calculations) the intrinsic reaction coordinate (IRC) method [3] is widely known and used. The techniques are related, and we plan to add IRC capability to the ConsistenForceField paclet in the future.
Ideally, we want to study geometries that are stationary, that is, the gradient norm is 0. For potential energy minima these geometries are found with unconstrained minimization of the energy. The
CFFOptimizeGeometry
function uses
FindMinimum
to accomplish this task. For the potential energy maxima minimization only works when the geometry is symmetric [
0
], and even then it sometimes fails. However, if one starts with a geometry that is sufficiently close to the stationary geometry, the gradient norm will be small and finding the root of the gradient vector will locate the precise stationary maximum energy geometry. The
CFFOptimizeGeometry
function uses
FindRoot
to accomplish this task.
In each of the next four subsections one of the selected conformations of n-butane from the potential energy scan is analyzed in detail. It is first constructed de novo, and then geometry-optimized to find the true stationary geometry. The function
CFFGradient
is used to confirm that the gradient norm is negligibly small (usually less than
-6
10
kcal/mol/Å). The symmetry of the conformation is computed with the
ResourceFunction["PointGroupSymbol"]
. The vibrational frequencies are computed with
CFFNormalModes
, and the number of imaginary frequencies is computed, which gives the number of dimensions with negative curvature of the stationary point on the potential energy hypersurface. Minima have 0, transition states have 1, and partial maxima have 2 or more. We then provide a summary as a way point on our adventure. Next, we use the transition vectors at each transition state to discover their adjacent minima, and produce an unconstrained energy profile. Lastly, we produce the conformation space map.

Stationary structures

anti

Construct the anti conformation of n-butane and optimize the geometry:
In[]:=
butaneAnti=Molecule["CCCC"]//MoleculeModify[#,{"SetTorsionAngle",{1,2,3,4}->180}]&//CFFOptimizeGeometry;​​MoleculePlot3D[butaneAnti,ImageSize->Tiny]
Out[]=
Compute the CFF energy and the gradient norm of the unconstrained geometry:
Compute the C1—C2—C3—C4 torsion angle and symmetry of the unconstrained geometry:
Count the number of imaginary vibrational frequencies:
The anti conformation is an energy minimum.

gauche

Construct the gauche conformation using the nominal torsion angle of 60° and optimize the geometry:
Compute the CFF energy and the gradient norm of the unconstrained geometry:
Compute the C1—C2—C3—C4 torsion angle and symmetry of the unconstrained geometry:
Count the number of imaginary vibrational frequencies:
The gauche conformation is also an energy minimum.

eclipsed

Compute the CFF energy and the gradient norm of the unconstrained geometry:
Compute the C1—C2—C3—C4 torsion angle and symmetry of the unconstrained geometry:
Count the number of imaginary vibrational frequencies:
The anti conformation is a 1-dimensional energy maximum, that is, a conformational transition state.

syn

Compute the CFF energy and the gradient norm of the unconstrained geometry:
Compute the C1—C2—C3—C4 torsion angle and symmetry of the unconstrained geometry:
Count the number of imaginary vibrational frequencies:
The syn conformation is also a 1-dimensional energy maximum, that is, a conformational transition state.

Stationary structures conformation map

The coordinate system used for a conformation map depends on the structural variation one is studying. For the torsion driving above, the natural coordinate is the torsion angle that was driven. When mapping unconstrained reaction paths more than one internal coordinate (bond length, bond angle, torsion angle, etc.) may undergo substantial variation. The nonrigid symmetry group of the molecule will allow one to construct a minimal set of symmetry distortion coordinates as linear combinations of several internal coordinates or atomic Cartesian coordinates. Often times, only one or a few of the full set of distortion coordinates are sufficient to describe most of the structural variation, and are therefore very convenient for mapping.
When dealing with molecular rotors, there is a periodic variation of the structure. That is, when the rotor undergoes a full rotation of 360° the structure is brought back to its original conformation. The conformational space is therefore also periodic and possesses translational symmetry. The best known examples are the crystallographic space groups. The mathematics of this kind of space is well studied, and and they can be constructed for any number of dimensions. Buergi was one of the first to recognise and apply space groups to conformational analysis [10].
To construct the conformation map we first place a marker at the appropriate position for each stationary conformation that we have have computed, that is, syn, gauche, eclipsed, and anti. The symmetry operations of the space, both point and translation, are used to make symmetry-related copies. Thus, we need include only one of the two mirror image gauche conformations and one of the two mirror image eclipsed conformations.
Now, generate list of symbols representing the stationary structures for the conformation map:

Unconstrained conformational reaction paths

As already noted, the eclipsed and syn conformations are conformational transition states. The normal modes corresponding to the imaginary frequencies can tell us something about the motion over barrier.

conformational energy profile

Now, joining the three segments into a continuous itinerary, reflecting a copy through the origin, and joining again will give the full itinerary.
This result is similar to the IRC, but there are some important differences that are beyond the scope of the present discussion.

Full conformation map

We can take the itineraries and make a set of paths for the conformation map, making use of the symmetry of the space.
The last element of the conformation map is the legend, showing the relative energy color scale and the labeled glyphs , which is defined in the initialization section at the end of the notebook.
Putting it all together:
Admittedly, the paths don’t add much for this one-dimensional problem, but it was worthwhile going through the steps so that they can be used for the following applications.
Back Matter

References

1
.
Longuet-Higgins, H. C. (1963). The symmetry groups of non-rigid molecules. Molecular Physics, 6(5), 445-460. https://doi.org/10.1080/00268976300100501
2
.
Eliel, E. L., Allinger, N. L., Angyal, S. J., & Morrison, G. A. (1965). Conformational analysis. New York, 200.
3
.
Fukui, K. (1970). Formulation of the reaction coordinate. The Journal of Physical Chemistry, 74(23), 4161-4163. https://doi.org/10.1021/j100717a029
4
.
Ruch, E., Hässelbarth, W., & Richter, B. (1970). Doppelnebenklassen als Klassenbegriff und Nomenklaturprinzip für Isomere und ihre Abzählung. Theoretica chimica acta, 19(3), 288-300. https://doi.org/10.1007/BF00532233
5
.
Hässelbarth, W., & Ruch, E. (1973). Classifications of rearrangement mechanisms by means of double cosets and counting formulas for the numbers of classes. Theoretica chimica acta, 29(3), 259-268. https://doi.org/10.1007/BF00529052
6
.
Ermer, O. (1975). Determination of molecular symmetries by force field calculations and evaluation of symmetric and nonsymmetric conformational transition states avoiding complete point-by-point mapping. Tetrahedron, 31(16), 1849-1854. https://doi.org/10.1016/0040-4020(75)87040-2
7
.
Murray-Rust, P., Bürgi, H. B., & Dunitz, J. D. (1978). Distortions of MX4 molecules from Td symmetry. I. Kernel, co-kernel and averaged configurations. Acta Crystallographica Section B: Structural Crystallography and Crystal Chemistry, 34(6), 1787-1793. https://doi.org/10.1107/S0567740878006718
8
.
Dunitz, J. D. (1979) X-ray Analysis and the Structure of Organic Molecules. Cornell University Press, Ithaca, NY.
9
.
Murray-Rust, P., Bürgi, H. B., & Dunitz, J. D. (1979). Description of molecular distortions in terms of symmetry coordinates. Acta Crystallographica Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography, 35(5), 703-713. https://doi.org/10.1107/S0567739479001686
10
.
Buergi, H. B. (1980) Space group-like symmetries in the conformational analysis of flexible molecules. MATCH, 9, 13-14. https://match.pmf.kg.ac.rs/electronic_versions/Match09/match9_13-14.pdf
11
.
Ermer, O., (1981) Aspekte von Kraftfeldrechnungen, Wolfgang Bauer Verlag, Muenchen.
12
.
Hahn, T., Shmueli, U., & Arthur, J. W. (Eds.). (1983). International tables for crystallography (Vol. A). Dordrecht: Reidel.
13
.
Dunitz, J. D., & Bürgi, H. B. (1994). Symmetry Aspects of Structure Correlation. in HB Bürgi & JD Dunitz Structure correlation. VCH Verlagsgesellschaft mbH.

Initialization

Paclet installation and loading

conformation map legends

CITE THIS NOTEBOOK

Consistent Force Field: alkylbenzenes, alkanes, alkenes. Conformational analysis of molecular rotors​
by Robert Nachbar​
Wolfram Community, STAFF PICKS, June 10, 2024
​https://community.wolfram.com/groups/-/m/t/3190718