Ab initio molecular dynamics study of overtone excitations in formic acid and its water complex

In this article, we present results from ab initio molecular dynamics simulation of overtone excitation in formic acid monomer and its water complex in the gas phase. For the monomer, a conformation change is observed employing both OH and CH vibrational excitations, which supports experimental findings. In the formic acid–water complex, interconversion also takes place, but it proceeds via hydrogen exchange rather than via intramolecular reaction. Simulations raise a question on effect of quantum and matrix effects to the results. Also, a brief test of different computation methods was done on the system.


Introduction
The high overtone induced processes are initiated by promoting a molecule or a complex to chemically relevant energies in high vibrational excited states through the thermal excitation or direct absorption of visible or near-infrared radiation.Unlike the thermally induced processes, the reactions initiated by vibrational overtonepumping of ground electronic state molecule through direct absorption of photon can be treated without the need to consider collisional activation processes.In general, it is required that the energy is deposited into the initially excited vibrational state and subsequently transferred by intramolecular vibrational redistribution to other modes of the molecule including the reaction coordinate.For example, it has been demonstrated that high overtone excitation of OH groups in carboxylic acids isolated in low temperature matrices leads to two dierent unimolecular processes: isomerization and/or decomposition of the molecules [1,2,3,4].Moreover, overtone pumping was used, for example, in the case of formic acid dimers to produce high-energy conformers within the hydrogen-bonded complexes [5,6,7].
In case of formic acid it has been shown to change conformation from trans to cis [1], see gure 1, and that formed cis conformer reorganises back to to trans form via a tunneling mechanism [8].But in formic acid water complex this conversion back from trans to cis does not happen [9].
The aim of this study was to try to simulate using ab-initio molecular dynamics these conformation changes that have been observed in experiments.In particular the aim was to try to see if 532 nm excitation energy could produce these results, as this has been the laser wavelength used to produce high-overtone interconversions, for example, for formic, oxalic and glyoxylic acids [2,10].

Method
Ab-initio molecular dynamics simulations were done with Car-Parrinello method [11,12], using CPMD v. 4.1 [13].Calculations were done using BLYP functional with D2 dispersion correction, using Troullier-Martins pseudo potentials [14] with 90 ry cuto, 3 au time step and 400 au ctitious electronic mass.Simulation box was cubic 10 Å for monomer and 14 Å for water com- plex.Periodic images were decoupled using Martyna-Tuckerman Poisson solver [15].Simulation box sizes were veried by increasing box sizes and doing geometry optimization, which had no dierence in energy or geometry.The plane wave cuto was chosen because it has been used on a similar system [16].
To generate trajectories system was rst propagated using velocity rescaling thermostat for ∼ 1 ps, in NVT conditions, followed by ∼ 5 ps using Nose-Hoover chain thermostat [12,17], in global mode.After which a total of 100 snapshots were saved in every ∼ 1 ps, while continuing with Nose-Hoover chain thermostat.
To simulate vibration excitations velocities in these snapshots were modied, so that kinetic energy is increased by a desired amount.According to hypothesis energy is added to OH stretching vibration.Also for formic acid monomer CH stretching excitations were done.Velocities were changed so that center of mass velocity does not change and that O and H atoms move away form each other.Due to mass dierence this means that hydrogen will take most of the added kinetic energy in relation to mass dierence.It also should be noted that because only O and H atoms, which form OH bond, have their velocities changed, there is also a little bit of energy added to COH bending and OCOH torsion motions.After adding velocities each snapshot was propagated forward in NVE conditions, thus creating a total of 100 trajectories.
For formic acid monomer the trajectories were propagated for about 1 ps.For water complex the trajectories had to be cut to about 200 fs, because complex dissociates due to added energy and fragments y out of the simulation box.Also water complex is not stable in 300 K, and thus only 15 K and 200 K trajectories were generated, while for monomer 15 K, 100 K, 200 K and 300 K trajectories were generated.These temperatures are arbitrary chosen, to see how temperature aects the system and such they have no specic meaning.Also when temperature is is addressed later in this article, it always refers to these values before excitation was applied.
In these trajectories we are interested in question is there is a a conformation change.This means that results from these trajectories are in form of probabilities to nd something, which can be considered to behave like a binomial-distribution.To assess these probabilities and errors Bayesian-statistics [18] was used.In Bayesian-statistics binomial-distribution has a Bayesianinterface via beta-distribution, from which a probability to nd a desired result can be taken as mean and error from associated cumulative distribution function.Even likelihood, β(1, 1)-distribution, was used as a prior.In presented results errors (dotted lines) represent 95% condence interval.

Monomer
Results from trajectories excited by 532 nm excitations to OH vibration are shown in gure 2. The most prominent feature is that conformation change is more probable in higher temperatures and it was not observed at all in 15 K.In case of CH excitation, gure 3, the conformation change is much more probable and temperature has only a minor eect.Conformation change was observed when excitation energy 1200 nm (8300 cm -1 ) or more, when CH bond was excited.We didn't nd any conformation changes with 1700 nm (6000 cm -1 ) excitation energy.However due to statistical nature of results, we can only say this in about 95 % certainty.
With OH excitation conformation change was observed when energy was 700 nm (14000 cm -1 ) or more and was not observed when energy was 800 nm (12500 cm -1 ).
To see how energy is transfered in molecule movement of certain parameters describing molecule was visualized see gures 4 and 5. Kinetic energy is related to square of change, so parameters chosen are square of change in distance and square of change in angle, although it should be noted that for e.g.bending motion kinetic energy depend on angular velocity and bond lengths, so care must be kept in mind when interpreting these gures.CH excitation data, gure 5, diers from OH excitation that OCH bending gains initially some energy and then starts to degrease.This points on possibility of this mode being part of energy channel and to test this a set of trajectories was made with a constraint on OCH motion.This lowered conformation change probability by ten fold of its initial value.From the gure can also be seen that COH bending motion is excited before OCOH torsion, so we could expect it to be part of energy channel, as it is in OH excitation.To test this a set of trajectories was again generated with a constraint on COH angle.This lowered conformation change probability to one third of its initial value.This suggest that there are two channels for energy transfer a one from CH stretch -> OCH bending -> OCOH torsion and an other from CH stretch -> COH bending -> OCOH torsion.However these two depend on each other, as the OCH constraint decreases movement in COH bending considerably.Which suggest that a third channel exists where energy transfers along path CH stretch -> OCH bending -> COH bending -> OCOH torsion.Also worth of mentioning is that energy added to COH bending and OCOH torsion motions due to our method of exciting molecule is small, as can be seen from the gures at time zero, and thus is not expected to aect results.

Water Complex
The most prominent feature of trajectories is hydrogen exchange with water molecule (see gure 6).This happens in all trajectories sampled from 15 K temperature (gure 7), while at 200 K the exchange probability diminishes due to increased rotation of the water to what is represented in gure 9.This would suggest that there are two dierent reaction paths and associating transition states.This however is not the case.
By inspecting trajectories more reveals that the rarer case is never present in rst wave of exchanges, which can be seen in gure 7 as the rst peak in around 25 fs.Also in the transition state the distance between water and formic acid is about 0.5 Å shorter than what it's in minimum energy geometry of the complex.So to get to the true transition state, water and formic acid would have to move closer.Here comes the fact that we have a lot of energy and that hydrogens move a lot of faster than oxygen.So when hydrogen is transfered to the water, the oxygen does not have enough time to move closer and instead one of the hydrogens will just y away from the the H 3 O + ion.In an analogy this could be thought of like playing billiard, when one ball hits hard a group of balls next to each other, or like Newton's cradle.
Energies calculated with dierent methods and basis for the transition state, gure 9, are listed in tables 1 and 2. These values dier a lot from method to method, which is inspected more later.This calculated energy represent the minimum energy path.But the true energy depends on how long the distance is between the two molecules, which can be explained by a potential energy between two ions increasing with distance.When the distance is around the minimum energy geometry of the complex transition state energy is a little bit over 12000 cm -1 , which can be estimated by a simple calculation with constraining the distance or by potential energy along a trajectory (see gure 10).
It should be noted too that H 3 O + and formate ions is not a local minimum energy conguration, when their distance is around the minimum energy distance of the complex.Only when the distance would be more than 0.5 Å longer would it be a local minimum.
All this sums up that depending of reaction path the exchange reaction can have over 7000 cm -1 dierence in potential energy barrier.

Comparing dierent calculation methods
The large spread of transition state energies between dierent methods, table 1, in raised a question on accuracy of dierent methods in general are describing geometries present in trajectories.To get more information in this matter a few trajectories where picked and potential energy was calculated with dierent methods among them.Some high lights of these calculations are depicted in gures 10, 11 and 12. Basis sets used in these  The most important information from these calculations is that ions are problematic.When ever formic acid water complex is in ionic form there is a large spread in potential energy among dierent methods.While all methods give in practice same results when Fig. 9 Transition state geometry for the hydrogen exchange reaction, DLPNO-CCSD(T)/aug-cc-pVTZ. Note that there is a symmetry plane in the middle.system is in form of formic acid and water molecules.
The energy dierence can be as high as 5000 cm -1 .
If we make a assumption that DLPNO-CCSD(T) would give the best results, we can then rank the other methods against it.This wound mean that BLYP is the worst method clearly.It does not give the correct potential energy curve shape and it gives energies that can be thousands of wave numbers o of values given by DLPNO-CCSD(T).B3LYP is a large improvement over BLYP, but it too fails when system is ionic in form.It does however give the correct shape of potential energy curve.MP2 is within a few exception better option than B3LYP.A special note also to fact that D2 performs better than D3BJ.
For formic acid monomer trajectories, gure 12, the dierence between dierent methods is a lot smaller.
Only BLYP gives considerably dierent values from the other methods.In practice this means that when OH or CH bond is long enough BLYP fails.B3LYP gives very good values for OH bond, but it too will fail a bit on CH bond.But the magnitude of dierence to DLPNO-CCSD(T) is of order for few hundred of wave numbers, while for BLYP it can be thousands.But these dierences are present only when energy is over 15000 cm -1 .MP2 and DLPNO-CCSD(T) on the other hand give in practice same values.

Monomer
Simulations show conformation change as is observed in experiments when the OH bond is excited, but not at the temperatures used in experiments.This raises a question on role of quantum eects in the process.
We already know that transition back from cis to trans happens with tunneling [8] and the whole process is about movement of hydrogen.So we can expect quantum eects, like tunneling and zero-point motion, to be important in trans to cis interconversion too.The Trajectory start from the point when energy is added.Its initial point is sampled from 300 K simulation.So initial potential energy represents thermal excitation at the beginning.Excitation in this trajectory is 532 nm on OH. fact that we can not see any trans to cis change in simulation when nucleons are treated classically, can be considered as an evidence that not only are quantum eects important, but they are actually crucial for the hole process.
Our simulations showed that increasing temperature increases the probability of the conformational interconversion.This was expected, as formic acid in it's minimum energy conguration is a planar and for symmetry reasons the OH stretch, that is plane symmetric, can not excite torsion, which is not plane symmetric.
In quantum nucleons case zero-point motion (and uncertainty principle) distorts molecule a little bit from plane.Which is somewhat similar to eect of increased temperature.
In the case of CH excitation the importance of zeropoint motion is not that clear, as in the simulations the temperature, where initial points of trajectories were generated, had only a minor eect to the nal outcome of the simulations.The same symmetry reasoning applies also for the CH stretching mode.This implicates that the CH stretching motion is more strongly coupled to the torsional motion, than the OH stretching mode.Which is also evident based on the nding that the conformational interconversion is much more probable, when excitation is applied to the CH stretching than to the OH.
Because our simulation used classical nuclei, we didn't have tunneling in our simulations.Thus it's expected to need more energy to achieve the conformational interconversion.In our simulation we needed 14000 cm -1 with excitation on the OH bond, which is about two and a half times more than what is observed experimentally.This dierence does sound to be too large eect to be caused by miss of tunneling.As a comparison, in a recent study by Tew and Mizukami [20].They calculated formic acid vibration overtones and showed that 6th and 7th excited states in torsion-mode tunnel trough barrier to cis conformer.Which in practice meant that barrier of about 4500 cm -1 height is eec- tively only about 3700 cm -1 .So we could expect tun- neling to lower barrier by about 1000 cm -1 .This could mean that our simulations are to short and running a longer time would show the conformational interconversion happening with less energy.This brings us to the question of role of matrix, which absorbs all the thermal energy [1].But in what time scale?However dynamics studies of matrix and steric eects of the systems studied here are outside the scope of the current paper.
Third possibility is that we have a wrong excitation.
It is assumed that the excitation happens on the OH vibration and that when it happens with 532 nm it excites 6th vibrational state of the OH stretching, based on experimental parameters [21].But we do not know for sure.Thus it would be good develop a better method to apply an excitation than a simple velocity rescaling.
In addition to repeating simulations with matrix atoms, it might be a good idea to do path integral simulations, to see how much quantum eects matter.These methods have some problems them selfs, but could still be good to see how they dier from purely classical nuclei.
For formic acid monomer the poor performance of BLYP is probably not and issue.Because we did get same results when excitation energy was smaller that the area where BLYP miss performs.But there could be some eect on probability and to how long it takes for conformational interconversion to happen.

Water complex.
Hydrogen exchange is well along with chemical intuition and it's already studied in many cases [22].The novelty here is that it can be induced with vibrational excitations and that it can be an intermediate step in a total reaction.But because formic acid is symmetric it cannot be proven experimentally without isotope substitution.
While this study can not provide reason to why in formic acid water complex cis conformer does not tunnel back to trans form.It does imply that mechanism for conformational interconversion very likely includes hydrogen exchange and that this may also be a part of the reason why cis is stable, as was observed by Marushkevich et.al.[9].
A big question for the water complex calculations come from accuracy of the used method BLYP+D2.
Normally it would be good idea to do the calculations with a better method.But to do B3LYP calculations with CPMD would need appropriate pseudopotentials, which we do not have.But in the future B3LYP or MP2 calculation should be minimum for the water complex calculations.
Then is the question of, how does the poorly performing method aect results from these simulations.
If we look at the gure 10.We see that in that trajectory BLYP does get the shape right, it's only lower.
Thus we could expect to get same kind of trajectory, if we just use more energy.So the hydrogen exchange reaction does most likely happen with other methods too.It's only about how much energy it needs and how probable it is.And in any case, this analysis of performance of dierent methods will prove to be useful in future works regarding the water complex.

Summary
The simulation results show that conformational interconversion takes place when OH and CH stretching vibrations are excited, which are in agreement with previous experimental observations.But the simulations also raise questions regarding a role of quantum and matrix eects.For excitations upon the OH vibrational mode, conformational interconversion does not happen at matrix isolation temperatures according to the simulations presented here.This points to the role of zero-point motion.Tunneling is also most likely important, because in simulations much more energy was needed than what has been shown to work in experiments.
Water complex simulation could not explain explicitly the observed stability of the cis-HCOOH water complex.To address this issue, most likely a simulation with added matrix atoms is needed to mimic the experimental conditions better.In water complex simulations hydrogen exchange reaction was found and it's expected that trans to cis transformation happens via it.Testing accuracy of used method BLYP+D2 raised issues when calculating system when it's in ionic from, which will most likely aect the magnitude of excitation energy needed and reaction probability.

Fig. 6 Fig. 7
Fig. 6 Hydrogen exchange reaction from left to right after exciting formic acid OH bond.

Fig. 8
Fig. 8 Chance to nd H 3 O + ion after exciting (532 nm) formic acid OH bond in water complex.

Fig. 10
Fig.10 Potential energy calculated along a trajectory with dierent methods in formic acid water complex hydrogen exchange reaction.Initial point of this trajectory is sampled from 15 K simulation and it is the point when energy has been added to the system.So this trajectory represents what is visualized in gure 6, as a back and forth reaction.So formic acid and water exist at 0 fs, 45-65 fs and around 130 fs.While the structure is ionic in between.

Fig. 11 Fig. 12
Fig.11Potential energy calculated along a trajectory with dierent methods in formic acid water complex hydrogen exchange reaction.In this trajectory snapshot initial point is near minimum energy transition state and them moves somewhat around it to form formic acid and water around 60 fs, which then goes back to formate and oxonium ions again.

Table 1
Transition state energy comparison for dierent methods.Dispersion correction (D3BJ) was used with DFT methods.