##################################### Bio 5476 Protocol for Lab Exercise #3 ##################################### (1) Below is a table containing the molecular weight, density, melting point and boiling point for several simple organic liquids. Your goal in this lab is to pick one of these liquids, perform a liquid molecular dynamics simulation on your liquid, and analyze the results to determine computed values for the density, heat of vaporization, and radial distribution functions, and to monitor approach to conformational equilibrium. Molecule MW Density MP BP Acetic Acid 60.05 1.049 46 118 Acetone 58.08 0.791 -94 56 Acetonitrile 41.05 0.786 -48 82 Ammonia 17.03 0.682 -78 -33 Aniline 93.13 1.022 -6 184 Benzene 78.11 0.874 6 80 Cyclohexane 84.16 0.779 7 81 Cyclopentene 68.12 0.774 -135 44 Dimethyl Formamide 73.10 0.944 -98 153 Ethanol 46.07 0.794 -114 78 Ethane Thiol 62.13 0.839 -145 35 Ethyl Acetate 88.11 0.902 -84 77 Ethyl Iodide 155.97 1.950 -108 71 Ethylene Glycol 62.07 1.113 -13 197 Formamide 45.04 1.134 3 210 Methanol 32.04 0.791 -98 65 Methyl Disulfide 94.20 1.046 -85 109 Methyl t-Butyl Ether 88.15 0.740 -109 52 Pentane 72.15 0.626 -130 36 Phenol 94.11 1.017 41 182 Nitromethane 61.04 1.127 -29 101 Trifluoroethanol 100.01 1.373 -44 79 Trimethyl Amine 59.11 0.636 -117 4 Water 18.02 1.000 0 100 (2) First, obtain a TINKER .xyz file containing a single molecule of your assigned substance. Files for isolated molecules, setup to use the OPLS-AA force are provided on the course website for this lab. You will also need to use the oplsaa.prm file in the above directory. This slightly modified version of the parameter file contains some extra parameters for some of the molecules in the list. (3) Examine the .xyz file and compare the atom types in the column just to the right of the x,y,z-coordinates to the OPLS-AA parameter file obtained in step (2) to verify that the correct atom types are used for your molecule. (4) Minimize your molecule using the "minimize" program to clean up any distortions in the original geometry. Then copy your molecule to a second file, edit this second file to translate the molecule to some new set of coordinates using any text editor, and finally use the TINKER "xyzedit" program to merge your two molecules into a single file containing a dimer. Use minimization to find the optimal energy for the dimer. What is the energy of interaction, ie, the energy of the dimer minus twice the energy of the monomer? Are you sure that you have found the best dimer structure? (5) Run the "xyzedit" program on your minimized monomer, using the option to create a cubic periodic box filled with solvent molecules. You should try to make a box about 20-25 Angstroms on a side, and you will need to calculate the number of molecules to place in the box based on molecular weight, density, etc. Note that 1 Angstrom is equal to 10(-10) meters, and Avogadro's number is 6.02214199 x 10(+23). (6) Make sure the box size is set as the appropriate keyword (a-axis) in your keyfile, then minimize the solvent box to again remove any bad or high energy contacts. You can also set the cutoff distance for the van der Waals interactions (keyword: VDW-CUTOFF) to be somewhat less than half the box dimension, but at least 9-10 Angstroms. For electrostatic interactions turn on Ewald summation using the EWALD keyword, and set the real-space Ewald cutoff to 7 Angstroms using the EWALD-CUTOFF keyword. Also turn on the pair neighbor lists via the NEIGHBOR-LIST keyword. (7) Run an MD simulation in the isothermal-iosbaric ensemble (constant temperature and pressure). Use 1 Atm as the target pressure. Use room temperature (298K = 25C) for the target temperature, as long as your substance is a liquid at that temperature. Otherwise use a target temperature midway between the melting and boiling points. (8) Collect your MD trajectory, saving coordinate snapshots every 0.1 ps or so, and using a 2 fs time step for integration. You will generally need to discard the first portion of the trajectory (perhaps 200 ps) as an "equilibration period". Following equilibration, you can use the remaining "production period" in your analysis to determine the density, heat of vaporization, and radial distributions for the OPLS-AA model of your liquid. Ideally, your production period should be at least a nanosecond or longer. (9) Average properties from your MD run can be determined by running the "mdavg" script on the log file from the dynamics calculation. This will provide the average temperature, pressure, density, potential energy, etc. The instantaneous values of these same quantities can be found by inspecting and "grepping" the log file. In particular, you should compare your computed density with the experimental value. (10) You should figure out how to compute the heat of vaporization from your MD results. See the Leach book, or other recommended course books for further information. It should also be possible to compute an estimate of heat capacity (try this!), although the value your derive will probably be uncertain at best (why?). (11) Radial distribution functions can be computed from your saved MD coordinates using the "radial" program. Try, for example, computing the distribution function between pairs of like polar atoms, if your liquid contains a polar functional group. The theory and form of radial distributions is described in section 6.2.5 of Leach. (12) Finally, if your liquid can exhibit alternative conformations, such as the "chair" and "boat" conformations of cyclohexane, determine the relative percentage of each conformation as a function of the simulation time. Does it look like your system has reached its conformational equilibrium state?