13850
J. Phys. Chem. A 2006, 110, 13850-13856
An All-Atom Force Field for Metallocenes
José N. Canongia Lopes,*,†, P. Cabral do Couto,§ and Manuel E. Minas da Piedade§
Centro de Quı́mica Estrutural, Instituto Superior Técnico, 1049-001 Lisboa, Portugal, and Departamento de
Quı́mica e Bioquı́mica, Faculdade de Ciências, UniVersidade de Lisboa, 1649-016 Lisboa, Portugal
ReceiVed: May 11, 2006; In Final Form: October 20, 2006
A new all-atom force field, for the molecular modeling of metallocenes was constructed. Quantum chemical
calculations were performed to obtain several force field terms not yet defined in the literature. The remainder
were transferred from the OPLS-AA/AMBER framework. The parametrization work included the obtention
of geometrical parameters, torsion energy profiles, and distributions of atomic charges that blend smoothly
with the OPLS-AA specification for a variety of organic molecular fragments. Validation was carried out by
comparing simulated and experimental data for five different ferrocene derivatives in the crystalline phase.
The present model can be regarded as a step toward a general force field for metallocenes, built in a systematic
way, easily integrated with OPLS-AA, and transferable between different metal-ligand combinations.
Introduction
The molecular models currently used in the area of statistical
mechanicssMonte Carlo (MC) and molecular dynamics (MD)
simulation methodssare generally developed from force-fields
containing both intra- and intermolecular sets of parameters.
The latter have a semiempirical nature and are fitted to
condensed-phase structural or thermodynamic data; the former
can be estimated from quantum mechanical calculations or gasphase spectroscopic data.
Databases containing force field parameters for a vast
selection of inorganic and organic molecules are now available,
for instance the OPLS-AA1 and AMBER2 force fields. These
databases include not only the semiempirical parameters defining
the intermolecular interactions between atoms or interacting
centerssgenerally divided into repulsion, dispersion, and electrostatic forcessbut also the parameters controlling the geometry
and energy of each individual molecule. A significant portion
of the existing databases is devoted to organic molecules that
are subunits of larger molecules of biological interest, for
instance amino acids, sugars or DNA bases. On the other hand,
the parametrization of metals or organometallic compounds is
comparatively less developed, especially in terms of a force field
comprising both inter- and intramolecular interactions. In fact,
the area of organometallic molecular modeling exhibits a
particularly pronounced difference between the large number
and scope of quantum chemical or molecular mechanics studies
and their statistical mechanics counterpart.3-5
In this paper we report the development of a simple but
transferable molecular force-field for MD simulations of metallocenes, to be used within the framework of statistical
mechanics. Metallocenes MCp2 (M ) transition metal, Cp )
C5H5) exhibit a sandwich-like structure with the metal coordinated between two cyclopentadienyl fragments (Figure 1) and
constitute a family of molecules whose most well-known
member is ferrocene (FeCp2). They are excellent models for
† To whom correspondence should be addressed. E-mail:jmlopes@
ist.utl.pt
‡ Centro de Quı́mica Estrutural, Instituto Superior Técnico.
§ Departamento de Quı́mica e Bioquı́mica, Faculdade de Ciências,
Universidade de Lisboa.
Figure 1. Labeling scheme of ferrocene and of the ferrocene derivatives
parametrized in this work. The labeling was based on that used in the
AMBER forcefield.2 The centroid of the cyclopentadienyl ring is
denoted by cpc.
systematic studies on how the nature of the metal and the ligand
determines the interactions that govern the intermolecular
potential since a variety of combinations of metals and
substituted Cp rings are possible.6-8
Transferability means that the force-field was built in an
additive way, starting with the parametrization of a simple forcefield for ferrocene, extending it to account for substitution on
the Cp rings in other molecules of the ferrocene family
(dimethylferrocene, decamethylferrocene, acetylferrocene, and
diacetylferrocene), and finally merging the parametrization of
the organic residues (methyl, acetyl, etc.) with the values given
by the widely used OPLS-AA force field (which borrows its
internal parametrization from the AMBER force field). This
strategy allows the future generalization of the force field to
other families of organometallic compounds in a straightforward
and consistent way.
Model Development and Testing
The molecular properties of an isolated ferrocene molecule
were taken as the starting point for the development of the
present force-field. These were obtained through quantum
chemistry calculations performed with the Gaussian 98 program.9 The choice of ferrocene was an obvious one, given the
enormous amount of experimental and theoretical data available
10.1021/jp062896l CCC: $33.50 © 2006 American Chemical Society
Published on Web 12/07/2006
Force Field for Metallocenes
J. Phys. Chem. A, Vol. 110, No. 51, 2006 13851
TABLE 1: Internal Parametrization of Ferrocenea
d(CA-HA)/pm
Tilt(HA)/deg
∆E/kJ‚ol-1
LANL2DZ
142.3 (-1.7)
145.3 (1.3)
143.9 (-0.1)
145.4 (1.4)
144.2 (0.2)
144.3 (0.3)
106.8 (-3.6)
108.8 (-1.6)
107.9 (-2.5)
108.8 (-1.6)
108.1 (-2.3)
108.1 (-2.3)
-1.0 (-4.7)
0.3 (-3.4)
0.5 (-3.2)
0.2 (-3.5)
0.4 (-3.3)
0.3 (-3.4)
0.6 (-3.2)
4.7 (0.9)
3.2 (-0.6)
4.4 (0.6)
3.4 (-0.4)
3.0 (-0.8)
188.1 (22.1)
165.4 (-0.6)
166.4 (0.4)
165.9 (-0.1)
166.5 (0.5)
170.1 (4.1)
SDDall
141.4 (-2.7)
144.4 (0.4)
143.0 (-1.0)
144.4 (0.4)
143.3 (-0.7)
143.4 (-0.6)
106.7 (-3.7)
108.6 (-1.8)
107.8 (-2.6)
108.6 (-1.8)
107.9 (-2.5)
107.9 (-2.5)
-0.9 (-4.6)
0.3 (-3.4)
0.4 (-3.3)
0.2 (-3.5)
0.4 (-3.3)
0.3 (-3.4)
-0.1 (-3.9)
3.8 (0.0)
2.6 (-1.2)
3.6 (-0.2)
2.8 (-1.0)
2.1 (-1.7)
223.5 (17.1)
207.0 (0.6)
207.3 (0.9)
207.4 (1.0)
207.4 (1.0)
210.1 (3.7)
187.9 (21.9)
166.1 (0.1)
167.3 (1.3)
166.5 (0.5)
167.3 (1.3)
170.5 (4.5)
SDD
142.3 (-1.7)
145.3 (1.3)
143.8 (-0.2)
145.4 (1.4)
144.1 (0.1)
144.3 (0.3)
106.8 (-3.6)
108.9 (-1.5)
107.9 (-2.5)
108.9 (-1.5)
108.1 (-2.3)
108.1 (-2.3)
-0.5 (-4.2)
0.4 (-3.3)
0.7 (-3.0)
0.4 (-3.3)
0.6 (-3.1)
0.4 (-3.3)
0.7 (-2.9)
5.8 (2.0)
4.1 (0.3)
5.4 (1.6)
3.4 (-0.4)
3.9 (0.1)
206.4 ( 0.3
166.0
d(Fe-CA)/pm
d(Fe-cpc)/pm
HF
PW91PW91
MPW1PW91
BPW91
B3PW91
B3LYP
225.7 (19.3)
208.5 (2.0)
209.1 (2.7)
209.0 (2.6)
209.1 (2.7)
212.0 (5.6)
190.5 (24.5)
167.8 (1.8)
169.5 (3.5)
168.4 (2.4)
169.4 (3.4)
172.8 (6.8)
HF
PW91PW91
MPW1PW91
BPW91
B3PW91
B3LYP
223.2 (16.8)
206.0 (-0.4)
206.1 (-0.3)
206.4 (0.0)
206.3 (-0.1)
209.4 (3.0)
HF
PW91PW91
MPW1PW91
BPW91
B3PW91
B3LYP
ferrocene
ferrocene
dimethylferrocene
decamethylferrocene
acetylferrocene
diacetylferrocene
Electron Diffraction Data19
144.0 ( 0.2
X-ray Diffraction Data20-24
141 ( 3
142 ( 2
142 ( 1
143 ( 2
142 ( 2
203 ( 3
204 ( 1
205 ( 1
205 ( 2
205 ( 1
206
d(CA-CA)/pm
166
Selected Force-Field Parameters
144
110.4 ( 0.6
3.4
103 ( 5
96 ( 1
1.0 ( 1.5
0.3 ( 0.1
101 ( 6
96 ( 3
6(4
1(1
108
0.0
3.8 ( 1.3
3.6
a
Comparison of DFT results obtained for ferrocene with X-ray and electron diffraction data for ferrocene and its derivatives. The last column
contains the energy differences, ∆E, between the staggered and eclipsed conformations of ferrocene in kJ mol-1. All distances are in picometers
and the HA tilt angle (deviation from the Cp ring plane toward the iron atom) is in degrees. The values in parentheses represent the differences
between the computed bond distances or conformational energies and the corresponding electron diffraction values.
for this particular molecule. In view of transferability, all the
approximations subsequently made in the development of the
force field avoided ferrocene-specific assumptions.
Quantum chemistry calculations are relevant for the construction of a consistent molecular force-field in three aspects: (i)
definition of molecular equilibrium geometries; (ii) parametrization of the internal force-field using conformational energy
profiles; (iii) allocation of partial charges on the interaction
centers (ESP charges), through a fit to the molecular electrostatic
potential. The strategies and approximations used to combine
the computed quantum chemical information with the already
available OPLS-AA parametrization in the construction of the
present force field are described in the next sections.
Equilibrium Geometry. DFT Calculations. It is a wellknown fact that metal-ligand bond distances are overestimated
when quantum chemistry calculations do not take electron
correlation into account.10 For this purpose, density functional
theory (DFT) has proven to be a reliable, economical, and
practical approach and was therefore used in the present work.3,11
In order to choose an adequate level of theory for molecular
structure prediction, the geometry of a single ferrocene molecule
with the two Cp rings in eclipsed (D5h) and staggered (D5d)
conformations was optimized by using PW91, BPW91, and
hybrid MPW1PW91, B3PW91, and B3LYP exchange-correlation functionals12-15 and LANL2DZ or SDD effective core
potentials and split valence basis sets.16-18
The iron-carbon, d(Fe-CA), iron-Cp ring centroid, d(Fecpc), carbon-carbon, d(CA-CA), and carbon-hydrogen d(CAHA) bond distances, and the tilt angle of the ring hydrogens,
Tilt(HA) (see Figure 1 for labeling scheme of the molecule,
consistent with the OPLS-AA/AMBER nomenclature) computed
by DFT are compared in Table 1 with the corresponding
Hartree-Fock theory (HF) predictions and with experimental
electron diffraction (ED)19 and X-ray diffraction (XRD) data.20-24
Also included in Table 1 (last column) are the energy differences, ∆E, between the staggered (D5d) and eclipsed (D5h)
conformations of ferrocene. As expected, important deviations
are found only at the HF level of theory, which overestimates
the Fe-CA or the Fe-cpc distances by ca. 10% and underestimates the experimental ∆E obtained by electron diffraction
by 3.0-3.8 kJ‚mol-1. The results of the DFT methods are rather
similar, with the B3LYP functional consistently giving slightly
larger differences on the bond distances than the remaining
functionals (ca. 1.5%). On average the best geometrical predictions seem to be those obtained by the BPW91/SDDall level of
theory (SDD effective core potentials and triple-ζ valence basis
sets on all heavy atoms). This theoretical approach was therefore
selected to define the set of structural parameters (bond distances
and tilt HA) given at the bottom of Table 1, which were
subsequently used to establish the equilibrium geometry of all
the molecules studied in this work.
Intramolecular Force-Field. In the parametrization of the
intramolecular force-field the Cp rings of ferrocene, dimethylferrocene, decamethylferrocene, acetylferrocene, and diacetylferrocene were taken as regular pentagons with the C-C and
C-H bond distances constrained to 144 and 108 pm, respectively (bottom of Table 1). This type of approximation is
common in other force-fields used in MD simulations (for
13852 J. Phys. Chem. A, Vol. 110, No. 51, 2006
Lopes et al.
Figure 2. Internal degrees of freedom (internal rotation of the cyclopentadienyl rings and torsion of methyl groups) parametrized in this work.
example, heterocyclic aromatic compounds such as pyrrole,
furan, diazole or oxazole are also modeled as rigid units in the
OPLS-AA force-field).25 In the present case it was further
supported by the experimental data in Table 1: (i) The size of
the ring and its distance from the metal are affected by the
environment of the moleculeseither isolated or in a crystal
lattice. This is apparent in the differences between the corresponding bond lengths obtained by electron diffraction, which
refer to the gas phase19 and by X-ray diffraction, which refer
to the crystalline state.20-24 However, the distortions are such
that the ring geometries in five completely different crystals
(ferrocene, dimethylferrocene, decamethylferrocene, acetylferrocene, and diacetylferrocene) are comparable and still within
the range of the values found for the gaseous ferrocene molecule.
(ii) The distortion of the ring caused by different substituent
groups is small, which means that the use of a symmetrical ring
geometrysa regular pentagonsconstitutes a good approximation. The parameters for the bond distances and angles of the
methyl and acetyl substituents were directly taken from OPLSAA.1
As a further simplification, the bond distances and angles
involving the iron atom and the centroids of the Cp rings were
held constant at 206 pm (Table 1), leaving the rotations of the
rings around their centroids as the only possible degrees of
freedom in the ferrocene moleculesthe so-called internal
rotations (Figure 2, cf. below). In the case of dimethyl and
decamethylferocene the motions associated with the torsion of
the methyl groups were also considered. In the acetyl substituents only the torsions of the methyl groups were taken into
account, since the carbonyl group is always coplanar with the
ring (CA-CA-C*-O dihedral angles of 0 or 180°) due to
resonance stabilization.
Internal Rotations. DFT Calculations. The torsion potentials
corresponding to the dihedral angles between the two Cp rings
of ferrocene and dimethylferrocene were obtained from BPW91/
SDDall potential energy scans, where a series of values for those
dihedral angles were held constant while the remaining internal
coordinates of the molecule were allowed to relax to energy
minima (Figure 3a).
In the case of ferrocene the torsion potential energy profile
obtained at the BPW91/SDDall level of theory was decomposed
Figure 3. Torsional energy profiles of ferrocene (dotted lines) and
dimethylferrocene (solid lines): (a) BPW91/SDDall results. (b) MD
results. The thin and thick solid lines refer to a dimethylferrocene
forcefield built using the ferrocene dihedral function plus nonbonded
contributions with different charge distributions in the Cp rings (see
text).
in 25 contributions from the 25 possible (and identical) CAcpc-cpc-CA dihedral angles and described by a simple cosine
function of the type Etorsion ) A[1 + cos(Bφ + C)], where φ is
the dihedral angle and A, B, and C are the parameters defining
the amplitude, frequency and phase of the profile, respectively.
The B and C parameters have values of 5 and 180°, respectively,
due to the 5-fold symmetry of the Cp rings and the fact that the
eclipsed conformations (where all dihedral angles are integer
multiples of 72°) correspond to energy minima. The A parameter
was fitted to a value of 0.072 kJ‚mol-1, corresponding to an
energy barrier of 3.6 kJ‚mol-1, cf. Figure 3a (dotted line) and
Table 1 bottom. In the case of ferrocene the contribution to the
torsion potential energy due to nonbonded interactions was
found to be negligible: a MD simulation with a single ferrocene
Force Field for Metallocenes
J. Phys. Chem. A, Vol. 110, No. 51, 2006 13853
TABLE 2: Point Charges and Van Der Waals Parameters for Ferrocene, Dimethylferrocene, and Acetylferrocene
Fe
DFT (ferrocene)
DFT (dimethylferrocene)a
DFT (acetylferrocene)a
OPLS charges
this force-field
CA
CT
C*
HA
Partial Atomic Point Charges (All in Percent of Atomic Charge Unit)
8.7((0.7)
-10.0((0.1)
9.2((0.1)
9.1((0.5)
-10((25)b
-25((2)
10.5((0.9)
8.2((0.6)
-9((5)
-20((2)
43((1)
9.9((0.6)
11.5
-18.0
47.0
11.5
+10
-10
-18
+47
+9
HC
O
6.9((0.6)
6((1)
6.0
+6
-46.2((0.1)
-47.0
-47
van der Waals Interactions
ǫ/kJ‚mol-1
σ/Å
2.016c
3.11c
OPLS Parametrization (Lennard-Jones (12-6) Potential)1
0.293
0.276
0.439
0.126
3.55
3.50
3.75
2.42
0.126
2.50
0.879
2.96
a Geometry computed at the BPW91/SDDall level of theory. b Large charge differences due to inductive and hyper-conjugation effects in the Cp
ring (see text). c Lennard-Jones parameters fitted to the corresponding Buckingham potential curve.
molecule yielded a null torsion energy profile when the only
contributions to that energy came from nonbonded interactions
(i.e., the A parameter is set to zero).
The dimethylferrocene torsion potential energy profile posed
a more complex problem. In this case different types of CAcpc-cpc-CA dihedral angles are present (the CA carbons on
the ring can be distinguished by their position relative to the
methyl group) and, in addition, the contribution of nonbonding
interactions to the torsion energy (steric effects involving the
methyl groups and different charge distribution around the Cp
rings) is no longer negligible. This means that the presence and
nature of the nonbonded interactions must also be taken into
account to describe the differences between the torsion potential
energy profiles of ferrocene and dimethylferrocene. In fact, in
the OPLS-AA force-field, atoms separated by three bonds within
the same molecule (a dihedral connection) also interact by
Lennard-Jones and electrostatic potentials scaled by a factor of
0.5.
In order to test the importance of the nonbonded interactions
in the definition of the torsion potential energy profile of
dimethylferrocene several MD simulation runs on single ferrocene and dimethylferrocene molecules were performed. In
these runs, the total configuration energy of each molecule was
computed for a series of rotation dihedral angles which were
kept constant throughout the simulation. The dihedral functions
used were those found for ferrocene by the BPW91/SDDall
method and the nonbonded interactions were described by the
sets of atomic partial charges and van der Waals parameters
reported in the following sections. Since other dihedral angles
as well as all bond distances and bond angles are constrained
in the present model, the profile of the configuration energy as
a function of the internal rotation dihedral angle represents the
sum of the dihedral function term plus the contribution from
the nonbonded interactions. The obtained MD results are shown
in Figure 3b, which also includes the BPW91/SDDall profile
(dotted line) obtained for ferrocene for comparison purposes.
It is clear from Figure 3, parts a and b, that the nonbonded
interactions can account for a large part of the differences
between the torsion potential energies of ferrocene (dotted line)
and dimethylferrocene (thick solid line): deviations between
the two profiles are always in the same direction as those found
in the DFT calculations. The main departure between the two
sets of results is in the energy minimum at 0°, which corresponds
to a higher value in the case of the MD calculation since the
molecule is not allowed to optimize its geometry by tilting the
two Cp rings. It should be noted that, since the two profiles are
obtained using different energy minima as reference, their
superimposition and comparison provides only a semiquantitative insight on the magnitude of the nonbonded interactions
contribution to the overall torsional energy.
If the force-field is slightly modified to account for the
different charge distribution among the carbon atoms due to
the inductive effect of the methyl group on the cyclopentadienyl
ring of dimethylferrocene, then the thin solid line in Figure 3b
is obtained. Here the differences are much larger, showing the
sensitivity of the torsion potential energy profile to the nature
and intensity of the nonbonded interactions.
It must be stressed at this point that the MD calculations on
dimethylferrocene were performed keeping the methyl groups
in each Cp ring in a fixed orientation that minimizes steric
repulsions (holding the HC-CT-CA-Fe dihedral angle at 180°).
This is also the lowest energy conformation found in the DFT
calculations. If the methyl groups are allowed to rotate and/or
other methyl or acetyl groups are added to the molecule, the
profiles will suffer further modifications due to the presence of
more complex and larger nonbonded contributions. In other
words, the deviations from the 5-fold symmetric torsional profile
of ferrocene are mainly motivated by nonbonded interactions
between ring substituents.
On the basis of the previous conclusion and in order to keep
the parametrization of the present force-field as simple and
general as possible, the parameters of the internal rotation
dihedrals of ferrocene obtained from DFT calculations were used
for all ferrocene derivatives and the differences between the
corresponding torsion potential energy profiles were assigned
exclusively to the existence of diverse and non-negligible
nonbonded interactions taken from the OPLS-AA force field
(see below).
Partial Charges. Atomic point charges (ESP charges)
presented in Table 2 were determined at the BPW91/6-311G(3df,3pd) level of theory, through a fit to the molecular
electrostatic potential, using the CHelpG procedure15 and the
BPW91/SDDall equilibrium geometry. The results for ferrocene
correspond to Boltzmann averages of eclipsed and staggered
equilibrium conformations (the amplitudes of charge variation
due to the differences in conformation are represented by the
numbers in parentheses in Table 2). In the case of dimethylferrocene and acetylferrocene the Boltzmann averages of the
charges corresponding to minima and maxima of conformational
energy (dimethylferrocene: CT-cpc-cpc-CT dihedrals of 0,
36, 72, 108, 144, and 180°; acetylferrocene: CA-cpc-cpcCA dihedrals of 0 and 36°), were first computed. The obtained
results were then averaged out over all CA or HA atoms, yielding
the values in Table 2 (in this case the values in parentheses
reflect the overall charge variation).
The partial charge distribution in the different ferrocene-like
molecules does not exhibit large fluctuations, the only exception
being the charge distribution between the aromatic carbons, CA,
in dimethylferrocene. The final parametrization of the atomic
point charges (Table 2, values in bold) was established taking
13854 J. Phys. Chem. A, Vol. 110, No. 51, 2006
Lopes et al.
TABLE 3: Density and Crystallographic Properties of Ferrocene Derivatives. Comparison between X-ray Diffraction (XRD)
and MD Simulation Data
ferrocene
dimethylferrocene
decamethylferrocene
space group
ref 20
14 (P21/a)
Crystallographic Properties
ref 21
ref 22
14 (P21/c)
64 (Cmca)
no. molecules
no. cells
cutoff distance/Å
T/K
280
4×5×7
16
298
Simulation Details
240
144
3×5×4
3×3×4
16
16
173
298
a/Å
b/Å
c/Å
β/deg
V/Å3
F/kg‚dm3
∆F/%
XRD
MD
XRD
MD
XRD
MD
XRD
MD
XRD
MD
XRD
MD
MD
10.53
10.63 ( 0.07
7.604
7.64 ( 0.02
5.921
5.94 ( 0.03
121.05
121.35 ( 0.07
406.166
410 ( 1
1.521
1.506 ( 0.004
-0.9
Crystallographic vs Simulated Data
12.19
15.21
12.18 ( 0.05
15.4 ( 0.1
7.466
11.887
7.56 ( 0.03
12.2 ( 0.1
10.839
9.968
11.04 + 0.04
9.52 ( 0.06
103.25
90
102.6 ( 0.2
87 ( 2
960.203
1802.227
986 ( 4
1780 ( 9
1.481
1.203
1.442 ( 0.006
1.217 ( 0.006
-2.6
1.2
into account (i) the average values between different conformations of the same molecule and between different ferrocenelike molecules, (ii) rounded values that do not alter the value
of the point charge more than 0.05 atomic charge units (acu),
and (iii) similarity with the partial charges used in the OPLSAA force-field for analogous organic residues.
In the case of dimethylferrocene the same CA values obtained
for the other ferrocene molecules and OPLS values for the
methyl carbon were selected. This represents a deviation of ca.
-0.27 acu in the aromatic carbon directly connected to the
methyl group, +0.09 acu in the two aromatic carbons adjacent
to it, and +0.07 acu in the aliphatic carbon. To take this fact
into account, a new set of charges for each new type of
substituted ferrocene would have to be calculated ab initio,
which would mean the loss of the transferability and robustness
of the model under development. Nevertheless, any future
refinements of the present model should address this problem,
taking into account polarization effects on the Cp rings or the
possibility of conjugation between the ring and its adjacent
groups.
Intermolecular Repulsion and Dispersion Forces. The
repulsion and dispersion terms of the intermolecular potential
were modeled via a Lennard-Jones (12-6) potential with the ǫ
and σ parameters taken from the OPLS-AA force-field. When
dealing with crystalline structures, the Buckingham potential
is usually preferred to the Lennard-Jones potential (as implemented in the OPLS-AA force-field), due to a better description
of repulsive forces using an exponential term. The LennardJones (12) term is generally considered too “soft” to correctly
describe the repulsion forces in tightly packed structures.4 On
the other hand the OPLS-AA force-field is more complete as
regards the parametrization of the van der Waals interactions
involving the organic residues present in the metallocene
molecules (the Lennard-Jones parameters in this force field were
obtained by fitting to condensed-phase thermodynamic properties of a wide variety of organic molecules). The main objective
of using the LJ potential in the present simulations was to check
if a liquid-oriented force-field (OPLS stands for optimized
potentials for liquid simulations) can also be used in the
simulation of crystalline phases within the scope of statistical
mechanics studies. The values of the Lennard-Jones (12-6) are
acetylferrocene
diacetylferrocene
ref 23
14 (P21/c)
ref 24
14 (P21/c)
192
2×6×2
16
298
216
6×3×3
16
298
20.595
20.8 ( 0.1
5.79
5.82 ( 0.03
18.84
19.2 ( 0.1
116.87
116.7 ( 0.1
2004.024
2061 ( 8
1.512
1.470 ( 0.006
-2.8
5.898
6.02 ( 0.04
13.036
13.1 ( 0.1
14.962
15 ( 0.1
90.68
91.77 ( 0.06
1150.292
1185 ( 7
1.560
1.514 ( 0.009
-2.9
given in Table 2. The cross interaction parameters (not shown
in the table) were calculated using the geometrical mean rule
for the parameters ǫ and σ.
Validation
The validation of the proposed force-field was based on
simulation results obtained using the molecular dynamics
technique, implemented with the DL_POLY software.26 Since
the most common condensed-phase experimental data for
metallocenes come from X-ray diffraction studies, the performance of the proposed force-field was tested by estimating the
unit cell parameters and density of ferrocene, dimethylferrocene,
decamethylferrocene, acetylferrocene, and diacetylferrocene.
The simulation boxes and initial configurations were set
taking into account the dimensions and occupancy of the unit
cells of each crystalline structure selected from the Cambridge
Structural Database (CSD).20-24 Since the dimensions of the
unit cells of these crystals are too small to accommodate a
sufficiently large cutoff distance, several cells were stacked
together to build a sufficiently large and well proportioned
simulation box, with cutoff distances of 16 Å. Long-range
corrections were applied beyond the cutoff distance and the
Ewald sum method (as implemented by the DLPOLY algorithm26 with a precision of 10-3) was implemented to take into
account the long-range character of the electrostatic interactions.
The initial position, orientation, and conformation of each
molecule within the simulation box were defined by the lattice
coordinates taken from the CSD. Since the overall size of the
simulation box is defined by the dimensions of the unit cell of
each crystal, simulations with different box sizes and cutoff
distances were run to check that the dimensions of the simulation
box/cutoff were sufficiently large to make negligible any finite
size effects (comparisons from results obtained from simulations
with cutoff distances of 10, 14, and 16 Å yielded lattice
parameters with precision within a few percent). The simulation
details for each crystal are summarized in Table 3 for the larger
cutoff distances.
The simulations were performed using a Nosé-Hoover
thermostat coupled with an anisotropic Hoover barostat that
Force Field for Metallocenes
J. Phys. Chem. A, Vol. 110, No. 51, 2006 13855
allowed the simulation box to change volume and shape under
N-p-T conditions. The temperature was fixed to match those
used during the crystallographic experiments while the pressure
was set to 1 bar. All runs were allowed to equilibrate for a period
of 100 ps, followed by production times of 200 ps. These
simulation times were found to be appropriate since the initial
configurations are close to the equilibrium structure and it is
observed that the relaxation is complete before the end of the
equilibration period.
It should be noted that the objective of these simulations is
not to test the ability of the present force-field to generate the
corresponding experimental crystal lattice but simply to check
if they are compatible. The stringency of the test was confirmed
with simulation runs where ad-hoc parameters were introduced
and large distortions of the unit cell parameters of the lattice
were observed. For instance, if the spacing between the two
Cp rings is increased by 10% the volume increase of the crystal
can be as high as 6% and the shifts in the cell parameters even
higher and with different signs. The use of an anisotropic
barostat coupled with the system allows for a relatively short
equilibration period due to the frequent rescaling of the position
of the particles. The rather short but effective relaxation time
was also confirmed by monitoring the length of the sides and
angles of the simulation box as the simulation proceeded through
the equilibration period.
Conclusion
The density results are shown in Table 3 and exhibit relative
deviations from the corresponding experimental values of the
order of a few percent (in most cases underestimating the density
by 2-3%). These deviations are of the same order of magnitude
as those obtained by other authors when comparing the
performance of a given force-field against experimental density
data.1 In our case we used the same simple force-field to
calculate the density of five different metallocenes, namely
ferrocene, dimethylferrocene, decamethylferrocene, acetylferrocene, and diacetylferrocene. The results are consistent within
this family of compounds. The agreement is very good
considering that the calculations are purely predictive: all
parameters used were either directly taken from the OPLS-AA
force-field or obtained from DFT calculations; none was
adjusted to match experimental data.
The structural properties of the monoclinic crystals considered
in the present study were also correctly predicted by the models.
After relaxation, the experimental unit cell dimensions were
reproduced with uncertainties of a few tenths of Å in the case
of the a, b, and c lengths and up to one degree for the β angle.
While the volumetric behavior of the ferrocene family is
correctly captured by the present model, the properties closely
related to its energetics (e.g., sublimation enthalpy) will be
analyzed in future work. Nevertheless, to confirm the validity
of the present model in terms of its energy-dependent parametrization, the standard molar enthalpy of sublimation of
ferrocene at 298.15 K was calculated as ∆subH°m (FeCp2) ) 76
( 3 kJ‚mol-1, from equation:
∆subH°m ) ∆subU°m + ∆(PV) ) -U°conf + RT
(1)
where U°conf represents the configurational internal energy of
the crystal at 1 bar and 298.15 K. This result is in very good
agreement with the value ∆subH°m (FeCp2) ) 73.4 ( 1.1
kJ·mol-1, at 298.15 K, currently recommended for the use of
ferrocene as a standard reference material for enthalpy of
sublimation measurements.27,28
The obtained results indicate that the present model can be
regarded as a step toward a general and simple force-field for
metallocenes, built in a coherent way, easily integrated with
the OPLS-AA force-field, and transferable within significantly
different members of the ferrocene family. The extension of
this proposed DFT/MD methodology to metallocenes of other
transition metals, in conjunction with accurate enthalpy of
sublimation measurements for validation of the energy-dependent parametrization is currently in progress.
Supporting Information Available: Table 1S, parameterization of substituted ferrocenes, and Table 2S, ferrocene FIELD
file (input of DLPOLY algorithm). This material is available
free of charge via the Internet at http://pubs.acs.org.
References and Notes
(1) (a) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem.
Soc. 1996, 118, 11225. (b) Kaminski, G.; Jorgensen, W. L. J. Phys. Chem.
1996, 100, 18010.
(2) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K.
M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman,
P. A. J. Am. Chem. Soc. 1995, 117, 5179. Parameters obtained from file
parm99.dat corresponding to AMBER versions 1999 and 2002.
(3) Fey, N. J. Chem. Technol. Biotechnol. 1999, 74, 852.
(4) (a) Theoretical Aspects and Computer Modelling of the Molecular
Solid State; Gavezzotti, A., Ed.; John Wiley: New York, 1997. (b)
Computational Organometallic Chemistry; Cundari, T. R., Ed.; Marcel
Dekker Inc.: New York, 2003.
(5) Shi, S.; Yan, L.; Yang, Y.; Fisher-Shaulsky, J.; Thacher, T. J.
Comput. Chem. 2003, 24, 1059.
(6) Cambridge Structural Database. Allen, F. H.; Kennard, O. Chem.
Design Automation News 1993, 8, 1 and 31-37.
(7) Metallocenes; Togni, A., Alterman, R. L., Eds.; Wiley: New York,
1998.
(8) ComprehensiVe Organometallic Chemistry III; Crabtree, R., Mingos,
M., Eds.; Elsevier: Amsterdam, 2006; Vol. 4-8.
(9) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb,
M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.;
Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A.
D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi,
M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.;
Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.;
Salvador, P.; Dannenberg, J. J.; Malick, D. K.; Rabuck, A. D.; Raghavachari,
K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. G.; Stefanov,
B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.;
Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.;
Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen,
W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; Head-Gordon, M.; Replogle,
E. S.; Pople, J. A. Gaussian 98 (ReVision A.10); Gaussian, Inc.: Pittsburgh,
PA, 2001.
(10) Jensen, F. Introduction to Computational Chemistry; John Wiley:
Chichester, U.K., 1999; pp 109-112.
(11) Ziegler, T. Can. J. Chem. 1995, 73, 743.
(12) Perdew, J. P.; Burke, K.; Wang, Y. Phys. ReV. B 1996, 54, 16533.
(13) Becke, A. D. Phys. ReV. A 1988, 38, 3098.
(14) Adamo, C.; Barone, V. J. Chem. Phys. 1998, 108, 664.
(15) Becke, A. D. J. Chem. Phys. 1993, 98, 5648.
(16) Dunning, T. H., Jr.; Hay, P. J. In Modern Theoretical Chemistry;
Schaefer, H. F., III, Ed.; Plenum: New York, 1976; Vol. 3; pp 1-28.
(17) (a) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 270. (b) Hay,
P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 284. (c) Hay, P. J.; Wadt, W.
R. J. Chem. Phys. 1985, 82, 299.
(18) Leininger, T.; Nicklass, A.; Stoll, H.; Dolg, M.; Schwerdtfeger, P.
J. Chem. Phys. 1996, 105, 1052.
(19) Haaland, A.; Nilson, J. E. Acta Chem. Scan. 1968, 23, 2653.
(20) (a) Takusagawa, F.; Koetzle, T. F. Acta Crystallogr. Sect. B: Struct.
Crystallogr. Cryst. Chem. 1979, 35, 1074. (b) Seiler, P.; Dunitz, J. D. Acta
Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1979, 35, 1068.
(c) Seiler, P.; Dunitz, J. D. Acta Crystallogr., Sect. B: Struct. Crystallogr.
Cryst. Chem. 1982, 38, 1741. (e) Brock, C. P.; Yigang, Fu. Acta Crystallogr.,
Sect. B: Struct. Sci. 1997, 52, 928 (CSD ref codes: FEROCE04-06., 13,
24, 27, 29, 31).
(21) Foucher, D. A.; Honeyman, C. H.; Lough, A. J.; Manners, I.;
Nelson, J. M. Acta Crystallogr., Sect. C: Cryst. Struct. Commun. 1995,
51, 1795 (CSD ref code: ZAYDUY).
13856 J. Phys. Chem. A, Vol. 110, No. 51, 2006
(22) Freyberg, D. P.; Robbins, J. L.; Raymond, K. N.; Smart, J. C. J.
Am. Chem. Soc. 1979, 101, 892 (CSD ref code: DMFERR01).
(23) Sato, K.; Konno, M.; Sano, H. Chem. Lett. 1982, 817 (CSD ref
code: BIWKEX).
(24) Palenik, G. J. Inorg. Chem. 1970, 9, 2424 (CSD ref code:
DACFER10).
(25) McDonald, N. A.; Jorgensen, W. L. J. Phys. Chem. B 1998, 102,
8049.
Lopes et al.
(26) Smith, W.; Forester, T. R. The DL_POLY Package of Molecular
Simulation Routines; version 2.12, The Council for The Central Laboratory of Research Councils, Daresbury Laboratory, Warrington, U.K.,
1999.
(27) Sabbah, R., Ed. Reference Materials for Calorimetry and Differential
Thermal Analysis. Thermochim. Acta 1999, 331, 93.
(28) Chickos, J. S.; Acree, W. E., Jr. J. Phys. Chem. Ref. Data 2002,
31, 537.