Overview

Review the Molecular Mechanics lecture supplement.

Objective

This week we will use ab initio and semi empirical molecular orbital calculations to obtain partial atomic charges for amino acids; for the most part we will use N-methyl acetamide (NMA) as a model. The electrostatic terms are quite significant in the calculation of the Ramachandran plot for a short peptide in water. The selection of solvation model and of dielectric constants is quite important, and fraught with subtlety. Equally important is the selection of "partial atomic charges" which are lumped charges describing the effect of the charges of the nucleus and electrons associated with each atom; this lumped partial charge is located at the nucleus.

Methods: Orbital Calculations

Partial atomic charges are often derived from quantum mechanical calculations; these calculations can be done at a variety of levels of accuracy. You will set up calculations through the Maestro interface at a "Semi-Empirical" level. You can also try to calculate these values at an ab initio level with Jaguar. The ab initio calculations take much longer than semi-empirical methods and are capable of giving more reliable answers for a much broader range of compounds. We will discuss the differences in these two methods in this excercise. The semi-empirical approach that you will use, (AM1), is capable of giving reasonable results for amides in a variety of results, and can run very quickly on these workstations. Using the semi-empirical approach you will look at trends in the resulting charges as a function of the structure: for example you will study the effects of solvation, hydrogen bonding, side chain, and conformation on the partial charges.

        Methods: Charge Fitting

We will use two methods for obtaining partial charges from the molecular orbitals. The first method, Mulliken charges, is a kind of "bookeeping" in which the electrons from each orbital are distributed between the parent atoms. In the second method, electrostatic fitting, the collection of charges that best reproduces the external field profile of the molecule is calculated. There is often a big discrepancy between the two kinds of charges. The electrostatic method is more consistent with force field calculations, and should be used here. When Kollman et al develop charges for use in the  AMBER force field, a restraint is used that minimizes the charges (yet produces a good agreement with the external field). For a discussion of these fitting methods see Kollman's page on amber and RESP . In this excercise we will develop charges that fit the field as well as possible, without any restraint, and compare their performance in the force field with the default charges developed with RESP.

        Testing the Charges: Using the Charges in Force Fields

The partial charges are a model for the charge distribution that is rather loosely associated with reality and cannot be measured experimentally. Nonetheless there are some important comparisons with experiment: the force field including these charges can be used to predict optimized geometry of the molecule (especially covalent bond lengths and bond angles) and the net electostatic dipole of the molecule. If you were developing force fields, you would also want to successfully predict solubilities, dimerization energies and other thermodynamic properties. Our criterion in this excersize will be the ability of the total force field to reproduce conformational preferences of the molecule-- especially the backbone dihedral or torsional preference of the protein structure.

Conformational data on these capped peptides in solution is of interest but in general is not available, mainly due to the fact that the conformation in solution is highly fluxional. We can compare with the ab initio conformational preferences and we can also compare with the experimental data on conformational preferences of the backbone residues in folded proteins, in structured pepetides and in solid state amides.

Of course the force fields have other parameterized physical attributes of the molecule as well, for example the equilibrium bond lengths and angles and the stretching and bending vibrational frequencies; which you could compare directly to experimental diffraction and vibrational data.

Comparing Force Field and Quantum Chemical Calculations of Relative Potential Energy

You could in principle do a detailed comparison of the Ramachandran plot as computed with QM and MM methods. Using AM1 calculations we will evaluate just the single point energies for the capped alanine at a few interesting locations on the Ramachandran plot and compare with single point energies obtained from the molecular mechanics/amber94 calculations: for example we will compare the helical and sheet regions and the region characteristic for a gamma turn. This comparison can be made for both the gas phase species and the water solvated species.

Calculating Properties of Capped Amino-acids

Force Field Calculations

• Building the structure

Build the amino-acid structure (using a clean workspace) using the build function (Edit -> Build or Alt-B) in maestro.

• Calculating Geometries

Optionally, you can minimize the structure and change the conformation using Macromodel.

Using the ECalc function in Macromodel (Current Energy), you can output some basic information and parameters for the force-field being used. Set the force-field to OPLS-AA and use Ecalc with minimal energy listing. Once the job has completed, open the output file (with mmo extension) and look at the atomic charges. Also take note of the other features of the MM force-field.

You can map the charges and coordinates of these charges to their respective atoms using the Display -> Atom labels function.


• Effects of Solvation

Repeat the above calculation ( without geometry optimization) for the amino-acid in solvent. Summarize the charges (using the jot editor or the emacs editor) and the molecular geometry as you proceed, so that you can make a comparison to the default charges in the force fields.
In the molecular geometry the two most important parameters are the N-C and the C-O bond lengths, which are good measures of the hybridization of the amide.


• Effects of Hydrogen Bonding

Build a dimer of your amino-acid with a linear hydrogen bonding geometry (that is the N-H...O angle and the N...O-C angle are 180 degrees, with the two amides nearly co-planar, and with a typical N...O distance of 2.9 A. To do this, bring two molecules to the screen at the same time, and position the molecule as best you can by hand. You select one molecule at one time by clicking the Advanced button under the Rotate and Translate buttons on the Maestro panel. Select pick to add the the cursor should change to an 'm' with a square to select the molecule. Once selected rotate/translate as you normally would. Deselect by clicking 'clear' in the advanced window.

Once the two molecules are fairly close together, conduct an constrained energy minimization. Under Macromodel -> Minimization, add the distance constraints mentioned above as well as the angle restraints. Conduct a progressive gradient minimization (using the OPLSAA FF) to get the final structure.


On the other hand, the structural coordinates of an isolated dimer from a naturally occuring beta sheet would be more accurate.

Compare the bond lengths and the charges for your dimer to the monomeric structure using ECalc. Compare the total energy for the dimer to the energy of a monomer to obtain an estimate of the hydrogen bond energy in the gas phase.


• Effects of the Non-planar Amide

Next build an amino-acid with a non-planar amide bond: one with a 45 degrees out-of-plane C-N rotation and one that is 90 degrees out-of-plane. Compare the charges, bond lengths and energies to those for the planar structures. If you minimize the structures after rotating them, make sure that you fix this degree of freedom..... it is far from a local minimum.

Calculating Charges for Capped Amino Acids and Peptide

Re-calculate these charges for a set of three capped amino acids in a beta sheet conformation.

Rotate the chi1 to observe the effects of charges on the capped Tyrosine. Chi1 has three stable conformations: trans, g+, and g-. For more of a discussion on chi1, please see the Birkbeck Course on Protein Structure.

Tabulate the energies, geometries and charges to compare them to molecular modelling results next week.

Discussion of More Advanced Calculations: Ab Initio 

You will use Jaguar to calculate point charges for these molecules repeat these exercises for all molecules, if time permits. Open Jaguar under the application menu item in Maestro. The nine buttons at the bottom are used in specifying the jobs to be calculated. You can use the geometry cleanup function to roughly minimize the structure - do not do this for the dimer because the constraints are not included in the minimzation.

First, select your electronic basis set under the basis set button and select the 6-31G** basis set.

Next we will set the electrostatic fitting parameters (ESP). Fit the potential to the atom centers and constrain the total charge. Use spherical grids and include the Mulliken charges. Now name your job and calculate the potentials. Once the job has completed, open the output file (.out extension) and read the point atomic charges from the electrostatic fit.

Compare and Discuss your Results

Summarize your results and discuss trends. Feel free to include calculations from structures of actual proteins - try to limit the number of atoms to ~20 in the QM calculations.

Go To Top Course Page