The minimum energy structures found during the conformational analysis represent structures that the molecule adopts at 0 K in the gas phase. Most interesting chemistry, however, happens at much warmer temperatures and in solutions. The set of methods that allow computational study of molecules at real temperatures in solutions is known as "molecular simulations". In essence, molecular simulations allow estimation of properties of neat liquids and solutions by averaging over a large number of solution configurations that are sampled according to the Boltzmann distribution. The Monte Carlo method provides one way to perform simulations of dissolved molecules. One versatile program that implements Monte Carlo molecular simulations in the BOSS, written by Prof. Jorgensen at Yale (http://zarbi.chem.yale.edu/) The following tutorial illustrates how BOSS can be used to learn about the structure of ethanol in aqueous solution.
BOSS is a complex program (see Manual) that uses several text files for input. It currently lacks a graphical user interface but the results of calculations can be easily visualized using standard molecular visualization (e.g. MOLDEN or gOpenMol) and plotting (e.g. Xmgrace) programs. The most important files are
Copy the files that you need for following the tutorial into a subdirectory in your directory. The command cp -r ../etoh_source ./etoh should create subdirectory etoh and copy the needed files. Change into this subdirectory and verify that you have six files.
Inspect the structure file, named here etoh.z. The molecular structure is described in the Z-matrix format. The Z-matrix format defines the structure via connectivity of atoms using the distance, angle, and torsion. For example, the sixth atom (C of CH3, atom type 135) is connected to the fifth atom (C of CH2, atom type 157) by a bond that is 1.53 Ang long, forms an angle of 108.9° with the first atom (O, atom type 154) and forms a torsional angle of 180° with the fourth atom (H, atom type 155). The z-matrix is a very important concept in computational chemistry and you need to understand how to construct z-matrices. MOLDEN provides a graphical interface that facilitates construction of z-matrices. BOSS z-matrices can be generated from other coordinate files using the program autozmat.
Inspect the bond stretching and angle bending parameter file. For the purpose of the tutorial, a simplified copy of the parameter file is provided by retaining only the parameters needed to describe ethanol. You can inspect this file by typing more etoh_oplsaa.sb into the UNIX shell. The first column defines the bonds and angles by atom classes, the second column gives the harmonic force constant and the third column gives the reference bond length or angle values. The optional last column is for comments. Because BOSS is written in Fortran, the position of fields in the parameter file should be kept as is.
Inspect the nonbonded and torsional parameter file. For the purpose of the tutorial, a simplified copy of parameter file is provided by retaining only the parameters needed to simulate ethanol in water. You can inspect this file by typing more etoh_oplsaa.par into the UNIX shell. This file has two parts. The first part defines nonbonded parameters--charge, sigma, and epsilon-- for each atom type. It also establishes connection between numerical atom types and letter-coded atom classes. For examples, both atom types 135 and 157 belong to the class of tetrahedral carbons (CT) but these two atoms carry different partial charges because type 135 is for the terminal methyl group while type 157 is for the methylene next to the oxygen.
The second part of this file contains torsional parameters. Torsional parameters are typically assigned based on atom classes. For example, the central CC-OH torsion in ethanol can be described by atom classes CT-CT-OH-HO. The program recognizes that the torsion defined in the z-matrix by atom types 135-157-154-155 should be described as the CT-CT-OH-HO torsion because of connection between types and classes that was established earlier. During the actual calculation of torsional energy, values -0.356, -0.174, and 0.492 will be substituted for V1, V2, and V3 in the torsional energy expression.
Inspect the job control file. A file suitable for the Monte Carlo simulation of dissolved ethanol in water is provided as mc_npt.head. Notice that we are doing a simulation of ethanol in a box of 260 TIP4P water molecules at 25 ° and 1 atm. The meaning of entries in this file is explained in the BOSS documentation.
Inspect the job command file. The command file that launches the Monte Carlo simulation is called etoh_boxcmd. The command file tells the BOSS program names of various input and output files and specifies how many simulation steps should be taken. Notice that the complete simulation is broken down into NVT equilibration phase, NPT equilibration phase, and the NPT averaging phase. The short NVT equilibration phase lets the molecules to to adopt low energy configurations but keeps the volume and temperature constant. The NPT equilibration phase allows the volume to change (recall that 50 mL ethanol + 50 mL of water gives less than 100 mL of booze) and prepares the system into a state near the equilibrium. The Boltzmann distribution of solute conformations is collected during the averaging period.
Launch the calculation by typing csh etoh_boxcmd into the UNIX shell. The progress of the calculation is shown in the shell while the calculation is running. When the NVT equilibration has finished, you'll see a message 'Running NPT Equilibration'. At this time, open another UNIX shell, and inspect the first structure file, named here MC_nvt.pdb, by typing go MC_nvt.pdb. This opens the PDB file in the program gOpenMol. If you prefer PyMol over gOpenMol, you can use pymol MC_nvt.pdb instead. Notice that the ethanol molecule is surrounded by a cubic box of water molecules. You can inspect the structure of the ethanol molecule by hiding water molecules. To hide water molecules in gOpenMol, go to View -> Atom Mask. First, hide everything by setting the Display state: Off and hitting Apply. Next, display the ethanol molecule by typing 1-9 into the Atom: field, setting the Display State: On and hitting Apply. Measure the CC-OH dihedral angle using View -> Measure Dialog. You need to fill in Atom: fields with numbers corresponding to atom order in the PDB file. Inspection of the text in the PDB file reveals that the atom numbers that define the CC-OH dihedral are 4, 3, 1, 2. Close gOpenMol when you have finished the visual inspection of the structure.
When the calculation has completed (this takes some time), extract the torsional angle probability distributions from the etoh_in_water.av file using the command bar etoh_in_water.av and visualize the probability profiles using the Xmgrace program. Inspection of the text in the MC_nvt.out file reveals that the first torsion (atoms 4- 1- 5- 6 in the Z-matrix) corresponds to the CC-OH torsion and torsion 4 corresponds to the rotation of the methyl group. Type xmgrace dihed1.xmgr in the UNIX shell and analyze the plot. Repeat the analysis for the methyl group rotation. Are these the expected results or do you spot signs that these results cannot be quite trustworthy?
The example job illustrated the essence of Monte Carlo simulations but the accuracy of results was affected seriously by our desire to finish the job quickly. Three most significant factors that affect the results are:
Level 1: Repeat the simulation of ethanol but use this time a larger box (512 waters) and collect 10 million (instead of earlier 2 million) configurations. Also, increase the temperature to 75° from the earlier 25°. Compare the dihedral distributions from this longer and warmer calculation to dihedral distributions from the short tutorial run.
Level 2: Isopropanol also has two low energy structures in the gas phase but the conformational preference in water is not well understood in this system. The /home/chem226/iproh_source directory contains iproh.z, iproh_oplsaa.sb, and iproh_oplsaa.par files as well as the iproh_boxcmd job command file. Perform the Monte Carlo simulation of isopropanol to learn about the torsional preference of the hydroxyl group in this molecule in water. Compare the observed probability distributions in solution with the gas phase energies of the two conformers of isopropanol with the OPLS2001 force field. You can easily obtain the gas phase energy of a conformation by using the xOPT script provided by BOSS. For example, if you have ipoh_cf1.z and ipoh_cf2.z use xOPT ipoh_cf1 to get the energy of the first conformer, and xOPT ipoh_cf2 to get the energy of the second conformer in the gas phase.
Level 3: Propanol has five low energy structures in the gas phase and virtually nothing sure is known about the conformational preference of this molecule in the aqueous solution. The /home/chem226/nproh_source directory contains proh.z structure file. However, the parameter files are these for ethanol and lack some parameters needed for 1-propanol. Inspect files oplsaa.sb and oplsaa.par in $BOSSdir (more $BOSSdir/oplsaa.sb) and append missing parameters to the force field files. Edit the header and command files appropriately to get reasonably converged results (at least 16 million configurations using box of 512 water molecules) and perform the simulation at 75°. Estimate the percentage of the trans-trans conformer (CC-OH around 180 ° and OC-CC around 180 °) among all the conformers from the solution simulation data. Compare this estimate to the relative energy of the trans-trans conformer in the gas phase with the OPLS2001 force field.