Quantum Chemical Tutorial


Gaussian03 files for this tutorial


Dynamo files for this tutorial



Determining the energetic profile of a model reaction in Glutathione S-Transferase (GSTs)

Before initiating a full-fledged QM/MM simulation of an enzyme reaction, it is often worthwhile to examine the gas phase energetics of a smaller representation of the enzyme reaction.  This exercise can then be used to judge the appropriateness of the QM method (especially when employing semiempirical QM methods in QM/MM simulations) and to compare the gas phase and condensed phase properties.

In this tutorial, we will examine the gas phase energetics of a model reaction that takes place within GSTs; the nucleophilic attack on an epoxide molecule by methylthiolate.  The epoxide molecule resembles a substrate for GST: 7,8-dihydroxy-9,10-oxy-7,8,9,10-tetrahydrobenzo[a]pyrene (BPDE) while methylthiolate is a good analogue of the reactive form of a cysteine residue.  Glutathione S-Transferases (GSTs) “detoxify” these type of compounds by catalyzing the addition of glutathione to the electrophilic compound.  For this tutorial, we will compare the results from RHF/6-31G(d) using Gaussian03 and from the AM1 method using DYNAMO.  In an actual application, researchers typically employ much higher levels of theory such as MP2, the DFT method B3LYP or possibly an accurate compound method such as CBS-Q.  These calculations are much more computationally intensive though.  The Gaussian03 calculations will be run on the Jonas platform through the training que while the DYNAMO calculations can be run interactively.

1. Optimize the reactants and non-bonded complex and determine their thermodynamic quantities.  The Gaussian input files are appended with “.dat” while the job submission scripts are appended with “.job”.  To submit the Gaussian job for the epoxide molecule on jonas:

qsub –q training epoxide-rhf.job


       Methylthiolate             +         epoxide (substrate for GST)

After the job is complete, look in the .oXXXX file where XXXX is a number.  This will indicate where the job was run and you can then examine and retrieve the results from there.  The reactant input files are labeled epoxide and methylthiolate.  The non-bonded complex is labeled complex-nb.

We would like to be able to compare the results from an ab initio QM method with our semiempirical QM method.  Most semiempirical methods return an enthalpy of formation (more on this in the lecture material).  Therefore the relevant thermodynamic quantity that we want to retrieve from our outputs is the enthalpy denoted in ***.

Zero-point correction= 0.037193 (Hartree/Particle)
Thermal correction to Energy=                                 0.040228
Thermal correction to Enthalpy=         0.041172
Thermal correction to Gibbs Free Energy=  0.013606
Sum of electronic and zero-point Energies= -437.372406
Sum of electronic and thermal Energies=  -437.369372
*** Sum of electronic and thermal Enthalpies=  -437.368427***
Sum of electronic and thermal Free Energies=    -437.395994


***This number is reported in Hartrees (627.5095 kcal = 1 Hartree)

To calculate the enthalpy of formation with AM1 in DYNAMO:

$ dynamo thiolate-am1 &

The result is returned in kJoules (4.184 kJ = 1 kcal)

2. Optimize the product structure and determine its thermodynamic quantities.  The product is labeled product.

        Product structure

3. Locate the transition state structure and verify by examining the one negative frequency.  These calculations can be notoriously difficult especially if one has a large molecule.  One way to discover a good guess is to perform potential energy scans along a defined simplistic reaction coordinate and take the maximum energy structure as your starting structure.  In this example, we will employ the STQN method in Gaussian indicated with the QST3 keyword.  In DYNAMO, we will use the Baker optimization routine to locate the transition state.  A transition state is found if one imaginary frequency is returned from a normal mode calculation.  The imaginary frequency must then be examined graphically to determine if it deforms the structure in such a manner as to connect the reactants and products.  A further test would be to run reaction path calculations.  The Gaussian output can be visualized with MOLDEN while the DYNAMO output can be visualized with VMD.

    Transition state structure

$ molden transition-rhf.log

To view the DYNAMO results in VMD, do the following:

$ vmd transition-state-am1.xyz

then read into this file the transition-state.dcd file.  Ask for help if needed on this step.
4. Calculating Basis Set Superposition Error (BSSE).

As the atoms of interacting molecules (or of different parts of the same molecule) or two molecules approach one another, their basis functions overlap. Each monomer "borrows" functions from other nearby components, effectively increasing its basis set and improving the calculation of derived properties such as energy.  We will use the counterpoise method to remove this error from the product structure energy.  This input (counterpoise-rhf.dat and .job) only serves as an example on how to run these calculations.  To complete the study, we would also have to run these on the non-bonded complex and the transition state.

5. Determine the gas-phase proton affinity of methylthiolate

The proton affinity is defined as the energy release when a proton is added to the system, computed as the energy difference between the system of interest and the same molecule with one additional proton.  We have already calculated the energy of methylthiolate.  To compute the proton affinity we must simply

Optimize methylthiol (sometime called methylmercaptan) and determine its thermodynamic properties using Gaussian03 and DYNAMO.