IR spectra of phosphate ions in aqueous solution: Predictions of a DFT/MM approach compared with observations

Article abstract of DOI:10.1021/jp048617g

Due to the progress of density functional theory (DFT) accurate computations of vibrational spectra of isolated molecules have become a standard task in computational chemistry. This is not yet the case for solution spectra. To contribute to the exploration of corresponding computational procedures, here we suggest a more efficient variant of the so-called instantaneous normal-mode analysis (INMA). This variant applies conventional molecular dynamics (MD) simulations, which are based on nonpolarizable molecular mechanics (MM) force fields, to the rapid generation of a large ensemble of different solvation shells for a solute molecule. Short hybrid simulations, in which the solute is treated by DFT and the aqueous solvent by MM, start from snapshots of the MM solute-solvent MD trajectory and yield a set of statistically independent hydration shells partially adjusted to the DFT/MM force field. Within INMA, these shells are kept fixed at their 300 K structures, line spectra are calculated from the DFT/MM Hessians of the solute, and its inhomogeneously broadened solution spectra are derived by second-order statistics. As our test application we have selected the phosphate ions HPO₄^2- and H₂PO₄ because sizable solvation effects are expected for the IR spectra of these strongly polarizable ions. The widths, intensities, and spectral positions of the calculated bands are compared with experimental IR spectra recorded by us for the purpose of checking the computational procedures. These comparisons provide insights into the merits and limitations of the available DFT/MM approach to the prediction of IR spectra in the condensed phase.

Full text of DOI:10.1021/jp048617g

6186
J. Phys. Chem. A 2004, 108, 6186-6194
IR Spectra of Phosphate Ions in Aqueous Solution: Predictions of a DFT/MM Approach
Compared with Observations
Marco Kla1hn,† Gerald Mathias,‡ Carsten Ko1tting,† Marco Nonella,‡ Ju1rgen Schlitter,†
Klaus Gerwert,† and Paul Tavan*,‡
Institut fu¨r Biophysik, ND04, Ruhr-UniVersita¨t, 44780 Bochum, Germany, and Lehrstuhl fu¨r Biomolekulare
Optik, Ludwig-Maximilians-UniVersita¨t, Oettingenstrasse 67, 80538 Mu¨nchen, Germany
ReceiVed: March 30, 2004; In Final Form: May 13, 2004
Due to the progress of density functional theory (DFT) accurate computations of vibrational spectra of isolated
molecules have become a standard task in computational chemistry. This is not yet the case for solution
spectra. To contribute to the exploration of corresponding computational procedures, here we suggest a more
efficient variant of the so-called instantaneous normal-mode analysis (INMA). This variant applies conventional
molecular dynamics (MD) simulations, which are based on nonpolarizable molecular mechanics (MM) force
fields, to the rapid generation of a large ensemble of different solvation shells for a solute molecule. Short
hybrid simulations, in which the solute is treated by DFT and the aqueous solvent by MM, start from snapshots
of the MM solute-solvent MD trajectory and yield a set of statistically independent hydration shells partially
adjusted to the DFT/MM force field. Within INMA, these shells are kept fixed at their 300 K structures, line
spectra are calculated from the DFT/MM Hessians of the solute, and its inhomogeneously broadened solution
spectra are derived by second-order statistics. As our test application we have selected the phosphate ions
2-
-
HPO4 and H2PO4 because sizable solvation effects are expected for the IR spectra of these strongly
polarizable ions. The widths, intensities, and spectral positions of the calculated bands are compared with
experimental IR spectra recorded by us for the purpose of checking the computational procedures. These
comparisons provide insights into the merits and limitations of the available DFT/MM approach to the prediction
of IR spectra in the condensed phase.
1. Introduction
Whereas for molecules isolated in the vacuum highly accurate
calculations of vibrational spectra have become simple and
computationally manageable tasks due to the development1,2 and
widespread accessibility3,4 of density functional theory (DFT)
(see, e.g., refs 5-7), the computation of solution spectra is by
no means simple nor standard. Here the question, which
computational procedures are adequate, is the subject of a
current debate.8-11
In principle, accurate descriptions of IR solution spectra have
become accessible with the recent progress of hybrid methods
(see ref 8 for an overview), whose key ideas were suggested
about three decades ago by Warshel and Levitt.12 Here, the
solute molecule is treated quantum mechanically by DFT and
is called the DFT fragment of the simulation system. The solvent
environment is microscopically modeled by a conventional
nonpolarizable molecular mechanics (MM) force field (see, e.g.,
ref 13) and, therefore, is called the MM fragment. The DFT/
MM hybrid force field approximately accounts for the interac-
tions between the two fragments. It includes the electrostatic
polarizability of the solute through DFT and neglects that of
the solvent beyond the mean field approximation, on which
nonpolarizable MM force fields generally rely (see refs 14-16
for critical discussions). For a phosphate ion in aqueous solution,
Figure 1 (right) shows a part of such a DFT/MM simulation
system and its partitioning into fragments.
2-
Figure 1. Left: Structures of the two phosphate ions HPO4 and
H2PO4- studied in this paper. Right: The ion HPO42- (dark gray, DFT
fragment) in aqueous solution (light gray, MM fragment).
Concerning the computation of vibrational spectra, the high
quality of the particular DFT/MM hybrid approach suggested
in ref 8 and implemented in the parallelized MD program EGO-
MMII17 has been previously demonstrated for DFT/MM models
of the water dimer and for a DFT water molecule dissolved in
MM water,8 for the quinone cofactors (DFT) in a reaction center
(MM) of bacterial photosynthesis,18 and for p-benzoquinone
(DFT) in MM water.11 Note, however, that all these applications
solely cover polar molecules in a similarly polar environment.
In the case of ions, whose electrostatic interactions with a polar
and polarizable environment are much stronger, the performance
of the approach remains yet to be checked.
* Corresponding author: tavan@physik.uni-muenchen.de.
† Ruhr-Universita¨t.
‡ Ludwig-Maximilians-Universita¨t.
10.1021/jp048617g CCC: $27.50 © 2004 American Chemical Society
Published on Web 06/25/2004

IR Spectra of Phosphate Ions in Aqueous Solution
J. Phys. Chem. A, Vol. 108, No. 29, 2004 6187
DFT/MM hybrid MD simulations of solvent-solute systems
enable the computation of vibrational spectra by Fourier
transformation of time-correlation functions (FTTCF).19 How-
ever, FTTCF requires extended MD trajectories and does not
give easy access to the calculation of isotope effects.11 Fur-
thermore, it does not give the normal modes of the solute
molecule, which are required for a detailed analysis of the force
field (see ref 11 for further discussions).
such as ATP or GTP play prominent parts in bioenergetics and
signal transduction and because FTIR techniques, when com-
bined with vibrational analyses, can help in understanding these
processes at a molecular level.21-24
We will complement our DFT/MM hybrid calculations on
the two monophosphates in solution by FTIR measurements
although the IR spectra of these ions have been repeatedly
measured (see ref 25 for an overview). The experimental
reinvestigation of the IR spectra serves to guarantee comparable
thermodynamic conditions (temperature, density, and ion con-
centration) in the two approaches, to verify the assignment of
the IR bands, and to enable quantitative checks with respect to
the accuracy of the computational results.
Deviations between calculated and observed IR spectra may
then be due to insufficiencies of the applied DFT method, to
the inadequacy of the applied harmonic approximation, or to
shortcomings of the MM model of the solvent. To distinguish
these sources of possible deviations, we will first have to get
an estimate on the quality of the harmonic DFT force field
employed by us, which is provided by the program CPMD.4
This program employs plane wave expansions for the construc-
tion of the Kohn-Sham MOs, implements the frozen core
approximation through the use of pseudopotentials, and is
integrated8 into the MD program package EGO-MMII17 to
enable DFT/MM calculations. The DFT quality estimate has to
be done in the conventional vacuum setting by comparisons
with results of standard DFT methods implemented in GAUSS-
IAN98,3 since the accuracy of the corresponding force fields is
well-established7 and since gas-phase spectra of the phosphate
ions, which could provide an experimental reference, are
unknown to us and hard to measure.
To gain access also to the normal modes of a solute molecule
in a DFT/MM setting, we will therefore apply the so-called
instantaneous normal-mode analysis (INMA).11 According to
INMA various statistically independent snapshots are taken from
an MD trajectory of the solute-solvent system, the Hessians
of the solute are calculated within the snapshot ensemble of
solvent cages, which are kept fixed at their 300 K structures,
and the spectral positions, inhomogeneous widths, and intensities
of the vibrational bands are derived from the average positions,
spectral fluctuations, and transition dipole moments of the
calculated vibrational lines. Finally, Gaussian models for the
band shapes are used for the construction of an IR spectrum.
This INMA protocol deviates from the earlier suggestions
of Eichinger9 and Cui and Karplus,10 according to which the
solute-solvent snapshot geometries should be jointly optimized
either by rapid freezing9 or by energy minimization.10 However,
one can expect that the joint optimization of the solute-solvent
structures strongly restricts the calculated ensemble of solvent
cages. This ensemble will then mimic a low temperature and,
therefore, should entail much too narrow bandwidths.
In ref 11 both FTTCF and the above INMA protocol have
been applied to compute the IR spectra of p-benzoquinone in
aqueous solution from a 17.5-ps DFT/MM hybrid MD trajectory.
The spectral positions of the IR bands calculated by the two
approaches turned out to agree quite well, demonstrating that
both approaches are valid. In contrast, the bandwidths and IR
intensities calculated by FTTCF and INMA showed marked
differences, indicating that the 17.5-ps trajectory, although
representing a huge computational effort, has been too short
for a statistically reliable computation of these observables. For
instance in the case of INMA, the requirement of the statistical
independence of the snapshots and the known relaxation time
of water at 300 K have enforced time delays of at least 2.5 ps
between subsequent snapshots, which has limited the number
of snapshots to only seven.20 Furthermore, due to the lack of
experimental data, a detailed evaluation of the merits and
deficiencies of the two methods had to be postponed.
2. Methods
2.1. Choice of DFT Procedures. As mentioned above, there
are no vibrational spectra available for the monophosphates
2-
-
HPO4 and H2PO4 isolated in the gas phase. Therefore, the
evaluation and calibration of the particular DFT method that
we want to apply to the solution spectra of these molecules has
to rely on vacuum calculations of well-established quality.
The very fact that Becke-Perdew DFT calculations26,27
(BP86) can yield unscaled harmonic force fields, which ac-
curately match the experimental, i.e., anharmonic vibrational
spectra, has first been demonstrated for p-benzoquinone5,28 and
for the main group homonuclear diatomics.29 These calculations
have been carried out with extended Gaussian basis sets3
comprising p and d polarization functions. Meanwhile the
excellent agreement of unscaled harmonic frequencies from
BP86 calculations with gas-phase vibrational spectra has been
confirmed for a large number of polyatomic molecules.7 Here,
the computation of anharmonic corrections has demonstrated
that the agreement is caused by a fortunate cancellation of errors.
This cancellation is inherent in BP86 and is lacking in other
DFT methods, like, e.g., in Beckes three-parameter Lee-Yang-
Parr (B3LYP) approach,30,31 which should be a superior method
according to theoretical considerations. Only after the inclusion
of anharmonic corrections do the B3LYP fundamentals agree
better with the observed spectra than the BP86 fundamentals.7
But as compared to the harmonic BP86 frequencies, this
improvement is very small. Therefore, BP86 is the DFT
functional of choice for the computation of vibrational spectra,
whenever one is forced to stick to the harmonic approximation
for reasons of computational feasibility.
Therefore, as to proceed in the establishment of methods, it
is the aim of this contribution to check whether INMA, when
applied to a sufficiently large statistical ensemble of snapshots,
can reliably predict the frequencies, inhomogeneous widths, and
intensities of the IR bands for a solute molecule in a DFT/MM
hybrid setting. To enable a rapid computation of a sufficiently
large snapshot ensemble we will modify the procedure. Instead
of using a costly DFT/MM hybrid MD simulation as in ref 11,
we will apply computationally much cheaper MM simulations
for the generation of statistically independent solvent shells.
Very short hybrid MD trajectories starting at different MM
configurations then yield the ensemble of solvent cages used
for INMA.
The two monophosphates depicted in Figure 1 should
represent sensitive and challenging test cases, because these
highly polarizable ions should exhibit strong solvation effects
in aqueous solution. Furthermore, from the standpoint of
biochemistry, it is an important task to establish a computational
approach to the IR spectra of phosphate ions, because derivatives
Having fixed the DFT functional to BP86, suitable basis sets
remain to be determined. Due to the strong polarizability of

6188 J. Phys. Chem. A, Vol. 108, No. 29, 2004
Kla¨hn et al.
the phosphate group we will use for the all-electron calculations
with Gaussian basis functions a most expensive 6-311++G**
basis set,3 which includes 3d and 1f functions on the heavy
atoms and 2p and 1d functions on the hydrogens, and adds
diffuse functions on both. The corresponding results will serve
as a reference for the evaluation of the quality of the basis sets
and of the pseudopotentials used by us in the plane-wave DFT
program CPMD.4
2.3. MD Simulations. To determine the solvation effects on
the structures and charge distributions of the phosphate ions
and to obtain solution-adapted partial charges for their MM
models, the freshly generated ionic solutions have been equili-
brated first for 100 ps in the MM setting (using the vacuum
estimates for the phosphate partial charges). Here, as in all other
MM simulations, a time step of 1 fs and an NVT ensemble
were used; the temperature was controlled through coupling to
an external heat bath36 (coupling time constant τ ) 0.1 ps, target
temperature T0 ) 300 K). In contrast, in the (HGH/BP)/MM
simulations the time step was reduced to 0.25 fs. In both settings
the long-range Coulomb interactions were treated by the multiple
moving-boundary reaction field approach of Mathias et al.37 with
a dielectric boundary toward a solvent continuum 23 Å distant
from each molecule. Within that boundary the electrostatics was
treated explicitly with use of the structure-adapted multipole
method.38,39 Here, the Coulomb sum was employed for the
atomic partial charges up to a distance of 10 Å. Beyond this
distance the electrostatics was computed by multipole expan-
sions of the partial charge distributions within the water
molecules and ions. The field, including the reaction field
correction, generated by the environment of a respective DFT
fragment was imported into the DFT Hamiltonian. For further
details on the DFT/MM interface we refer to ref 8.
INMA requires the consideration of many solvent cages,
within which the energies of the phosphates have to be
minimized and their Hessians have to be calculated. To keep
these computations manageable, the size of the plane wave basis
set should be at the lower limit compatible with the necessary
accuracy. In this respect a series of test calculations has
demonstrated (data not shown) that in combination with the
BP86 functional the pseudopotentials suggested by Hartwigsen
et al.32 yield a satisfactory accuracy at sufficiently small basis
sets. The corresponding plane-wave cutoff and the extension
of the box, in which the charge density is calculated on a grid,
have been determined by gradual decreases of the cutoff and
box sizes until the energies of the phosphates started to sizably
change. The result was a cutoff at 80 Ry corresponding to 22 716
plane waves and a minimal distance between any DFT atom
and the surface of the box at 2.7 Å. We denote this DFT method
by HGH/BP and the corresponding DFT/MM hybrid approach
by (HGH/BP)/MM.
After the initial MM equilibrations the simulations were
continued for another 100 ps serving to generate different initial
conditions for a series of 0.1-ps (HGH/BP)/MM simulations.
In each of these simulations one of the phosphate ions was
selected as the DFT fragment. The calculated ESP charges were
averaged over all hybrid trajectories and over equivalent atoms,
i.e., the atoms of the different types P, O, OH, and H, respec-
tively, occurring in the phosphates (cf. Appendix). The thus
obtained values have been taken as refined estimates for the
MM partial charges of the phosphate ions in solution. Table 1
in Section 3 lists the values of the corresponding solvent-adapted
partial charges. The hybrid simulations were also used for
statistics over other structural parameters of the phosphates, for
instance for the computation of the average bond lengths, which
will be mentioned in Section 3.
2.2. Simulation System and MM Force Field. For the
computation of the IR solution spectra of the monophosphates
two periodic boxes (46.6 Å3) were filled with 3375 water
2-
-
molecules and eight HPO4 or H2PO4 ions, respectively.
Appropriate numbers of Na+ ions were added to balance the
negative charges. These compositions represent concentrations
of 130 mM, which are close to the phosphate concentrations
used in the FTIR experiments. At 300 K the densities of the
solutions correspond to the experimental density of water at
this temperature.
The water molecules have been treated by MM and described
by a variant of the simple TIP3P (three-point transferable
intermolecular potential) force field, in which, in contrast to
the original TIP3P model,33 the water molecules are completely
flexible. The required force constants, partial charges, and
Lennard-Jones parameters have been taken from the CHARMM19
force field.34,35
By using the MM force field with the solvent-adapted partial
charges for the phosphates the simulation systems were equili-
brated at 300 K for another 100 ps. Subsequently, MD
simulations with durations of 300 ps were carried out to produce
the desired ensembles of independent solvent cages. Starting
for each ionic species from three snapshots (at 100, 200, and
300 ps) of the corresponding simulations system, short 75-fs
(HGH/BP)/MM simulations were executed for each of the eight
solvated phosphates yielding a total of 24 different solvent
configurations per species surrounding the respective DFT
fragments.
To ensure a close matching between the descriptions of the
phosphates in the envisaged MM and (HGH/BP)/MM simula-
tions, respectively, a CHARMM22-type13 MM force field has
been constructed for these ions as follows. First, vacuum HGH/
BP calculations were executed for the two ions, through which
equilibrium geometries and Hessian matrices were obtained. The
parameters of the harmonic contributions to the force field (cf.
Appendix) were derived from the equilibrium structures and
through transformations of the Hessians into a basis of the
relevant internal degrees of freedom. The parameters of the
torsional potentials were derived from fits to the corresponding
HGH/BP potential curves. Further details and the results are
given in the Appendix, which also addresses the choice of
Lennard-Jones parameters required both in the MM and in the
(HGH/BP)/MM settings. The vacuum HGH/BP calculations
yielded a set of ESP charges (see, e.g., ref 8 for details), which
were used as initial guesses for the MM partial charges in setting
up the simulation system. Improved estimates were then
obtained by MD simulations based on the (HGH/BP)/MM
Hamiltonian as described below.
2.4. INMA Calculations of IR Spectra. For each of the two
ionic species and each of the 24 corresponding snapshots, the
solvent cages were now kept fixed and the (HGH/BP)/MM
Hamiltonian was minimized by gradient descent with respect
to the coordinates of the DFT fragment until a local minimum
was reached as detected by energy gradients smaller than 8.4
× 10-5 au. Next, the Hessian and the dipole gradients of the
phosphates were calculated in the (HGH/BP)/MM setting. Here,
analytical first derivatives and numerical second derivatives were
used. These data served as inputs for normal-mode analyses
executed in the usual way yielding normal modes, corresponding
frequencies, and IR intensities for each snapshot. The normal
modes were classified by visual inspection, and, separately in

IR Spectra of Phosphate Ions in Aqueous Solution
J. Phys. Chem. A, Vol. 108, No. 29, 2004 6189
vibrational modes involving motions of the heavy atoms have
frequencies below this limit.25
All DFT methods agree that the frequencies of eight of the
twelve normal modes of HPO42- should be located in this range.
These modes cover four different bond stretches ν between the
heavy atoms P and O and four deformations δ of the angles
between these atoms. For a more detailed characterization of
2-
the normal modes two properties of HPO4 can be used. (i)
There are two different oxygen bonding patterns, i.e., three
oxygens are exclusively bound to P and one oxygen is
additionally bound to the hydrogen (cf. Figure 1, top left).
Therefore, the oxygens can be distinguished by introducing the
atom types O and OH, respectively. (ii) Due to the presence of
three equivalent oxygens, the normal modes will comprise in-
phase and out-of-phase combinations of the internal coordinates
involving these oxygens. Therefore, they can be distinguished
by the labels s and a, respectively. In Figure 2 the thus defined
nomenclature has been used to classify the energy levels in the
depicted term scheme. For instance, νs(PO) characterizes the
normal mode, in which all three P-O bonds are simultaneously
elongated or compressed.
Figure 2. Frequencies of normal modes calculated by various
2-
DFT methods for HPO4 isolated in a vacuum. Results of the HGH/
BP approach provided by CPMD4 are compared with results obtained
from the BP86, BLYP, and B3LYP functionals implemented in
GAUSSIAN983 (see Section 2 for computational details and the text
for discussion).
each class, second order statistics was applied to the frequencies
and intensities. With use of the standard deviations and average
frequencies as parameters of normal distributions, the inhomo-
geneously broadened shapes and spectral positions of the various
IR bands were obtained. After weighting these normalized
Gaussians with the respective average intensities, they were
added up to give the desired INMA IR spectra that can be
directly compared with the observed IR spectra after proper
adjustment of the intensity scale.
Figure 2 shows that HGH/BP predicts frequencies that are
upshifted by 23 cm-1 in a close to uniform fashion as compared
to the BP86 reference. Therefore, the respective spectral
orderings are identical. Identical orderings are also predicted
by BLYP and B3LYP. Compared to BP86, these spectra show
either an average red-shift by 11 cm-1 (BLYP) or an average
blue-shift by 39 cm-1 (B3LYP). Thus, the harmonic HGH/BP
force field is a little stiffer than that of the reference method
BP86 (but softer than that of B3LYP) demonstrating that the
use of pseudopotentials and of different basis sets is of minor
importance. The consideration of our corresponding results for
2.5. Experimental IR Spectra. The IR spectra of solvated
K2HPO4 and KH2PO4 and of deuterated and 18O-labeled
phosphates have been measured at 300 K via FTIR. The FTIR
measurements have been performed on a Bruker IFS 66/v
spectrometer, using a DTGS detector and ZnSe windows.
Spectra have been recorded with a spectral resolution of 4 cm-1
in the frequency range between 480 and 1300 cm-1 and have
been apodized with the Blackman-Harris three-term function.40
A zero-filling factor of 2 has been used. All samples consisted
of ∼1 µm thick films of 100 mM solutions. Absorbance spectra
have been calculated by means of a reference spectrum of pure
H2O (D2O). K2HPO4 (99%) and KH2PO4 (99%) have been
obtained from Merck and Baker, respectively, and have been
measured without further purification. K2DPO4 and KD2PO4
have been prepared by drying the solution three times and
rehydrating it with D2O (99.9%, Deutero). 18O-labeled phos-
phates have been synthesized by hydrolysis of PCl5 (98%, Fluka)
with H218O (99.1%, Isotec). Excess H218O has been removed
by drying and rehydrating with H2O. The pH value has been
-
the other phosphate ion H2PO4 provides no further insights
and, therefore, is omitted. Note, in addition, that refs 18 and 11
arrived at similar conclusions for a different class of molecules,
i.e., the quinones, concerning the comparable qualities of
Becke-Perdew DFT descriptions using GAUSSIAN98 and
CPMD, respectively.
In view of the known quality of the harmonic BP86 force
field in describing gas-phase spectra,7 we thus conclude that
the predictions of the closely related HGH/BP method also can
be considered as high-quality estimates for the unknown gas-
phase spectra of phosphate ions. Therefore, the effects of
solvation on the structure and vibrational spectra of phosphate
ions can be judged by comparison of HGH/BP vacuum results
with experimental solution data.
3.2. Solvation Effects on the Phosphate Structures. In our
MM simulations of the ionic solutions, the ions stayed away
from each other as to maximize their hydration shells. The
average distance of a sodium ion to the next phosphate was
about 8 Å and never got smaller than 4 Å. Due to the strong
dielectric shielding of the Na+ charges by the surrounding water,
the solvation effects, which are predicted by our (HGH/BP)/
MM simulations for the phosphate ions, should be mainly due
to their hydration shells, whereas the influence of the sodium
ions should be small.
2-
-
calibrated by HCl and NaOH. The HPO4 and H2PO4 ions
have been measured at pH 9.2 and 5.2. At pH 7.2 the two ions
are found at a 1:1 ratio.
3. Results and Discussion
3.1. Checks of the HGH/BP Approach. We start with the
required estimate on the quality of the DFT approach used by
us, that is of HGH/BP. Figure 2 compares the gas-phase
vibrational spectrum predicted for the doubly charged phosphate
We first turn to the hydration-induced changes of the charge
distributions within the phosphate ions. As observables we use
the electrostatic potential (ESP) partial charges derived from
the HGH/BP vacuum calculations and from the (HGH/BP)/MM
solvent simulations described in Section 2. These charges are
optimized so as to match the electrostatic potential generated
by a molecule at its van der Waals surface. Table 1 lists the
resulting partial charges for the two phosphate ions under
2-
anion HPO4 by the HGH/BP approach with results of DFT
calculations obtained with the GAUSSIAN98 standard func-
tionals BP86, BLYP, and B3LYP, respectively. For this
comparison, as for all other comparisons below, we have chosen
the spectral range above 400 cm-1, because the transmission
properties of the ZnSe windows set a corresponding lower limit
to our FTIR spectra, and below 1300 cm-1, because all

6190 J. Phys. Chem. A, Vol. 108, No. 29, 2004
Kla¨hn et al.
TABLE 1: ESP Partial Charges q Given in Units of the
2-
-
Elementary Charge e for HPO4 and H2PO4 in a Vacuum
and in Solution, Respectively, As Derived from HGH/BP and
(HGH/BP)/MM Calculationsa
2-
-
HPO4
vacuum
H2PO4
vacuum
atom type
solution
solution
P
O
OH
H
+1.71
-1.05
-0.87
+0.31
+1.43
-1.02
-0.76
+0.39
+1.26
-0.79
-0.68
+0.34
+1.20
-0.86
-0.64
+0.40
a The solution values represent averages over a series of short (HGH/
BP)/MM trajectories (see Section 2 for details and the Appendix for
the definition of the atom types).
consideration. The partial charges in solution represent averages
over the ensembles of 24 short (HGH/BP)/MM simulations
described in Section 2. The fluctuations of the partial charges
are very small as documented by the average standard deviations
of 0.03 e. For both ions an analysis of the data reveals that the
hydration shifts about 30% of the negative charge, which in
vacuo is located at the OH-H groups, toward the PO3 and PO2
moieties of HPO42- and H2PO4-, respectively. These substantial
charge shifts are convincing witnesses for the strong polariz-
abilities of the phosphate ions.
The shifts of charge density are accompanied by hydration-
induced decreases of the P-OH bond lengths (from 1.78 Å to
an average value of 1.67 Å for HPO42- and from 1.67 Å to an
average value of 1.63 Å for H2PO4-). According to the usual
arguments, the shortening of a bond length points toward an
increase of the corresponding bond stretching force constant.
Therefore, one expects strong blue-shifts for the frequencies of
all normal modes, to which the P-OH bond stretch gives an
important contribution. For the case of the HPO42- ion, whose
vacuum spectrum has been discussed above in connection with
Figure 2, we therefore expect that the modes ν(POH) and
δ(POHH) will be strongly blue-shifted by hydration.
2-
Figure 3. Comparison of the IR spectra obtained for HPO4 (a, b)
-
and H2PO4 (c, d) by FTIR (a, c) and INMA (b, d), respectively. For
computational and experimental details see Section 2 and for a
discussion see the text.
remaining candidate. Note that the latter suggestion agrees with
an earlier assignment.25
2-
The above assignment of the FTIR spectrum of HPO4 in
aqueous solution has been based on the intensities, the band
shapes, and the spectral ordering of the sharp bands predicted
by the INMA calculations. Somewhat surprisingly it implies
that the (HGH/BP)/MM approach has severe deficiencies
concerning the description of band frequencies, because the
spectral positions of several INMA bands in Figure 3b sizably
deviate from those of their FTIR counterparts in Figure 3a.
Therefore we should check as to whether similar deficiencies
3.3. Assignment of FTIR Bands by INMA Results. Figure
3 compares the IR spectra of HPO42- (a, b) and H2PO4- (c, d)
as measured by FTIR (a, c) with the results of the INMA (b, d)
calculations.
The FTIR spectrum of HPO42- depicted in Figure 3a exhibits
four rather sharp peaks a-d at 1080, 990, 855, and 541 cm-1
,
respectively, as well as a very weak and broad shoulder s. Baril
et al.25 have attributed this shoulder, which was also visible in
their IR spectra, to a band at 1197 cm-1. Likewise the
corresponding INMA spectrum shown in Figure 3b predicts four
sharp bands, one of them being generated by two out-of-phase
P-O stretching modes νa(PO), and a very broad band, which
is due to the P-OH-H bending mode δ(POHH).
If one takes the calculated shapes and relative intensities of
the INMA bands as a guideline and if one assumes that our
(HGH/BP)/MM approach correctly reproduces the spectral
ordering of the sharp peaks, then the visual comparison of the
FTIR and INMA spectra in parts a and b of Figure 3
immediately suggests an assignment of the FTIR peaks. In
particular, the most intense and high-frequency peak a has to
be assigned to a superposition of the two intense INMA bands
νa(PO), the less intense peak b to the equally less intense INMA
band νs(PO) associated with the in-phase P-O stretch, the
broader peak c to the broader INMA band ν(POH), and peak d
to the three neighboring INMA bands δ(OPO). Thus, in the
FTIR spectrum of Figure 3a only the weak shoulder s at 1197
cm-1 remains to be assigned and, here, the very broad INMA
band δ(POHH) of the P-OH-H bending mode is the only
appear also in the related case of the singly charged anion
-
H2PO4
.
The FTIR spectrum of H2PO4- depicted in Figure 3c exhibits
five sharp bands a-e at 1156, 1077, 944, 879, and 521 cm-1
,
respectively, and a broad shoulder s at 1213 cm-1. Also the
corresponding INMA spectrum in Figure 3d shows one broad
and five sharp bands. Applying once again the criteria explained
above, the FTIR peak a is assigned to νa(PO), peak b to νs-
(PO), peaks c and d to the two P-OH stretches νa(POH) and
νs(POH), respectively, and peak e to the assembly of O-P-O
bending modes surrounding the sharp δ(OPO) band predicted
by INMA at 483 cm-1. In agreement with an earlier sugges-
tion,25 finally the broad FTIR shoulder at 1213 cm-1 has to be
assigned to the two P-OH-H bending modes generating the
broad INMA bands δs(POHH) and δa(POHH).
Having tentatively assigned the IR bands also for H2PO4-, a
comparison of the spectral positions of corresponding bands in
Figure 3c,d now strengthens the suspicion that the (HGH/BP)/
MM approach cannot predict the vibrational frequencies of small
ions in aqueous solution at the same level of accuracy at which
pure DFT treatments such as HGH/BP usually describe the gas-

IR Spectra of Phosphate Ions in Aqueous Solution
J. Phys. Chem. A, Vol. 108, No. 29, 2004 6191
phase spectra of a wide class of molecules.7 Also for H2PO4
-
some of the FTIR frequencies strongly deviate from the INMA
results.
Because this conclusion depends on the correctness of the
assignments, we will now shortly discuss isotope effects for
clarification. To save space we will not show the corresponding
FTIR and INMA spectra (see, however, the Supporting Infor-
mation for some of these spectra). Instead we will merely infer
the key results.
3.4. Corroboration of Assignment by Isotope Labeling.
The most intense FTIR band a at 1080 cm-1 of HPO4 (cf.
2-
Figure 3a) splits in D2O into two bands at 1097 and 1079 cm-1
,
respectively. This observation proves two aspects of our
assignment. The splitting demonstrates that the HPO42- band a
is actually due to two nearly degenerate P-O stretching modes
νa(PO) described by INMA. The 17-cm-1 blue-shift induced
by deuteration at one of the νa(PO) modes shows that another
mode, which in H2O is coupled to this νa(PO) mode and is
located at a higher frequency, must have been strongly red-
shifted and correspondingly decoupled upon deuteration. Ac-
cording to our assignment this higher frequency mode should
be the δ(POHH) vibration, which we have associated to the very
weak shoulder at 1197 cm-1
.
2-
The INMA results on DPO4 corroborate this analysis. As
main effects of deuteration they predict a large 269-cm-1 red-
shift for δ(POHH), which decouples this local mode from one
of the nearby νa(PO) modes, and a sizable 18-cm-1 red-shift
for this νa(PO) mode, proving the decoupling. In this context it
is irrelevant that INMA assigns the wrong sign to the decoupling
shift of νa(PO). Because INMA has erroneously located the
δ(POHH) frequency below that of νa(PO), the wrong sign of
the calculated decoupling shift of νa(PO) is a necessary
2-
Figure 4. Hydration effects on the vibrational spectra (a) of HPO4
and (b) of H2PO4-. The gas-phase frequencies calculated by HGH/BP
are compared with the average INMA frequencies obtained by (HGH/
BP)/MM for the solution case and with the band positions measured
by FTIR. For a discussion see the text.
of band d is nicely reproduced by INMA predicting a 22-cm-1
red-shift of the associated ν(POH) band.
2-
consequence. The FTIR spectrum of DPO4 furthermore
reveals that the peak c of the HPO42- spectrum in Figure 3a is
red-shifted by 20 cm-1. INMA predicts for the corresponding
Apart from corroborating for both ions the assignments of
the sharp FTIR peaks in the 800-1300 cm-1 range to P-O
and P-OH stretching vibrations, the substitution of the oxygens
by 18O isotopes does not add any further insights. On the average
these peaks experience 33-cm-1 red-shifts. In nice agreement
INMA predicts an average red-shift of 35 cm-1 for the
corresponding modes.
3.5. Hydration Effects on Vibrational Frequencies. Having
fixed the assignment of all IR bands in the 400-1300 cm-1
spectral range we are now prepared to answer the questions (i)
as to how the spectral locations of the various vibrational modes
are modified by hydration and (ii) to what extent these effects
are covered by our (HGH/BP)/MM approach.The answers are
provided by the two-term schemes in Figure 4, which compare
for the two phosphate ions the gas-phase spectra estimated by
vacuum HGH/BP calculations (first columns) with the spectral
locations of the IR bands observed in solution (last columns)
and with their INMA descriptions (central columns). Because
of the established quality of our DFT method (cf. Section 3.1),
the HGH/BP estimates of the unknown gas-phase spectra should
be accurate within a few wavenumbers. Therefore, the com-
parisons of these estimates with the FTIR solution spectra can
provide insights into the effects of hydration.
ν(POH) mode a similar deuteration shift of 14 cm-1
.
2-
For HPO4 the decoupling argument is our key indication
on the spectral location of δ(POHH) at 1197 cm-1. As is apparent
from Figure 3a, the weak intensity in the 1200-cm-1 spectral
range prevents the clear-cut and direct identification of this
δ(POHH) band, which according to INMA should be very broad.
Due to baseline problems in the D2O spectrum, we were also
unable to detect the expected disappearance of the weak shoulder
from the 1197-cm-1 spectral location.
For H2PO4-, however, the experimental evidence for the
location of the corresponding δ(POHH) bands at the blue edge
of the spectral range is direct. Here, the broad shoulder at 1213
cm-1 is more intense (cf. Figure 3c) and this shoulder clearly
2-
disappears upon deuteration. Quite like in the case of HPO4
,
the FTIR spectrum of D2PO4- also reveals a deuteration-induced
decoupling of the νa(PO) mode from the higher frequency
δ(POHH) modes as witnessed by a 21-cm-1 blue-shift of the
νa(PO) band a in D2O. Because of the additional direct
spectroscopic evidence on the high-frequency location of the
δ(POHH) modes, the decoupling argument can thus be consid-
2-
ered as hard evidence also for HPO4 and finally proves the
assignment of the very weak shoulder in Figure 3a to δ(POHH).
Already a first glance at the two-term schemes in Figure 4
reveals that hydration causes very large shifts of the vibrational
frequencies by at least 25 cm-1 up to as much as 288 cm-1
(compare the first and last columns). With the sole exception
of νa(PO) in H2PO4-, all frequencies are shifted to the blue.
The hydration effects are particularly pronounced for all modes,
which involve motions of the protonated oxygen atoms OH.
Also for H2PO4- the decoupling is described by INMA. The
deuteration shifts the frequency of the νa(PO) mode by 20 cm-1
to the red, because the δ(POHH) modes are shifted by about
300 cm-1 to much lower frequencies. Again the sign of the 20
cm-1 decoupling shift is wrong, because the δ(POHH) frequen-
cies are underestimated. Finally, the observed 17-cm-1 red-shift

6192 J. Phys. Chem. A, Vol. 108, No. 29, 2004
Kla¨hn et al.
Here, the P-OH stretches experience blue-shifts on the average
3.6. Inhomogeneous Bandwidths and Intensities. In Section
3.3 we have assigned the FTIR bands to normal modes by
visually comparing the INMA band shapes and intensities with
corresponding features of the FTIR spectra. Apart from the
assumption that INMA predicts a correct spectral ordering for
the sharp bands, we have largely ignored the detailed spectral
positions of the bands in this process. Its apparent success
demonstrates that INMA performs quite well on the shapes and
relative intensities of the IR bands.
Unfortunately the insufficient quality of our FTIR spectra
prevents for many bands a quantitative comparison of their
widths and intensities with the INMA descriptions. Solely for
the three high-frequency FTIR bands s, a, and b of H2PO4- (cf.
Figure 3c), a reasonable Gaussian fit was possible. This fit
yielded bandwidths of 114, 46, and 26 cm-1, respectively. The
INMA values are 99, 33, and 24 cm-1, respectively, for the
corresponding bands associated with the two δ(POHH), the
νa(PO) and νs(PO) modes. Thus, on the average, INMA under-
estimates the FTIR bandwidths only by 17% and we may
conclude that the INMA procedure of sampling different
hydration cages can reasonably account for the fluctuations of
the electric field, by which the strong electrostatic dipoles of
the thermally moving water molecules polarize the phosphate
ions.
by 220 cm-1 and the P-OH-H bending vibrations by 165 cm-1
.
Recall that our analysis of the hydration effects on the phosphate
structures (cf. Section 3.2) had led us to expect the former blue-
shifts. In line with these arguments, the unusual magnitude of
these shifts is a signature of the large polarizability of the
phosphate ions and of the strong polarity of the aqueous solvent.
As a result, the vibrational spectra of the phosphate ions
represent a particularly hard testing case for any method, which
is supposed to describe the vibrational spectra of molecules in
polar solvents.
To apply this test to our (HGH/BP)/MM hybrid method, we
compare the central columns of the term schemes in Figure 4
with the gas-phase and FTIR solution spectra in the neighboring
columns. Apart from two apparent exceptions to be discussed
further below, the inclusion of the MM solvent into the DFT
description shifts all vibrational levels from their gas-phase
positions toward the FTIR solution values. However, as seen
most clearly at the blue-shift of the P-OH stretching mode of
HPO42-, this approach toward the FTIR spectrum covers only
about half the distance: The (HGH/BP)/MM treatment shifts
the frequency of ν(POH) only by 161 cm-1 to the blue of the
HGH/BP vacuum result instead of the 288 cm-1 required by
FTIR.
As quantitatively determined here and as is apparent from
the INMA results in Figure 3b,d, the δ(POHH) bands are much
broader than the other bands. This large broadening is readily
understood by considering the fact that the P-OH-H bending
modes involve the motion of a single proton, whereas all other
modes exclusively comprise heavy nuclei. Therefore, the charge
density localized by the weak Coulomb potential of the proton
and associated with the OH-H bond is much more strongly
perturbed by the fluctuations of the external electric field than
the charge densities of the P-O and P-OH bonds. In Section
3.2, we have already noted another signature of this effect, that
is the hydration-induced charge transfer from the peripheral
OH-H group toward the core of the molecule.
The only exceptions from this general trend apparent in the
term schemes are the large blue-shifts of the two νa(PO) modes
of HPO42-, which exceed the required 25-cm-1 shifts on the
average by about 40 cm-1. But these exceptions are easily
understood in terms of the coupling between the corresponding
local νa(PO) modes and the local δ(POHH) mode analyzed in
Section 3.4. The (HGH/BP)/MM modeling of the hydration
apparently does not suffice to shift the intrinsic frequency of
the local δ(POHH) mode above those of the local νa(PO) modes
as required by FTIR. Therefore, in INMA the coupling between
these local modes induces a spectral splitting, which shifts the
νa(PO) frequencies too far to the blue and diminishes the blue-
shift of δ(POHH) by the same amount.
Note here that the absolute value predicted by (HGH/BP)/
MM for this spectral splitting is about as large (≈100 cm) as
the one found by FTIR. When taken together with the fact that
the deuteration effects have revealed similarly sized couplings
in both cases, the comparable magnitudes and the opposite signs
of the splittings show that the intrinsic frequency of the local
δ(POHH) mode must be located experimentally just above and,
according to the (HGH/BP)/MM description, just below those
of the two local νa(PO) modes. Similar coupling effects between
the local νa(PO) mode and the two local δ(POHH) modes can
be deduced from Figure 4b also for H2PO4-. In this case,
however, they do not disturb the homogeneity by which the
vibrational levels approach the FTIR spectra upon inclusion of
an MM solvent model into the DFT description of the phosphate.
Concerning the accuracy of the calculated relative intensities
the above fit renders only two numbers and, therefore, one
cannot derive any further conclusions. In this respect we have
to be content with the positive judgment on the corresponding
performance of INMA stated above.
4. Discussion and Outlook
Apart from small underestimates, our computationally more
efficient INMA protocol turned out to reproduce the widths of
the phosphate IR bands quite well. This result proves that the
sampled ensemble of hydration shells appropriately covers the
fluctuations of the electric field, by which the surrounding water
dipoles polarize the phosphate ions. Also the IR intensities
calculated by INMA seem to be reasonable. In line with the
arguments in ref 11 this finding demonstrates that the neglect
of the kinetic couplings between the solute and the solvent,
which is a characteristic of the INMA approach, hardly affects
the intensities. However, the hydration effects on the spectral
locations of the phosphate IR bands predicted by our DFT/MM
Hamiltonian turned out to be underestimated by about a factor
of 2, and this result requires an explanation.
In Section 3.1 we have demonstrated that the DFT part of
the hybrid Hamiltonian, that is the HGH/BP method provided
by CPMD, is as accurate as standard high-quality DFT
methods: Despite the applied harmonic approximation the
calculated IR spectra should therefore deviate from correspond-
ing observations only by a few wavenumbers per band.
In summary, the (HG/BP)/MM treatment of the two phos-
phate ions in aqueous solution can qualitatively account for the
effects of hydration on the vibrational spectra, but it underes-
timates these effects by about a factor of 2. The hydration effects
are particularly large for all modes involving the motions of a
protonated oxygen atom, that is for the P-OH stretching and
P-OH-H bending vibrations. Therefore, the deviations of the
(HG/BP)/MM frequencies from the FTIR observations are
particularly large for these modes. In the case of the P-OH-H
bending vibrations, hydration shifts the frequencies according
to FTIR to the blue edge of the investigated spectral range. Here,
the (HG/BP)/MM underestimate of hydration effects even causes
a wrong spectral ordering.

IR Spectra of Phosphate Ions in Aqueous Solution
J. Phys. Chem. A, Vol. 108, No. 29, 2004 6193
TABLE 2: Parameters of the Bonded Force Field 1 for the
Phosphate Ionsa
predicted for these molecules by the applied DFT functional.
Such a parametrization will then enable a smooth transition from
an initial MMpol to a subsequent DFT description of the water
molecules in the first hydration shells of the DFT phosphate
ions. Finally, a comparison of the phosphate vibrational spectra
resulting from these two sets of additional calculations will
provide estimates on the relative sizes, by which the two effects
neglected by us (polarizability vs charge transfer) can correct
the descriptions. But whatever these future studies may show,
our results have clearly demonstrated that the quantitatively
correct DFT/MM description of the phosphate vibrational levels
represents a hard and sensitive test for the quality of the chosen
computational method and, here, particularly for the quality and
transferability of polarizable MM force fields.
ij
fij
r0,ij, Å
O-P
OH-P
OH-H
395
98
519
1.61
1.87
0.98
ijk
fijk
æ
0,ijk, deg
O-P-O
O-P-OH
P-OH-H
73
73
45
116
105
105
ijkl
kijkl
0.6
δ
ijkl, deg
180
O-P-OH-H
a The force constants fij of the bonds are given in kcal/(mol‚Å2),
those of the angles (fijk) in kcal/(mol‚deg2), and the torsional barriers
kijkl in kcal/mol. The phosphate atom types P, O, OH, and H are
explained in the text.
Acknowledgment. G.M. and M.N. acknowledge support
by the Volkswagen-Stiftung (project I/73 224) and P.T. a grant
by the Deutsche Forschungsgemeinschaft (SFB533/C3). The
Bochum group is grateful for support by the Deutsche Fors-
chungsgemeinschaft and to B. Portmann, C. Wegener, and F.
Wennmohs for setting up the Myrinet linked PC cluster used
for the calculations. The authors thank Matthias Schmitz for
his careful proofreading of the manuscript.
Correspondingly, the detected underestimate of the hydration
effects must be due to insufficiencies of the MM part and, here,
the most likely cause is the neglect of the electronic polariz-
ability in the MM force field. The following arguments support
this explanation.
Nonpolarizable MM force fields for water, like the one by
which we try to cover the effects of hydration in our hybrid
calculations, have been designed as mean field models for water
molecules in a polar environment, that is in pure water (cf. the
recent review in ref 14). In particular, the dipole moments of
the MM water models are fixed at a value which exceeds that
of an isolated water molecule by about 40%. As demonstrated
by DFT and DFT/MM simulations,8,41 in liquid water the dipole
moments of the water molecules increase by about this
percentage due to their polarization by surrounding water
dipoles. As a result, the transferability of these MM water
models into a different environment, like into the hydration shell
of a small ion, is not guaranteed.14
A small ion generates a much stronger electric field than a
water molecule. Therefore, the induced dipole moments of the
water molecules in its hydration shell will be larger than those
of water molecules, which are in the environment of pure liquid
water. In particular, the resulting total dipole moments of the
water molecules hydrating the ion will be larger than the fixed
values provided by the MM force field. For a small ion, these
arguments lead to the general expectation that the hydration
effects, which are calculated from such a force field in a DFT/
MM setting, will be sizably underestimated. Our results for the
phosphate ions, in particular the nearly homogeneous but
incomplete approach of the DFT/MM vibrational levels toward
the FTIR solution values, support this expectation. Concurrently
the above arguments also explain why this DFT/MM setting
performs well8,11 for polar molecules solvated in water. Here,
the effects of polarizability resemble more closely those in pure
water and the nonpolarizable MM force fields are more
appropriate.
Supporting Information Available: FTIR spectra of H2PO4-
and of the isotopomers D2PO4 and H2P18O4 (Figure S1) as
well as a band fit to the FTIR spectrum of H2PO4- (Figure S2).
This material is available free of charge via the Internet at http://
pubs.acs.org.
-
-
Appendix: MM Phosphate Force Field
During the MM simulations the so-called “bonded” interac-
tions within the phosphate ions were described by a model
potential
EMM,bond
)
f (r - r )2 + f (æ - æ )2 +
∑
∑
ijk
0
0
ij
(k/2)(1 + cos(3θ + δ)) (1)
∑
ijkl
which is based on the CHARMM22 concepts.13 The indices
run over all atoms, which are covalently connected through at
most three chemical bonds. The variables r, æ, and θ character-
ize the associated internal coordinates (for simplicity the required
indices have been omitted). The first term assigns a harmonic
bond stretching potential to all covalently connected atom pairs.
The harmonic potentials of the second term serve to keep the
bond angles close to their equilibrium values, and the third term
parametrizes the torsional potentials. Identical values are chosen
for the force constants f, the torsional energy barriers k, the
equilibrium bond lengths r0, and angles æ0, as well as for the
phase shifts δ, if the atoms involved have identical atom types.
In CHARMM22 atom types are assigned to all atoms of a
molecule according to their respective covalent bonding patterns.
In the phosphate ions, four atom types are required for the
classification of the atoms. These are the phosphorus P, the
hydroxyl oxygens OH, the remaining oxygens O, and the
hydroxyl hydrogens H. Table 2 collects the parameters of the
phosphate model potential 1, which have been determined by
vacuum HGH/BP calculations as sketched in Section 2.
However, as correctly argued by both referees, also a charge
transfer between the phosphate ions and the hydrating water
molecules, which is unaccounted for by our simple DFT/MM
approach confining the DFT description to the phosphates, could
contribute to the noted underestimate of the hydration effects.
For a check of these hypotheses one should recalculate the
vibrational spectra of the phosphate ions using a polarizable
MM force field for the water and, in addition, more extended
DFT fragments covering the first hydration shells of the
phosphates. Here, the chosen MMpol force field for the water
molecules should closely model the properties, which are
The Lennard-Jones parameters, which cover the van der
Waals interactions in the simulation system both in the pure
MM and in the hybrid DFT/MM settings, have been taken from
the CHARMM22 file parallh22x.nuc by using the CHARMM22
atom types P for P, ON4 for OH, ON3 for O, and HN4 for H.

6194 J. Phys. Chem. A, Vol. 108, No. 29, 2004
Kla¨hn et al.
(14) Guillot, B. J. Mol. Liquids 2002, 101, 219-260.
Solution-adapted partial charges qi calculated by (HGH/BP)/
MM simulations (cf. Section 2) and required for pure MM
simulations, are collected in Table 1.
(15) Ponder, J. W.; Case, D. A. AdV. Protein Chem. 2003, 66, 27-85.
(16) Tavan, P.; Carstens, H.; Mathias, G. Molecular dynamics simula-
tions of proteins and peptides: Problems, achievements, and perspectives.
In Handbook of Protein Folding; Buchner, J., Kiefhaber, T., Eds.; Wiley-
VCH: Weinheim, Germany, 2004, in press.
The quality of the thus defined MM force field has been
checked by the computation of the vacuum vibrational spectra.
In the case of HPO42-, for instance, we have obtained MM
frequencies for the modes shown in Figure 2, which exhibit
rms deviations of 55 cm-1 from the HGH/BP results and of 44
cm-1 from BP86, respectively, while preserving the energetic
ordering of the modes. Since vibrational spectra are very
sensitive to the details of a force field, this quite good agreement
demonstrates that the MM parameters in Table 2 represent a
reasonable choice.
(17) Mathias, G.; Eichinger, M.; Carstens, H.; Stork, M.; Weiss, A.;
Heller, H.; Grubmu¨ller, H.; Egwolf, B.; Niedermeier, C.; Tavan, P. EGO-
MMII users guide; Lehrstuhl fu¨r BioMolekulare Optik, Ludwig Maximilian
Universita¨t Mu¨nchen, Oettingenstrasse 67, D-80538 Mu¨nchen, unpublished.
(18) Nonella, M.; Mathias, G.; Eichinger, M.; Tavan, P. J. Phys. Chem.
B 2003, 107, 316-322.
(19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids;
Oxford University Press: Oxford, UK, 1996.
(20) Computational limitations have enforced a similarly bad statistic
(10 snapshots) in the INMA calculations of Cui and Karplus.10
(21) Gerwert, K. Molecular Reaction Mechanism of Proteins Monitored
by Time-resolved FT-IR Difference Spectroscopy. In Handbook of Vibra-
tional Spectroscopy; Chalmers, J. M., Griffiths, P. R., Eds.; John Wiley:
Chichester, UK, 2001; Vol. 5.
References and Notes
(1) Hohenberg, P.; Kohn, W. Phys. ReV. B 1964, 136, 864-870.
(2) Kohn, W.; Sham, L. J. Phys. ReV. A 1965, 140, 1133-1138.
(3) 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.; Malick,
D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.;
Ortiz, J. V.; 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.; Gonzalez, C.; Challacombe, M.; Gill, P. M.
W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon,
M.; Replogle, E. S.; Pople, J. A. Gaussian 98, revision A.5; Gaussian, Inc.:
Pittsburgh, PA, 1999.
(4) Hutter, J.; Alavi, A.; Deutsch, T.; Bernasconi, M.; Goedecker, S.;
Marx, D.; Tuckermann, T.; Parrinello, M. CPMD, Version 3.4.0; MPI fu¨r
Festko¨rperforschung and IBM Zurich Research Laboratory.
(5) Nonella, M.; Tavan, P. Chem. Phys. 1995, 199, 19-32.
(6) Zhou, X.; Mole, S. J.; Liu, R. Vib. Spectrosc. 1996, 12, 73-79.
(7) Neugebauer, J.; Hess, B. A. J. Chem. Phys. 2003, 118, 7215-
7225.
(22) Cepus, V.; Goody, R. S.; Scheidig, A. J.; Gerwert, K. Biochemistry
1998, 37, 10263-10271.
(23) Allin, C.; Gerwert, K. Biochemistry 2001, 40, 3037-3046.
(24) Allin, C.; Ahmadian, M. R.; Wittinghofer, A.; Gerwert, K. Proc.
Natl. Acad. Sci. 2001, 98, 7754-7759.
(25) Baril, J.; Max, J.-J.; Chapados, C. Can. J. Chem. 2000, 78, 490-
507.
(26) Becke, A. D. Phys. ReV. A 1988, 38, 3098-3100.
(27) Perdew, J. P. Phys. ReV. B 1986, 33, 8822-8824.
(28) Boesch, S. E.; Wheeler, R. A. J. Phys. Chem. 1995, 99, 8125-
8134.
(29) Hertwig, R. H.; Koch, W. J. Comput. Chem. 1995, 16, 576-585.
(30) Becke, A. D. J. Chem. Phys. 1993, 98, 1372-1377, 5648-5652.
(31) Lee, C.; Yang, W.; Parr, R. C. Phys. ReV. B 1988, 37, 785-789.
(32) Hartwigsen, C.; Goedecker, S.; Hutter, J. Phys. ReV. B 1998, 58,
3641-3662.
(33) Jorgensen, W. L.; et al. J. Chem. Phys. 1983, 79, 926-935.
(34) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.;
Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217.
(35) Neria, E.; Fischer, S.; Karplus, M. J. Chem. Phys. 1996, 105, 1902-
1921.
(36) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola,
A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690.
(37) Mathias, G.; Egwolf, B.; Nonella, M.; Tavan, P. J. Chem. Phys.
2003, 118, 10847-10860.
(38) Niedermeier, C.; Tavan, P. J. Chem. Phys. 1994, 101, 734-748.
(39) Niedermeier, C.; Tavan, P. Mol. Simul. 1996, 17, 57-66.
(40) Herres, W.; Gronholz, J. Comput. Appl. Lab. 1984, 2, 216-220.
(41) Laasonen, K.; Sprik, M.; Parrinello, M.; Car, R. J. Chem. Phys.
1993, 99, 9080-9089.
(8) Eichinger, M.; Tavan, P.; Hutter, J.; Parrinello, M. J. Chem. Phys.
1999, 110, 10452-10467.
(9) Eichinger, M. Thesis, Ludwig-Maximilians Universita¨t Mu¨nchen,
Germany, 1999.
(10) Cui, Q.; Karplus, M. J. Chem. Phys. 2000, 112, 1133-1149.
(11) Nonella, M.; Mathias, G.; Tavan, P. J. Phys. Chem. A 2003, 107,
8638-8647.
(12) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227-249.
(13) MacKerell, A. D.; et al. J. Phys. Chem. 1998, 102, 3586-3616.