Tutorial - Quantum Chemistry - pKa Estimations using Gaussian

Determining the pKa of Simple Molecules Using Gaussian 2016

Edition 1, September 16, 2019. Updated/reviewed December 12, 2023.

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. ΔG_s for H₃O⁺ and G_s for H₂O + G_x for phenol and phenoxide are found by subtracting the Gibbs free energy for its gaseous phase from its Gibbs free energy for its solvation. ΔG_g is calculated following the equation from the thermodynamic cycle 1 (Chem. Phys. Lett., 2003): 

Add ΔG_s for H₃O⁺ and G_s for phenoxide and subtract ΔG_s for H₂O + G_s for phenol. Convert to kcal/mole by multiplying by 627.509. This is now the value for ΔG_g. Now, all the values to find ΔG_aq have been acquired: 

Now ΔG_sol 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: 

pK_a (corrected) = pK_a (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 + H₂O –> A^– + H₃O⁺

This presents the equation of ΔG_sol then as the following: 

G_sol = G_g + G_solv(A^–) + G_solv(H₃O⁺) – G_solv(HA) – G_solv(H₂O)

The solved G_sol 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: 

pK_a (corrected) = pK_a (corrected) – 4.54