example/thermodynamic_calculation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thermodynamics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`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.\n",
    "\n",
    "The source jupyter notebook can be downloaded from [here](https://github.com/isayevlab/Auto3D_pkg/blob/main/example/thermodynamic_calculation.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "root = os.path.dirname(os.path.dirname(os.path.abspath(\"__file__\")))\n",
    "sys.path.append(root)\n",
    "\n",
    "import Auto3D\n",
    "from Auto3D.auto3D import options, main"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.2.10\n"
     ]
    }
   ],
   "source": [
    "#Always ensure that you have the latest version\n",
    "print(Auto3D.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate low-energy conformers with Auto3D\n",
    "\n",
    "Thermodynamic calculations requires proper 3D structures as the input. This time, we will use `ANI2x` as the optimizing engine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Checking input file...\n",
      "\tThere are 4 SMILES in the input file /home/jack/Auto3D_pkg/example/files/smiles.smi. \n",
      "\tAll SMILES and IDs are valid.\n",
      "Suggestions for choosing isomer_engine and optimizing_engine: \n",
      "\tIsomer engine options: RDKit and Omega.\n",
      "\tOptimizing engine options: ANI2x, ANI2xt and AIMNET.\n",
      "The available memory is 251 GB.\n",
      "The task will be divided into 1 jobs.\n",
      "Job1, number of inputs: 4\n",
      "\n",
      "\n",
      "Isomer generation for job1\n",
      "Enumerating cis/tran isomers for unspecified double bonds...\n",
      "Enumerating R/S isomers for unspecified atomic centers...\n",
      "Removing enantiomers...\n",
      "Enumerating conformers/rotamers, removing duplicates...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 4/4 [00:00<00:00,  9.38it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Optimizing on job1\n",
      "Preparing for parallel optimizing... (Max optimization steps: 5000)\n",
      "Total 3D conformers: 259\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 10%|█         | 501/5000 [01:51<06:51, 10.94it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total 3D structures: 259  Converged: 205   Dropped(Oscillating): 0    Active: 54\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 20%|█▉        | 999/5000 [02:15<01:48, 36.92it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total 3D structures: 259  Converged: 257   Dropped(Oscillating): 0    Active: 2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 27%|██▋       | 1356/5000 [02:22<06:23,  9.51it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization finished at step 1357:   Total 3D structures: 259  Converged: 259   Dropped(Oscillating): 0    Active: 0\n",
      "Begin to select structures that satisfy the requirements...\n",
      "Energy unit: Hartree if implicit.\n",
      "Program running time: 3 minute(s)\n",
      "Output path: /home/jack/Auto3D_pkg/example/files/smiles_20240329-151158-762973/smiles_out.sdf\n",
      "/home/jack/Auto3D_pkg/example/files/smiles_20240329-151158-762973/smiles_out.sdf\n"
     ]
    }
   ],
   "source": [
    "if __name__ == \"__main__\":\n",
    "    path = os.path.join(root, \"example/files/smiles.smi\")\n",
    "    args = options(path, k=1, optimizing_engine=\"ANI2x\", use_gpu=False)   #args specify the parameters for Auto3D \n",
    "    out = main(args)            #main acceps the parameters and run Auto3D\n",
    "    print(out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate thermodynamic properties with the 3D structures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Auto3D.ASE.thermo import calc_thermo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 0/4 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "smi1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/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.)\n",
      "  cell = torch.tensor(self.atoms.get_cell(complete=True),\n",
      " 25%|██▌       | 1/4 [00:01<00:05,  1.78s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enthalpy components at T = 298.00 K:\n",
      "===============================\n",
      "E_pot             -13045.663 eV\n",
      "E_ZPE                  5.113 eV\n",
      "Cv_trans (0->T)        0.039 eV\n",
      "Cv_rot (0->T)          0.039 eV\n",
      "Cv_vib (0->T)          0.197 eV\n",
      "(C_v -> C_p)           0.026 eV\n",
      "-------------------------------\n",
      "H                 -13040.251 eV\n",
      "===============================\n",
      "Entropy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017751 eV/K        0.529 eV\n",
      "S_rot              0.0013167 eV/K        0.392 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0012313 eV/K        0.367 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0043220 eV/K        1.288 eV\n",
      "=================================================\n",
      "Enthalpy components at T = 298.00 K:\n",
      "===============================\n",
      "E_pot             -13045.663 eV\n",
      "E_ZPE                  5.113 eV\n",
      "Cv_trans (0->T)        0.039 eV\n",
      "Cv_rot (0->T)          0.039 eV\n",
      "Cv_vib (0->T)          0.197 eV\n",
      "(C_v -> C_p)           0.026 eV\n",
      "-------------------------------\n",
      "H                 -13040.251 eV\n",
      "===============================\n",
      "\n",
      "Entropy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017751 eV/K        0.529 eV\n",
      "S_rot              0.0013167 eV/K        0.392 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0012313 eV/K        0.367 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0043220 eV/K        1.288 eV\n",
      "=================================================\n",
      "\n",
      "Free energy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=======================\n",
      "    H     -13040.251 eV\n",
      " -T*S         -1.288 eV\n",
      "-----------------------\n",
      "    G     -13041.539 eV\n",
      "=======================\n",
      "smi2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 50%|█████     | 2/4 [00:02<00:02,  1.14s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enthalpy components at T = 298.00 K:\n",
      "===============================\n",
      "E_pot              -6324.139 eV\n",
      "E_ZPE                  3.132 eV\n",
      "Cv_trans (0->T)        0.039 eV\n",
      "Cv_rot (0->T)          0.039 eV\n",
      "Cv_vib (0->T)          0.096 eV\n",
      "(C_v -> C_p)           0.026 eV\n",
      "-------------------------------\n",
      "H                  -6320.808 eV\n",
      "===============================\n",
      "Entropy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0016811 eV/K        0.501 eV\n",
      "S_rot              0.0011114 eV/K        0.331 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0006091 eV/K        0.182 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0034004 eV/K        1.013 eV\n",
      "=================================================\n",
      "Enthalpy components at T = 298.00 K:\n",
      "===============================\n",
      "E_pot              -6324.139 eV\n",
      "E_ZPE                  3.132 eV\n",
      "Cv_trans (0->T)        0.039 eV\n",
      "Cv_rot (0->T)          0.039 eV\n",
      "Cv_vib (0->T)          0.096 eV\n",
      "(C_v -> C_p)           0.026 eV\n",
      "-------------------------------\n",
      "H                  -6320.808 eV\n",
      "===============================\n",
      "\n",
      "Entropy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0016811 eV/K        0.501 eV\n",
      "S_rot              0.0011114 eV/K        0.331 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0006091 eV/K        0.182 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0034004 eV/K        1.013 eV\n",
      "=================================================\n",
      "\n",
      "Free energy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=======================\n",
      "    H      -6320.808 eV\n",
      " -T*S         -1.013 eV\n",
      "-----------------------\n",
      "    G      -6321.821 eV\n",
      "=======================\n",
      "smi3\n",
      "optiimize the input geometry\n",
      "      Step     Time          Energy         fmax\n",
      "BFGS:    0 15:16:07    -9406.938834        0.0107\n",
      "BFGS:    1 15:16:07    -9406.938833        0.0152\n",
      "BFGS:    2 15:16:07    -9406.938835        0.0024\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 75%|███████▌  | 3/4 [00:03<00:01,  1.00s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enthalpy components at T = 298.00 K:\n",
      "===============================\n",
      "E_pot              -9406.939 eV\n",
      "E_ZPE                  3.435 eV\n",
      "Cv_trans (0->T)        0.039 eV\n",
      "Cv_rot (0->T)          0.039 eV\n",
      "Cv_vib (0->T)          0.140 eV\n",
      "(C_v -> C_p)           0.026 eV\n",
      "-------------------------------\n",
      "H                  -9403.261 eV\n",
      "===============================\n",
      "Entropy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017235 eV/K        0.514 eV\n",
      "S_rot              0.0011976 eV/K        0.357 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0008407 eV/K        0.251 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0037606 eV/K        1.121 eV\n",
      "=================================================\n",
      "Enthalpy components at T = 298.00 K:\n",
      "===============================\n",
      "E_pot              -9406.939 eV\n",
      "E_ZPE                  3.435 eV\n",
      "Cv_trans (0->T)        0.039 eV\n",
      "Cv_rot (0->T)          0.039 eV\n",
      "Cv_vib (0->T)          0.140 eV\n",
      "(C_v -> C_p)           0.026 eV\n",
      "-------------------------------\n",
      "H                  -9403.261 eV\n",
      "===============================\n",
      "\n",
      "Entropy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017235 eV/K        0.514 eV\n",
      "S_rot              0.0011976 eV/K        0.357 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0008407 eV/K        0.251 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0037606 eV/K        1.121 eV\n",
      "=================================================\n",
      "\n",
      "Free energy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=======================\n",
      "    H      -9403.261 eV\n",
      " -T*S         -1.121 eV\n",
      "-----------------------\n",
      "    G      -9404.382 eV\n",
      "=======================\n",
      "smi4\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 4/4 [00:05<00:00,  1.29s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enthalpy components at T = 298.00 K:\n",
      "===============================\n",
      "E_pot              -9499.111 eV\n",
      "E_ZPE                  4.839 eV\n",
      "Cv_trans (0->T)        0.039 eV\n",
      "Cv_rot (0->T)          0.039 eV\n",
      "Cv_vib (0->T)          0.195 eV\n",
      "(C_v -> C_p)           0.026 eV\n",
      "-------------------------------\n",
      "H                  -9493.974 eV\n",
      "===============================\n",
      "Entropy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017382 eV/K        0.518 eV\n",
      "S_rot              0.0012564 eV/K        0.374 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0013657 eV/K        0.407 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0043591 eV/K        1.299 eV\n",
      "=================================================\n",
      "Enthalpy components at T = 298.00 K:\n",
      "===============================\n",
      "E_pot              -9499.111 eV\n",
      "E_ZPE                  4.839 eV\n",
      "Cv_trans (0->T)        0.039 eV\n",
      "Cv_rot (0->T)          0.039 eV\n",
      "Cv_vib (0->T)          0.195 eV\n",
      "(C_v -> C_p)           0.026 eV\n",
      "-------------------------------\n",
      "H                  -9493.974 eV\n",
      "===============================\n",
      "\n",
      "Entropy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017382 eV/K        0.518 eV\n",
      "S_rot              0.0012564 eV/K        0.374 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0013657 eV/K        0.407 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0043591 eV/K        1.299 eV\n",
      "=================================================\n",
      "\n",
      "Free energy components at T = 298.00 K and P = 101325.0 Pa:\n",
      "=======================\n",
      "    H      -9493.974 eV\n",
      " -T*S         -1.299 eV\n",
      "-----------------------\n",
      "    G      -9495.273 eV\n",
      "=======================\n",
      "Number of failed thermo calculations:  0\n",
      "Number of successful thermo calculations:  4\n",
      "/home/jack/Auto3D_pkg/example/files/smiles_20240329-151158-762973/smiles_out_ANI2x_G.sdf\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "If the thermodynamic properties are calculated in 298 K, it's straightforward to get the thermodynamic properties.\n",
    "\"\"\"\n",
    "\n",
    "out_thermo = calc_thermo(out, \"ANI2x\", opt_tol=0.003)\n",
    "print(out_thermo)  #enthalpy, entropy and Gibbs free energy are stored in the SDF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 0/4 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "smi1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 25%|██▌       | 1/4 [00:01<00:04,  1.35s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enthalpy components at T = 273.00 K:\n",
      "===============================\n",
      "E_pot             -13045.663 eV\n",
      "E_ZPE                  5.113 eV\n",
      "Cv_trans (0->T)        0.035 eV\n",
      "Cv_rot (0->T)          0.035 eV\n",
      "Cv_vib (0->T)          0.164 eV\n",
      "(C_v -> C_p)           0.024 eV\n",
      "-------------------------------\n",
      "H                 -13040.292 eV\n",
      "===============================\n",
      "Entropy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017562 eV/K        0.479 eV\n",
      "S_rot              0.0013054 eV/K        0.356 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0011171 eV/K        0.305 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0041776 eV/K        1.140 eV\n",
      "=================================================\n",
      "Enthalpy components at T = 273.00 K:\n",
      "===============================\n",
      "E_pot             -13045.663 eV\n",
      "E_ZPE                  5.113 eV\n",
      "Cv_trans (0->T)        0.035 eV\n",
      "Cv_rot (0->T)          0.035 eV\n",
      "Cv_vib (0->T)          0.164 eV\n",
      "(C_v -> C_p)           0.024 eV\n",
      "-------------------------------\n",
      "H                 -13040.292 eV\n",
      "===============================\n",
      "\n",
      "Entropy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017562 eV/K        0.479 eV\n",
      "S_rot              0.0013054 eV/K        0.356 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0011171 eV/K        0.305 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0041776 eV/K        1.140 eV\n",
      "=================================================\n",
      "\n",
      "Free energy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=======================\n",
      "    H     -13040.292 eV\n",
      " -T*S         -1.140 eV\n",
      "-----------------------\n",
      "    G     -13041.432 eV\n",
      "=======================\n",
      "smi2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 50%|█████     | 2/4 [00:02<00:01,  1.05it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enthalpy components at T = 273.00 K:\n",
      "===============================\n",
      "E_pot              -6324.139 eV\n",
      "E_ZPE                  3.132 eV\n",
      "Cv_trans (0->T)        0.035 eV\n",
      "Cv_rot (0->T)          0.035 eV\n",
      "Cv_vib (0->T)          0.081 eV\n",
      "(C_v -> C_p)           0.024 eV\n",
      "-------------------------------\n",
      "H                  -6320.832 eV\n",
      "===============================\n",
      "Entropy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0016622 eV/K        0.454 eV\n",
      "S_rot              0.0011000 eV/K        0.300 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0005548 eV/K        0.151 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0033159 eV/K        0.905 eV\n",
      "=================================================\n",
      "Enthalpy components at T = 273.00 K:\n",
      "===============================\n",
      "E_pot              -6324.139 eV\n",
      "E_ZPE                  3.132 eV\n",
      "Cv_trans (0->T)        0.035 eV\n",
      "Cv_rot (0->T)          0.035 eV\n",
      "Cv_vib (0->T)          0.081 eV\n",
      "(C_v -> C_p)           0.024 eV\n",
      "-------------------------------\n",
      "H                  -6320.832 eV\n",
      "===============================\n",
      "\n",
      "Entropy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0016622 eV/K        0.454 eV\n",
      "S_rot              0.0011000 eV/K        0.300 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0005548 eV/K        0.151 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0033159 eV/K        0.905 eV\n",
      "=================================================\n",
      "\n",
      "Free energy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=======================\n",
      "    H      -6320.832 eV\n",
      " -T*S         -0.905 eV\n",
      "-----------------------\n",
      "    G      -6321.737 eV\n",
      "=======================\n",
      "smi3\n",
      "optiimize the input geometry\n",
      "      Step     Time          Energy         fmax\n",
      "BFGS:    0 15:17:15    -9406.938834        0.0107\n",
      "BFGS:    1 15:17:16    -9406.938833        0.0152\n",
      "BFGS:    2 15:17:16    -9406.938835        0.0024\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 75%|███████▌  | 3/4 [00:02<00:00,  1.12it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enthalpy components at T = 273.00 K:\n",
      "===============================\n",
      "E_pot              -9406.939 eV\n",
      "E_ZPE                  3.435 eV\n",
      "Cv_trans (0->T)        0.035 eV\n",
      "Cv_rot (0->T)          0.035 eV\n",
      "Cv_vib (0->T)          0.118 eV\n",
      "(C_v -> C_p)           0.024 eV\n",
      "-------------------------------\n",
      "H                  -9403.291 eV\n",
      "===============================\n",
      "Entropy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017046 eV/K        0.465 eV\n",
      "S_rot              0.0011862 eV/K        0.324 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0007638 eV/K        0.209 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0036535 eV/K        0.997 eV\n",
      "=================================================\n",
      "Enthalpy components at T = 273.00 K:\n",
      "===============================\n",
      "E_pot              -9406.939 eV\n",
      "E_ZPE                  3.435 eV\n",
      "Cv_trans (0->T)        0.035 eV\n",
      "Cv_rot (0->T)          0.035 eV\n",
      "Cv_vib (0->T)          0.118 eV\n",
      "(C_v -> C_p)           0.024 eV\n",
      "-------------------------------\n",
      "H                  -9403.291 eV\n",
      "===============================\n",
      "\n",
      "Entropy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017046 eV/K        0.465 eV\n",
      "S_rot              0.0011862 eV/K        0.324 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0007638 eV/K        0.209 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0036535 eV/K        0.997 eV\n",
      "=================================================\n",
      "\n",
      "Free energy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=======================\n",
      "    H      -9403.291 eV\n",
      " -T*S         -0.997 eV\n",
      "-----------------------\n",
      "    G      -9404.289 eV\n",
      "=======================\n",
      "smi4\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 4/4 [00:03<00:00,  1.06it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enthalpy components at T = 273.00 K:\n",
      "===============================\n",
      "E_pot              -9499.111 eV\n",
      "E_ZPE                  4.839 eV\n",
      "Cv_trans (0->T)        0.035 eV\n",
      "Cv_rot (0->T)          0.035 eV\n",
      "Cv_vib (0->T)          0.166 eV\n",
      "(C_v -> C_p)           0.024 eV\n",
      "-------------------------------\n",
      "H                  -9494.012 eV\n",
      "===============================\n",
      "Entropy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017193 eV/K        0.469 eV\n",
      "S_rot              0.0012450 eV/K        0.340 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0012634 eV/K        0.345 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0042266 eV/K        1.154 eV\n",
      "=================================================\n",
      "Enthalpy components at T = 273.00 K:\n",
      "===============================\n",
      "E_pot              -9499.111 eV\n",
      "E_ZPE                  4.839 eV\n",
      "Cv_trans (0->T)        0.035 eV\n",
      "Cv_rot (0->T)          0.035 eV\n",
      "Cv_vib (0->T)          0.166 eV\n",
      "(C_v -> C_p)           0.024 eV\n",
      "-------------------------------\n",
      "H                  -9494.012 eV\n",
      "===============================\n",
      "\n",
      "Entropy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=================================================\n",
      "                           S               T*S\n",
      "S_trans (1 bar)    0.0017193 eV/K        0.469 eV\n",
      "S_rot              0.0012450 eV/K        0.340 eV\n",
      "S_elec             0.0000000 eV/K        0.000 eV\n",
      "S_vib              0.0012634 eV/K        0.345 eV\n",
      "S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV\n",
      "-------------------------------------------------\n",
      "S                  0.0042266 eV/K        1.154 eV\n",
      "=================================================\n",
      "\n",
      "Free energy components at T = 273.00 K and P = 101325.0 Pa:\n",
      "=======================\n",
      "    H      -9494.012 eV\n",
      " -T*S         -1.154 eV\n",
      "-----------------------\n",
      "    G      -9495.166 eV\n",
      "=======================\n",
      "Number of failed thermo calculations:  0\n",
      "Number of successful thermo calculations:  4\n",
      "/home/jack/Auto3D_pkg/example/files/smiles_20240329-151158-762973/smiles_out_ANI2x_G.sdf\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "If the thermodynamic properties are NOT calculated in 298 K,\n",
    "you need to define a function that gets the unique ID and temperature for each molecules in the input path,\n",
    "then pass the custom function as the value of argument get_mol_idx_t.\n",
    "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.\n",
    "You can customize the function for other special needs.\n",
    "\"\"\"\n",
    "def custom_func(mol):\n",
    "    '''The mol object is an RDKit Molecule object.'''\n",
    "    id = mol.GetProp(\"_Name\")\n",
    "    t = 273\n",
    "    return (id, t)\n",
    "\n",
    "\n",
    "out_thermo = calc_thermo(out, \"ANI2x\", get_mol_idx_t=custom_func, opt_tol=0.003)\n",
    "print(out_thermo)  #enthalpy, entropy and Gibbs free energy are stored in the SDF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function calc_thermo in module Auto3D.ASE.thermo:\n",
      "\n",
      "calc_thermo(path: str, model_name: str, get_mol_idx_t=None, gpu_idx=0, opt_tol=0.0002, opt_steps=5000)\n",
      "    ASE interface for calculation thermo properties using ANI2x, ANI2xt or AIMNET.\n",
      "    \n",
      "    :param path: Input sdf file\n",
      "    :type path: str\n",
      "    :param model_name: ANI2x, ANI2xt or AIMNET\n",
      "    :type model_name: str\n",
      "    :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\n",
      "    :type get_mol_idx_t: function, optional\n",
      "    :param gpu_idx: GPU cuda index, defaults to 0\n",
      "    :type gpu_idx: int, optional\n",
      "    :param opt_tol: Convergence_threshold for geometry optimization, defaults to 0.0002\n",
      "    :type opt_tol: float, optional\n",
      "    :param opt_steps: Maximum geometry optimization steps, defaults to 5000\n",
      "    :type opt_steps: int, optional\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(calc_thermo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 ('py39')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.17"
  },
  "vscode": {
   "interpreter": {
    "hash": "61b4062b24dfb1010f420dad5aa3bd73a4d2af47d0ec44eafec465a35a9d7239"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}