New Theoretical Insight into the Interactions and Properties of Formic Acid: Development of a Quantum-based Pair Potential for Formic Acid

We performed ab initio quantum-chemical studies for the development of intra-and intermolecular interaction potentials for formic acid for use in molecular-dynamics simulations of formic acid molecular crystal. The formic acid structures considered in the ab initio studies include both the cis and trans monomers which are the conformers that have been postulated as part of chains constituting liquid and crystal phases under extreme conditions. Although the cis to trans transformation is not energetically favored, the trans isomer was found as a component of stable gas-phase species. Our decomposition scheme for the interaction energy indicates that the hydrogen-bonded complexes are dominated by the Hartree-Fock forces while parallel clusters are stabilized by the electron correlation energy. The calculated three-body and higher interactions are found to be negligible, thus rationalizing the development of an atom-atom pair potential for formic acid based on high-level ab initio calculations of small formic acid clusters. Here we present an atom-atom pair potential that includes both intra-and inter molecular degrees of freedom for formic acid. The newly developed pair potential is used to examine formic acid in the condensed phase via molecular-dynamics simulations. The isothermal compression under hydrostatic pressure obtained from molecular-dynamics simulations is in good agreement with experiment. Further, the calculated equilibrium melting temperature is found to be in good agreement with experiment.


INTRODUCTION
The simplest carboxylic acid, formic acid ͑FA͒, has been the subject of numerous studies and still provokes considerable interest, even for small gas-phase complexes. 1 The two geometrical isomers of formic acid, namely, cis and trans of the hydroxy group relative to the carbonyl group, exhibit distinct electrostatic properties that may play an important role in the processes of molecular cluster aggregation. 2 Despite its simplicity, formic acid monomer offers a number of possible arrangements to form complex aggregates, which are stabilized by hydrogen bonding ͑HB͒. FA molecules in the gas phase aggregate to form mainly cyclic dimers, consisting of either eight-membered rings with two ͓͑C v ͒O...H-O͔ hydrogen bonds ͑FAD1 in Fig. 1͒ ͑although other isomers are also present͒, 3  Chocholousova et al. 3 have found dimers of formic acid in the gas phase that contain cis monomers as components.
To the contrary, low-temperature crystalline FA, as determined by x-ray and neutron powder diffraction, 4-6 contains infinite planar chains of hydrogen-bonded cis formic acid molecules. The neighboring pairs on the same chain possess the same structure as the FAD2 isomer ͑Fig. 1͒ found in the gas phase. On the basis of the IR spectra, it has been proposed 7 that mixed cis and trans polymorphic structures exist at elevated temperatures, confirming the coexistence of FAD2 and FAD4 ͑Fig. 2͒ dimers, whereas FAD4 is not observed at lower temperatures. The FAD4 structure contains a seven-membered ring of trans monomers connected through the ͓͑C v ͒O...H-O͔ and ͓C-H…O͑H͔͒ hydrogen bonds. The high-pressure Raman study of Schimizu 8 revealed a change in the slopes of O-D and C v O stretch bands indicating a phase transition between the low-pressure planar chains of cis-cis dimers ͑FAD2͒ to a high-pressure transtrans ͑FAD4͒ organization. 8 Furthermore a high-temperature and high-pressure study by Allan and Clark 9 using the combination of single-crystal x-ray-diffraction and ab initio density-functional calculations, yielded a different crystal structure from that found at ambient pressure and low temperatures. It was concluded that in the high-temperaturehigh-pressure crystal, molecules align on planes with alternating cis and trans monomer conformations ͓Fig. 2͑c͔͒. A more recent x-ray-diffraction and IR-absorption study 10 confirms the softening of the O-D and C v O stretch bands, but finds the crystal structure to be identical to that at the low temperature. 8 The results of the high-pressure studies are ambiguous, [8][9][10] thus necessitating further studies.
Liquid formic acid has been studied by the use of x-ray-diffraction, 11 neutron-scattering, 12,13 infrared, and Raman 14 techniques. Monte Carlo ͑MC͒ simulations have also been applied. 15, 16 Jedlovszky and Turi 15,17 ͑JT͒ have developed pairwise intermolecular potentials based on LJ ͑12-6͒ potentials using density-functional theory ͑DFT͒ computations of the dimer with empirically fit charges. However, they did not include the intramolecular terms, nor did they include the effect of trans isomers or higher-order many-body interactions beyond the dimer. Experimental studies have deduced that preferred configurations are chains, cyclic dimers, or other more complicated structures. The MC study revealed that about 7% of the molecules form cyclic dimers. These simulations have also suggested that the structures in the liquid and crystalline phases exhibit dramatic contrasts. The common feature of both phases is the simultaneous presence of the O-H…O and C-H…O interactions. 13 Experimental studies have not yielded much information regarding the cis or trans forms of monomers.
FIG. 1. The most abundant cis-cis dimers in the gas phase. Distances of hydrogen-bonded bridges were determined at the MP2/aug-cc-pVDZ, MP2/ aug-pVTZ ͑in parentheses͒ levels, and experimental ͑Refs. 5 and 44͒ values ͑square brackets͒ for the gas phase and crystal values are given for FAD1 and FAD2, respectively. All bond distances are in angstroms. Only the cis monomer was considered in the previous MC simulation. 15,16 Analogous to the crystal phase, the trans conformer could play an important role in the liquid at hightemperature and -pressure conditions.
The available experimental data for condensed phases of formic acid is rather restricted and its interpretation leads to many ambiguities. It is especially true for materials under extreme conditions. It is evident from recent discussions 12 that there is a clear and compelling need for further experimental and theoretical studies. The available theoretical studies on condensed phases are restricted to Monte Carlo and molecular-dynamics ͑MD͒ simulations of the liquid under mild conditions, where only cis isomers were considered. 18 The intermolecular atom-atom potentials were developed exclusively for cis FA isomers. 15,18,19 The need for the extension of theoretical approaches to include trans conformers is explored in this paper, using a substantially larger basis set, and thus, this work is computationally more intensive than previous studies. 15 The quantum-chemical calculations for the most important cis and trans conformers of monomers, dimers, and trimers are presented here. Moreover, we have applied an interaction energy decomposition scheme to provide new insight into the nature of the bonding in condensed phases. The quantum-chemical results are adopted for the extension of our modified JT intermolecular potentials to account for the possible cis-trans transformation. We have employed MD simulations, using our newly developed FA interaction potential to reproduce the ͑V , P͒ isotherm of crystalline FA over a wide range of pressure, and have compared our results with the experimental results of Goncharov et al. 20

Quantum-chemical methods
The quantum-chemical calculations were carried out using the MP2 method. 21 The geometry optimizations were carried out using the Berny algorithm with redundant internal coordinates and the location of proper minima and transition states were confirmed. The calculations were performed using aug-cc-pVDZ and aug-cc-pVTZ basis sets. 22,23 Properties of the formic acid monomer were also studied using a very extended aug-cc-pVQZ basis set. 23,24 The electronic density distribution was studied by applying the Mertz-Kollman ͑MK͒ population analysis scheme. 25 The MK charges were shown to properly reproduce the electrostatic interactions in the molecular crystal of triaminotrinitrobenzene ͑TATB͒. 26 The ab initio MK charges thus derived were used in the MD simulations.
The intermolecular interactions were studied using a hybrid variational-perturbational interaction energy decomposition scheme. 27 In this approach the self-consistent-field ͑SCF͒ interaction energy is partitioned into first-order electrostatic el ͑10͒ , Heitler-London exchange ex HL , and higherorder delocalization ⌬E del HF energy terms, The ⌬E del HF component accounts for the induction and exchange-induction effects associated with the relaxation of the electronic densities of monomer interactions, whereas the Heitler-London terms ͑ el ͑10͒ and ex HL ͒ represent the electrostatic interactions and exchange repulsion between the subsystems for which the electron-density distributions are unperturbed by the presence of a neighboring molecule. The electron correlation effects are taken into account by means of the MP2 method. The MP ͑2͒ interaction energy term, which includes the dispersion contribution and correlation corrections to the Hartree-Fock components, is calculated in the supermolecular approach as a difference between the appropriate second-order Møller-Plesset perturbation theory ͑MPPT͒ energy corrections, All interaction energy terms are calculated consistently in the dimer centered basis set and are therefore free from the basis set superposition error ͑BSSE͒ due to a full counterpoise correction scheme. 28 Such a scheme can be extended to more complex molecular systems comprising those beyond dimers, for example, trimers and tetramers. The total interaction energy of a system consisting of three subunits can be expressed as a sum of the two-and three-body interaction energies. 29 We can also generalize this scheme to encompass higher-order decompositions of exchange, delocalization, and correlation components into two-and three-body contributions 30 by taking into account the additive nature of electrostatic interactions and the applicability of the Löwdin formula 31 for the Heitler-London interaction energy to any n-body system. 32 The intermolecular interaction energy corrected for the energy of structural deformation of monomers constitutes the dissociation energy of the complex ͑D e ͒. The further correction for the zero-point vibrational energy leads to the spectroscopic dissociation energy ͑D 0 ͒. All quantum-chemical calculations have been performed using the GAUSSIAN 98 suite of codes. 33 The interaction energy decomposition scheme was implemented 34 in the GAMESS program. 35

Molecular-dynamics method
All MD simulations were performed on condensed phase FA using the new formic acid atom-atom pair potentials presented below. The FA condensed phase structures were simulated using three-dimensional ͑3D͒ cubic periodic boundary conditions. Simulations of ensembles of purely crystalline FA and ensembles of coexisting crystalline solid and liquid FA phases were performed, where the initial FA crystal structures used in these simulations were taken from Albinati et al. 5 The periodic MD simulation cell of crystalline FA consisted of 648 formic acid molecules, which were obtained by replicating the experimental unit cell three, nine, and six times in the a, b, and c crystallographic axes, respectively. The periodic MD simulation cell of the coexisting two-phase crystalline-liquid FA ensembles consisted of the crystalline ensemble ͑same as above͒, which was brought into contact, along the ͑100͒ crystallographic plane, with 648 randomly distributed "liquid" FA molecules.
The crystalline FA ensemble was initially energy minimized to relax the ensembles during which the cell dimen-sions of the simulation cell were allowed to adjust. Following this step, as a pretreatment for overcoming local minima by imposing thermal energy, constant particle number, pressure, and temperature ͑NPT͒ MD simulations were carried out for 100 ps at 250 K to equilibrate the system and allow the simulation cell to dynamically adjust. Thus, polymorphic transitions of the FA unit cell are possible. After this step, a final data collection NPT MD equilibration simulation was performed for another 100 ps at 250 K and a hydrostatic pressure of 0 GPa. Using this equilibrated ensemble as a starting point, the pressure-volume relation was determined for crystalline FA by increasing the hydrostatic pressure during NPT runs in increments of 1 GPa ͑0 -10-GPa pressure range͒ or 5 GPa ͑10-70-GPa pressure range͒ to a maximum pressure of 70 GPa ͑23 total pressures͒. The model FA crystal was allowed to dynamically evolve at each pressure for a total of 200 ps ͑pressures and volumes were time-averaged over a period of 100 ps after an equilibration time of 100 ps͒ before the hydrostatic pressure was increased. All computations were carried out using large-scale atomic/molecular massively parallel simulator ͑LAMMPS͒ code. 36 The equations of motion were integrated using the Verlet algorithm 37 with a time step of 1.0 fs. A Nosé-Hoover-type thermostat 38 with a relaxation time of 0.1 ps was used to control the temperature, and the pressure was controlled isotropically. 39 The particleparticle particle-mesh Ewald ͑PPPM͒ method 40 ͑the accuracy criterion was set to 1 ϫ 10 −5 and the near-field cutoff to 10.0 Å͒ was used for the long-range correction of electrostatic interactions.

Potential functions
Our newly developed FA force field consists of valence, nonbonded van der Waals, and Coulombic interaction terms. The functional form of the potential employs harmonic bond stretching and bond angle bending and a two-term Fourier expansion for dihedrals, where all valence degrees of freedom were explicitly treated and unconstrained. The atomic charges for the nonbonded Coulombic interaction are parametrized using the MK method ͑described above͒. The van der Waals interaction parameters are those of Jedlovszky et al., 15,17 implemented using the well-known Lennard-Jones ͑LJ͒ potential, however, we implemented the LJ potential using an inverse ninth-power term for the repulsive part rather than the 12th-power term as used by Jedlovszky and Turi. 15,17 The LJ ͑9-6͒ potential was implemented due to the softer repulsive core, which ultimately provided better agreement to experiment under pressure ͑see below͒.
The nonbonded LJ interactions were treated by truncating atom pairs with interatomic distances greater than 10 Å. Furthermore, to avoid steric overlap, nonbonded intramolecular interactions were not included, since all intramolecular interactions are accounted for explicitly in the valence interaction terms.

Formic acid
The molecule possesses two stable geometrical isomers in the gas phase. The ground-state minimum is represented by a planar cis structure with respect to the hydroxy and carbonyl groups. The corresponding trans isomer is higher in energy by 4.24 kcal/ mol at the MP2/aug-cc-pVQZ level. The calculated energy difference reproduces the measured data of 3.90± 0.09 kcal/ mol. 4 More importantly the transition state for the cis-trans rearrangement lies 12.97 kcal/ mol above the lower minimum. The calculated structural parameters of the monomer are within the error bars for the gasphase measurements 41,42 ͑Table I͒. The subtle differences between the cis and trans isomers are also well reproduced. Interestingly, the gas phase and crystal structures of formic acid are quite close. The structural changes corresponding to different points on the potential-energy surface for the cistrans rotation are negligible ͑Table I͒. The potential-energy curves representing a constrained rotation of the hydroxy hydrogen and the rotation with the rest of molecule fully relaxed are practically indistinguishable ͑see Fig. 3͒. The dipole moment to a large extent depends on the position of hydroxy hydrogen. Available experimental dipole moments 2 for the cis ͑1.420 D͒ and trans ͑3.790 D͒ structures are well

FA dimers
The interactions among dimer fragments play a major role in the determination of the overall interaction energy and the structures of complex aggregates. The first group of species studied here consists of the most stable dimers in the gas phase originating from the cis and trans monomers ͑see Figs. 1 and 2͒. The most stable FAD1 dimer has also been detected in the liquid state while the most stable trans-trans and cis-trans formic acid complexes in the gas phase are considered to be the precursors to high-pressure crystals. 7,8 Two open structure dimers that we have considered are the building blocks of potential arrangements in crystals ͑see Fig. 4͒. A third group consists of dimers that are parts of different chains found in the crystalline phase. This set of geometrical arrangements that we call a parallel arrangement exists only in the crystal environment and the gas-phase optimization of these structures results in hydrogen-bonded dimers that have already been considered ͑see Fig. 5͒.
As can be seen from Table III, the MP2/aug-cc-pVTZ level of theory reproduces the available experimental gasphase structure for the cyclic dimer of formic acid ͑FAD1͒ ͑Ref. 44͒ quite well. More importantly the reduction of the size of the basis set to aug-cc-pVDZ only slightly alters the  geometry. The experimental dissociation energy 45,46 is between 13 and 15.5 kcal/ mol considering the error bars. As can be seen from Table III, our calculations yield a dissociation energy within this range of values. The presented calculations also agree well with the best theoretical estimates available in the literature 47,48 validating the level of accuracy of our calculations for other larger species and other isomers.
The second most abundant FA dimer in the gas phase, FAD2, ͑see Fig. 1͒ is a precursor for the crystal structure under mild conditions. Its properties have not yet been measured in the gas phase. The available crystal geometry, however, is well reproduced by our optimized structure ͑Fig. 1͒, indicating that even structures that contain fragments sensitive to intermolecular interactions are computed quite accurately. The inspection of cyclic dimers indicates that parameters of hydrogen bond fragments ͓͑C v ͒O...H-O͔ and ͓͑͑H͒O…H-O͔͒ are to a large extent preserved in different dimers.
A hypothetical open cis-cis dimer proposed as a possible arrangement of the chain in the polymorphic crystal 6 when optimized in the gas phase transforms to the energetically lower six-membered cyclic complex ͓FAD7, Fig. 4͑a͔͒. The structure resulting from the optimization was found to be the more abundant isomer in the gas phase. 3 The hydrogen bonding characteristics of this structure is similar to those found in seven-membered ring isomers. The gas-phase optimized open cis-trans dimer ͓FAD8, Fig. 4͑b͔͒ is similar to the experimentally determined fragment of the crystal under high pressure. 9 This indicates a higher propensity of molecules in the crystal to form dimers under very high pressures.
As can be seen from Table IV, the interactions of hydrogen-bonded formic acid dimers are predominantly determined by the Hartree-Fock forces and the electron correlation contributions are not significant. The interactions of the studied complexes are mostly governed by the ͓O-H…O͔ hydrogen bonds. The group of isomers with ͓C v O...H͑O͔͒ hydrogen bonds ͑FAD1, FAD2, FAD4, and FAD6͒ exhibit interaction energies of −8.2, −9.1, −7.6, and −10.1 kcal/ mol, respectively. Note, that there are two hydrogen bonds in the FAD1 dimer complex, and thus the interaction energy listed in Table IV is twice the hydrogen bond interaction energy of −8.2 kcal/ mol. In the cases of FAD3 and FAD5 isomers, which possess hydrogen bonding of the ͓O-H…O͑H͔͒ type, interaction energies are significantly lower ͑Table IV͒. The decompositions of interaction energies are similar in all studied hydrogen-bonded complexes. To a large extent the attractive electrostatic interactions ͑ el ͑10͒ ͒ are neutralized by exchange repulsion interactions ͑ ex HL ͒. The major contribution to Hartree-Fock interactions originates from the delocalization energy term. The interactions in the cis-trans complex ͑FAD6͒ are 1 kcal/ mol higher than those in the cis-cis isomer ͑FAD2͒. The increase in the total energy of the molecular system due to the cis-trans transformation in one of the monomers is partially compensated by the formation of more stable hydrogen bonds. The dissociation energies are close to the interaction energies indicating a limited influence of the hydrogen bonding on the structures of monomers that comprise the dimers.
The parallel complexes are found in low-temperature crystals 5

FA trimers
We have studied a number of formic acid trimer structures that included the optimized hydrogen-bonded complexes shown in Fig. 6, as well as representative structures selected from the crystallographic data 5 shown in Fig. 5. The two complexes shown in Fig. 6 were selected from a large number of possible isomers representing hydrogen-bonded species. The linear chain of cis molecules found in lowtemperature crystals, 5 when optimized in the gas phase, transforms into a structure with two seven-membered rings analogous to the one found in the FAD2 dimer. The structure we label t-t-t ͑trans-trans-trans isomer, shown in Fig. 6͒, which is a fragment of chains in the high-pressure crystal, 8 was also optimized. The t-t-t trimer is stable in the gas phase and forms two seven-membered rings derived from its transtrans dimer, FAD4 parent isomer. Both optimized trimers represent the most compact structures characterized by the highest interaction energy of the dimers. The hydrogen bond-ing in these complexes is very similar to that observed in the corresponding dimers. As can be seen from Table V, for these trimers the interaction energies are dominated by the Hartree-Fock term with very little electron correlation energy contribution. The three-body contributions are below 10%, justifying the use of atom-atom pair potentials used in classical molecular-dynamics studies. The contribution to threebody interactions originates mostly from the delocalization energy. In the case of the t-t-t trimer, the three-body energy contribution to the delocalization energy amounts to 15% of that of the two-body energy. This energy contribution is twice the interaction energy of the FAD2 dimer ͑the component of the cis trimer͒ of 18.30 kcal/ mol, and almost perfectly matches with the magnitude of two-body interaction energy part of the trimer ͑18.62 kcal/ mol͒ which is consistent with the local character of structural changes.
The complexes selected from the crystal that belong to different molecular chains ͑cryst-1 and cryst-2 in Fig. 5͒ are influenced by the van der Waals interaction with the main contribution from the electron correlation energy ͑Table V͒. The three-body interactions make only minor contributions to the total interaction energy. The structure of the mixed trimer shown in Fig. 5͑c͒, which exhibits hydrogen bonding and van der Waals interacting parts, is consistent with the additive nature of interactions in the FA crystal. As concluded before, three-body interactions contribute less than TABLE IV. Interaction energy components and the dissociation energies for studied formic acid dimers. All values are in kcal/mol. The D 0 value includes the zero-point vibrational energy correction while the D e value is the absolute dissociation energy without zero-point energy correction.

FUNCTIONAL REPRESENTATION OF INTERMOLECULAR INTERACTIONS FOR THE FORMIC ACID CRYSTAL
Our ab initio quantum-chemical studies support the assumption of the additivity of intermolecular interactions. Thus, the interactions can be decomposed into an atom-atom pair potential, where the nonbonded interaction energy between two molecules is calculated according to the following equation: where E el represents the Coulombic electrostatic term and E LJ is a LJ 9-6 term. The expression for the Coulombic electrostatic interaction is given as where q i is the charge on atom i, 0 is the vacuum permittivity, and R ij is the distance between atoms i and j. The charges were calculated for the cis monomer using the MK population analysis for the MP2/aug-cc-pVQZ electronic density. For the remaining term in the nonbonded interaction potential, namely, the Lennard-Jones portion of the pair potential, we use the following combination rules: ij = ͱ i j and r 0 ij = 1 2 ͑r 0 i + r 0 j ͒. The parameters ij and r 0 ij were adopted from Jedlovszky and Turi. 17 The valence terms, developed here, used in FA interaction potential, consisted of harmonic bond stretch ͓E s =1/2k b ͑b − b 0 ͒ 2 ; two-body term͔, bond-angle bend ͓E =1/2k ͑ − 0 ͒ 2 ; three-body term͔, and dihedral angle torsion terms ͓E = k ͑1 − cos͑n͒͒; four-body term͔. The two-term torsional potential was fit to the quantumchemically derived rotational curve for the cis to trans conformational transformation ͑see Fig. 3͒. The atom-atom potential-energy parameters thus derived are given in Table  VI. No attempts have been made here to incorporate proton transfer in our newly developed atomistic pair potential presented here.

Force field evaluation and validation
The MD simulations of FA were carried out on a variety of condensed phase ensembles using our newly developed atomistic pair potential. All atoms in the FA molecules were treated explicitly.
The resulting pressure-volume isotherm curve of the bulk orthorhombic FA crystal was obtained by using the potential parameters in Table VI, and is shown in Fig. 7. The pressure-volume isotherm ͑at a constant temperature of 250 K͒ varies over a hydrostatic pressures range of 0 -70 GPa. The simulation results are compared with the experimental crystal data of Goncharov et al. 20 The isotherm calculated using the potential parameters of Jedlovszky and Turi 16 ͓using the LJ ͑12-6͒ form in conjunction with our  is the well depth in kcal/mol; r 0 the distance at the minimum of the pair potential in Å; charges ͑computed in this work͒ in electron units, 332.0637 converts energies to kcal/mol. and r 0 parameters taken from Ref. 17. newly developed valence terms͔ are also shown. As seen from Fig. 7, our computed results, using the LJ ͑9-6͒ potential, agree well with the experimental data over the entire pressure range. However, the results using the LJ ͑12-6͒ potential form resulted in less satisfactory overall agreement. Furthermore, no cis to trans conformational changes were identified over the entire pressure range studied. Thus, the simulations performed here do not support the proposed pressure-induced phase transition suggested by Allan and Clark. 9 To further validate our FA potential, we focus on determining the melting temperature for our model FA at a constant hydrostatic pressure of 1 atm using the approach described in Ref. 50. To this end, an NPT simulation of 648 crystalline FA molecules was performed near the estimated melting temperature to achieve an approximate equilibrated initial condition. Next, a liquid, consisting of 648 FA molecules, randomly placed into a simulation cell of dimensions equal to the dimensions of the equilibrated solid ͑from above͒, was created. The liquid simulation cell was fixed in the y and z dimensions, but allowed to relax in the x dimension at 450 K until equilibrium was reached. The liquid was then cooled to 250 K and again equilibrated. The equilibrated crystalline and liquid simulation cells were then brought into contact by joining the two simulation cell along the x axis ͓͑100͒ crystallographic plane͔, thus creating a supercell ͑1296 total FA molecules͒ with two interfaces between the solid and liquid systems. The two-phase supercell ensemble was then equilibrated for 1 ns at 250 K under NPT conditions, where all three cell dimensions were allowed to relax independently ͑anisotropic relaxation͒. The melting temperature was determined by slowly heating the supercell and locating the point at which the density peaks at the solidliquid interface in the density profile begin to diminish, which is indicative of melting. The crystalline structure in the density profiles noticeably diminishes at ϳ265 K ͑our computed melting temperature͒, where the experimental melting temperature has been determined to be 281 K ͑Ref. 10͒ ͑giving an error of less than 6%͒.

CONCLUSIONS
Formic acid both in gas-phase and condensed forms exhibits a variety of structures and poses considerable challenges for theoretical and experimental studies. Despite its simplicity, even the monomer of formic acid exhibits a number of possible orientations to form complex aggregates of clusters, thus exemplifying the variety and complexity of formic acid forms. The stable chain complexes are formed due to hydrogen-bond interactions. The short chains are observed in the liquid phase 13 while the crystalline phase is built from infinite planar hydrogen-bonded chains which are weakly connected through van der Waals interactions. 8 The composition of FA chains in the high-pressure crystals continues to be a topic of considerable controversy. 10 An important feature that we have discussed is the structure of the monomer ͑cis or trans͒ that forms crystal chains. The less stable trans isomer of formic acid is separated from the cis isomer by an energy barrier of 8.7 kcal/ mol. It is interesting to note that Wei et al. 51 have computed the barrier for proton transfer to be 9.8 kcal/ mol, which is quite comparable to that of the barrier for cis-trans conversion. Moreover the barrier for internal rotation is expected to be lowered in liquid state due to the polar nature of the trans isomer. Since the barrier for proton transfer is comparable to that of the internal rotation, we may expect, and have shown, 20 that proton transfer can play a role in the organization of crystal structure and liquid state especially at high temperatures and pressures. The cis-cis dimer with a single ͓͑C=͒O...H-O͔ hydrogen bond ͑FAD2͒ is a precursor for structures found in lowtemperature crystals of formic acid. It was inferred that highpressure crystals are composed of trans-trans ͑FAD4͒ and cis-trans ͑FAD6͒ dimers. The hydrogen bonding in the complexes studied here is similar in nature regardless of their compositions. The hydrogen-bonded complexes are predominantly influenced by the Hartree-Fock forces, while van der Waals forces play an important role for the parallel complexes. The energetically unfavorable transformation of FA from cis to trans form is partially compensated by the increased stabilization of hydrogen bonding. The three-body and higher interactions are reasonably small, thus making the development of an atom-atom pair potential feasible. The success of the proposed FA model presented here in reproducing the experimental pressure-volume isotherm as well as the melting temperature for FA illustrates the feasibility of developing intermolecular potential functions using an ab initio approach. Finally, we find no evidence of a pressure-induced phase transition in crystalline FA as proposed by Allan and Clark. 9

ACKNOWLEDGMENTS
The work at LLNL was performed in part under the auspices of the US Department of Energy by the University of California, LLNL under Contract No. W-7405-Eng-48. The   FIG. 7. Pressure-volume isotherm for crystalline formic acid, as determined from MD. All MD results are performed using the formic acid intermolecular Lennard-Jones parameters of Jedlovszky and Turi ͑Ref. 15͒, coupled with formic acid intramolecular, Coulombic parameters, and MK charges presented in this paper. Experimental values ͑Ref. 20͒ are at room temperature, shown as symbols ‫,͒ء͑‬ while simulation results are shown as curves; the Lennard-Jones 12-6 intermolecular formic acid potential at 250 K ͑dotted curve͒; the Lennard Jones 9-6 intermolecular formic acid potential at 250 K ͑solid curve͒.