Thermodynamics
Auto3D is mainly designed for generating low-energy 3D structures from the SMILES. It aslo provides a wrapper function calc_thermo to get enthalpy, entropy and Gibbs free energy from the 3D structures. Please note that the thermodynamic calculations are in vaccum.
The source jupyter notebook can be downloaded from here
[ ]:
import os, sys
root = os.path.dirname(os.path.dirname(os.path.abspath("__file__")))
sys.path.append(root)
import Auto3D
from Auto3D import Auto3DOptions, main
[2]:
#Always ensure that you have the latest version
print(Auto3D.__version__)
2.2.10
Generate low-energy conformers with Auto3D
Thermodynamic calculations requires proper 3D structures as the input. This time, we will use ANI2x as the optimizing engine.
[ ]:
if __name__ == "__main__":
path = os.path.join(root, "example/files/smiles.smi")
config = Auto3DOptions(path=path, k=1, optimizing_engine="ANI2x", use_gpu=False) # Configure Auto3D
out = main(config) # Run Auto3D
print(out)
Calculate thermodynamic properties with the 3D structures
[4]:
from Auto3D.ASE.thermo import calc_thermo
[5]:
"""
If the thermodynamic properties are calculated in 298 K, it's straightforward to get the thermodynamic properties.
"""
out_thermo = calc_thermo(out, "ANI2x", opt_tol=0.003)
print(out_thermo) #enthalpy, entropy and Gibbs free energy are stored in the SDF file
0%| | 0/4 [00:00<?, ?it/s]
smi1
/home/jack/miniconda3/envs/py39/lib/python3.9/site-packages/torchani/ase.py:50: UserWarning: Creating a tensor from a list of numpy.ndarrays is extremely slow. Please consider converting the list to a single numpy.ndarray with numpy.array() before converting to a tensor. (Triggered internally at /opt/conda/conda-bld/pytorch_1682343964576/work/torch/csrc/utils/tensor_new.cpp:245.)
cell = torch.tensor(self.atoms.get_cell(complete=True),
25%|██▌ | 1/4 [00:01<00:05, 1.78s/it]
Enthalpy components at T = 298.00 K:
===============================
E_pot -13045.663 eV
E_ZPE 5.113 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.039 eV
Cv_vib (0->T) 0.197 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -13040.251 eV
===============================
Entropy components at T = 298.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017751 eV/K 0.529 eV
S_rot 0.0013167 eV/K 0.392 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0012313 eV/K 0.367 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0043220 eV/K 1.288 eV
=================================================
Enthalpy components at T = 298.00 K:
===============================
E_pot -13045.663 eV
E_ZPE 5.113 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.039 eV
Cv_vib (0->T) 0.197 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -13040.251 eV
===============================
Entropy components at T = 298.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017751 eV/K 0.529 eV
S_rot 0.0013167 eV/K 0.392 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0012313 eV/K 0.367 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0043220 eV/K 1.288 eV
=================================================
Free energy components at T = 298.00 K and P = 101325.0 Pa:
=======================
H -13040.251 eV
-T*S -1.288 eV
-----------------------
G -13041.539 eV
=======================
smi2
50%|█████ | 2/4 [00:02<00:02, 1.14s/it]
Enthalpy components at T = 298.00 K:
===============================
E_pot -6324.139 eV
E_ZPE 3.132 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.039 eV
Cv_vib (0->T) 0.096 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -6320.808 eV
===============================
Entropy components at T = 298.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0016811 eV/K 0.501 eV
S_rot 0.0011114 eV/K 0.331 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0006091 eV/K 0.182 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0034004 eV/K 1.013 eV
=================================================
Enthalpy components at T = 298.00 K:
===============================
E_pot -6324.139 eV
E_ZPE 3.132 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.039 eV
Cv_vib (0->T) 0.096 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -6320.808 eV
===============================
Entropy components at T = 298.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0016811 eV/K 0.501 eV
S_rot 0.0011114 eV/K 0.331 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0006091 eV/K 0.182 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0034004 eV/K 1.013 eV
=================================================
Free energy components at T = 298.00 K and P = 101325.0 Pa:
=======================
H -6320.808 eV
-T*S -1.013 eV
-----------------------
G -6321.821 eV
=======================
smi3
optiimize the input geometry
Step Time Energy fmax
BFGS: 0 15:16:07 -9406.938834 0.0107
BFGS: 1 15:16:07 -9406.938833 0.0152
BFGS: 2 15:16:07 -9406.938835 0.0024
75%|███████▌ | 3/4 [00:03<00:01, 1.00s/it]
Enthalpy components at T = 298.00 K:
===============================
E_pot -9406.939 eV
E_ZPE 3.435 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.039 eV
Cv_vib (0->T) 0.140 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -9403.261 eV
===============================
Entropy components at T = 298.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017235 eV/K 0.514 eV
S_rot 0.0011976 eV/K 0.357 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0008407 eV/K 0.251 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0037606 eV/K 1.121 eV
=================================================
Enthalpy components at T = 298.00 K:
===============================
E_pot -9406.939 eV
E_ZPE 3.435 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.039 eV
Cv_vib (0->T) 0.140 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -9403.261 eV
===============================
Entropy components at T = 298.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017235 eV/K 0.514 eV
S_rot 0.0011976 eV/K 0.357 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0008407 eV/K 0.251 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0037606 eV/K 1.121 eV
=================================================
Free energy components at T = 298.00 K and P = 101325.0 Pa:
=======================
H -9403.261 eV
-T*S -1.121 eV
-----------------------
G -9404.382 eV
=======================
smi4
100%|██████████| 4/4 [00:05<00:00, 1.29s/it]
Enthalpy components at T = 298.00 K:
===============================
E_pot -9499.111 eV
E_ZPE 4.839 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.039 eV
Cv_vib (0->T) 0.195 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -9493.974 eV
===============================
Entropy components at T = 298.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017382 eV/K 0.518 eV
S_rot 0.0012564 eV/K 0.374 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0013657 eV/K 0.407 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0043591 eV/K 1.299 eV
=================================================
Enthalpy components at T = 298.00 K:
===============================
E_pot -9499.111 eV
E_ZPE 4.839 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.039 eV
Cv_vib (0->T) 0.195 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -9493.974 eV
===============================
Entropy components at T = 298.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017382 eV/K 0.518 eV
S_rot 0.0012564 eV/K 0.374 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0013657 eV/K 0.407 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0043591 eV/K 1.299 eV
=================================================
Free energy components at T = 298.00 K and P = 101325.0 Pa:
=======================
H -9493.974 eV
-T*S -1.299 eV
-----------------------
G -9495.273 eV
=======================
Number of failed thermo calculations: 0
Number of successful thermo calculations: 4
/home/jack/Auto3D_pkg/example/files/smiles_20240329-151158-762973/smiles_out_ANI2x_G.sdf
[6]:
"""
If the thermodynamic properties are NOT calculated in 298 K,
you need to define a function that gets the unique ID and temperature for each molecules in the input path,
then pass the custom function as the value of argument get_mol_idx_t.
For example, the following functions works for Auto3D output. It sets the thermodynamic calculation temperature at 273 K and gets the molecular ID. You can set the temperature to any value you want for a given molecule.
You can customize the function for other special needs.
"""
def custom_func(mol):
'''The mol object is an RDKit Molecule object.'''
id = mol.GetProp("_Name")
t = 273
return (id, t)
out_thermo = calc_thermo(out, "ANI2x", get_mol_idx_t=custom_func, opt_tol=0.003)
print(out_thermo) #enthalpy, entropy and Gibbs free energy are stored in the SDF file
0%| | 0/4 [00:00<?, ?it/s]
smi1
25%|██▌ | 1/4 [00:01<00:04, 1.35s/it]
Enthalpy components at T = 273.00 K:
===============================
E_pot -13045.663 eV
E_ZPE 5.113 eV
Cv_trans (0->T) 0.035 eV
Cv_rot (0->T) 0.035 eV
Cv_vib (0->T) 0.164 eV
(C_v -> C_p) 0.024 eV
-------------------------------
H -13040.292 eV
===============================
Entropy components at T = 273.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017562 eV/K 0.479 eV
S_rot 0.0013054 eV/K 0.356 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0011171 eV/K 0.305 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0041776 eV/K 1.140 eV
=================================================
Enthalpy components at T = 273.00 K:
===============================
E_pot -13045.663 eV
E_ZPE 5.113 eV
Cv_trans (0->T) 0.035 eV
Cv_rot (0->T) 0.035 eV
Cv_vib (0->T) 0.164 eV
(C_v -> C_p) 0.024 eV
-------------------------------
H -13040.292 eV
===============================
Entropy components at T = 273.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017562 eV/K 0.479 eV
S_rot 0.0013054 eV/K 0.356 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0011171 eV/K 0.305 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0041776 eV/K 1.140 eV
=================================================
Free energy components at T = 273.00 K and P = 101325.0 Pa:
=======================
H -13040.292 eV
-T*S -1.140 eV
-----------------------
G -13041.432 eV
=======================
smi2
50%|█████ | 2/4 [00:02<00:01, 1.05it/s]
Enthalpy components at T = 273.00 K:
===============================
E_pot -6324.139 eV
E_ZPE 3.132 eV
Cv_trans (0->T) 0.035 eV
Cv_rot (0->T) 0.035 eV
Cv_vib (0->T) 0.081 eV
(C_v -> C_p) 0.024 eV
-------------------------------
H -6320.832 eV
===============================
Entropy components at T = 273.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0016622 eV/K 0.454 eV
S_rot 0.0011000 eV/K 0.300 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0005548 eV/K 0.151 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0033159 eV/K 0.905 eV
=================================================
Enthalpy components at T = 273.00 K:
===============================
E_pot -6324.139 eV
E_ZPE 3.132 eV
Cv_trans (0->T) 0.035 eV
Cv_rot (0->T) 0.035 eV
Cv_vib (0->T) 0.081 eV
(C_v -> C_p) 0.024 eV
-------------------------------
H -6320.832 eV
===============================
Entropy components at T = 273.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0016622 eV/K 0.454 eV
S_rot 0.0011000 eV/K 0.300 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0005548 eV/K 0.151 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0033159 eV/K 0.905 eV
=================================================
Free energy components at T = 273.00 K and P = 101325.0 Pa:
=======================
H -6320.832 eV
-T*S -0.905 eV
-----------------------
G -6321.737 eV
=======================
smi3
optiimize the input geometry
Step Time Energy fmax
BFGS: 0 15:17:15 -9406.938834 0.0107
BFGS: 1 15:17:16 -9406.938833 0.0152
BFGS: 2 15:17:16 -9406.938835 0.0024
75%|███████▌ | 3/4 [00:02<00:00, 1.12it/s]
Enthalpy components at T = 273.00 K:
===============================
E_pot -9406.939 eV
E_ZPE 3.435 eV
Cv_trans (0->T) 0.035 eV
Cv_rot (0->T) 0.035 eV
Cv_vib (0->T) 0.118 eV
(C_v -> C_p) 0.024 eV
-------------------------------
H -9403.291 eV
===============================
Entropy components at T = 273.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017046 eV/K 0.465 eV
S_rot 0.0011862 eV/K 0.324 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0007638 eV/K 0.209 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0036535 eV/K 0.997 eV
=================================================
Enthalpy components at T = 273.00 K:
===============================
E_pot -9406.939 eV
E_ZPE 3.435 eV
Cv_trans (0->T) 0.035 eV
Cv_rot (0->T) 0.035 eV
Cv_vib (0->T) 0.118 eV
(C_v -> C_p) 0.024 eV
-------------------------------
H -9403.291 eV
===============================
Entropy components at T = 273.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017046 eV/K 0.465 eV
S_rot 0.0011862 eV/K 0.324 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0007638 eV/K 0.209 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0036535 eV/K 0.997 eV
=================================================
Free energy components at T = 273.00 K and P = 101325.0 Pa:
=======================
H -9403.291 eV
-T*S -0.997 eV
-----------------------
G -9404.289 eV
=======================
smi4
100%|██████████| 4/4 [00:03<00:00, 1.06it/s]
Enthalpy components at T = 273.00 K:
===============================
E_pot -9499.111 eV
E_ZPE 4.839 eV
Cv_trans (0->T) 0.035 eV
Cv_rot (0->T) 0.035 eV
Cv_vib (0->T) 0.166 eV
(C_v -> C_p) 0.024 eV
-------------------------------
H -9494.012 eV
===============================
Entropy components at T = 273.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017193 eV/K 0.469 eV
S_rot 0.0012450 eV/K 0.340 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0012634 eV/K 0.345 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0042266 eV/K 1.154 eV
=================================================
Enthalpy components at T = 273.00 K:
===============================
E_pot -9499.111 eV
E_ZPE 4.839 eV
Cv_trans (0->T) 0.035 eV
Cv_rot (0->T) 0.035 eV
Cv_vib (0->T) 0.166 eV
(C_v -> C_p) 0.024 eV
-------------------------------
H -9494.012 eV
===============================
Entropy components at T = 273.00 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0017193 eV/K 0.469 eV
S_rot 0.0012450 eV/K 0.340 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0012634 eV/K 0.345 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0042266 eV/K 1.154 eV
=================================================
Free energy components at T = 273.00 K and P = 101325.0 Pa:
=======================
H -9494.012 eV
-T*S -1.154 eV
-----------------------
G -9495.166 eV
=======================
Number of failed thermo calculations: 0
Number of successful thermo calculations: 4
/home/jack/Auto3D_pkg/example/files/smiles_20240329-151158-762973/smiles_out_ANI2x_G.sdf
[7]:
help(calc_thermo)
Help on function calc_thermo in module Auto3D.ASE.thermo:
calc_thermo(path: str, model_name: str, get_mol_idx_t=None, gpu_idx=0, opt_tol=0.0002, opt_steps=5000)
ASE interface for calculation thermo properties using ANI2x, ANI2xt or AIMNET.
:param path: Input sdf file
:type path: str
:param model_name: ANI2x, ANI2xt or AIMNET
:type model_name: str
:param get_mol_idx_t: A function that returns (idx, T) from a pybel mol object, by default using the 298 K temperature, defaults to None
:type get_mol_idx_t: function, optional
:param gpu_idx: GPU cuda index, defaults to 0
:type gpu_idx: int, optional
:param opt_tol: Convergence_threshold for geometry optimization, defaults to 0.0002
:type opt_tol: float, optional
:param opt_steps: Maximum geometry optimization steps, defaults to 5000
:type opt_steps: int, optional
[ ]: