Parameter generation using SPARTAN (Optional)

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.

Part 1: Starting Spartan

To start Spartan, open a terminal and type spartan.

Part 2: Building a new molecule

When you launch Spartan, the first thing you see is a large green screen. This is the main Spartan window. If you look at the top of your display, you will see the main Spartan menu.

You can build a new molecule in Spartan by choosing (from the main Spartan menu):

1. File $\rightarrow$ New

Figure 5: The Spartan Program

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: $\mathrm{CH_3-S-CO-CH_3}$. 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 $\rightarrow$ Undo. To build the truncated molecule:

1. Select a methyl group by clicking on it the sp3 carbon option. The sp3 carbon is the carbon with 4 single bonds sticking out of it.
   
2. Now click anywhere inside the black window, and see the carbon appear with 4 open valences.
   
3. Return to the side bar and select the sulfur with two open valances (i.e. the sulfur with two single bonds sticking out of it).
   
4. Click on one of the open valences of the sp3 carbon. The sulfur group should appear!
   
5. Again, return to the side bar and in the Groups dropdown menu, select Carbonyl. A picture of a carbonyl group, i.e. a carbon with a double bond to an oxygen, should appear in the Model Kit window.
   
6. Click on the open valence of the sulfur. The carbonyl group should now be attached.
   
7. Finally, select the another sp3 carbon group
   
8. Now add this last carbon to your existing molecule by clicking on the open valence coming off of the carbonyl group.
   

Save the molecule using the menu File $\rightarrow$ 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.

Figure 6: Building CYG in Spartan

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.

Figure 7: CYG after minimization in Spartan

Part 3: Geometry Optimization and Vibrational Modes

Everything is ready now for a semiempirical quantum mechanics calculation. Now you will optimize the geometry of the structure, calculate electrostatic potential (ESP) charges and vibrational frequencies.

To do this, from the main Spartan menu, choose Setup $\rightarrow$ Calculations and set the following options:

• Calculate: Equilibrium Geometry at Ground state
• with: Semi-Empirical PM3
• Start from: Initial Geometry
• Compute: Frequencies, Electrostatic charges
• Print: Atomic Charges, Vibrational Modes
• Options: check Converge

It should look identical to Fig. 8. When it does, click Save.

Figure 8: Setting up the energy calculation.

You can submit the job and view its output by choosing:

1. Click on Setup $\rightarrow$ Submit. Choose a name for the job save the job files to the directory you're working in.
2. The job will take a couple of minutes to complete. A window will appear indicating the completion of the job
3. Then go to Display$\rightarrow$ Output

The output file lists details about the method used, and lists the Cartesian coordinates of the atoms. You can measure atomic distances and bond angles by using the Geometry pulldown.

Now, we are going to display the vibrational modes. The modes can be viewed through the animation part of the program.

Click: Display $\rightarrow$ 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!

Figure 9: Vibrations Window

Figure: Thioester linkage found in 1ODV.pdb, the photoactive yellow protein.

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.

Part 4: Atomic Charges

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:

1. Select Model $\rightarrow$ Ball and Wire
2. Select Model $\rightarrow$ Configure Labels In the Configure window, select Mulliken first and later, Electrostatic Charges.
3. Select Model $\rightarrow$ Labels

Generally you will see that the partial charge values derived from an ESP calculation are higher than the Mulliken charges. ESP charges are calculated to reproduce the electrostatic potential around an atom, whereas the Mulliken charges are derived from the electron occupancy of orbitals.

When you are finished investigating the charges, close the Configure window.

You can display a potential surface by:

1. Clicking on Setup $\rightarrow$ Surfaces.
2. Then in the new "Add Surface" window, select Surface $\rightarrow$ density and Property $\rightarrow$ potential and click Save.
3. In the new window, select Add.
4. Submit the surface calculation with Setup $\rightarrow$ Submit.
5. When the job completes, you can view the surface with Display $\rightarrow$ Surfaces. In the new "Surfaces" window, you can toggle the surface on and off by checking the box next to density.

Figure 11: Surface of molecule with ESP charges

Be sure to toggle the surface off to proceed to the next part of the exercise and close the Surface window.

Methodologies for Developing the Partial Charges

Now is a good time to discuss the determination of the partial charges appearing in the topology file. Two widely used force fields are CHARMM and AMBER - and they both use different protocols to develop the partial atomic charges.

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.

Part 5: Coordinate driving energy profile

The next exercise will be to calculate the energy of rotation around the H1-C1-S1-C2 dihedral. In Spartan, this is called an Energy Profile calculation (there are other names for the same task, such as conformational sampling, coordinate drive, or dihedral search).

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:

1. Change the labels back to the atom name so you are sure you pick the correct atoms by selecting Model $\rightarrow$ Configure $\rightarrow$ Labels from the main Spartan menu
2. From the main Spartan menu, select Geometry $\rightarrow$ Measure Dihedral.
3. Now return back to the main Spartan window and click on the four atoms of interest in this order: H1-C1-S1-C2

Figure 12: Dihedral selection

1. From the main Spartan menu, select Display $\rightarrow$ Properties. A Molecule Properties window will appear.
2. Go to Build $\rightarrow$ Define Profile menu. Click on Drive dihedral botton, and select the dihedral angle we want to investigate by clicking on the following atoms in the following order: H1, C1, S1, C2.
3. In the window that appears, window, change the values to:
   • Value: 180 To: 540
   • Steps: 12

4. Save the changes, and close the drive dihedral window

5. Select: Setup $\rightarrow$ Calculations from the main Spartan menu.
6. In the Calculations window, select Calculate: Energy Profile at Ground state, with: Semi-empirical PM3, turn off Symmetry, Print: Atomic Charges, check Converge. All of the options under Compute: UV/VIS, IR, and NMR should NOT be checked. Then hit Save.
7. Select Setup $\rightarrow$ Submit from the main Spartan menu to submit the job.

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.

1. In the main Spartan menu, select Display $\rightarrow$ Spreadsheet.
2. In the new window that appears, click on Column $\rightarrow$ E$_{gas}$.
3. Close the Spreadsheet window.
4. In the main Spartan menu, select Display $\rightarrow$ Plots $\rightarrow$ Create. In the Add Column window, keep the X axis set to Molecule, and select the Y axis to be E (kcal/mol), then click OK.

Figure 13: Dihedral energy profile results.

Spartan automatically superimposes the plot of the energy profile on top of the molecule. If you click on the ''play forward" button on the bottom bar in the main window, you can watch the torsion angle change and simultaneously see the corresponding point on the plot. See Fig. 13.

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.