This section relies on the program Spartan for semi-empirical parameter generation.
If you do not have Spartan available to you, read the section in its entirety, but without
doing any of the specific steps. You will not need anything produced in this section to
complete the tutorial.
As you know by now, the glutamine-cysteine linkage is defined as
a new residue in the CHARMM27 topology file with the name CYG. To actually
determine the optimum parameter values for equilibrium bond lengths, angles, and to determine
the energetic barriers associated with them, we need to use quantum chemistry. Here, we will use
SPARTAN to generate these values. For the purposes of this tutorial, we will calculate the values
semi-empirically; for a true parameterization in accordance with the CHARMM22/CHARMM27 force fields,
we would need to do full ab-initio calculations with the Hartree-Fock 6-31G* basis set. These types
of calculations would take too long to complete in one afternoon. To compute more rigorous quantum
chemistry calculations on small molecules, one often uses parallelizable quantum chemistry packages,
e.g. GAMESS (freely available, highly recommended) or Gaussian.
For this exercise, we also need to use a highly truncated form of the residue because of computational restrictions on the number of atoms one can feasibly calculate with quantum chemistry packages. Presently, the maximum number of atoms one can expect to realistically compute with high level quantum chemistry is approximately 120 atoms. You will likely find that attempting to do larger systems will take a prohibitively long time to compute. When developing the new parameters for this nonstandard residue, we will neglect any solvent effects and consider the molecule to be in the gas phase. The truncated form of CYG that we will be working with has the following formula: CH3-S-CO-CH3 and its IUPAC name is methyl thioethanoate.
You can build a new molecule in Spartan by choosing (from the main Spartan menu):
This command changes the active environment from the main Spartan window to the Model Kit window where you will build and minimize your molecule. Remember that you will need to save your molecule before exiting the Model Kit window, which will take you back to the main Spartan window.
In the Model Kit window, you can build
molecules by selecting atom (and bond) types from the right side bar, and placing
them by in the black window. Atoms can be
connected by clicking on the open valences. Now we will start building the truncated form of CYG:
. All the bonds are single except for
the double bond between the carbon and oxygen (a carbonyl group!). If you make a mistake, you
can easily fix it by selecting Edit
Undo. To build the truncated molecule:
Save the molecule using the menu File Save. Use a name
of your choice for the file.
Now minimize the structure by selecting the Minimize botton. Once it's done, you should see something identical to Fig. 7.
Save the molecule again.
Turn on atom labels by choosing Model: Labels from the main Spartan menu at the top
of your display. When you have done it correctly, you should have
something very similar to Fig. 6.
Close the Model Kit window by clicking on the red button in the top
right corner of the window. Once you close it, you will see explicit
hydrogens appear on any open valences.
To do this, from the main Spartan menu, choose Setup
Calculations and set the following options:
You can submit the job and view its output by choosing:
Now, we are going to display the vibrational modes. The
modes can be viewed through the animation part of the program.
Click: Display Vibrations; when you check the boxes you
can watch individual vibrational modes. Try increasing the amplitude to see what happens; you will
find that the motion of the vibrations become more exaggerated as you increase their amplitudes. Note that
the highest frequencies belong to the hydrogens coming off the methyl groups!
A systematic calculation of the force constants for the
bond stretching and angle bending motion can be obtained from the
ab initio calculations [10], but is beyond the scope of
this exercise.
To complete the exercise, be sure you toggle off any selected vibrations. In order to proceed to the next steps, the molecule should no longer be moving! When this is the case, close the Vibration window.
Now we are going to investigate the partial charges of the atoms. First, we are going to display a surface around the atoms, then we will map the ESP charges of the molecule onto the surface we generated. But before we do this, let's check out the numerical values of the Mulliken charges and ESP charges and compare the difference between them. This can be achieved by doing the following:
When you are finished investigating the charges, close the Configure window.
You can display a potential surface by:
CHARMM starts with the Mullikin charges and then strategically places a water molecule near any
polar group. CHARMM's charges are refined so that they best reproduce the interaction energy
between the water and the polar group. To recreate these values, the partial charges of the
polar group are then altered until this interaction energy is accurately replicated with the MD
simulations. Several concrete examples of this charge-fitting strategy are worked out in Ref. [7].
AMBER, on the other hand, begins with the ESP charges and then uses a program they distribute called
RESP (for Restricted ESP). This charge-fitting scheme replicates the electrostatic potential at a distance out
from the nucleus (much like the surface you built in the previous exercise) with only the atom centers having
charge. Of course, the electrostatic surface generated using this method is often degenerate; meaning there could be more
than one set of charges that would reproduce the potential surface. To resolve this issue, it is best
to further restrict your potential surface by supplying it an initial set of charges for the program's iterations, or
further restrictions derived from chemical insight.
However they are generated, the importance of a good set of partial charges cannot be overstated. Much of the energy (and indirectly then, the behavior) of the system will be calculated using these values.
We recommend that you quit Spartan at this point, and start a new session. Build the molecule again using the steps described in Part 2, Minimize it, Save it using a different name, close the Model Kit window, and continue with the following steps in the main Spartan window.
First we need to identify which dihedral angle we want to investigate. To do this:
You will see popup windows informing you of the start and termination of your job. Click OK when you receive thiese messages. After the job is finished, the program prompts you that the molecular system is being promoted to a list molecules (since you have now multiple conformations of the same molecule). Click OK.
You can animate the molecule through its calculated conformations by pushing the Play botton at the bottom of the window, or you can step through the animation one step at a time using the step forward and step backward bottons.
You can plot the Energy Profile.
If you were not successful in creating the energy drive plot, you may need to close Spartan completely and rebuild
the truncated molecule; then repeat all the steps again carefully.
When you are finished with the exercise, please quit Spartan.