Topics Map > SCS Computing > Software
Topics Map > SCS Computing > How-To
Tutorial - Quantum Chemistry - pKa Estimations using Gaussian
Determining the pKa of Simple Molecules Using Gaussian 2016
Edition 1, September 16, 2019. Updated/reviewed May 17, 2024.
Sydney A. McKee & Taras V. Pogorelov
Department of Chemistry, School of Chemical Sciences,
University of Illinois Urbana-Champaign
This is a tutorial designed to introduce users to calculating pKa values of simple molecules utilizing Gaussian 16. Users are expected to have a foundation in Unix/Linux (see Linux primer tutorial).
Expected tutorial completion time: 2 hrs
Outline: we will use the SCS computer cluster, Lop, to perform Gaussian calculations and GaussView to construct phenol and calculate its pKa in water. Please see the appendix for further information regarding the basis of these calculations, which are based on Chem. Phys. Lett. 367, 145-149, 2003.
1. Software. Lop: Gaussian, GaussView. Personal computer: xserver (MobaXTerm for a Windows) or Terminal (Mac).
2. Starting up. Open Windows PowerShell (or Terminal) window. Connect to Lop. If connecting outside of campus, connect to the VPN first. See SCS Clusters
>ssh –YC netID@lop.scs.illinois.edu
Make a directory for this tutorial:
>mkdir calculate-pka-tutorial
Change your directory to the new directory you just made:
>cd calculate-pka-tutorial
Start up a session of Gaussian with the following command:
>module load gaussian/g16
Open GaussView with the following command:
>gv
Note: you will likely get an error that says, "failed to locate Gaussian scratch directory." Ignore this error by simply hitting the "X". GaussView should load.
3. Constructing a molecule. Construct your molecule of interest in the blue window (which should initially be blank):
Right click anywhere on the blue window and select "Builder." Click "Ring fragment." Right click, select "Builder" again, and then click "Fragment" and select benzene.
Now, right click again and select "Builder". Select "Element Fragment" which should pull up a periodic table. Select the oxygen fragment on the far right. Click on any hydrogen atom of the ring to substitute one of the hydrogens for a hydroxyl group.
With the structure now built:
Select "Edit" and "Clean" to clean up the structure – this will perform a short geometry optimization.
4. Optimization and Frequency Calculations in the gas phase. Now set up a calculation with the built phenol. Right click and select "Calculate" then "Gaussian Calculation Setup..."
Under "Job Type" select "Opt+Freq" to optimize the geometry of the molecule and then calculate energy values of interest. It is important to always optimize first prior to running a frequency calculation.
Under "Method", select "Ground State" "DFT..." "Default Spin" "B3LYP." This is the level of theory for the calculation. For basis set, select "6-31G++d,p". The basis set chosen for this calculation is pretty robust. Note: The lowest level of theory and basis sets necessary should be utilized when running calculations to reduce the expense of them. In this case, each calculation will take approximately 10 minutes to run with the level of theory and basis sets chosen.
This will be phenol without solvation (gaseous phase), so now hit "Submit." Save the Gaussian input file phenol.com. Submit the file to Gaussian through GaussView –you will immediately get a window describing an error with the input file. Ignore this–we will be submitting through the cluster queueing software.
5. Submission of the Gaussian input file. Go back to your terminal window with Triton signed in. In order to submit a job, you must fill GaussView. Execute this with "crt+C". Submit the calculation to Gaussian with the following command:
>submit–g16 –n # –q ibl phenol.com
Replace the "#" in the command with the number of cores desired to be utilized for the calculation. Do not exceed 16.
After submission, check the progress of the calculation with
>qstat -u myscs-login
Once your SCS username no longer appears, the calculation has completed.
6. Calculation of phenoxide in the gas phase. Now, repeat the same process for phenoxide as was done for phenol. The only difference is that the hydrogen from oxygen needs to be deleted, and the charge of the molecule needs to be changed to -1. Delete a hydrogen atom by right clicking and selecting "Builder" then "Select Atoms by Clicking". Click on the hydrogen and then press delete on the keyboard or the X button on the main window. Change the charge in the Method tab for the Gaussian calculation set up:
7. Calculating for solvated molecules. Now the calculations have been completed for phenol and phenoxide in their gaseous phases. Now we need to obtain the same calculations but for their solvated counterparts. Construct the molecules and set up the calculations just as before. Select the "Solvation" tab, then select the IPCM solvation model. Choose water as the solvent.
8. Submission of solvated molecules. Submit the solvated molecules, just as it was done for the gaseous counterparts.
9. Submission for water. Based on using thermodynamics cycle 1 (below), the energies of water and hydronium also need to be calculated. Construct and submit the jobs in the same way as before, ensuring that the charge denoted for hydronium is now +1.
10. Retrieving relevant values. After all submitted jobs completed, extraction of necessary information is done with the Linux "grep" command. Our values of interest are the sum of thermal and electronic energies (the Gibbs free energy). Input the following command:
>grep thermal energies *log
This command will search every log file (using the wildcard '*') in the current directory for the keywords "thermal energies" in the Gaussian calculation output.
The values of interest are the sum of electronic and thermal energies.
11. Calculation of pKa. Open Microsoft Excel and format the sheet in the following fashion:
The sum of electron and thermal energies are the values extracted utilizing the grep command. ΔGs for H3O+ and Gs for H2O + Gx for phenol and phenoxide are found by subtracting the Gibbs free energy for its gaseous phase from its Gibbs free energy for its solvation. ΔGg is calculated following the equation from the thermodynamic cycle 1 (Chem. Phys. Lett., 2003):
Add ΔGs for H3O+ and Gs for phenoxide and subtract ΔGs for H2O + Gs for phenol. Convert to kcal/mole by multiplying by 627.509. This is now the value for ΔGg. Now, all the values to find ΔGaq have been acquired:
Now ΔGsol can be plugged into the following pKa equation to get the uncorrected value:
The concentration of bulk water is 55.49 M.
To find the corrected pKa, simply subtract 4.54 (see Appendix) from the initial estimate to get the final pKa value:
pKa (corrected) = pKa (corrected) – 4.54
The final pKa for phenol was found to be 9.81. The experimental value of phenol's pKa in water is reported to be 9.95. Not bad!
Appendix
Calculating pKas in Gaussian is not straightforward in the sense that it will predict the pKa with a one simple command. Rather, calculations for the Gibbs free energy for the acid, its conjugate base, and water/hydronium need to be obtained in both gaseous phases and solution (solvation). The pKa can then be determined by employing the following thermodynamic cycle:
The equilibrium between the acid and water:
HA + H2O –> A– + H3O+
This presents the equation of ΔGsol then as the following:
Gsol = Gg + Gsolv(A–) + Gsolv(H3O+) – Gsolv(HA) – Gsolv(H2O)
The solved Gsol then is plugged into the following pKa equation, where the bulk concentration of water is 55.49 M:
A more detailed derivation for the pKa equation can be found in Chemical Physics Letters 367, 145-149, 2003. A source of error in the thermodynamic cycle arises with the solvation free energy of hydronium. To mitigate this, the following equation is used:
pKa (corrected) = pKa (corrected) – 4.54