DIADORIM : a Monte Carlo Program for Liquid Simulations including Quantum Mechanics and Molecular Mechanics ( QM / MM ) Facilities : Applications to Liquid Ethanol

Um programa computacional para a simulação de líquidos e soluções utilizando o Método de Monte Carlo com algoritmo de Metropolis (MMC) é apresentado. Neste programa a energia total de interação é dividida em clássica e quântica, QM/MM: a contribuição clássica é obtida utilizando campo de força e a contribuição quântica é calculada utilizando métodos de química quântica. As seguintes facilidades foram implementadas: (i) o programa MOPAC 6 foi incorporado como uma sub-rotina do programa principal; (ii) alternativamente, um banco de dados de energias previamente obtidas pode ser utilizado; (iii) Simulações nos ensembles isotérmico e isobárico, NpT, e canônico, NVT, podem ser realizadas; (iv) a trajetória obtida ao longo do processo de MMC pode ser armazenada em disco rígido, utilizando o formato xtc, para análise e visualização utilizando outros programas. O programa foi utilizado para calcular propriedades estruturais e termodinâmicas do líquido etanol no ensemble NpT a 1,0 atm e 298 K. Distâncias e ângulos de ligação das moléculas de etanol foram mantidos constantes, mas amostragem do ângulo diedral definido pelos átomos H-O-C-C foi efetuada. Para calcular a energia interna em função deste ângulo foram utilizados os métodos AM1, PM3 e MNDO implementados no programa MOPAC 6. Cálculos também foram efetuados utilizando energias obtidas com HF/6-31g*, MP2/6-31g* e B3LYP/6-31g*. Diferenças significativas foram observadas na distribuição dos ângulos diedrais na fase líquida. Entretanto, distribuições de energia para as interações de pares são muito semelhantes, o mesmo acontecendo com as distribuições radiais de pares. Resultados obtidos para a densidade e calor de vaporização do líquido estão em ótimo acordo com dados experimentais.


Introduction
Computer simulations methodologies such as Monte Carlo 1 and Molecular Dynamics 2 have been largely used to study thermodynamic and structural properties of condensed phase systems. 3,4As a general trend in such simulations classical potential functions are used to calculate the potential energy. 5,6A very large collection of achievements have been reported in the literature, showing the success of these methodologies to obtain insights on the behavior of complex systems at the molecular level. 7,8Nevertheless, it is well known that various important chemical issues such as bond breaking, bond making and molecular rearrangements cannot be describe by classical force fields and the inclusion of quantum mechanics is needed.To fulfill this need, methodologies including particle interactions via quantum chemistry calculations have been reported. 9,10evertheless, computer simulation of large systems using potential energy surface calculated with quantum mechanics is very time consuming.Therefore, hybrid methodologies based in the partitioning of the system phase space into two regions, the quantum motif and their surroundings, which is treated classically, are very useful. 11,124][15] In some methodologies the quantummotif is not included directly in the trajectory production and the QM calculations are performed a posteriori using molecular arrangements chosen from the classical force field trajectory. 168][19] Therefore, differently conceived quantum chemistry methodologies can be merged to build new computer simulation programs aimed to treat large chemical systems in its complexity.In this paper a Metropolis Monte Carlo simulation program incorporating the MOPAC 20 package as a routine is presented.Coupling to an external data bank containing energies previously calculated with other quantum chemistry software is also possible.Next, the main features of this program, which was named DIADORIM, will be described and some results obtained for the pure liquid ethanol (ET) will be presented.

Methodology
The DIADORIM QM/MM program was developed to study liquids and solutions using the Metropolis Monte Carlo algorithm. 1,4,21Some features and capability of this program are: (i)The configuration space is sampled using standard Metropolis Monte Carlo procedures in both canonical (NVT) and isothermic-isobaric (NpT) ensembles; (ii) intramolecular energies can be calculated using either a classical force field, the MOPAC program or read in from an external data bank; (iii) intermolecular interactions are calculated using classical force field; (iv) periodic boundary conditions, minimum image convention and cut-off ranges are used; (v) corrections for Lennard-Jones potential interactions beyond the cut-off range are taken into account using the methodology presented elsewhere. 4(vi) Long-range corrections for coulomb interactions beyond the cut-off range are considered using the reaction field formalism. 22(vii) Free energy differences can be calculated using the statistical perturbation theory (FEP). 23(viii) The trajectory generated along the Monte Carlo Metropolis sampling can be saved using the xtc file format; 24 (ix) interfaces with molecular visualization programs such as vmd, 25 molden 26 and the plotting tool xmgrace 27 are available.
Earlier versions of this program, not including QM/MM facilities, were extensively used to study pure liquids [28][29][30] and solutions [31][32][33] using classical force fields.The present QM/MM version was used to study liquid ethanol and the results are presented in the next sections.

The potential function
The quality of the potential energy surface describing the particle interactions is of major importance to obtain a reliable modeling of a given system.6][7] Assuming effective pair wise potentials, the total potential energy U total is partitioned as: The energy terms in equation 1 are discussed in the next sections.

The intermolecular potential function
Following usual procedures in force field calculations the molecules are modeled by collections of interacting sites and the intermolecular interaction potential is represented by a sum of Coulomb and Lennard-Jones potentials centered on the sites. 6Therefore, the energy E ab between molecules a and b is represented by a pairwise sum of Lennard-Jones and Coulomb potentials centered on the sites: Where r ij is the distance between site i in a and site j in b and q i and q j are point charges located on these sites.For a given site k, the parameters A kk and B kk are given by A kk = 4 ε k σ k 12 and B kk = 4 ε k σ k 6 , where ε k and σ k are the Lennard-Jones parameters for this site.Parameters A ij and B ij for a non-diagonal interaction [i,j] are obtained using the geometric mean combining rules A ij = (A ii A jj ) 1/2 and B ij = (B ii B jj ) 1/2 .Therefore, the U intermolecular energy contribution in equation 1 is obtained as sum of the E ab terms.

The intramolecular energy
It is well know the failure of force fields to calculate the exact dependence of internal molecular energy as a function of geometry. 6This difficulty is expected as force fields are classical representations of atomistic interactions and cannot be intended to substitute in full the formalism of quantum mechanics.Consequently, the dependence of molecular energy with internal degrees of freedom will never be adequately achieved by classical models.0][11][12][13][14][15][16][17][18] Nevertheless, despite the extraordinary achievements in computer power computational cost is still a limitation against the use of full ab initio quantum mechanical methods to study large system.To accomplish this need, semiempirical methodologies have been developed to obtain quantum mechanical information of large systems at a feasible computational cost. 34Among others, model Hamiltonians implemented in the computer programs MOPAC, 20 ZINDO, 35,36 AMPAC 37,38 have been successfully used to study a diversity of molecular systems and chemical processes.New parameterizations of MOPAC Hamiltonians are available, extending its capabilities and usefulness. 39,40onte Carlo QM/MM studies using MOPAC have been reported and the results obtained are very encouraging. 41,42n the present QM/MM version of DIADORIM program the semiempirical program MOPAC 6 was introduced as a subroutine of the main Monte Carlo code.Therefore, all the Hamiltonian models and other facilities implemented in the MOPAC program can be used to calculate the energy and other information for a given molecule or cluster along the phase space sampling.The data and keywords needed to obtain MOPAC results are passed from the main program to MOPAC module using standard subroutine call defined in the FORTRAN language.The energy calculated in MOPAC is returned to the main program and used on the fly to perform Metropolis sampling.Alternatively, quantum chemistry data obtained with other methodologies can be used to calculate the intramolecular energy.

The interaction between quantum and classical sub-systems
As described above, the MOPAC package is used to obtain the internal energy of a previously defined quantum sub-system.Therefore, as this sub-system may be solvated by others molecules, it is wise to include the solvent effects in the calculation of the wave function and properties.Two procedures are delineated to accomplish this task: (i) implicit solvent models such as the ones implemented by Thrular and Cramer; 43 (ii) using sparkles atoms included in the MOPAC Hamiltonian, that is, replacing nearby solvent atoms by point charges whose values can be calculated using the charge models proposed by Truhlar and Cramer. 44hese alternatives are being implemented but were not used in the present calculation.It is interesting to note that point charges obtained with the models proposed by Truhlar and Cramer 44 can be used to compute the Coulomb interaction energy between particles in the quantum-motif and solvent molecules.This procedure have been successfully used in the QM/MM method implemented by Jorgensen. 17,45,46herefore, for an arbitrary configuration, the total energy can be obtained as follows: (i) Intermolecular energies are calculated using classical force field parameters; (ii) For each molecule the internal molecular energy can be calculated either using a chosen semiempirical Hamiltonian implemented in MOPAC 6 or using a classical force field or read in from an external data bank previously defined by the user.
As a general representation, the models obtained with the hybrid methodology will be called QM-FF, with QM indicating the quantum chemistry level of calculation and FF the force field.
Using the definitions above, the total configurational energy can be calculated and straightforwardly used with standard Metropolis Monte Carlo formalism to sample the phase space.Results obtained for liquid ethanol using different QM/MM approaches to calculate the internal energy are presented in the next sections.

Molecular models and Monte Carlo simulation protocol
The QM/MM models and the protocols used in the present calculation are presented below.

Ethanol molecules
The intermolecular interactions energy between ethanol molecules were calculated using the OPLS-AA force field parameters. 8Bond distances and angles optimized for the trans conformation of ethanol molecule with MP2/6-31g* level of theory 47 were used and kept constant during the simulation.Rotation along the dihedral angle defined be H-O-C-C sequence of bonded atoms was probed at every 100 configurations.After an arbitrary dihedral angle rotation the new molecular geometry was transferred to MOPAC and the contribution of this molecule to the intramolecular energy needed in equation 1 was calculated.In the present work, the heat of formation calculated in MOPAC was used to obtain the dependence of the internal molecular energy with the dihedral angle.Calculations with AM1, PM3 and MNDO semiempirical Hamiltonians implemented in MOPAC were performed.Considering the dihedral angle defined by the H-O-C-C atoms, simulations were also performed using rotational energies obtained with HF, MP2 and B3LYP and 6-31g* basis set. 47In this case, for each level of theory the rotational energy barrier was calculated at dihedral angles increments of df = 3 o and the other geometry parameters were not optimized.Each discrete energy curve obtained was fitted to a four-term Fourier expansion of the dihedral angle f using standard mathematical procedures implemented in the xmgrace 25 program.

Monte Carlo simulation protocol
Monte Carlo Metropolis simulations were carried out in the NpT ensemble at 1.0 atm and 298 K. Cubic box boundary conditions and minimum image convention were used.In the calculation of the configurational energy using equation 1 a full intermolecular interaction energy Eab between molecules a and b was considered whenever any of the site-to-site distances r ij fell bellow the cut-off radius r c = 10 Å. Lennard-Jones potential contributions due to interactions beyond the cutoff radius were calculated using the formalism presented by Allen and Tildesley. 4he generalized reaction field approach 48 was used to account for Coulomb interactions beyond the cut-off range.The experimental value of the dielectric constant of ethanol, 24.3, was used.The simulation box was filled with 256 ethanol molecules and the initial box dimensions were calculated to reproduce the experimental density.Starting from an arbitrary distribution of molecules on the simulation box a new configuration was generated by randomly translating and rotating the Cartesian coordinates of a randomly chosen molecule.Rotations along the dihedral angle defined by atoms H-O-C-C were attempted at every 100 configurations.As the calculations were carried out in the NpT ensemble new configurations were also generated by probing the density of the liquid with volume changes.After a volume change all the centre of mass coordinates of all molecules in the reference box were scaled in the usual way. 4 Ranges for monomers translations, rotations and volume changes were adjusted to yield an acceptance/trial ratio of about 0.45 for new configurations.Each calculation started with an equilibration phase of 2 × 10 6 configurations and the Monte Carlo averages were then obtained after a new sampling segment with 2 × 10 7 configurations.Statistical uncertainties were calculated from separate averages over blocks of 2 × 10 5 configurations.To accomplish the present study more than 10 8 energy calculations were performed using MOPAC.All calculations were carried out on a dual processors PC running the FreeBSD 4.2 operational system. 49

Results and Discussion
Dihedral angle and radial distributions obtained in the present investigation are presented bellow.These data are useful to reveal the preferential molecular conformations in the liquid state.Due to the presence of -O-H groups significant contributions from hydrogen bonding interactions are expected.

Dihedral angle distributions
The dihedral angle distributions obtained for ethanol molecules using different methodologies are presented in Figures 1 and 2.
Differences in the dihedral angle distributions presented in Figures 1 and 2 are clear.
Using the ab initio energy barriers, the amplitude of peaks related to the trans population (dihedral f = 180 o ) decreased in the order HF > MP2 > DFT.The dihedral angle distributions obtained with the semiempirical Hamiltonians also show noticeable differences.The MNDO parameterization gives large contribution of trans conformer population, contrasting with the predominance of gauche conformers population obtained with AM1 and PM3 Hamiltonians.The peak positions for gauche and trans conformations obtained by all methodologies but MNDO-OPLS are in reasonable agreement.Nevertheless, the conformer populations are very different indicating disagreements between the methodologies used.OPLS-AA results were compared to B3LYP-OPLS and the agreement is very good.As the OPLS-AA rotational energy barrier was developed fitting to experimental data, 8 one concludes that among the hybrid models the B3LYP-OPLS has a better agreement with experiment.Therefore, the overall analysis of the dihedral angle distributions obtained indicates the need of adjusting the semiempirical parameterizations to reproduce the dependence of the internal energy with the rotational degree of freedom.Further details are provided in the electronic supplementary information section of this article.

Radial distribution functions
The radial distributions functions obtained with different potential models are very similar, as can be seen inspecting the set of O-O radial distributions shown in Figure 3. Therefore, only distributions obtained with AM1-OPLS and MP2-OPLS models are presented in full.The curves exhibit the characteristic features found in hydrogen bonded liquids. 8,28,50(i) the first peak position in the O-O radial distribution is near 2.8 Å; (ii) the first peak position in the H-O radial distribution is near 1.8 Å.The peaks positions are in good agreement with results obtained with the pure OPLS force field published elsewhere. 8The similarities observed in these results indicate a weak influence of the dihedral angle distributions on the population of hydrogen bonded dimers.
To obtain more details about molecular interaction the lower energy dimer found along a 5000 steps MMC   trajectory for each QM-FF model is displayed in Figure 4.One observes in Figure 4 differences in the geometries but the intermolecular interaction energies are very close.In Figure 4 the values obtained for the O --H distances and the H --O-H angles are in the intervals generally accepted to characterize hydrogen bonding: (H --O) distance in the interval ca.1.2-2.2Å and the (H --O-H ) angle in the range ca.130 o -180 o . 50It is worth to note the agreement between the H---O distances in Figure 4 and the peak position in the corresponding radial distribution functions shown in Figure 3.The dimers structures in Figure 4 can be interpreted as a supermolecule with weak energy of torsion along the O-H ---O segment.This rotation like effect contributes to increase the liquid entropy with low influence in the cohesion energy.To further investigate the similarity between dimer interaction, pair energy distributions are shown in Figure 5.As the curves obtained with all QM-FF models are very similar, only two examples are shown.As discussed elsewhere, the bimodal shape of these curves is characteristic of hydrogen bonded liquids. 8,50n Figure 6 the average values of <cos (theta)>, where theta is the angle between the dipole moments of the two ethanol molecules being considered, are presented.The agreement between the three QM-OPLS models in the distance interval [2, 3.7] is clear but differences are found beyond 3.7 Å.There is a remarkable agreement between the positions of features in Figure 6 with the ones in the O-O radial distribution presented on Figure 3. Therefore, the behavior of <cos(theta)> in the distance interval [2,  3.7] can be also explained by the occurrences of hydrogen bonding.As shown in Table 1 there is remarkable agreement between the values of theta obtained for low energy dimers searched along the MMC trajectory.At distances greater than 3.7 Å one observes differences between the QM-OPLS results obtained for <cos(theta)>.Therefore, properties related to dipole moment correlations, such as the dielectric permittivity, DP, are expected to have dependences with the dihedral angle distribution.As discussed elsewhere, the DP is proportional to <M 2 > where M is total dipole moment considering molecules inside a sphere of radius Ro. 51 Therefore, differences are expected on the values of DP calculable with these QM-FF models and, consequently, on the pure liquid and solvation properties which are related to the dielectric permittivity.

Densities and enthalpies of vaporization
Theoretical values obtained for density and enthalpy of vaporization of liquid ethanol at 1.0 atm and 298 K are compared with experimental data in Table 1.One observes in Table 1 a good agreement between theoretical results and experimental values.Therefore, the agreement Table 1.Average dipole-dipole angle, oxygen-oxygen distances and energies obtained for the lowest energy dimers searched along a 5000 steps MC trajectory.The magnitude of this angle is mainly determined by hydrogen bonding interaction instead of dipole-dipole orientation.The oxygen-oxygen distances are in good agreement with the corresponding peak position in the radial distribution functions presented in Figure 3 Model   between the theoretical results obtained for density and enthalpy of vaporization using different QM-FF methods is in contrast with the discrepancies found in the dihedral angle distributions, as shown in Figures 1 and 2. Apparently, the differences between gauche and trans populations obtained using different QM-FF methods are not a problem for obtaining density and enthalpy of vaporization in close agreement with experimental data.This finding is also in agreement with the dimer energies and geometry parameters shown if Figure 4 and also with the pair energy distributions presented in Figure 5.

Conclusions
A Monte Carlo Metropolis simulation program including QM/MM facilities was presented.Following usual procedures the total energy of the system is partitioned into force field and quantum mechanical contributions.Classical contributions are calculable using standard force field and the main MMC program includes MOPAC 6 is as a sub-routine to perform quantum calculations using the facilities implemented in this program.The program was used to study liquid ethanol in the NpT ensemble at 1 atm and 298 K.The OPLS-AA force field was used to obtain intermolecular interactions energies.Internal energy for rotations along de dihedral angle defined by the H-O-C-C atoms were obtained using AM1, PM3 and MNDO parameterizations implemented in MOPAC.Rotational energies calculated with HF, MP2 and B3LYP using 6-31g* basis set were also used.The results obtained for the dihedral angle distributions shows noticeable differences regarding the populations of trans and gauche conformers.These results are further discussed in the electronic supplementary information section of this article.Features observed in the Oxygen-Oxygen and Hydrogen-Oxygen radial distributions functions show strong influence of hydrogen bonded dimmers.These distributions are mainly unaffected by the differences in conformers population.The dipole-dipole correlation, as a function of the distance, shows that the average angle between dipole moments is in the interval 50 o -70 o .This result is very different if compared to the preferential head-to-tail (0 o ) and anti-parallel (180 o ) dipoledipole orientations expected for dipolar liquids.Therefore, the hydrogen bond drives the dipole-dipole orientation to a different pattern when compared to the ones expected for non-hydrogen bonded dipolar molecules.This behavior was also observed for other hydrogen bonded liquids. 31herefore, a competition between configurations favorable to hydrogen bonding and to dipole-dipole correlations is a destabilizing electrostatic force towards dimers with different configurations, increasing the entropy.Pair energy distributions obtained with different QM/MM models are very similar.As a consequence, the values obtained for liquid density and enthalpy of vaporization are also in very good accord.The agreement with experimental and pure OPLS-AA force field calculation is also very good.As concluding remark, the possibility of using MOPAC to calculate the internal energy as well as other quantum mechanical information is a considerable improvement to study complex systems using computer simulations.Other applications are under progress.

Figure 1 .
Figure 1.Dihedral angle distributions for rotation along the dihedral angle defined by the H-O-C-C atoms obtained with HF-OPLS, MP2-OPLS and B3LYP-OPLS models.

Figure 4 .
Figure 4. Intermolecular energies (kcal mol -1 ) and geometries of ethanol dimers.A dimer structure can be interpreted as a supermolecule with weak energy of torsion along the [O-H ---O] segment.This rotation like effect contributes to increase the liquid entropy but keeping the cohesion energy almost unaffected.

Figure 6 .
Figure 6.Average dipole-dipole correlation as a function of oxygenoxygen distance.Theta is the angle between the dipole moments of the two ethanol molecules being considered.

Table 2 .
Results obtained for density and enthalpy of vaporization of liquid ethanol at 298 K and 1 atm