## Journal of the Brazilian Chemical Society

*Print version* ISSN 0103-5053

### J. Braz. Chem. Soc. vol.20 no.8 São Paulo 2009

#### http://dx.doi.org/10.1590/S0103-50532009000800022

**ARTICLE **

**DIADORIM: a Monte Carlo Program for liquid simulations including quantum mechanics and molecular mechanics (QM/MM) facilities: applications to liquid ethanol**

**Luiz Carlos Gomide Freitas ^{*}**

Departamento de Química, Universidade Federal de São Carlos, 13565-905 São Carlos-SP, Brazil

**ABSTRACT**

This paper presents a Metropolis Monte Carlo (MMC) computer program developed to simulate liquids and solutions including QM/MM facilities: the energy from intermolecular interactions is calculated with classical force field functions and the internal molecular energies are calculated using Quantum Chemistry methods. The following facilities were implemented: (*i*) the semiempirical MOPAC 6 quantum chemistry package was included as a subroutine of the main MMC simulation program; (*ii*) alternatively, an interface with an external data bank was developed to allow the use of energies previously obtained; (*iii*) calculation in NpT and NVT ensembles are available; (*iv*) the trajectory generated along the MMC sampling can be saved using standard xtc file format allowing trajectory visualization and data analysis using external programs. The program was used to calculate thermodynamical and structural properties of liquids ethanol (ET) in the NpT ensemble at 1.0 atm and 298 K. Bond angles and bond distances of ethanol molecules were kept constant but torsions along the dihedral angle defined by the H-O-C-C atoms were sampled in the simulation. QM/MM calculations were performed using the MOPAC AM1, PM3 and MNDO Hamiltonians to calculate the internal molecular energy. Calculations using internal energies obtained with OPLS-AA force field, *ab initio* HF/6-31g*, MP2/6-31g* and B3LYP/6-31g* calculations were performed. The dihedral angles distributions obtained with different methodologies reveal extraordinary qualitative and quantitative differences. However, the pair energies and radial distributions functions obtained with the different methodologies are in good agreement. The values calculated for the density and enthalpy of vaporization of liquid ethanol are in good agreement with experimental data.

**Keywords:** Monte Carlo method, QM/MM, computer simulation, liquid ethanol, liquids

**RESUMO**

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,4}As a general trend in such simulations classical potential functions are used to calculate the potential energy.^{5,6} A 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,8 }Nevertheless, 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,10} Nevertheless, 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,12} These methodologies, classified by the general denomination of QM/MM (Quantum Mechanics-Molecular Mechanics) have been progressing and the literature reports are very encouraging.^{13-15 }In some methodologies the quantum-motif 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.^{16} Due to the development of quantum chemical methods, different approaches have been used to describe the quantum-motif.^{17-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,21 }Some 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-30} and solutions^{31-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. To fulfill this need adequately parameterized potential functions have been developed to study molecular liquids and biological molecules.^{5-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.^{6} Therefore, 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.^{6} This 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. Therefore, the development of methodologies to incorporate quantum mechanics corrections in molecular simulation is a subject that has been pursued by many research groups.^{9-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.^{34} Among 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,40} Monte Carlo QM/MM studies using MOPAC have been reported and the results obtained are very encouraging.^{41,42}

In 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.^{44} These 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,46}

Therefore, 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.^{8} Bond 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.^{47} In this case, for each level of theory the rotational energy barrier was calculated at dihedral angles increments of δφ = 3º 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 φ 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.^{4 }The 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 φ = 180º) 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.^{8} The 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º-180º.^{50} It 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,50}

In 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 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º-70º. This result is very different if compared to the preferential head-to-tail (0º) and anti-parallel (180º) dipole-dipole 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.^{31} Therefore, 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.

**Acknowledgments**

This work was supported by FAPESP, Fundação de Amparo à Pesquisa do Estado de São Paulo, under Proc. 07/54381-4. The author is grateful to his former students and collaborators for support and encouragement.

**Supplementary Information**

Supplementary data are available free of charge at http://jbcs.sbq.org.br, as PDF file.

**References**

1. Heermann, D. H.; *Computer Simulation Methods in Theoretical Physics*; Spring-Verlag: London, 1986. [ Links ]

2. Rappaport, D. C.; *The Art of Molecular Dynamics Simulation*, Cambridge Press: Cambridge, 1995. [ Links ]

3. Frenkel, D.; Smit, B.; *Understanding Molecular Simulation*, Academic Press: San Diego, 2002. [ Links ]

4. Allen, T. D.; Tildesley, D. J.; *Computer Simulation of Liquids*, Oxford University Press: Oxford, 2002. [ Links ]

5. Dykstra, C. E.; *Chem. Rev*. **1993**, *93*, 2339. [ Links ]

6. Rappé, A. K.; Casewit, C. J.; *Molecular Mechanics Across Chemistry*, University Science Books: Sausalito, 1997. [ Links ]

7. Jorgensen, W. L. In *Encyclopedia of Computational Chemistry*; Schleyer, P. v. R., ed., Wiley: 1998, 5, 3281-3285. [ Links ]

8. Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rieves, J.; *J. Am. Chem. Soc*. **1996**, *118*, 11225. [ Links ]

9. Vreven, T.; Byun, K. S.; Komaromi, I., Dapprich, S.; Montgomery, J. A., Jr.; Morokuma, K.; Frisch, M. J.; *J. Chem. Theory Comput. ***2006**, *2*, 815. [ Links ]

10. Tse, T. J.; *Annu. Rev. Phys. Chem*. **2002**, *53*, 249. [ Links ]

11. Aqvist, J.; Warshel, A.; *Chem. Rev*. **1993**, *93*, 2523. [ Links ]

12. Warshel, A.; *Annu. Rev. Biophys. Biomol. Struct.* **2003**, *32*, 425 [ Links ]

13. Kapral, R.; *Annu. Rev. Phys. Chem*. **2006**, *57*, 129. [ Links ]

14. Bruice, T. C.; *Chem. Rev*. **2006**, *106*, 3119. [ Links ]

15. Prado, C. E. R.; Pópolo, M. G. Del; Youngs, T. G. A.; Kohanoff, J.; Lynden-Bell, R. M.; *Mol. Phys*. **2006**, *104*, 2477. [ Links ]

16. Fileti, E. E.; Georg, H. C.; Coutinho, K.; Canuto, S.; *J. Braz. Chem. Soc*. **2007**, *18*, 74. [ Links ]

17. Kaminski, G. A.; Jorgensen, W. L.; *J. Phys. Chem. B* **1998**, *102*, 1787. [ Links ]

18. Biswas, P. K.; Gogonea, V.; *J. Chem. Phys*. **2005**, *123*, 1. [ Links ]

19. Tishchenko, O.; Truhlar, D. G.; *J. Phys. Chem. A* **2006**, *110*, 13530; [ Links ] *J. Chem. Theory Comput. ***2007**, *3*, 938.

20. Coolidge, M. B.; Stewart, J. J. P.;* MOPAC Manual*, Frank J. Seiler Res. Lab., United States Air Force Academy, CO 80840 1990. [ Links ]

21. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E.; *J. Chem. Phys*. **1953**, *21*, 1087. [ Links ]

22. Onsager, L.; *J. Am. Chem. Soc*. **1936**, *58*, 1486; [ Links ] Marenich, A. V.; Olson, R. M.; Kelly, C. P.; Cramer, C. J.; Truhlar, D. G.; *J. Chem. Theory Comput*. **2007**, *3*, 2011. [ Links ]

23. Zwanzig, R. W.; *J. Phys. Chem*. **1954**, *22*, 1420; [ Links ] Bennett, C. A.; *J. Comput. Phys*. **1976**, *22*, 245. [ Links ]

24. *XTC* Format, http://www.gromacs.org/documentation/reference/online/xtc.html, acessed in January 2009. [ Links ]

25. *Visual Molecular Dynamics, VMD*, http://www.ks.uiuc.edu/Research/vmd/, accessed in January 2009. [ Links ]

26. Schaftenaar, G.; *MOLDEN Program*; ttp://www.cmbi.ru.nl/molden/molden.html, accessed in January 2009. [ Links ]

27. *Xmgrace Program*, http://plasma-gate.weizmann.ac.il/Grace/, acessed in January 2009. [ Links ]

28. Freitas, L.C.G.; *THEOCHEM* **1993**, *101*, 153. [ Links ]

29. Silva, L. B.; Freitas, L. C. G.; *THEOCHEM* **2007**, *806*, 23. [ Links ]

30. Bertran, C. A.; Cirino, J. J. V.; Freitas, L. C. G.; *J. Braz. Chem.Soc*. **2002**, *13*, 238. [ Links ]

31. Sinoti, A. L. L.; Politi J. R. S.; Freitas, L. C. G.; *J. Braz. Chem. Soc*. **1996**, *7*, 133. [ Links ]

32. Freitas, L. C. G.; Cordeiro J. M. M.; Garbujo F. L. L.; *J. Mol. Liq*. **1999**, *79*, 1. [ Links ]

33. Freitas, L. C. G.; Cordeiro, J. M. M.; *THEOCHEM* **1995**, *335*,189. [ Links ]

34. Stewart, J. J. P. In *Encyclopedia of Computational Chemistry*; Schleyer, P. v. R., ed., Wiley: 1998, 4, 2574-2578. [ Links ]

35. Stavrev, K. K.; Zerner, M. C.; *Chem. Phys. Lett. ***1996**, *253*, 667. [ Links ]

36. Canuto, S.; Coutinho, K.; *Adv. Quantum Chem*. **2002**, *41*, 161. [ Links ]

37. Dewer, M. J. S.; Stewart, J. J. P.; *AMPAC, Quantum Chemistry Program Exchange Bulletin* **1986**, *6*, 24. [ Links ]

38. Thompson, J. D.; Cramer, C. J.; Truhlar, D. G.; *J. Phys. Chem. A* **2004**, *108*, 6532. [ Links ]

39. Chandrasekhar, J.; Jorgensen, W. L.; *J. Comput. Chem.* **2002**, *23*, 1601 [ Links ]

40. Rocha, G. B.; Freire, R. O.; Simas, A. M.; Stewart, J. J. P.; *J. Comput. Chem*. **2006**,* 27*, 1101. [ Links ]

41. Gao, J. L.; Freindorf, M.; *J. Phys. Chem. A* **1997**, *101*, 3182. [ Links ]

42. Chandrasekhar, J.; Shariffskul, S.; Jorgensen, W. L.; *J. Phys. Chem. B* **2002**, *106*, 8078 [ Links ]

43. Dolney, D. M.; Hawkins, G. D.; Winget, P.; Liotard, D. A.; Cramer, C. J.; Truhlar, D. G.; *J. Comput. Chem*. **2000**, *21*, 340. [ Links ]

44. Storer, J. W.; Giesen, D. J.; Cramer, C. J.; Truhlar, D. G.; *J. Comput.-Aided Mol. Des. ***1995**, *9*, 87. [ Links ]

45. Chandrasekhar, J.; Shariffskul, S.; Jorgensen, W. L.; *J. Phys. Chem. B* **2002**, *106*, 8078 [ Links ]

46. Acevedo, O.; Jorgensen, W. L.; *J. Chem. Theory Comput.* **2007**, *3*, 1412. [ Links ]

47. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Peng, C. Y.; Ayala, P. A.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; Head-Gordon, M.; Gonzalez, C.; Pople, J. A.; *Gaussian 94 (Revision A.1)*, Gaussian, Inc.: Pittsburgh, PA, 1995. [ Links ]

48. Tironi, I. G.; Smith, P. E.; van Gunsteren, W. F.; *J. Chem. Phys*. **1995**, *102*, 5451. [ Links ]

49. FreeBSD Operational System, http://www.freebsd.org, accessed in January 2009. [ Links ]

50. Jeffrey, G. A.; *An Introduction to Hydrogen Bonding*, Oxford University Press: Oxford, England, 1997. [ Links ]

51. Neumann, M.; *J. Chem. Phys.***1985**, *82*, 5663. [ Links ]

52. Marcus, Y.; *Ion Solvation*, John Wiley: New York, 1985. [ Links ]

53. Riddick, J. A.; Bunger, W. B.; *Organic Solvents*, Wiley-Interscience: New York, 1970. [ Links ]

Received: September 10, 2008

Web Release Date: August 31, 2009

**FAPESP helped in meeting the publication costs of this article.**

* e-mail: gomide@dq.ufscar.br

**Supplementary Information**

*Comments*

The rotational energies profiles obtained with the semiempirical methods are qualitatively different from the ones calculated with Hartree-Fock, second order Moller-Plesset (MP2) theory and B3LYP functional using 6-31g* basis set. The qualitative disagreement between the semiempirical and the others methodologies are clear. These differences are mainly related to the small values of the barriers to internal rotation as well as to the inadequacy of semiempirical methods in reproducing conformational analysis.^{34,37} For each rotational energy profile the corresponding dihedral angle distribution obtained in the Monte Carlo simulation is also shown. As expected from the Boltzmann distribution, minima in the potential energy curve are correlated with maxima in the corresponding distribution. This correspondence explains the predominance of *gauche *populations obtained with AM1-OPLS and PM3-OPLS models as well as of the *trans* population observed in the MNDO-OPLS results.