Calculating properties#
The compechem.functions submodule contains a series of methods to “manually” calculate some physical properties of the system under study. The user is expected to provide the exact intended states for the system(s) under study, as the functions will simply take the inputs and generate the output without carrying out any checks. For example, if the user calculates the pKa of a molecule in an unoptimised geometry, the resulting pKa will be that of the unoptimised geometry.
The compechem.functions functions can be imported via the following syntax:
from compechem.functions import calculate_pka
Empirical corrections
The functions to calculate pKa and reduction potentials take into account the self-energy of proton and electron for calculations carried out with GFN2-xTB:
electron self energy = \(111.75\, \mathrm{kcal/mol}\)
proton self energy = \(164.22\, \mathrm{kcal/mol}\)
pKa#
The compechem.functions.pka submodule provides an interface for the computation of the pKa or a given species. The module is, at this time, composed by two functions:
calculate_pka: Computes the pKa of a molecular system and its deprotomer. The user must provide both streucture in the form ofSystemobjects with an already defined electronic energy and possibly a vibronic one.auto_calculate_pka: Computes the pKa of a given molecule by automatically searching the lowest-energy deprotomer using CREST. Once the proper deprotomer has been identified the function take care of the geometry optimization of both structures, the calculation of electronic energies and frequencies.
Note
Both functions return the computed pKa value and set the pKa as a property (system.properties.pka) of the protonated system. The auto_calculate_pka also returns the deprotonated system.
In general terms both functions calculate the pKa of a molecule \(HA\) considering a reaction of the type:
The equilibrium constant of the reaction is calculated as
where \(G_{A^{-}}\) and \(G_{HA}\) are calculated summing the electronic + vibronic energies at the selected/provided level of theory, \(G_{H^{+}} = -270.29 kcal/mol\), \(R = 1.987 \cdot 10^{-3} kcal/(mol \cdot K)\) and \(T = 298.15 K\).
The calculate_pka() function:#
The calculate_pka function takes as arguments the following elements:
protonated(System): molecule in its protonated formdeprotonated(System): molecule in its deprotonated form
Please notice how both the protonated and deprotonated molecules must already be optimized (in water) and must posses a valid electronic energy value. If the vibronic energy is provided, its contribution is taken into account during the calculation.
An example script that can be used to compute the pKa of a molecule is provided in what follows:
from compechem.engines.xtb import XtbInput
from compechem.systems import System
from compechem.functions.pka import calculate_pka
protonated = System("protonated.xyz", charge=0, spin=1)
deprotonated = System("deprotonated.xyz", charge=-1, spin=1)
xtb = XtbInput(solvent="water")
xtb.opt(protonated, inplace=True)
xtb.opt(deprotonated, inplace=True)
pka = calculate_pka(protonated, deprotonated)
The auto_calculate_pka() function:#
The auto_calculate_pka function takes as main argument the protonated molecule structure (in the form of a System object). The molecule is sequentially deprotonated using the CREST deprotomer search routine until the lowest energy deprotomer is identified. Once the deprotomer search has been completed, the structure of both molecules is optimized using the specified level of theory and both electronic and vibronic energies are computed at the user defined level of theory. The routine takes as arguments the following elements:
protonated(System): The protonated molecule for which the pKa must be computed.method_el(Engine): The computational engine to be used in the electronic level of theory calculations.method_vib(Engine): The computational engine to be used in the vibronic level of theory calculations. (optional)method_opt(Engine): The computational engine to be used in the geometry optimization of the protonated molecule and its deprotomers. (optional)ncores(int): The number of cores to be used in the calculations. (optional)maxcore(int): For the engines that supprots it, the memory assigned to each core used in the computation. (optional)
An example script that can be used to compute the pKa of a molecule is provided in what follows:
from compechem.engines.xtb import XtbInput
from compechem.systems import System
from compechem.functions.pka import calculate_pka
protonated = System(f"protonated.xyz", charge=0, spin=1)
xtb = XtbInput(solvent="water")
pka, deprotonated = auto_calculate_pka(
protonated,
method_el=xtb,
method_vib=xtb,
method_opt=xtb,
)
Please notice how the optimized structure of the deprotonated system is also returned together with the pKa value.
1-el redox potential#
Calculates the one-electron reduction potential of a molecule \(MH_{n}\), considering a generic reaction of the type:
provided the following arguments:
oxidised(System): molecule in its oxidised statereduced(System): molecule in its reduced statepH(float, default:7.0): pH at which the reduction potential is calculated
and returns the reduction potential of the molecule considering the provided states at the provided pH, including eventual PCET mechanisms, calculated as:
where \(G_{M^{\cdot(n-1)-}}\) and \(G_{MH_{n}}\) are calculated summing the electronic + vibronic energies at the selected level of theory, \(G_{H^{+}} = -270.29 kcal/mol\), \(F = 23.061 kcal/volt–gram-equivalent\), and \(E_{SHE} = 4.28 V\).
Fukui functions#
The compechem.functions.calculate_fukui function calculates the Fukui functions \(f^+(r)\), \(f^-(r)\) and \(f^0(r)\) associated with a given molecular geometry. The Fukui functions are computed according to the definitions:
Where, given a molecule with \(N\) electrons, \(\rho_{N}(r)\) represents its electronic density while \(\rho_{N\pm1}(r)\) represents the electronic density of the molecule, in the same nuclear configuration, when one electron is either added (\(+1\)) or removed (\(-1\)).
The Fukui functions are both computed as volumetric quantities and saved in a Gaussian Cube compatible format in the output_density folder and as condensed values saved in the System object properties attribute in the form of a dictionary. The condensed Fukui functions are computed by applying the \(f^+\), \(f^-\) and \(f^0\) definitions replacing the charge density with either the Mulliken charges or the Hirshfeld charges (changing the sign accordingly given that a localized electronic density represents an accumulation of electrons hence of negative charge). Please notice how the Hirshfeld charges are supported only by the OrcaInput engine.
Important
Please notice how the Fukui cubes contain the localized Mulliken-charge-based Fukui values in place of the atomic charges. This is explained in the first comment line of each cube file and, for sake of clarity, all the files are saved using the extension .fukui.cube.
The function can be called with the following minimal arguments:
molecule(System): The molecular structure to be used in the computationengine(OrcaInputorXtbInput): The engine defining the level of theory to be used in the calculation.
The function assumes that the molecule supports only singlet and doublet states and switches the spin multeplicity according to the number of electrons. If different spin states needs to be considered the spins_states option can be used to provide the spin multeplicity values as a list.
An example code snippet is provided in what follows:
from compechem.systems import System
from compechem.engines.orca import OrcaInput
from compechem.functions.fukui import calculate_fukui
mol = System("./acetaldehyde.xyz")
orca = OrcaInput(method="PBE", basis_set="def2-SVP")
orca.opt(mol, inplace=True)
calculate_fukui(mol, orca)
print(mol)
That for the acetaldehyde molecule returns the following result:
=========================================================
SYSTEM: acetaldehyde
=========================================================
Number of atoms: 7
Charge: 0
Spin multiplicity: 1
********************** GEOMETRY *************************
Total system mass: 44.0526 amu
----------------------------------------------
index atom x (Å) y (Å) z (Å)
----------------------------------------------
0 C -3.00790 0.17495 0.13810
1 C -3.66800 1.53034 0.11359
2 H -3.22195 2.19896 0.87341
3 H -4.75895 1.41874 0.29379
4 H -3.56793 1.97830 -0.89850
5 O -2.11052 -0.15105 0.88722
6 H -3.42225 -0.55833 -0.62459
----------------------------------------------
********************** PROPERTIES *************************
Geometry level of theory: OrcaInput || method: PBE | basis: def2-SVP | solvent: None
Electronic level of theory: OrcaInput || method: PBE | basis: def2-SVP | solvent: None
Vibronic level of theory: OrcaInput || method: PBE | basis: def2-SVP | solvent: None
Electronic energy: -153.528266323306 Eh
Vibronic energy: 0.02861407 Eh
Helmholtz free energy: None Eh
Gibbs free energy: -153.49965226 Eh
pKa: None
MULLIKEN ANALYSIS
----------------------------------------------
index atom charge spin
----------------------------------------------
0 C 0.05436 0.00000
1 C -0.04823 0.00000
2 H 0.05176 0.00000
3 H 0.05695 0.00000
4 H 0.05677 0.00000
5 O -0.17156 0.00000
6 H -0.00005 0.00000
CONDENSED FUKUI - MULLIKEN
----------------------------------------------
index atom f+ f- f0
----------------------------------------------
0 C 0.23610 0.08077 0.15843
1 C -0.11416 -0.02731 -0.07074
2 H 0.12882 0.11125 0.12003
3 H 0.15563 0.12121 0.13842
4 H 0.15492 0.12121 0.13806
5 O 0.22073 0.35571 0.28822
6 H 0.21797 0.23717 0.22757
HIRSHFELD ANALYSIS
----------------------------------------------
index atom charge spin
----------------------------------------------
0 C 0.14702 0.00000
1 C -0.07279 0.00000
2 H 0.04088 0.00000
3 H 0.04575 0.00000
4 H 0.04567 0.00000
5 O -0.22623 0.00000
6 H 0.01970 0.00000
CONDENSED FUKUI - HIRSHFELD
----------------------------------------------
index atom f+ f- f0
----------------------------------------------
0 C 0.28818 0.15946 0.22382
1 C 0.07190 0.09734 0.08462
2 H 0.05820 0.06063 0.05941
3 H 0.09907 0.07355 0.08631
4 H 0.09850 0.07352 0.08601
5 O 0.25823 0.37205 0.31514
6 H 0.12592 0.16345 0.14468
The volumetric fukui functions can then be plotted using the built in vmd based rendering tool. As an example the following code can be used to render the the \(f^+(r)\) Fukui function.
from compechem.tools.vmdtools import render_fukui_cube
render_fukui_cube(
"./output_densities/acetaldehyde_Fukui_plus.fukui.cube",
include_negative=True,
isovalue=0.02,
)
The rendered volume is saved in a .bmp bitmap image format. In the case of the acetaldehyde molecule considered in the example, the following image is obtained:
