Article pubs.acs.org/jchemeduc
Computational Study of Environmental Effects on Torsional Free Energy Surface of N‑Acetyl‑N′‑methyl‑L‑alanylamide Dipeptide Silvia Carlotto* and Mirco Zerbetto Dipartimento di Scienze Chimiche, Università degli Studi di Padova, via Marzolo 1, 35131, Padova, Italy S Supporting Information *
ABSTRACT: We propose an articulated computational experiment in which both quantum mechanics (QM) and molecular mechanics (MM) methods are employed to investigate environment effects on the free energy surface for the backbone dihedral angles rotation of the small dipeptide N-acetyl-N′-methyl-L-alanylamide. This computation exercise is appropriate for an upper-level physical chemistry course. The purposes are (i) to show the importance of solvent effects on the free energy surface of molecules in solution and (ii) to compare different computational chemistry tools in terms of computation complexity versus detail of information obtained from calculus. The QM section is divided in three parts: in vacuo, implicit solvent, and explicit solvent calculations. Similarly, also the MM section is divided in three parts: molecular mechanics, molecular dynamics, and adaptive biasing force molecular dynamics. QM and MM sections can be proposed by the instructor as different and complementary parts of the experiment or as independent QM or MM class lessons. Together, the two sections can be used to set up a computational chemistry laboratory exercise targeted to the study of solvent effects on molecular properties. KEYWORDS: Upper-Division Undergraduate, Physical Chemistry, Computer-Based Learning, Computational Chemistry, Molecular Mechanics/Dynamics, Molecular Modeling, Quantum Chemistry
E
properties using as a case study the short peptide N-acetyl-N′methyl-L-alanylamide (NANMA). In particular, the solvent effects on the backbone conformational energetics of NANMA will be explored. This dipeptide is the ideal prototype system for computational studies of biomolecules because its reduced dimensions allow us to conduct (in a reasonable time) a full computational study. The relevant coordinates are clearly identified as the two backbone dihedral angles (see Figure 1).
nvironmental effects can play an important and not negligible role on the free energy of solute molecules in the condensed phase. Both the internal energy contribution, originating from specific and nonspecific interactions with solvent molecules, and the entropic contribution, arising from statistical averaging over solvent configurations, should generally be taken into account in the evaluation of the free energy surface in a system. Important variations in both the physical and chemical properties of a molecule can reflect changes of the energy surface caused by the presence of the solvent. For these reasons, the inclusion of environmental effects in the calculation of molecular properties is mandatory for the accurate comparison of the results of computations with experimental outcomes. Recently, reports of a number of numerical experiments have been published in this Journal, underscoring the importance of calculus in predicting the properties of chemical systems (from single molecule properties, to bulk thermodynamics). A review of these studies reveals that the solvent is usually introduced in the calculations.1−3 Several important considerations arise when determining how the solvent should be accounted for in numerical calculations. One important decision that must be made is the level of description that should be used for the solvent in a numerical experiment. Here we present a comprehensive computational study aimed at comparing different levels of theory of solvent description in computations, suggesting a systematic computational study of solvent effects on molecular © 2013 American Chemical Society and Division of Chemical Education, Inc.
Figure 1. Schematic representation of N-acetyl-N′-methyl-L-alanylamide (NANMA). The two backbone torsion angles φ and ψ are highlighted. Numeric labels (in blue) have been assigned to heavy atoms and hydrogen atoms not belonging to methyl groups.
The potential energy surface (PES) profile, which is a function of the two angles, is influenced by both nonspecific and specific interactions with the solvent. In vacuo, the formation of an intramolecular hydrogen bond is observed,4 which leads to an important stabilization of the corresponding backbone configuration. In contrast, it is known in polar Published: November 22, 2013 96
dx.doi.org/10.1021/ed400346v | J. Chem. Educ. 2014, 91, 96−102
Journal of Chemical Education
Article
solvents that the specific interaction with solvent molecules may prevent the formation of the intramolecular hydrogen bond. The computational experiment is articulated in the application of both quantum mechanical (QM) and molecular mechanical (MM) methods to investigate the environmental effects on backbone energetics. The QM section is divided into three parts: in vacuo, implicit solvent, and explicit solvent calculations. Similarly, the MM section is divided into three parts: molecular mechanics, molecular dynamics (MD), and adaptive biasing force (ABF) molecular dynamics calculations. The QM and MM sections can be proposed by the instructor as different and complementary parts of the experiment, or as independent QM or MM class lessons. The two sections can be jointly used to set up a computational chemistry laboratory exercise targeted to the study of solvent effects on molecular properties. The rationale of the full computational experiment is presented in Figure 2.
Packages section) and carry out QM and MM calculations. Through this work, students also begin to understand that different computational approaches allow the derivation or extraction of different kinds of information, leading to an awareness that is important when determining the best protocol for approaching a problem in chemistry by means of computation. The computation exercise is appropriate for an upper-level physical chemistry course. To supplement the text, we provide tentative estimations of the teaching time for all portions of the experiment. Note that these times include the preparatory part of the computations (discussion of the objectives of the experiment, explanation of the basics of the computational methods, preparation of the input files) and the discussion of the results. We do not mention pure computation times because these times are strongly dependent on the available hardware for computation.
■
SOFTWARE PACKAGES The full set of software packages used in conducting the proposed experiment includes Molden 4.05 and VMD 1.96 to model the system, Gaussian 097 for the QM calculations, NAMD8 for the MM and MD calculations, and Octave9 for the analysis. Clearly, any other software packages performing the same tasks may be used, based on their availability in the laboratory. Here, we assume that the instructor is familiar with the programs. However, to simplify the setup of the experiment, we provide, as Supporting Information to this paper, the full set of input files for the QM and MM portions of the computational exercise.
■
QUANTUM MECHANICAL CALCULATIONS
System Modeling
In QM calculations, the first requirement is to specify a set of Cartesian or internal coordinates (Z matrix) for the atoms of the NANMA molecule. To this end, we suggest the use of Molden 4.0, which provides students with the capability to build many types of systems from their inception and to learn the use of template databases for the rapid building of peptides. A set of non-optimized, internal coordinates of NANMA produced by Molden is available in the Supporting Information. We chose the internal coordinates instead of Cartesian coordinates because a number of conformations using different values for the backbone torsion angles should be produced. This operation is easier with internal coordinates; for other QM
Figure 2. The rationale of the full computational experiment, divided into six different calculations, each with a different method (as described in the main text). The top half of the schematic summarizes the three quantum mechanical approaches (A, B, C), whereas the bottom half reports the classical approaches: (D) molecular mechanics, MM; (E) molecular dynamics, MD; (F) adaptive biasing force, ABF, molecular dynamics.
At the end of this experiment, the students learn how to use a set of computational chemistry tools that are frequently employed in the modeling of molecular systems (see Software
Figure 3. Backbone angles potential energy surfaces (PES) obtained from quantum mechanical calculations (A) in vacuo, and (B) in explicit solvent. Energy units are kJ/mol. 97
dx.doi.org/10.1021/ed400346v | J. Chem. Educ. 2014, 91, 96−102
Journal of Chemical Education
Article
software packages, the choice could be different. Molden permits the user to set and save the coordinates of the system for common QM software packages (eg, Gaussian 09). These coordinates will be used in the input files for QM in vacuo and for implicit solvent calculations. The estimated time for comprehension of the QM exercise, use of Molden, and student work preparing the Z matrix is approximately 1 h.
molecules to account for the most important specific solute− solvent interactions. We suggest repeating the PES calculations for the water solvent using two different solvation models. In the first model, the continuum solvent model is employed to account for bulk effects; in the second model, an explicit water molecule is introduced to account for the specific solute−solvent interaction.
QM in Vacuo Calculations (Part A)
Implicit Solvent QM Calculations (Part B)
Before starting the scan over torsion angles, the initial geometry of NANMA should be fully optimized to reduce the computational time required to calculate the relaxed potential energy surface. We suggest optimization in vacuo at an HF/631G level of theory using the Gaussian 09 software package. A sample input file can be found in the Supporting Information. This input file is built by merging the non-optimized coordinates obtained by Molden software and the Gaussian optimization directive: #HF/6-31G Opt. For simplicity we use the standard optimization algorithm of the Gaussian 09 software package. After geometry optimization, the PES can be determined as a function of the two dihedral angles φ and ψ. A two-dimensional scan featuring 30 deg steps for both angles is employed, producing a grid of 169 conformations. For each conformation, a constrained geometry optimization at the HF/6-31G level of theory is performed by keeping fixed the values of the two dihedral angles. The Gaussian input file for the scan as a function of the two dihedral angles, together with the Octave script to plot the surface, are provided in the Supporting Information. The resulting PES is shown in Figure 3A. As we are comparing the depth of the extrema of the energy surfaces, all of them are shifted such that the energy minimum is located at 0 kJ/mol. Students are encouraged to note that a number of energy minima are observed (dark blue areas in the figure), and that these minima are compatible with those found using other approaches.4,10 In particular, the most stable regions are found around the six couples of (φ,ψ) values: (−90, −60), (−140, 160) (as found by Drozdov et al.11), (70, −60), (80, 60), (−145, −145), and (−180, 180) deg. By analyzing the molecular geometry corresponding to the different minima, an intramolecular hydrogen bond is found in the (70, −60) deg configuration with a 1.91 Å distance between NANMA O3 and H12 atoms. The estimated time for both input preparation and output analysis is approximately 1 h.
The relaxed PES calculation in implicit solvent using the polarizable continuum model (PCM)12 can be performed in Gaussian 09 using a slightly modified command line with respect to that used in part A, that is, #HF/6-31G Opt=ZMatrix SCRF=(pcm, solvent=water). The first keyword in the SCRF options specifies the implicit solvent model to be used (here PCM), whereas the second one specifies the solvent (and thus its dielectric constant and radius). In this procedure, the solute molecule is embedded into a cavity opened in the dielectric medium (water). A standard cavity is used, and a relative dielectric constant of water, at room temperature, of 78.35 is employed. As in the case of the in vacuo calculations, the two angles are varied by steps of 30 deg. Because the PCM, as with the other implicit solvent approaches, does not account for specific solute−solvent interactions, the resulting PES (data not shown) simply shows a general global stabilization of the system (energy shift) with a slightly more pronounced energy stabilization of conformations with higher dipole moment. This is an expected behavior due to the high value of the dielectric constant of water. For both input preparation and output analysis, the estimated time is 1 h. Explicit Solvent QM Calculations (Part C)
For our explicit solvent QM calculation, we introduced a water molecule to the system to show that the situation (with respect to calculation) radically changes even if only one explicit water molecule is introduced. We suggest using the Molden software and positioning the molecule in such a way that one of its hydrogen atoms is near the NANMA O3 atom. For the sake of clarity, we show in Figure 4 how we chose the initial position of
QM Solvent Calculations
The next step is to understand solvent effects on the obtained PES. Usually, in QM calculations, solvent effects are taken into account by means of two different approaches: the implicit solvent approach (or continuum) and the explicit solvent approach (or molecular) approach. In the continuum model, a homogeneous medium is employed to represent the solvent, which influences the molecular properties in an averaged way. Solvents are identified by dielectric constants and solvent radii. The molecular model includes an atomistic representation of the solvent molecules, permitting the introduction of specific interactions between solute and solvent. Clearly, this second model provides a more realistic description of the system, but poses limitations due to the large number of solvent molecules that would be required to reproduce the bulk properties. Thus, the typical practice is to introduce a small number of solvent
Figure 4. Non-optimized structure of NANMA with one explicit water molecule used as initial configuration in the explicit-solvent QM calculations.
the water molecule with respect to NANMA before the geometry optimization. Among other initial configurations, the one reported in Figure 4 has been chosen as it leads to a convergence of the optimization algorithm in which the solute−solvent interaction is taken into consideration. Attention should be paid in preparing the initial configuration 98
dx.doi.org/10.1021/ed400346v | J. Chem. Educ. 2014, 91, 96−102
Journal of Chemical Education
Article
Figure 5. Backbone angles potential energy surfaces (PES) obtained from (A) classical MM calculations in vacuo and (B) ABF calculations in (explicit) water. Energy units are kJ/mol.
box is sufficient. To solve the problem, it is possible to invoke the time-scale separation between the dynamics of electrons and those of nuclei, which leads to the molecular mechanics (MM) approach.
because placement of the water molecule too far from the peptide runs the risk of describing a noninteracting system. Calculations can be performed at the HF/6-31G level, by varying the two backbone angles by steps of 30 degrees (the Gaussian command line is HF/6-31G Opt=Z-Matrix). In the presence of one explicit water molecule, we noticed some instability of the Gaussian program in performing the 2D scan. If this happens, we suggest performing different runs in which the angle φ is kept fixed to a number of increasing values and a 1D scan along ψ is performed. The resulting energy surface is shown in Figure 3B. An important change in the minima is observed with respect to the in vacuo calculations (Figure 3A). It is also noteworthy that the presence of just one explicit water molecule greatly increases the barriers among the minima (up to 4−5 times) in comparison to the barriers in vacuo. This leads to higher stability in the minimum energy configurations in solvent with respect to those in vacuo. The deeper minima are located around (−180, −120), (−180, 160), (90, −90), and (90, 90) deg, whereas the maximum is in (0, 0) deg. The in vacuo calculations, instead, show the maximum around (0, ±180) deg. One of the most evident changes is the larger energy stabilization around (90, −90) deg, which is observed in the presence of the explicit water molecule. Following our suggestion to split the PES 2D scan into 13 different 1D scans, the estimated time for the preparation of input files and output analysis is 2 h.
System Modeling
MM and MD simulations can be carried out with NAMD, which requires both a protein data bank (PDB) file and a protein structure file (PSF). The Molden and VMD software packages can be used to produce the required files. Since the preparation of these files can be cumbersome for students performing MM and MD calculations for the first time, we provide the full set of ready-to-use files for NAMD in the Supporting Information. For the in vacuo MM calculations, a total of 169 PDB files with different values for the (φ,ψ) angles (discretized by steps of 30 deg) were prepared starting from the Z Matrix created in part A of the QM calculations. Only a single PSF file is required as the topology of the molecule remains the same; it can be generated using the psfgen utility of NAMD. The psfgen script is reported in the Supporting Information. For the standard and biased MD simulations, a water box can be added to the system using the solvation tool of VMD. In this particular case, we use a cubic box with a side length of 24 Å. The PDB and PSF files for the solvated system are available in the Supporting Information. The preparation of the PDB, PSF, and configuration files for the MM and MD calculations may require approximately 2 h if students are guided by an instructor. However, we suggest providing the needed files for the students; this reduces the time needed to about 30 min, during which the instructor would explain the concepts of force fields in MD calculations and of molecular topology. This process allows students to remain focused on the objectives of the experiment without getting lost into the technicalities of the specific software package (NAMD in this case), which are available in the software user guide for those who take an interest.
■
MOLECULAR MECHANICS AND MOLECULAR DYNAMICS CALCULATIONS To gather information on the energetics of the torsion around the two φ and ψ angles one should in general consider both the specific interactions with the solvent (contribution to internal energy) and the statistical distribution of the solvent (contribution to entropy). The former contribution has been taken into account in the QM calculations described above, in which one explicit solvent molecule has been added to the system. The second effect is of statistical nature and requires the averaging over all possible configurations of the solvent molecules. To this end, one needs to sample the statistical distribution of the solvent. In the ergodic hypothesis, this can be done by performing a molecular dynamics simulation. Because a large number of molecules are required, it is still impossible to perform this kind of calculation at the QM level, even for molecules such as NANMA for which a small solvation
MM Calculations (Part D)
Before dealing with the calculations of the NANMA molecule in solvent, we suggest performing an MM calculation of the backbone potential energy surface in vacuo. The main purpose of this calculation is to show students that when reproducing a previous observation using coarse-grained models (if it is possible to define the MM as a coarse graining of QM), it is important to verify (at least qualitatively) that the picture is the 99
dx.doi.org/10.1021/ed400346v | J. Chem. Educ. 2014, 91, 96−102
Journal of Chemical Education
Article
same. Also the qualitative results should not be too different from those of the previous calculation. Students may observe that the available force field (FF) is able to satisfactorily reproduce the position and depth of minima of the PES as obtained in the analogous QM calculations in vacuo. This reproducibility ensures that the two simulations are compatible and that the QM and MM calculations can also be compared quantitatively in the presence of a water solvent. A simultaneous scan of the two torsion angles, as described in the QM calculations, that evaluates the total energy of the molecule for each conformation at the MM level can be conducted using the NAMD software package. The CHARMM27 force field with CMAP correction is employed in this experiment. CMAP introduces energy terms coupling adjacent (φ,ψ) couples to better describe backbone energy in peptides. Other parameters of the simulation are a 12 Å cutoff for electrostatic interactions (with function switching at 10 Å); a pair list distance of 13.5 Å (pair list update every 10 steps); scaling of 1−4 interactions set to 1.0; and particle mesh Ewald (PME) treatment of electrostatics, which provides a fast method to calculate long-range (Coulomb) interactions. The full set of the suggested parameters is considered a standard choice in MD calculations of peptides in explicit water, as it assures a good balance between accuracy of the simulation and computation time. The NAMD configuration file is provided in the Supporting Information. The resulting PES is shown in Figure 5A. The PES surface compares satisfactorily with Figure 3A (PES in vacuo at QM level). Thus, quantitative comparison of the effects of the environment on the backbone free energy of NANMA calculated at QM or MM level is possible. Energy calculations can be performed using the files and scripts provided in the Supporting Information. They take a fraction of minute to run, whereas about 1 h is suggested for the analysis and discussion of the results.
taking the natural logarithm of the equilibrium (Boltzmann) probability density Peq(φ , ψ ) = exp[− V (φ , ψ )/kBT ]/⟨exp[− V (φ , ψ )/kBT ]⟩ (1)
where T is the absolute temperature, kB is the Boltzmann constant, and ⟨...⟩ means averaging over torsion angles phase space. The probability density can be obtained from a histogram analysis of the (φ,ψ) time series. In the Supporting Information we provide the required NAMD configuration file and the NAMD colvar module input file, which handles the calculation of the probability density. The obtained PES is shown in Figure 6. It is worthy to note that the 10 ns simulation time is insufficient to allow the system
Figure 6. Backbone angles potential energy surface (PES) obtained from unbiased molecular dynamics simulations. Colored in gray is the unexplored area. Energy units are kJ/mol.
to explore the complete surface (gray area in the figure). This is due to the very large maxima present in the PES; the unbiased trajectory is only able to sample the two most stable minima. The presence of high energy barriers may prevent the full sampling of the phase-space, irrespective of the trajectory length. This is due to the fact that large changes of energy are not allowed in MD simulations. The possibility of being able to build only an incomplete PES from unbiased MD calculations represents a failure of the experiment, as it is impossible to recover the same kind of information previously obtained by QM calculations. Thus, a better approach should be followed. In particular, biased MD simulations should be used. These simulations should involve coordinates that are chosen to explore the conformational phase space; this is the aim of the latter portion of the MM and MD part of the experiment. The full MD calculation is divided into four parts: energy minimization, heating, equilibration, and production. The proper NAMD configuration files are provided in the Supporting Information. The preparation of the four NAMD input files requires approximately 2 h, including explanation on simulation parameters. However, we suggest providing all needed input files to the students. Assuming this procedure is followed, 30 min would be required for submission of jobs. Finally, we estimate 30 min for analysis of the results.
MD Simulations (Part E)
In the second part of the classical mechanics study, a tentative determination of the PES is carried out by means of a standard unbiased MD simulation using NAMD. As before, the CHARMM27 force field with CMAP correction is employed for NANMA and the water molecules are modeled with the TIP3P force field. TIP3P is a rigid water model in which the Lennard−Jones site is located on the oxygen atom while bare charge sites are placed on the hydrogen atoms.13 As anticipated in the system modeling section, the molecule is immersed in a cubic box with a side length of 24 Å, together with a total of 401 water molecules. A 10 ns trajectory in canonical (constant number of particles, pressure and temperatureNPT) conditions can be easily produced with 2 fs of integration step, saving coordinates every 50 fs. Simulation parameters are 12 Å cutoff for electrostatic interactions (with function switching at 10 Å); pair list distance of 13.5 Å (pair list update every 10 steps); scaling of 1−4 interactions set to 1.0; and particle mesh Ewald (PME) treatment of electrostatics. A Langevin thermostat at 300 K and a Langevin piston (target pressure 1 atm, piston period 200 fs, piston decay 100 fs) are coupled to the system to ensure maintenance of the canonical ensemble properties. Also, the SHAKE algorithm is used to fix all stretches of bonds with hydrogen atoms. We use the colvar module of NAMD to automatically extract the time series of the two backbone angles and produce the potential energy surface V(φ,ψ). The latter is obtained by
ABF Simulations (Part F)
To overcome problems related to passing barriers, one possible approach is to use biased MD simulations. We applied a bias (potential or force) to the coordinates for which the PES has to 100
dx.doi.org/10.1021/ed400346v | J. Chem. Educ. 2014, 91, 96−102
Journal of Chemical Education
Article
solvent) calculations. In particular, a difference that is worthy to note is that in vacuo, and in implicit solvent (i.e., when no specific interactions are explicitly taken into account), the preferred configurations of the peptide are those in which the formation of an intramolecular hydrogen bond is possible, while with explicit water (that is, the explicit solvent) other minima are found. From the classical point of view, a student can, in the first step of this experiment, validate the classical force field (MM calculations; part D of the experiment) before running a standard MD simulation as the second step. Since the latter approach is destined to fail, the student faces the necessity of looking for more sophisticated approaches, such as ABF, to efficiently investigate the entire phase-space of molecular conformations, which in general cannot be sampled using standard unbiased MD calculations. Students may also learn from this observation that, if gaining information on the complete PES for a “simple” molecule as NANMA is so difficult, the sampling of the potential surface of large systems (e.g., proteins) is a very delicate and complex process during which much attention should be devoted to the avoidance of trivial errors. The full experiment described herein should take a maximum of approximately 10 h, including learning and discussion of results. Computation times for QM, MD, and ABF calculations should be added in the estimation of the effective duration of the experiment (depending on available hardware). Thus, the full experiment shall be split into at least 4 or 5 lessons (2 h per lesson). Another option is to divide the students in two groups (one dedicated to the QM exercise and the other one to the MM and MD calculations) or in six groups (each carrying out one single part of the experiment). At the end of the simulations, the groups shall present their results to each other, and the instructor can lead the comparative discussion.
be constructed. For this experiment, we decided to employ the adaptive biasing force method,14 in which a biasing force is applied to the selected reaction coordinates. The strength of the biasing force was adapted to keep a zero average force on the reaction coordinates over time. Due to this adaptation, the profile of the biasing force along the reaction coordinates (the torsion angles, in our case) corresponds to the opposite of the gradient of the free energy surface. The latter is recovered by means of numerical integration. Instructors can follow different approaches, if preferred. In the Supporting Information we provide the input files for the ABF simulations. Parameters of the simulation are the same as before with the inclusion of the following biasing protocol: four independent runs were conducted to explore the four quadrants of the Ramachandran plot. The colvar module of NAMD is used to define both the two backbone angles to be driven and the kind of bias. ABF is simultaneously applied to the two torsion angles with the f ullSamples parameter set to 100, a bin width of 5.0 deg and lower, and upper wall constants set to 0.2 kcal/mol. The resulting PES is shown in Figure 5B, which features minima and maxima comparable to those found with similar approaches.15,16 In this part, the only input file to be modified is the colvar.inp (see the Supporting Information). If part E has just been carried out, it is possible to restart from the equilibration to produce the biased trajectory. Job preparation and submission would require approximately 45 min. Otherwise, if the instructor decides to propose only part F of the experiment, the complete protocol (energy minimization, heating, equilibration, and production trajectory) should be followed (2 h). For the analysis and discussion of results, we estimate 1 h. The inclusion of the solvent produces the major effects that are also shown in the explicit solvent QM calculations. The part of the surface corresponding to negative φ values remains practically unchanged, the maximum of the PES moves from (0, ± 180) to (0, 0) deg, and a minimum energy area appears around (90, −90) deg. This last change represents the most relevant solvent effect on the PES. Some differences are, however, observed between the two surfaces. The most important discrepancy is found around the (90, 90) deg configuration, which is a minimum in the QM PES, but not in the MD/ABF surface. Inspecting the output files, it is found that the energy of the system at (90, 90) deg is 72 kJ/mol from QM, while it is 62 kJ/mol from MD/ABF (nearly 4 kBT smaller than QM PES, at 300 K). From this observation, the hypothesis that can be drawn is that the barriers in the MD/ABF surface are much lower than those in the QM calculation due to entropic effects, which are taken into account in the MD/ABF surface by averaging over all possible solvent molecular configurations. In the explicit solvent QM calculations, in contrast, only one water molecule is present; also, no averaging is carried out over the water molecule rototranslational degrees of freedom.
■
ASSOCIATED CONTENT
S Supporting Information *
A complete ready-to-use archive containing all the input files. This material is available via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS All the simulations for the setup and test of the experiment have been run on the LICC (“Laboratorio Interdipartimentale di Chimica Computazionale”) HPC facility of the Dipartimento di Scienze Chimiche of the Università degli Studi di Padova. This work is supported by Ministero dell’Istruzione, Università e Ricerca (MIUR), grant FIRB 2012.
■
■
DISCUSSION AND CONCLUSIONS In this article, we present a general overview on the employment of quantum mechanical and classical methods to take into account solvent effects on molecular properties. From a quantum mechanical point of view, a student can verify the influence of the solvent on the stabilization of the conformations. It should be observed that the presence of just one water molecule (part C) stabilizes different configurations with respect to the in vacuo (and implicit
REFERENCES
(1) Simpson, S.; Autschbach, J.; Zurek, E. Computational modeling of the optical rotation of amino acids: An ‘in silico’ experiment for physical chemistry. J. Chem. Educ. 2013, 90, 656−660. (2) Simeon, T.; Aikens, C. M.; Tejerina, B.; Schatz, G. Northwestern University initiative for teaching nanosciences (NUITNS): An approach for teaching computational chemistry to engineering undergraduate students. J. Chem. Educ. 2011, 88, 1079−1084.
101
dx.doi.org/10.1021/ed400346v | J. Chem. Educ. 2014, 91, 96−102
Journal of Chemical Education
Article
(3) Box, V. G. S. Using molecular modeling to understand some of the more subtle aspects of aromaticity and antiaromaticity. J. Chem. Educ. 2011, 88, 898−906. (4) Tobias, D. J.; Brooks, C. L. Conformational equilibrium in the alanine dipeptide in the gas phase and aqueous solution: A comparison of theoretical results. J. Phys. Chem. 1992, 96, 3864−3870. (5) Schaftenaar, G.; Noordik, J. H. Molden: A pre- and postprocessing program for molecular and electronic structures. J. Comput.Aided Mol. Des. 2000, 14, 123−134; http://www.cmbi.ru.nl/molden/ (accessed Nov 2013). (6) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38; http://www.ks.uiuc.edu/ Research/vmd/vmd-1.9/ (accessed Nov 2013). (7) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision A.1; Gaussian, Inc.: Wallingford, CT, 2009. (8) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (9) Gnu Octave Web site: http://www.gnu.org/software/octave/ (accessed Nov 2013). (10) Jang, H.; Woolf, T. B. Multiple pathways in conformational transitions of the alanine dipeptide: An application of dynamic importance sampling. J. Comput. Chem. 2006, 27, 1136−1141. (11) Drozdov, A. N.; Grossfield, A.; Pappu, R. V. Role of solvent in determining conformational preferences of alanine dipeptide in water. J. Am. Chem. Soc. 2004, 126, 2574−2581. (12) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999−3093. (13) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (14) Darve, E.; Pohorille, A. Calculating free energies using average force. J. Chem. Phys. 2001, 115, 9169−9183. (15) Smith, P. E. The alanine dipeptide free energy surface in solution. J. Chem. Phys. 1999, 111, 5568−5579. (16) Chipot, C.; Pohorille, A. Conformational equilibria of terminally blocked single amino acids at the water-hexane interface. A molecular dynamics study. J. Phys. Chem. B 1998, 102, 281−290.
102
dx.doi.org/10.1021/ed400346v | J. Chem. Educ. 2014, 91, 96−102