Repository: isayevlab/Auto3D_pkg (language: python)

Issue:
pull from the fork



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 9741c81..c59c27a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,6 @@
+example/files/2022*
 .vscode/*
+*.sdf
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]
diff --git a/.gitmodules b/.gitmodules
deleted file mode 100644
index 4084f67..0000000
--- a/.gitmodules
+++ /dev/null
@@ -1,3 +0,0 @@
-[submodule "src/Auto3D"]
-	path = src/Auto3D
-	url = https://github.com/LiuCMU/Auto3D
diff --git a/README.md b/README.md
index 106dd34..0162994 100644
--- a/README.md
+++ b/README.md
@@ -1,47 +1,65 @@
-<!-- <img src="https://github.com/isayevlab/workflow/blob/master/auto3d_logo.png" width="300"> -->
-
 # **Auto3D**
+<a href="https://pypi.org/project/Auto3D/" target="_blank"><img src="https://img.shields.io/badge/pypi-url-informational" alt="pypi_link"></a>
+![PyPI](https://img.shields.io/pypi/v/Auto3D)
+![PyPI - Downloads](https://img.shields.io/pypi/dm/Auto3D)
+![PyPI - License](https://img.shields.io/pypi/l/Auto3D)
 
 # Introduction
-**Auto3D** automatically find the lowest-energy structures for the input SMILES. The user can get optimal 3D structures from plain SMIES files with a single line of command. All the processes, like isomer enumeration, duplicate and enantiomer filtering, optimization and ranking, are all taken care of by our package. The user can also try out different isomer enumeration programs and evaluation programs based on their demands.
+**Auto3D** automatically find the lowest-energy structures for the input SMILES. The user can get optimal 3D structures from plain SMIES files within 6 lines of code. All the processes, like isomer enumeration, duplicate and enantiomer filtering, optimization and ranking, are all taken care of by our package. The user can also try out different isomer enumeration programs and evaluation programs based on their demands.
 
-## Major Dependencies
-1. Python == 3.7
-2. One of the following packages for isomer enumeration step:
-- [RDKit](https://www.rdkit.org/docs/Install.html) >= 2020.09.01
-- [Omega](https://anaconda.org/openeye/openeye-toolkits) from [OpenEye Software](https://www.eyesopen.com/omega)
-3. The following cheminformatical tools for ranking and optimization step:
-- [AIMNET](https://github.com/aiqm/aimnet)
-- [TorchANI](https://github.com/roitberg-group/torchani_sandbox/tree/repulsion_calculator) with repulsion calculator fueatures
-- [OpenBabel](https://open-babel.readthedocs.io/en/latest/index.html)
-- [PyTorch](https://pytorch.org/get-started/locally/)
 
 # Installatioin
-1. Cloning the repository into a folder, all dependencies can be installed by entering the following command in your terminal:
+
+## Minimum Dependencies Installatioin
+1. Python >= 3.7
+2. [RDKit](https://www.rdkit.org/docs/Install.html) >= 2020.09.01 (For isomer engine)
+3. [OpenBabel](https://open-babel.readthedocs.io/en/latest/index.html) >= 3.1.1 (For molecular file processing)
+4. [PyTorch](https://pytorch.org/get-started/locally/) (For optimizing engine)
+
+If you have an environment with the above dependencies, Auto3D can be installed by
+```{bash}
+pip install Auto3D
+```
+Otherwise, you can create an environment and install Auto3D. In a terminal, excute the following code will create a environment named `auto3D` with the minimum dependencies installed, and install Auto3D.
 ```{bash}
+git clone https://github.com/isayevlab/Auto3D_pkg.git
+cd Auto3D_pkg
 conda env create --file auto3D.yml --name auto3D
 conda activate auto3D
+pip install Auto3D
 ```
-2. Install [TorchANI](https://github.com/roitberg-group/torchani_sandbox/tree/repulsion_calculator) with repulsion calculator fueatures.
-
-
-The [open eye toolkit](https://anaconda.org/openeye/openeye-toolkits) needs to be installed if you want to use Omega. Omega is a commercial software, so you will need a license to run it. However, you can still use all features about `Auto3D` with the free RDKit, because we implemented alternatives for all Omega functions that are used in `Auto3D`.
-
-# Usage
-A `.smi` file that stores your chemical structures is needed as the input for the program. You can find some example `.smi` files in the `input` folder. Basically, an `.smi ` file contains SMILES and their IDs.  Running the following command in the terminal will give you the 3-dimensional structures, which are stored in a file that has the same name as your input file, but is appended with `_3d.sdf`.
+## Optional Denpendencies Installation
+By installing Auto3D with the above minimum dependencies, you can use Auto3D with RDKit and [AIMNET](https://github.com/aiqm/aimnet) as the isomer engine and optimizing engine, respectively.
+Two additional optimizing engines are available: ANI-2x and ANI-2xt, which can be installed by:
 ```{bash}
-python auto3D.py input_SMILES_file_path --k=1
+conda activate auto3D
+conda install -c conda-forge torchani
 ```
-The above command runs the program and keeps 1 loweest-energy structure for each SMILES in your input file. If you want to keep n structures for each SMILES, simply set `--k=n `. You can also keep structures that are within 0.0048 Hatree compared with the lowest-energy structure for each SMILES if you run the following command:
+One additional isomer engine is availabel: OpenEye toolkit. It's a commercial software from [OpenEye Software](https://www.eyesopen.com/omega). It can be iinstalled by
 ```{bash}
-python workflow.py your_smiles_file_path --isomer_program=rdkit --window=0.0048
+conda activate auto3D
+conda install -c openeye openeye-toolkits
 ```
-
-The documentaion for different arguments are accessible via:
+To calculate thermodynamical properties (such as Gibbs free energy, enthalpy, entropy, geometry optimization) with Auto3D, [ASE](https://wiki.fysik.dtu.dk/ase/) needs to be installed:
 ```{bash}
-python workflow.py --h
+conda activate auto3D
+conda install -c conda-forge ase
 ```
 
+# Usage
+A `.smi` file that stores your chemical structures is needed as the input for the program. You can find some example `.smi` files in the `examplesinput` folder. Basically, an `.smi ` file contains SMILES and their IDs.  **ID can contain anything like numbers or letters, but not "_", the underscore.**
+Running the following command in the terminal will give you the 3-dimensional structures, which are stored in a file that has the same name as your input file, but is appended with `_3d.sdf`.
+```{Python}
+from Auto3D.auto3D import options, main
+
+if __name__ == "__main__":
+    path = "example/files/smiles.smi"
+    args = options(path, k=1)   #args specify the parameters for Auto3D 
+    out = main(args)            #main acceps the parameters and run Auto3D
+```
+The above command runs the Auto3D and keeps 1 loweest-energy structure for each SMILES in your input file. It uses RDKit as the isomer engine and AIMNET as the optimizng engine by default. `out` will be a path that stores the optimized 3D structures. If you want to keep n structures for each SMILES, simply set `k=n `. You can also keep structures that are within x kcal/mol compared with the lowest-energy structure for each SMILES if you replace `k=1` with `window=x`. More options are available by type `help(f)` where `f` is the function name.
+
+
 # Tunable parameters in Auto3D
 
 |State|Type|Name|Explanation|
@@ -49,6 +67,7 @@ python workflow.py --h
 |       |required argument|path   |a path of `.smi` file to store all SMILES and IDs|
 |ranking|required argument|--k    |Outputs the top-k structures for each SMILES. Only one of `--k` and `--window` need to be specified. |
 |ranking|required argument|--window|Outputs the structures whose energies are within a window (Hatree) from the lowest energy. Only one of `--k` and `--window` need to be specified. |
+|job segmentation|optional argument|--capacity|By default, 42. This is the number of SMILES that each small job will contain incase a large input file is given.|
 |isomer enumeration|optional argument|--enumerate_tautomer|By default, False. When True, enumerate tautomers for the input|
 |isomer enumeration|optional argument|--taut_program|Programs to enumerate tautomers, either 'rdkit' or 'oechem'. This argument only works when `--enumerate_tautomer=True`|
 |isomer enumeration|optional argument|--isomer_engine|By default, rdkit. The program for generating 3D isomers for each SMILES. This parameter is either rdkit or omega. RFKit is free for everyone, while Omega reuqires a license.))|
diff --git a/auto3D.yml b/auto3D.yml
new file mode 100644
index 0000000..8f1ebf3
--- /dev/null
+++ b/auto3D.yml
@@ -0,0 +1,15 @@
+name: auto3D
+channels:
+  - rmg
+  - conda-forge
+  - rdkit
+  - pytorch
+  - defaults
+dependencies:
+  - python>=3.7
+  - numpy
+  - pytorch
+  - rdkit[version='>=2021.09.02']
+  - openbabel[version='>=3.1.1']
+  - pandas
+prefix: /home/jack/miniconda3/envs/auto3D
diff --git a/example/files/smiles.smi b/example/files/smiles.smi
new file mode 100644
index 0000000..eba6749
--- /dev/null
+++ b/example/files/smiles.smi
@@ -0,0 +1,4 @@
+CNCC(C1=CC=CC=C1)=O smi1
+CC(CC)=O smi2
+COC(/C=C/C)=O smi3
+CC=CCC(CC)=O smi4
\ No newline at end of file
diff --git a/example/tutorial.ipynb b/example/tutorial.ipynb
new file mode 100644
index 0000000..abd6a8c
--- /dev/null
+++ b/example/tutorial.ipynb
@@ -0,0 +1,189 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import Auto3D\n",
+    "from Auto3D.auto3D import options, main"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "0.1.2\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(Auto3D.__version__)"
+   ]
+  },
+  {
+   "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 /Users/liu5/Documents/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, ANI_2xt and AIMNET.\n",
+      "There are 4 SMILES, available memory is 16.0 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",
+      "\n",
+      "\n",
+      "Optimizing on job1\n",
+      "Preparing for parallel optimizing... (Max optimization steps: 10000)\n",
+      "Total 3D conformers: 54\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      " 37%|███▋      | 3682/10000 [02:08<03:41, 28.56it/s] \n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Total 3D structures: 54  Converged: 38   Active: 16\n",
+      "Total 3D structures: 54  Converged: 53   Active: 1\n",
+      "Total 3D structures: 54  Converged: 53   Active: 1\n",
+      "All structures converged at step 3683:   Total 3D structures: 54  Converged: 54   Active: 0\n",
+      "Beggin to slelect structures that satisfy the requirements...\n",
+      "Energy unit: Hartree if implicit.\n",
+      "Program running time: 132.81213808059692 seconds\n",
+      "/Users/liu5/Documents/Auto3D_pkg/example/files/20220428-215702_smiles/smiles_out.sdf\n"
+     ]
+    }
+   ],
+   "source": [
+    "if __name__ == \"__main__\":\n",
+    "    path = \"/Users/liu5/Documents/Auto3D_pkg/example/files/smiles.smi\"\n",
+    "    args = options(path, k=1, 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": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Help on function options in module Auto3D.auto3D:\n",
+      "\n",
+      "options(path, k=False, window=False, verbose=False, job_name='', enumerate_tautomer=False, tauto_engine='rdkit', isomer_engine='rdkit', cis_trans=True, mode_oe='classic', mpi_np=4, max_confs=1000, use_gpu=True, gpu_idx=0, capacity=42, optimizing_engine='AIMNET', opt_steps=10000, convergence_threshold=0.003, threshold=0.3)\n",
+      "    Arguments for Auto3D main program\n",
+      "    path: A input.smi containing SMILES and IDs. Examples are listed in the input folder\n",
+      "    k: Outputs the top-k structures for each SMILES.\n",
+      "    window: Outputs the structures whose energies are within window (kcal/mol) from the lowest energy\n",
+      "    verbose: When True, save all meta data while running.\n",
+      "    job_name: A folder name to save all meta data.\n",
+      "    \n",
+      "    enumerate_tautomer: When True, enumerate tautomers for the input\n",
+      "    taut_engine: Programs to enumerate tautomers, either 'rdkit' or 'oechem\n",
+      "    isomer_engine: The program for generating 3D isomers for each SMILES. This parameter is either rdkit or omega.\n",
+      "    cis_trans: When True, cis/trans and r/s isomers are enumerated.\n",
+      "    mode_oe: \"The mode that omega program will take. \"It can be either 'classic', 'macrocycle' \"\"'rocs', 'pose' or 'dense'. \"\n",
+      "                            \"By default, the 'classic' mode is used. \"\n",
+      "                            \"For detailed information about each mode, see \"\n",
+      "                            \"https://docs.eyesopen.com/applications/omega/omega/omega_overview.html\"\n",
+      "    mpi_np: Number of CPU cores for the isomer generation step.\n",
+      "    max_confs: Maximum number of isomers for each SMILES\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "help(options)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Help on function main in module Auto3D.auto3D:\n",
+      "\n",
+      "main(args: dict)\n",
+      "    Take the arguments from options and run Auto3D\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "help(main)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
+  },
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.10.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/setup.cfg b/setup.cfg
index c85c089..892e711 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = Auto3D
-version = 0.0.4
+version = 0.1.2
 author = 'Zhen Liu, Tetiana Zubatiuk, Olexandr Isayev'
 author_email = liu5@andrew.cmu.edu, tetiana@zubatyuk.com, olexandr@olexandrisayev.com
 description = Generating the Low-energy 3D Conformers from SMILES with Isomer Enumeration and Blitz Geometry Optimization
@@ -20,6 +20,10 @@ package_dir =
 packages = find:
 python_requires = >=3.7
 include_package_data = True
+#...
+install_requires =
+    tqdm
+    psutil
 
 [options.packages.find]
 where = src
\ No newline at end of file
diff --git a/src/Auto3D b/src/Auto3D
deleted file mode 160000
index 09e68b6..0000000
--- a/src/Auto3D
+++ /dev/null
@@ -1,1 +0,0 @@
-Subproject commit 09e68b65f0dc51dbe134d2f06a6d94bc1517ded0
diff --git a/src/Auto3D/.gitignore b/src/Auto3D/.gitignore
new file mode 100644
index 0000000..c0df954
--- /dev/null
+++ b/src/Auto3D/.gitignore
@@ -0,0 +1,18 @@
+*.sdf
+oe_license.txt
+oeomega_*
+__MACOSX
+*.zip
+*.xyz
+flipper_*
+Emin
+sdf_files
+xyz_file
+*.tar.gz
+.ipynb_checkpoints
+__pycache__
+smiles_enumerated_hashed.smi
+2021*
+input/2022*
+*.txt
+.vscode/*
diff --git a/src/Auto3D/ASE/__init__.py b/src/Auto3D/ASE/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/Auto3D/ASE/geometry.py b/src/Auto3D/ASE/geometry.py
new file mode 100644
index 0000000..9f265db
--- /dev/null
+++ b/src/Auto3D/ASE/geometry.py
@@ -0,0 +1,241 @@
+#!/usr/bin/env python
+"""
+Geometry optimization with ANI2xt, AIMNET or ANI2x
+"""
+import sys
+import os
+root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+sys.path.append(root)
+
+import torch
+from ase import Atoms
+from ase.optimize import BFGS
+import ase.calculators.calculator
+try:
+    from ..batch_opt.ANI2xt import ANI2xt
+except:
+    pass
+import torchani
+from openbabel import pybel
+
+
+class EnForce_ANI(torch.nn.Module):
+    """Takes in an torch model, then defines forward functions for it.
+    Especially useful for AIMNET (class torch.jit)
+    Arguments:
+        name: ['ANI2xt', 'AIMNET']
+        model_parameters: path to the state dictionary
+    Returns:
+        the energies and forces for the input molecules.
+    """
+    def __init__(self, name, model_parameters=None, device=torch.device("cpu")):
+        super().__init__()
+        self.name = name
+        self.model_parameters = model_parameters
+        if self.name == 'ANI2xt':
+            model = ANI2xt(device, model_parameters)
+        elif self.name == "AIMNET":
+            model = torch.jit.load(model_parameters, map_location=device)
+        self.model = model
+        self.device = device
+
+    def forward(self, coord, numbers, charge=0):
+        """Calculate the energies and forces for input molecules. Called by self.forward_batched
+        
+        Arguments:
+            coord: coordinates for all input structures. size (B, N, 3), where
+                  B is the number of structures in coord, N is the number of
+                  atoms in each structure, 3 represents xyz dimensions.
+            numbers: the atomic numbers
+            
+        Returns:
+            energies
+            forces
+        """
+
+        if self.name == "AIMNET":
+            charge = torch.tensor(charge, dtype=torch.float, device=self.device)
+            d = self.model(dict(coord=coord, numbers=numbers, charge=charge))
+            e = (d['energy'] + d['disp_energy']).to(torch.double)
+            g = torch.autograd.grad([e.sum()], [coord])[0]
+            assert g is not None
+            f = -g
+        elif self.name == "ANI2xt":
+            # d = {1:0, 6:1, 7:2, 8:3, 16:4, 9:5, 17:6}
+            d = {1:0, 6:1, 7:2, 8:3, 9:4, 16:5, 17:6}
+            numbers2 = numbers.to('cpu').apply_(d.get).to(self.device)
+            e, f = self.model(numbers2, coord)
+
+        return e, f
+
+class Calculator(ase.calculators.calculator.Calculator):
+    """ASE calculator interface for AIMNET and ANI2xt"""
+    implemented_properties = ['energy', 'forces']
+    def __init__(self, model, charge=0):
+        super().__init__()
+        self.charge = charge
+        self.species = {'H':1, 'C':6, 'N':7, 'O':8, 'F':9, 'Si':14, 'P':15,
+                        'S':16, 'Cl':17, 'As':33, 'Se':34, 'Br':35, 'I':53,
+                        'B':5}
+        self.model = model 
+        for p in self.model.parameters():
+            p.requires_grad_(False)
+        a_parameter = next(self.model.parameters())
+        self.device = a_parameter.device
+        self.dtype = a_parameter.dtype
+
+    def calculate(self, atoms=None, properties=['energy'],
+                  system_changes=ase.calculators.calculator.all_changes):
+        super().calculate(atoms, properties, system_changes)
+
+        species = torch.tensor([self.species[symbol] for symbol in self.atoms.get_chemical_symbols()],
+                               dtype=torch.long, device=self.device)
+        coordinates = torch.tensor(self.atoms.get_positions()).to(self.device).to(self.dtype)
+        coordinates = coordinates.requires_grad_(True)
+
+        species = species.unsqueeze(0)
+        coordinates = coordinates.unsqueeze(0)
+        
+        energy, forces = self.model(coordinates, species, self.charge)
+        self.results['energy'] = energy.item()
+        self.results['forces'] = forces.squeeze(0).to('cpu').numpy()
+
+
+def opt_geometry(path: str, model_name:str, gpu_idx=0, opt_tol=0.003, opt_steps=10000):
+    """Geometry optimization interface with Auto3D
+    path: Input sdf file
+    model_name: ANI2xt or AIMNET
+    gpu_idx: GPU cuda index
+    opt_tol: Convergence_threshold for geometry optimization (eV/A)
+    opt_steps: Maximum geometry optimizaiton steps
+    """
+    ev2hatree = 1/27.211385
+    #create output path that is in the same directory as the input file
+    dir = os.path.dirname(path)
+    basename = os.path.basename(path).split(".")[0] + f"_{model_name}_opt.sdf"
+    outpath = os.path.join(dir, basename)
+    out_mols = []
+    # mols_failed = []
+    species2numbers = {'H':1, 'C':6, 'N':7, 'O':8, 'F':9, 'Si':14, 'P':15,
+                       'S':16, 'Cl':17, 'As':33, 'Se':34, 'Br':35, 'I':53, 'B':5}
+    numbers2species = dict([(val, key) for (key, val) in species2numbers.items()])
+
+    if torch.cuda.is_available():
+        device = torch.device(f"cuda:{gpu_idx}")
+    else:
+        device = torch.device("cpu")
+    if model_name == "ANI2xt":
+        dict_path = os.path.join(root, "models/ani2xt_seed0.pt")
+        model = EnForce_ANI('ANI2xt', dict_path, device=device)
+    elif model_name == "AIMNET":
+        dict_path = os.path.join(root, "models/aimnet2nqed_pc14iall_b97m_sae.jpt")
+        model = EnForce_ANI('AIMNET', dict_path, device=device)
+    elif model_name == "ANI2x":
+        calculator = torchani.models.ANI2x().to(device).ase()
+    else:
+        raise ValueError("model has to be 'ANI2x', 'ANI2xt' or 'AIMNET'")
+
+    mols = pybel.readfile("sdf", path)
+    for mol in mols:
+        coord = [a.coords for a in mol.atoms]
+        charge = mol.charge
+        species = [numbers2species[a.atomicnum] for a in mol.atoms]
+        atoms = Atoms(species, coord)
+        
+        if model_name != "ANI2x":
+            calculator = Calculator(model, charge)
+        atoms.set_calculator(calculator)
+
+        opt = BFGS(atoms)
+        opt.run(fmax=opt_tol, steps=opt_steps)
+        e = atoms.get_potential_energy()
+        mol.data['E_hatree'] = e * ev2hatree
+
+        #Updating ASE atoms coordinates into pybel mol
+        coord = atoms.get_positions()
+        for atom, c in zip(mol.atoms, coord):
+            atom.OBAtom.SetVector(*c)
+        out_mols.append(mol)
+
+    print(len(out_mols))
+    with open(outpath, 'w+') as f:
+        for mol in out_mols:
+            f.write(mol.write('sdf'))
+    return outpath
+
+
+if __name__ == "__main__":
+    import os
+    from openbabel import pybel
+    import argparse
+
+    
+    ev2hatree = 1/27.211385  #1 Hartree = 27.211385 eV
+    parser = argparse.ArgumentParser(description="ASE interface for calculation G for ANI2xt and AIMNET")
+    parser.add_argument("path", type=str, help="Input sdf file")
+    parser.add_argument("model", type=str, help="ANI2xt or AIMNET")
+    parser.add_argument("--gpu_idx", default=0, type=int, help="GPU cuda index")
+    parser.add_argument("--convergence_threshold", default=0.003, type=float, 
+                        help="Convergence_threshold for geometry optimization")
+    parser.add_argument("--opt_steps", default=10000, type=int,
+                        help="Maximum geometry optimizaiton steps")
+
+    args = parser.parse_args()
+    path = args.path
+    model_name = args.model
+    gpu_idx = args.gpu_idx
+    # idx_t = args.idx_t
+    opt_steps = args.opt_steps
+    opt_tol = args.convergence_threshold
+
+    dir = os.path.dirname(path)
+    basename = os.path.basename(path).split(".")[0] + f"_{model_name}_opt.sdf"
+    outpath = os.path.join(dir, basename)
+    out_mols = []
+    mols_failed = []
+    species2numbers = {'H':1, 'C':6, 'N':7, 'O':8, 'F':9, 'Si':14, 'P':15,
+                       'S':16, 'Cl':17, 'As':33, 'Se':34, 'Br':35, 'I':53, 'B':5}
+    numbers2species = dict([(val, key) for (key, val) in species2numbers.items()])
+
+    if torch.cuda.is_available():
+        device = torch.device(f"cuda:{gpu_idx}")
+    else:
+        device = torch.device("cpu")
+    if model_name == "ANI2xt":
+        dict_path = os.path.join(root, "models/ani2xt_seed0.pt")
+        model = EnForce_ANI('ANI2xt', dict_path, device=device)
+    elif model_name == "AIMNET":
+        dict_path = os.path.join(root, "models/aimnet2nqed_pc14iall_b97m_sae.jpt")
+        model = EnForce_ANI('AIMNET', dict_path, device=device)
+    elif model_name == "ANI2x":
+        calculator = torchani.models.ANI2x().to(device).ase()
+    else:
+        raise ValueError("model has to be 'ANI2x', 'ANI2xt' or 'AIMNET'")
+
+    mols = pybel.readfile("sdf", path)
+    for mol in mols:
+        coord = [a.coords for a in mol.atoms]
+        charge = mol.charge
+        species = [numbers2species[a.atomicnum] for a in mol.atoms]
+        atoms = Atoms(species, coord)
+        
+        if model_name != "ANI2x":
+            calculator = Calculator(model, charge)
+        atoms.set_calculator(calculator)
+
+
+        opt = BFGS(atoms)
+        opt.run(fmax=opt_tol, steps=opt_steps)
+        e = atoms.get_potential_energy()
+        mol.data['E_hatree'] = e * ev2hatree
+
+        #Updating ASE atoms coordinates into pybel mol
+        coord = atoms.get_positions()
+        for atom, c in zip(mol.atoms, coord):
+            atom.OBAtom.SetVector(*c)
+        out_mols.append(mol)
+
+    print(len(out_mols))
+    with open(outpath, 'w+') as f:
+        for mol in out_mols:
+            f.write(mol.write('sdf'))
diff --git a/src/Auto3D/ASE/thermo.py b/src/Auto3D/ASE/thermo.py
new file mode 100644
index 0000000..dc4ad3c
--- /dev/null
+++ b/src/Auto3D/ASE/thermo.py
@@ -0,0 +1,364 @@
+#!/usr/bin/env python
+"""
+Calculating thermodynamic perperties using Auto3D output
+"""
+import sys
+import os
+root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+sys.path.append(root)
+import warnings
+import torch
+from ase import Atoms
+from ase.optimize import BFGS
+from openbabel import pybel
+from ase.vibrations import Vibrations
+from ase.thermochemistry import IdealGasThermo
+import ase.calculators.calculator
+try:
+    from ..batch_opt.ANI2xt import ANI2xt
+except:
+    pass
+import torchani
+
+
+ev2hatree = 1/27.211385  #1 Hartree = 27.211385 eV
+class EnForce_ANI(torch.nn.Module):
+    """Takes in an torch model, then defines forward functions for it.
+    Especially useful for AIMNET (class torch.jit)
+    Arguments:
+        name: ['ANI2xt', 'AIMNET']
+        model_parameters: path to the state dictionary
+    Returns:
+        the energies and forces for the input molecules.
+    """
+    def __init__(self, name, model_parameters=None, device=torch.device("cpu")):
+        super().__init__()
+        self.name = name
+        self.model_parameters = model_parameters
+        if self.name == 'ANI2xt':
+            model = ANI2xt(device, model_parameters)
+        elif self.name == "AIMNET":
+            model = torch.jit.load(model_parameters, map_location=device)
+        self.model = model
+        self.device = device
+
+    def forward(self, coord, numbers, charge=0):
+        """Calculate the energies and forces for input molecules. Called by self.forward_batched
+        
+        Arguments:
+            coord: coordinates for all input structures. size (B, N, 3), where
+                  B is the number of structures in coord, N is the number of
+                  atoms in each structure, 3 represents xyz dimensions.
+            numbers: the atomic numbers
+            
+        Returns:
+            energies
+            forces
+        """
+
+        if self.name == "AIMNET":
+            charge = torch.tensor(charge, dtype=torch.float, device=self.device)
+            d = self.model(dict(coord=coord, numbers=numbers, charge=charge))
+            e = (d['energy'] + d['disp_energy']).to(torch.double)
+            g = torch.autograd.grad([e.sum()], [coord])[0]
+            assert g is not None
+            f = -g
+        elif self.name == "ANI2xt":
+            # d = {1:0, 6:1, 7:2, 8:3, 16:4, 9:5, 17:6}
+            d = {1:0, 6:1, 7:2, 8:3, 9:4, 16:5, 17:6}
+            numbers2 = numbers.to('cpu').apply_(d.get).to(self.device)
+            e, f = self.model(numbers2, coord)
+
+        return e, f
+
+class Calculator(ase.calculators.calculator.Calculator):
+    """ASE calculator interface for AIMNET and ANI2xt"""
+    implemented_properties = ['energy', 'forces']
+    def __init__(self, model, charge=0):
+        super().__init__()
+        self.charge = charge
+        self.species = {'H':1, 'C':6, 'N':7, 'O':8, 'F':9, 'Si':14, 'P':15,
+                        'S':16, 'Cl':17, 'As':33, 'Se':34, 'Br':35, 'I':53,
+                        'B':5}
+        self.model = model 
+        for p in self.model.parameters():
+            p.requires_grad_(False)
+        a_parameter = next(self.model.parameters())
+        self.device = a_parameter.device
+        self.dtype = a_parameter.dtype
+
+    def calculate(self, atoms=None, properties=['energy'],
+                  system_changes=ase.calculators.calculator.all_changes):
+        super().calculate(atoms, properties, system_changes)
+
+        species = torch.tensor([self.species[symbol] for symbol in self.atoms.get_chemical_symbols()],
+                               dtype=torch.long, device=self.device)
+        coordinates = torch.tensor(self.atoms.get_positions()).to(self.device).to(self.dtype)
+        coordinates = coordinates.requires_grad_(True)
+
+        species = species.unsqueeze(0)
+        coordinates = coordinates.unsqueeze(0)
+        
+        energy, forces = self.model(coordinates, species, self.charge)
+        self.results['energy'] = energy.item()
+        self.results['forces'] = forces.squeeze(0).to('cpu').numpy()
+
+
+def get_mol_idx_t1(mol):
+    """Get idx and temperature from openbabel molecule
+    example: /Users/liu5/Documents/tautomer/dG/dG3/DFT_output_analyze2/817-2-473.xyz
+    """
+
+    idx = str(mol).split()[1].strip().split("/")[-1].strip().split(".")[0]
+    T = int(idx.split("-")[-1])
+    return (idx, T)
+
+def get_mol_idx_t2(mol):
+    """Get idx and temperature from openbabel molecule"""
+    idx = mol.data['ID']
+    ref_idx_t = idx.split("_")[0].strip()
+    T = round(float(ref_idx_t.split("-")[-1]))
+    return (idx, T)
+
+def get_mol_idx_t3(mol):
+    "Setting default index and temperature"
+    idx = ""
+    T = 298
+    return (idx, T)
+
+def get_mol_idx_t4(mol):
+    idx = mol.title.strip()
+    T = 298
+    return (idx, T)
+
+def get_mol_idx_t5(mol):
+    idx = mol.data["ID"].strip()
+    T = 298
+    return (idx, T)
+
+def get_mol_idx_t6(mol):
+    idx = mol.data["ID"].strip()
+    T = float(mol.data["T"])
+    return (idx, T)
+
+
+def calc_thermo(path: str, model_name: str, get_mol_idx_t:callable, gpu_idx=0, opt_tol=0.005, opt_steps=10000):
+    """ASE interface for calculation thermo properties using ANI2xt and AIMNET
+    path: Input sdf file
+    model_name: ANI2x, ANI2xt or AIMNET
+    get_mol_idx_t: a functioin that returns (idx, T) from a pybel mol object
+    gpu_idx: GPU cuda index
+    opt_tol: Convergence_threshold for geometry optimization
+    opt_steps: Maximum geometry optimizaiton steps"""
+    #Prepare output name
+    dir = os.path.dirname(path)
+    basename = os.path.basename(path).split(".")[0] + f"_{model_name}_G.sdf"
+    outpath = os.path.join(dir, basename)
+
+    #internal parameters
+    out_mols = []
+    mols_failed = []
+    species2numbers = {'H':1, 'C':6, 'N':7, 'O':8, 'F':9, 'Si':14, 'P':15,
+                       'S':16, 'Cl':17, 'As':33, 'Se':34, 'Br':35, 'I':53, 'B':5}
+    numbers2species = dict([(val, key) for (key, val) in species2numbers.items()])
+    if torch.cuda.is_available():
+        device = torch.device(f"cuda:{gpu_idx}")
+    else:
+        device = torch.device("cpu")
+    if model_name == "ANI2xt":
+        dict_path = os.path.join(root, "models/ani2xt_seed0.pt")
+        model = EnForce_ANI('ANI2xt', dict_path, device=device)
+    elif model_name == "AIMNET":
+        dict_path = os.path.join(root, "models/aimnet2nqed_pc14iall_b97m_sae.jpt")
+        model = EnForce_ANI('AIMNET', dict_path, device=device)
+    elif model_name == "ANI2x":
+        calculator = torchani.models.ANI2x().to(device).ase()
+    else:
+        raise ValueError("model has to be 'ANI2x', 'ANI2xt' or 'AIMNET'")
+
+    mols = pybel.readfile("sdf", path)
+    for mol in mols:
+        coord = [a.coords for a in mol.atoms]
+        charge = mol.charge
+        species = [numbers2species[a.atomicnum] for a in mol.atoms]
+        atoms = Atoms(species, coord)
+        
+        if model_name != "ANI2x":
+            calculator = Calculator(model, charge)
+        atoms.set_calculator(calculator)
+
+
+        opt = BFGS(atoms)
+        opt.run(fmax=opt_tol, steps=opt_steps)
+        e = atoms.get_potential_energy()
+
+        idx, T = get_mol_idx_t(mol)
+
+        vib = Vibrations(atoms)
+        try:
+            vib.clean()
+            vib.run()
+            vib_e = vib.get_energies()
+
+            thermo = IdealGasThermo(vib_energies=vib_e,
+                                    potentialenergy=e,
+                                    atoms=atoms,
+                                    geometry='nonlinear',
+                                    symmetrynumber=1, spin=0)
+            H = thermo.get_enthalpy(temperature=T) * ev2hatree
+            S = thermo.get_entropy(temperature=T, pressure=101325) * ev2hatree
+            G = thermo.get_gibbs_energy(temperature=T, pressure=101325) * ev2hatree
+            vib.clean()
+            print(G)
+
+            mol.data['H_hatree'] = H
+            mol.data['S_hatree'] = S
+            mol.data['T_K'] = T
+            mol.data['G_hatree'] = G
+            mol.data['E_hatree'] = e * ev2hatree
+            out_mols.append(mol)
+
+            #Updating ASE atoms coordinates into pybel mol
+            coord = atoms.get_positions()
+            for atom, c in zip(mol.atoms, coord):
+                atom.OBAtom.SetVector(*c)
+        
+        except:
+            print("Failed: ", idx)
+            vib.clean()
+            mols_failed.append(mol)
+
+    print(len(mols_failed))
+    print(len(out_mols))
+    with open(outpath, 'w+') as f:
+        all_mols = out_mols + mols_failed
+        for mol in all_mols:
+            f.write(mol.write('sdf'))
+    return outpath
+
+
+if __name__ == "__main__":
+    import os
+    from openbabel import pybel
+    from ase.vibrations import Vibrations
+    from ase.thermochemistry import IdealGasThermo
+    import argparse
+
+    
+    ev2hatree = 1/27.211385  #1 Hartree = 27.211385 eV
+    parser = argparse.ArgumentParser(description="ASE interface for calculation G for ANI2xt and AIMNET")
+    parser.add_argument("path", type=str, help="Input sdf file")
+    parser.add_argument("model", type=str, help="ANI2x, ANI2xt or AIMNET")
+    parser.add_argument("idx_t", type=int, help="Need to define a function to get index and T from openbabel molecule")
+    parser.add_argument("--gpu_idx", default=0, type=int, help="GPU cuda index")
+    parser.add_argument("--convergence_threshold", default=0.0005, type=float, 
+                        help="Convergence_threshold for geometry optimization")
+    parser.add_argument("--opt_steps", default=10000, type=int,
+                        help="Maximum geometry optimizaiton steps")
+
+    args = parser.parse_args()
+    path = args.path
+    model_name = args.model
+    gpu_idx = args.gpu_idx
+    idx_t = args.idx_t
+    opt_steps = args.opt_steps
+    opt_tol = args.convergence_threshold
+
+    dir = os.path.dirname(path)
+    basename = os.path.basename(path).split(".")[0] + f"_{model_name}_G.sdf"
+    outpath = os.path.join(dir, basename)
+    out_mols = []
+    mols_failed = []
+    species2numbers = {'H':1, 'C':6, 'N':7, 'O':8, 'F':9, 'Si':14, 'P':15,
+                       'S':16, 'Cl':17, 'As':33, 'Se':34, 'Br':35, 'I':53, 'B':5}
+    numbers2species = dict([(val, key) for (key, val) in species2numbers.items()])
+
+    if torch.cuda.is_available():
+        device = torch.device(f"cuda:{gpu_idx}")
+    else:
+        device = torch.device("cpu")
+    if model_name == "ANI2xt":
+        dict_path = os.path.join(root, "models/ani2xt_seed0.pt")
+        model = EnForce_ANI('ANI2xt', dict_path, device=device)
+    elif model_name == "AIMNET":
+        dict_path = os.path.join(root, "models/aimnet2nqed_pc14iall_b97m_sae.jpt")
+        model = EnForce_ANI('AIMNET', dict_path, device=device)
+    elif model_name == "ANI2x":
+        calculator = torchani.models.ANI2x().to(device).ase()
+    else:
+        raise ValueError("model has to be 'ANI2x', 'ANI2xt' or 'AIMNET'")
+
+
+    if int(idx_t) == 1:
+        get_mol_idx_t = get_mol_idx_t1
+    elif int(idx_t) == 2:
+        get_mol_idx_t = get_mol_idx_t2
+    elif int(idx_t) == 3:
+        get_mol_idx_t = get_mol_idx_t3
+    elif int(idx_t) == 4:
+        get_mol_idx_t = get_mol_idx_t4
+    elif int(idx_t) == 5:
+        get_mol_idx_t = get_mol_idx_t5
+    elif int(idx_t) == 6:
+        get_mol_idx_t = get_mol_idx_t6
+    else:
+        raise ValueError("Has to define a function for get index and temperature for OB mol object.")
+
+    mols = pybel.readfile("sdf", path)
+    for mol in mols:
+        coord = [a.coords for a in mol.atoms]
+        charge = mol.charge
+        species = [numbers2species[a.atomicnum] for a in mol.atoms]
+        atoms = Atoms(species, coord)
+        
+        if model_name != "ANI2x":
+            calculator = Calculator(model, charge)
+        atoms.set_calculator(calculator)
+
+
+        opt = BFGS(atoms)
+        opt.run(fmax=opt_tol, steps=opt_steps)
+        e = atoms.get_potential_energy()
+
+        idx, T = get_mol_idx_t(mol)
+
+        vib = Vibrations(atoms)
+        try:
+            vib.clean()
+            vib.run()
+            vib_e = vib.get_energies()
+
+            thermo = IdealGasThermo(vib_energies=vib_e,
+                                    potentialenergy=e,
+                                    atoms=atoms,
+                                    geometry='nonlinear',
+                                    symmetrynumber=1, spin=0)
+            H = thermo.get_enthalpy(temperature=T) * ev2hatree
+            S = thermo.get_entropy(temperature=T, pressure=101325) * ev2hatree
+            G = thermo.get_gibbs_energy(temperature=T, pressure=101325) * ev2hatree
+            vib.clean()
+            print(G)
+
+            mol.data['H_hatree'] = H
+            mol.data['S_hatree'] = S
+            mol.data['T_K'] = T
+            mol.data['G_hatree'] = G
+            mol.data['E_hatree'] = e * ev2hatree
+            out_mols.append(mol)
+
+            #Updating ASE atoms coordinates into pybel mol
+            coord = atoms.get_positions()
+            for atom, c in zip(mol.atoms, coord):
+                atom.OBAtom.SetVector(*c)
+        
+        except:
+            print("Failed: ", idx)
+            vib.clean()
+            mols_failed.append(mol)
+
+    print(len(mols_failed))
+    print(len(out_mols))
+    with open(outpath, 'w+') as f:
+        all_mols = out_mols + mols_failed
+        for mol in all_mols:
+            f.write(mol.write('sdf'))
diff --git a/src/Auto3D/LICENSE b/src/Auto3D/LICENSE
new file mode 100644
index 0000000..15f0186
--- /dev/null
+++ b/src/Auto3D/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021 isayevlab
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/src/Auto3D/README.md b/src/Auto3D/README.md
new file mode 100644
index 0000000..106dd34
--- /dev/null
+++ b/src/Auto3D/README.md
@@ -0,0 +1,66 @@
+<!-- <img src="https://github.com/isayevlab/workflow/blob/master/auto3d_logo.png" width="300"> -->
+
+# **Auto3D**
+
+# Introduction
+**Auto3D** automatically find the lowest-energy structures for the input SMILES. The user can get optimal 3D structures from plain SMIES files with a single line of command. All the processes, like isomer enumeration, duplicate and enantiomer filtering, optimization and ranking, are all taken care of by our package. The user can also try out different isomer enumeration programs and evaluation programs based on their demands.
+
+## Major Dependencies
+1. Python == 3.7
+2. One of the following packages for isomer enumeration step:
+- [RDKit](https://www.rdkit.org/docs/Install.html) >= 2020.09.01
+- [Omega](https://anaconda.org/openeye/openeye-toolkits) from [OpenEye Software](https://www.eyesopen.com/omega)
+3. The following cheminformatical tools for ranking and optimization step:
+- [AIMNET](https://github.com/aiqm/aimnet)
+- [TorchANI](https://github.com/roitberg-group/torchani_sandbox/tree/repulsion_calculator) with repulsion calculator fueatures
+- [OpenBabel](https://open-babel.readthedocs.io/en/latest/index.html)
+- [PyTorch](https://pytorch.org/get-started/locally/)
+
+# Installatioin
+1. Cloning the repository into a folder, all dependencies can be installed by entering the following command in your terminal:
+```{bash}
+conda env create --file auto3D.yml --name auto3D
+conda activate auto3D
+```
+2. Install [TorchANI](https://github.com/roitberg-group/torchani_sandbox/tree/repulsion_calculator) with repulsion calculator fueatures.
+
+
+The [open eye toolkit](https://anaconda.org/openeye/openeye-toolkits) needs to be installed if you want to use Omega. Omega is a commercial software, so you will need a license to run it. However, you can still use all features about `Auto3D` with the free RDKit, because we implemented alternatives for all Omega functions that are used in `Auto3D`.
+
+# Usage
+A `.smi` file that stores your chemical structures is needed as the input for the program. You can find some example `.smi` files in the `input` folder. Basically, an `.smi ` file contains SMILES and their IDs.  Running the following command in the terminal will give you the 3-dimensional structures, which are stored in a file that has the same name as your input file, but is appended with `_3d.sdf`.
+```{bash}
+python auto3D.py input_SMILES_file_path --k=1
+```
+The above command runs the program and keeps 1 loweest-energy structure for each SMILES in your input file. If you want to keep n structures for each SMILES, simply set `--k=n `. You can also keep structures that are within 0.0048 Hatree compared with the lowest-energy structure for each SMILES if you run the following command:
+```{bash}
+python workflow.py your_smiles_file_path --isomer_program=rdkit --window=0.0048
+```
+
+The documentaion for different arguments are accessible via:
+```{bash}
+python workflow.py --h
+```
+
+# Tunable parameters in Auto3D
+
+|State|Type|Name|Explanation|
+|---|---|---|---|
+|       |required argument|path   |a path of `.smi` file to store all SMILES and IDs|
+|ranking|required argument|--k    |Outputs the top-k structures for each SMILES. Only one of `--k` and `--window` need to be specified. |
+|ranking|required argument|--window|Outputs the structures whose energies are within a window (Hatree) from the lowest energy. Only one of `--k` and `--window` need to be specified. |
+|isomer enumeration|optional argument|--enumerate_tautomer|By default, False. When True, enumerate tautomers for the input|
+|isomer enumeration|optional argument|--taut_program|Programs to enumerate tautomers, either 'rdkit' or 'oechem'. This argument only works when `--enumerate_tautomer=True`|
+|isomer enumeration|optional argument|--isomer_engine|By default, rdkit. The program for generating 3D isomers for each SMILES. This parameter is either rdkit or omega. RFKit is free for everyone, while Omega reuqires a license.))|
+|isomer enumeration|optional argument|--max_confs|Maximum number of isomers for each configuration of the SMILES|
+|isomer enumeration|optional argument|--cis_trans|By default, True. When True, cis/trans and r/s isomers are enumerated|
+|isomer enumeration|optional argument|--mode_oe|By default, classic. "The mode that omega program will take. It can be either 'classic' or 'macrocycle'. Only works when `--isomer_engine=omega`|
+|isomer enumeration|optional argument|--mpi_np|Number of CPU cores for the isomer generation step. Only works when `--isomer_engine=omega`|
+|optimization|optional argument|--optimizing_engine|By default, AIMNET. Choose either 'ANI2x', 'ANI2xt', or 'AIMNET' for energy calculation and geometry optimization.|
+|optimization|optional argument|--use_gpu|By deafult, True. If True, the program will use GPU|
+|optimization|optional argument|--gpu_idx|"GPU index. It only works when --use_gpu=True|
+|optimization|optional argument|--opt_steps|By deafult, 10000. Maximum optimization steps for each structure|
+|optimization|optional argument|--convergence_threshold|By deafult, 0.003. Optimization is considered as converged if maximum force is below this threshold|
+|duplicate removing|optional argument|--threshold|By default, 0.3. If the RMSD between two conformers are within threhold, they are considered as duplicates. One of them will be removed. Duplicate removing are excuted after conformer enumeration and geometry optimization|
+|  housekeeping     |optional argument| --verbose |By default, True. When True, save all meta data while running|
+
diff --git a/src/Auto3D/SPE.py b/src/Auto3D/SPE.py
new file mode 100644
index 0000000..8b67484
--- /dev/null
+++ b/src/Auto3D/SPE.py
@@ -0,0 +1,214 @@
+#!/usr/bin/env python
+"""Calculating single point energy using ANI2xt, ANI2x or AIMNET"""
+import sys
+import os
+root = os.path.dirname(os.path.abspath(__file__))
+sys.path.append(root)
+
+import torch
+import warnings
+from ase import Atoms
+import ase.calculators.calculator
+try:
+    from .batch_opt.ANI2xt import ANI2xt
+except:
+    pass
+import torchani
+from openbabel import pybel
+from tqdm import tqdm
+
+
+ev2hatree = 1/27.211385
+class EnForce_ANI(torch.nn.Module):
+    """Takes in an torch model, then defines forward functions for it.
+    Especially useful for AIMNET (class torch.jit)
+    Arguments:
+        name: ['ANI2xt', 'AIMNET']
+        model_parameters: path to the state dictionary
+    Returns:
+        the energies and forces for the input molecules.
+    """
+    def __init__(self, name, model_parameters=None, device=torch.device("cpu")):
+        super().__init__()
+        self.name = name
+        self.model_parameters = model_parameters
+        if self.name == 'ANI2xt':
+            model = ANI2xt(device, model_parameters)
+        elif self.name == "AIMNET":
+            model = torch.jit.load(model_parameters, map_location=device)
+        self.model = model
+        self.device = device
+
+    def forward(self, coord, numbers, charge=0):
+        """Calculate the energies and forces for input molecules. Called by self.forward_batched
+        
+        Arguments:
+            coord: coordinates for all input structures. size (B, N, 3), where
+                  B is the number of structures in coord, N is the number of
+                  atoms in each structure, 3 represents xyz dimensions.
+            numbers: the atomic numbers
+            
+        Returns:
+            energies
+            forces
+        """
+
+        if self.name == "AIMNET":
+            charge = torch.tensor(charge, dtype=torch.float, device=self.device)
+            d = self.model(dict(coord=coord, numbers=numbers, charge=charge))
+            e = (d['energy'] + d['disp_energy']).to(torch.double)
+            g = torch.autograd.grad([e.sum()], [coord])[0]
+            assert g is not None
+            f = -g
+        elif self.name == "ANI2xt":
+            # d = {1:0, 6:1, 7:2, 8:3, 16:4, 9:5, 17:6}
+            d = {1:0, 6:1, 7:2, 8:3, 9:4, 16:5, 17:6}
+            numbers2 = numbers.to('cpu').apply_(d.get).to(self.device)
+            e, f = self.model(numbers2, coord)
+
+        return e, f
+
+class Calculator(ase.calculators.calculator.Calculator):
+    """ASE calculator interface for AIMNET and ANI2xt"""
+    implemented_properties = ['energy', 'forces']
+    def __init__(self, model, charge=0):
+        super().__init__()
+        self.charge = charge
+        self.species = {'H':1, 'C':6, 'N':7, 'O':8, 'F':9, 'Si':14, 'P':15,
+                        'S':16, 'Cl':17, 'As':33, 'Se':34, 'Br':35, 'I':53,
+                        'B':5}
+        self.model = model 
+        for p in self.model.parameters():
+            p.requires_grad_(False)
+        a_parameter = next(self.model.parameters())
+        self.device = a_parameter.device
+        self.dtype = a_parameter.dtype
+
+    def calculate(self, atoms=None, properties=['energy'],
+                  system_changes=ase.calculators.calculator.all_changes):
+        super().calculate(atoms, properties, system_changes)
+
+        species = torch.tensor([self.species[symbol] for symbol in self.atoms.get_chemical_symbols()],
+                               dtype=torch.long, device=self.device)
+        coordinates = torch.tensor(self.atoms.get_positions()).to(self.device).to(self.dtype)
+        coordinates = coordinates.requires_grad_(True)
+
+        species = species.unsqueeze(0)
+        coordinates = coordinates.unsqueeze(0)
+        
+        energy, forces = self.model(coordinates, species, self.charge)
+        self.results['energy'] = energy.item()
+        self.results['forces'] = forces.squeeze(0).to('cpu').numpy()
+
+
+def calc_spe(path:str, model_name:str, gpu_idx=0):
+    """Calculating single point energy
+    path: Input sdf file
+    model_name: ANI2x, ANI2xt or AIMNET
+    gpu_idx: GPU cuda index"""
+    #Create a output path that is the in the same directory as the input
+    dir = os.path.dirname(path)
+    basename = os.path.basename(path).split(".")[0] + f"_{model_name}_E.sdf"
+    outpath = os.path.join(dir, basename)
+
+    out_mols = []
+    species2numbers = {'H':1, 'C':6, 'N':7, 'O':8, 'F':9, 'Si':14, 'P':15,
+                       'S':16, 'Cl':17, 'As':33, 'Se':34, 'Br':35, 'I':53, 'B':5}
+    numbers2species = dict([(val, key) for (key, val) in species2numbers.items()])
+
+    if torch.cuda.is_available():
+        device = torch.device(f"cuda:{gpu_idx}")
+    else:
+        device = torch.device("cpu")
+
+    if model_name == "ANI2xt":
+        dict_path = os.path.join(root, "models/ani2xt_seed0.pt")
+        model = EnForce_ANI('ANI2xt', dict_path, device=device)
+    elif model_name == "AIMNET":
+        dict_path = os.path.join(root, "models/aimnet2nqed_pc14iall_b97m_sae.jpt")
+        model = EnForce_ANI('AIMNET', dict_path, device=device)
+    elif model_name == "ANI2x":
+        calculator = torchani.models.ANI2x().to(device).ase()
+    else:
+        raise ValueError("model has to be 'ANI2x', 'ANI2xt' or 'AIMNET'")
+
+    mols = pybel.readfile("sdf", path)
+    for mol in tqdm(mols):
+        coord = [a.coords for a in mol.atoms]
+        charge = mol.charge
+        species = [numbers2species[a.atomicnum] for a in mol.atoms]
+        atoms = Atoms(species, coord)
+        
+        if model_name != "ANI2x":
+            calculator = Calculator(model, charge)
+        atoms.set_calculator(calculator)
+
+        e = atoms.get_potential_energy()
+        mol.data['E_hatree'] = e * ev2hatree
+        out_mols.append(mol)
+
+    with open(outpath, 'w+') as f:
+        for mol in out_mols:
+            f.write(mol.write('sdf'))
+    return outpath
+
+
+if __name__ == "__main__":
+    import os
+    from openbabel import pybel
+    import argparse
+
+    
+    ev2hatree = 1/27.211385  #1 Hartree = 27.211385 eV
+    parser = argparse.ArgumentParser(description="ASE interface for calculation G for ANI2xt and AIMNET")
+    parser.add_argument("path", type=str, help="Input sdf file")
+    parser.add_argument("model", type=str, help="ANI2x, ANI2xt or AIMNET")
+    parser.add_argument("--gpu_idx", default=0, type=int, help="GPU cuda index")
+
+    args = parser.parse_args()
+    path = args.path
+    model_name = args.model
+    gpu_idx = args.gpu_idx
+
+    dir = os.path.dirname(path)
+    basename = os.path.basename(path).split(".")[0] + f"_{model_name}_E.sdf"
+    outpath = os.path.join(dir, basename)
+    out_mols = []
+    species2numbers = {'H':1, 'C':6, 'N':7, 'O':8, 'F':9, 'Si':14, 'P':15,
+                       'S':16, 'Cl':17, 'As':33, 'Se':34, 'Br':35, 'I':53, 'B':5}
+    numbers2species = dict([(val, key) for (key, val) in species2numbers.items()])
+
+    if torch.cuda.is_available():
+        device = torch.device(f"cuda:{gpu_idx}")
+    else:
+        device = torch.device("cpu")
+
+    if model_name == "ANI2xt":
+        dict_path = os.path.join(root, "models/ani2xt_seed0.pt")
+        model = EnForce_ANI('ANI2xt', dict_path, device=device)
+    elif model_name == "AIMNET":
+        dict_path = os.path.join(root, "models/aimnet2nqed_pc14iall_b97m_sae.jpt")
+        model = EnForce_ANI('AIMNET', dict_path, device=device)
+    elif model_name == "ANI2x":
+        calculator = torchani.models.ANI2x().to(device).ase()
+    else:
+        raise ValueError("model has to be 'ANI2x', 'ANI2xt' or 'AIMNET'")
+
+    mols = pybel.readfile("sdf", path)
+    for mol in tqdm(mols):
+        coord = [a.coords for a in mol.atoms]
+        charge = mol.charge
+        species = [numbers2species[a.atomicnum] for a in mol.atoms]
+        atoms = Atoms(species, coord)
+        
+        if model_name != "ANI2x":
+            calculator = Calculator(model, charge)
+        atoms.set_calculator(calculator)
+
+        e = atoms.get_potential_energy()
+        mol.data['E_hatree'] = e * ev2hatree
+        out_mols.append(mol)
+
+    with open(outpath, 'w+') as f:
+        for mol in out_mols:
+            f.write(mol.write('sdf'))
diff --git a/src/Auto3D/__init__.py b/src/Auto3D/__init__.py
new file mode 100644
index 0000000..2c70b78
--- /dev/null
+++ b/src/Auto3D/__init__.py
@@ -0,0 +1,26 @@
+import warnings
+from pkg_resources import get_distribution, DistributionNotFound
+
+try:
+    __version__ = get_distribution(__name__).version
+except DistributionNotFound:
+    pass
+
+with warnings.catch_warnings():
+    warnings.simplefilter("ignore")
+    try:
+        from openeye import oechem
+        from openeye import oequacpac
+        from openeye import oeomega
+    except:
+        warnings.warn("Openeye software is not available.")
+
+    try:
+        import torchani
+    except:
+        warnings.warn("TorchANI is not installed")
+
+    try:
+        from .batch_opt.ANI2xt import ANI2xt
+    except:
+        warnings.warn("ANI2xt model is not available")
diff --git a/src/Auto3D/auto3D.py b/src/Auto3D/auto3D.py
new file mode 100644
index 0000000..1c38898
--- /dev/null
+++ b/src/Auto3D/auto3D.py
@@ -0,0 +1,456 @@
+#!/usr/bin/env python
+"""
+Generating 3-dimensional structures from SMILES based on user's demands.
+"""
+
+import argparse
+import os
+import shutil
+import sys
+import time
+import torch
+import math
+import psutil, tarfile
+import glob
+import pandas as pd
+import multiprocessing as mp
+from Auto3D.isomer_engine import rd_isomer, tautomer_engine
+from Auto3D.isomer_engine import oe_isomer
+from Auto3D.ranking import ranking
+from Auto3D.utils import housekeeping, check_input
+from Auto3D.utils import hash_taut_smi,  my_name_space
+from Auto3D.batch_opt.batchopt import optimizing
+from send2trash import send2trash
+try:
+    mp.set_start_method('spawn')
+except:
+    pass
+
+
+def create_chunk_meta_names(path, dir):
+    """Output name is based on chunk input path and directory
+    path: chunck input smi path
+    dir: chunck job folder
+    """
+    dct = {}
+    output_name = os.path.basename(path).split('.')[0].strip() + '_3d.sdf'
+    output = os.path.join(dir, output_name)
+    optimized_og = output.split('.')[0] + '0.sdf'
+
+    output_taut = os.path.join(dir, 'smi_taut.smi')
+    smiles_enumerated = os.path.join(dir, 'smiles_enumerated.smi')
+    smiles_reduced = smiles_enumerated.split('.')[0] + '_reduced.smi'
+    smiles_hashed = os.path.join(dir, 'smiles_enumerated_hashed.smi')
+    enumerated_sdf = os.path.join(dir, 'smiles_enumerated.sdf')
+    sorted_sdf = os.path.join(dir, 'enumerated_sorted.sdf')
+    housekeeping_folder = os.path.join(dir, 'verbose')
+    # dct["output_name"] = output_name
+    dct["output"] = output
+    dct["optimized_og"] = optimized_og
+    dct["output_taut"] = output_taut
+    dct["smiles_enumerated"] = smiles_enumerated
+    dct["smiles_reduced"] = smiles_reduced
+    dct["smiles_hashed"] = smiles_hashed
+    dct["enumerated_sdf"] = enumerated_sdf
+    dct["sorted_sdf"] = sorted_sdf
+    dct["housekeeping_folder"] = housekeeping_folder
+    dct["path"] = path
+    dct["dir"] = dir
+    return dct
+
+def isomer_wraper(chunk_info, args, queue):
+    """
+    chunk_info: (path, dir) tuple for the chunk
+    args: auto3D arguments
+    queue: mp.queue
+    """
+    for i, path_dir in enumerate(chunk_info):
+        print(f"\n\nIsomer generation for job{i+1}")
+        path, dir = path_dir
+        meta = create_chunk_meta_names(path, dir)
+
+        # Tautomer enumeratioin
+        if args.enumerate_tautomer:
+            output_taut = meta["output_taut"]
+            taut_mode = args.taut_engine
+            print("Enumerating tautomers for the input...", end='')
+            taut_engine = tautomer_engine(taut_mode, path, output_taut)
+            taut_engine.run()
+            hash_taut_smi(output_taut, output_taut)
+            path = output_taut
+            print(f"Tautomers are saved in {output_taut}")
+
+        smiles_enumerated = meta["smiles_enumerated"]
+        smiles_reduced = meta["smiles_reduced"]
+        smiles_hashed = meta["smiles_hashed"]
+        enumerated_sdf = meta["enumerated_sdf"]
+        max_confs = args.max_confs
+        duplicate_threshold = args.threshold
+        mpi_np = args.mpi_np
+        cis_trans = args.cis_trans
+        isomer_program = args.isomer_engine
+        # Isomer enumeration step
+        if isomer_program == 'omega':
+            mode_oe = args.mode_oe
+            oe_isomer(mode_oe, path, smiles_enumerated, smiles_reduced, smiles_hashed,
+                    enumerated_sdf, max_confs, duplicate_threshold, cis_trans)
+        elif isomer_program == 'rdkit':
+            engine = rd_isomer(path, smiles_enumerated, smiles_reduced, smiles_hashed, 
+                            enumerated_sdf, dir, max_confs, duplicate_threshold, mpi_np, cis_trans)
+            engine.run()
+        else: 
+            raise ValueError('The isomer enumeration engine must be "omega" or "rdkit", '
+                            f'but {args.isomer_engine} was parsed. '
+                            'You can set the parameter by appending the following:'
+                            '--isomer_engine=rdkit')
+
+        queue.put((enumerated_sdf, path, dir))
+    queue.put("Done")
+
+
+def optim_rank_wrapper(args, queue):
+    job = 1
+    while True:
+        sdf_path_dir = queue.get()
+        if sdf_path_dir == "Done":
+            break
+        print(f"\n\nOptimizing on job{job}")
+        enumerated_sdf, path, dir = sdf_path_dir
+        meta = create_chunk_meta_names(path, dir)
+
+        # Optimizing step
+        opt_steps = args.opt_steps
+        opt_tol = args.convergence_threshold
+        config = {"opt_steps": opt_steps, "opttol": opt_tol}
+        optimized_og = meta["optimized_og"]
+        optimizing_engine = args.optimizing_engine
+        if args.use_gpu:
+            idx = args.gpu_idx
+            device = torch.device(f"cuda:{idx}")
+        else:
+            device = torch.device("cpu")
+        optimizer = optimizing(enumerated_sdf, optimized_og,
+                               optimizing_engine, device, config)
+        optimizer.run()
+
+        # Ranking step
+        output = meta["output"]
+        duplicate_threshold = args.threshold
+        k = args.k
+        window = args.window
+        rank_engine = ranking(optimized_og,
+                              output, duplicate_threshold, k=k, window=window)
+        rank_engine.run()
+
+        # Housekeeping
+        housekeeping_folder = meta["housekeeping_folder"]
+        os.mkdir(housekeeping_folder)
+        housekeeping(dir, housekeeping_folder, output)
+        #Conpress verbose folder
+        housekeeping_folder_gz = housekeeping_folder + ".tar.gz"
+        with tarfile.open(housekeeping_folder_gz, "w:gz") as tar:
+            tar.add(housekeeping_folder, arcname=os.path.basename(housekeeping_folder))
+        shutil.rmtree(housekeeping_folder)
+        if not args.verbose:
+            send2trash(housekeeping_folder_gz)
+        job += 1
+
+
+def options(path, k=False, window=False, verbose=False, job_name="",
+    enumerate_tautomer=False, tauto_engine="rdkit",
+    isomer_engine="rdkit", cis_trans=True, mode_oe="classic", mpi_np=4, max_confs=1000,
+    use_gpu=True, gpu_idx=0, capacity=42, optimizing_engine="AIMNET",
+    opt_steps=10000, convergence_threshold=0.003, threshold=0.3):
+    """Arguments for Auto3D main program
+    path: A input.smi containing SMILES and IDs. Examples are listed in the example/files folder
+    k: Outputs the top-k structures for each SMILES.
+    window: Outputs the structures whose energies are within x (kcal/mol) from the lowest energy conformer
+    verbose: When True, save all meta data while running.
+    job_name: A folder name to save all meta data.
+    
+    enumerate_tautomer: When True, enumerate tautomers for the input
+    taut_engine: Programs to enumerate tautomers, either 'rdkit' or 'oechem'
+    isomer_engine: The program for generating 3D isomers for each SMILES. This parameter is either rdkit or omega.
+    cis_trans: When True, cis/trans and r/s isomers are enumerated.
+    mode_oe: The mode that omega program will take. It can be either 'classic' or 'macrocycle'. By default, the 'classic' mode is used. For detailed information about each mode, see https://docs.eyesopen.com/applications/omega/omega/omega_overview.html
+    mpi_np: Number of CPU cores for the isomer generation engine.
+    max_confs: Maximum number of isomers for each SMILES.
+    use_gpu: If True, the program will use GPU when available
+    gpu_idx: GPU index. It only works when --use_gpu=True
+    capacity: Number of SMILES that the model will handle for 1 G memory
+    optimizing_engine: Choose either 'ANI2x', 'ANI2xt', or 'AIMNET' for energy calculation and geometry optimization.
+    opt_steps: Maximum optimization steps for each structure.
+    convergence_threshold: Optimization is considered as converged if maximum force is below this threshold.
+    threshold: If the RMSD between two conformers are within threhold, they are considered as duplicates. One of them will be removed.
+    """
+    d = {}
+    args = my_name_space(d)
+    args['path'] = path
+    args['k'] = k
+    args['window'] = window
+    args['verbose'] = verbose
+    args['job_name'] = job_name
+    args["enumerate_tautomer"] = enumerate_tautomer
+    args["tauto_engine"] = tauto_engine
+    args["isomer_engine"] = isomer_engine
+    args["cis_trans"] = cis_trans
+    args["mode_oe"] = mode_oe
+    args["mpi_np"] = mpi_np
+    args["max_confs"] = max_confs
+    args["use_gpu"] = use_gpu
+    args["capacity"] = capacity
+    args["gpu_idx"] = gpu_idx
+    args["optimizing_engine"] = optimizing_engine
+    args["opt_steps"] = opt_steps
+    args["convergence_threshold"] = convergence_threshold
+    args["threshold"] = threshold
+    return args
+
+def main(args:dict):
+    """Take the arguments from options and run Auto3D"""
+
+
+    chunk_line = mp.Manager().Queue(1)   #A queue managing two wrappers
+
+
+    start = time.time()
+    job_name = time.strftime('%Y%m%d-%H%M%S')
+
+    path = args.path
+    k = args.k
+    window = args.window
+    if (not k) and (not window):
+        sys.exit("Either k or window needs to be specified. "
+                "Usually, setting '--k=1' satisfies most needs.")
+    if args.job_name == "":
+        args.job_name = job_name
+    job_name = args.job_name
+
+    basename = os.path.basename(path)
+    # initialiazation
+    dir = os.path.dirname(os.path.abspath(path))
+    job_name = job_name + "_" + basename.split('.')[0].strip()
+    job_name = os.path.join(dir, job_name)
+    os.mkdir(job_name)
+    check_input(args)
+
+
+    ## Devide jobs based on memory
+    smiles_per_G = args.capacity  #Allow 40 SMILES per GB memory
+    if args.use_gpu:
+        gpu_idx = int(args.gpu_idx)
+        t = int(math.ceil(torch.cuda.get_device_properties(gpu_idx).total_memory/(1024**3)))
+    else:
+        t = psutil.virtual_memory().total/(1024**3)
+    chunk_size = t * smiles_per_G
+
+    #Get indexes for each chunk
+    df = pd.read_csv(path, sep=" ", header=None)
+    data_size = len(df)
+    num_chunks = int(data_size // chunk_size + 1)
+    print(f"There are {len(df)} SMILES, available memory is {t} GB.")
+    print(f"The task will be divided into {num_chunks} jobs.")
+    chunk_idxes = [[] for _ in range(num_chunks)]
+    for i in range(num_chunks):
+        idx = i
+        while idx < data_size:
+            chunk_idxes[i].append(idx)
+            idx += num_chunks
+
+    #Save each chunk as smi
+    chunk_info = []
+    basename = os.path.basename(path).split(".")[0].strip()
+    for i in range(num_chunks):
+        dir = os.path.join(job_name, f"job{i+1}")
+        os.mkdir(dir)
+        new_basename = basename + "_" + str(i+1) + ".smi"
+        new_name = os.path.join(dir, new_basename)
+        df_i = df.iloc[chunk_idxes[i], :]
+        df_i.to_csv(new_name, header=None, index=None, sep=" ")
+        path = new_name
+
+        print(f"Job{i+1}, number of inputs: {len(df_i)}")
+        chunk_info.append((path, dir))
+
+    # if __name__ == "__main__":
+    p1 = mp.Process(target=isomer_wraper, args=(chunk_info, args, chunk_line))
+    p2 = mp.Process(target=optim_rank_wrapper, args=(args, chunk_line,))
+    p1.start()
+    p2.start()
+    p1.join()
+    p2.join()
+
+    #Combine jobs into a single sdf
+    data = []
+    paths = os.path.join(job_name, "job*/*.sdf")
+    files = glob.glob(paths)
+    for file in files:
+        with open(file, "r") as f:
+            data_i = f.readlines()
+        data += data_i
+    combined_basename = basename + "_out.sdf"
+    path_combined = os.path.join(job_name, combined_basename)
+    with open(path_combined, "w+") as f:
+        for line in data:
+            f.write(line)
+
+    # Program ends
+    end = time.time()
+    print("Energy unit: Hartree if implicit.")
+    print(f'Program running time: {end - start} seconds')
+    return path_combined
+
+
+if __name__ == "__main__":
+    chunk_line = mp.Manager().Queue(1)   #A queue managing two wrappers
+    print("""
+        _              _             _____   ____  
+       / \     _   _  | |_    ___   |___ /  |  _ \ 
+      / _ \   | | | | | __|  / _ \    |_ \  | | | |
+     / ___ \  | |_| | | |_  | (_) |  ___) | | |_| |
+    /_/   \_\  \__,_|  \__|  \___/  |____/  |____/ 
+        // Automatically generating the lowest-energy 3D structure                                       
+    """)
+
+    start = time.time()
+    job_name = time.strftime('%Y%m%d-%H%M%S')
+
+    parser = argparse.ArgumentParser(prog='auto3D',
+                        description=('This program takes in a SMILES(s)'
+                                    ' and returns the optimal 3D structure(s)'),
+                        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    parser.add_argument('path', type=str,
+                        help=('A input.smi containing SMILES and IDs. '
+                            'Examples are listed in the input folder'))
+    parser.add_argument('--k', type=int, default=False,
+                        help='Outputs the top-k structures for each SMILES.')
+    parser.add_argument('--window', type=float, default=False,
+                        help=('Outputs the structures whose energies are within '
+                            'window (kcal/mol) from the lowest energy'))
+    parser.add_argument('--verbose', default=True, type=lambda x: (str(x).lower() == 'true'),
+                        help='When True, save all meta data while running.')
+    parser.add_argument('--job_name', type=str, default=job_name,
+                        help="A folder name to save all meta data.")
+
+    # parameters for isomer enumeration step
+    parser.add_argument('--enumerate_tautomer', default=False, type=lambda x: (str(x).lower() == 'true'),
+                        help="When True, enumerate tautomers for the input")
+    parser.add_argument('--taut_engine', type=str, default='rdkit',
+                        help="Programs to enumerate tautomers, either 'rdkit' or 'oechem'")
+    parser.add_argument('--isomer_engine', type=str, default='rdkit',
+                        help=('The program for generating 3D isomers for each '
+                            'SMILES. This parameter is either '
+                            'rdkit or omega. rdkit is free for everyone, '
+                            'while omega reuqires a license.'))
+    parser.add_argument('--cis_trans', default=True, type=lambda x: (str(x).lower() == 'true'),
+                        help='When True, cis/trans and r/s isomers are enumerated.')
+    parser.add_argument('--mode_oe', type=str, default='classic',
+                        help=("The mode that omega program will take. "
+                            "It can be either 'classic', 'macrocycle' "
+                            "'rocs', 'pose' or 'dense'. "
+                            "By default, the 'classic' mode is used. "
+                            "For detailed information about each mode, see "
+                            "https://docs.eyesopen.com/applications/omega/omega/omega_overview.html"))
+    parser.add_argument('--mpi_np', type=int, default=4,
+                        help="Number of CPU cores for the isomer generation step.")
+    parser.add_argument('--max_confs', type=int, default=1000,
+                        help="Maximum number of isomers for each SMILES")
+    parser.add_argument('--use_gpu', default=True, type=lambda x: (str(x).lower() == 'true'),
+                        help="If True, the program will use GPU.")
+    parser.add_argument('--gpu_idx', default=0, type=int, 
+                        help="GPU index. It only works when --use_gpu=True")
+    parser.add_argument('--capacity', default=42, type=int,
+                        help="Number of SMILES that the model will handle for 1 G memory")
+    parser.add_argument('--optimizing_engine', type=str, default='AIMNET',
+                        help=("Choose either 'ANI2x', 'ANI2xt', or 'AIMNET' for energy "
+                            "calculation and geometry optimization."))
+    parser.add_argument('--opt_steps', type=int, default=10000,
+                        help="Maximum optimization steps for each structure.")
+    parser.add_argument('--convergence_threshold', type=float, default=0.003,
+                        help="Optimization is considered as converged if maximum force is below this threshold.")
+    parser.add_argument('--threshold', type=float, default=0.3,
+                        help=("If the RMSD between two conformers are within threhold, "
+                            "they are considered as duplicates. One of them will be removed."))
+
+    args = parser.parse_args()
+    path = args.path
+    k = args.k
+    window = args.window
+    if (not k) and (not window):
+        sys.exit("Either k or window needs to be specified. "
+                "Usually, setting '--k=1' satisfies most needs.")
+    job_name = args.job_name
+
+
+
+    basename = os.path.basename(path)
+    # initialiazation
+    dir = os.path.dirname(os.path.abspath(path))
+    job_name = job_name + "_" + basename.split('.')[0].strip()
+    job_name = os.path.join(dir, job_name)
+    os.mkdir(job_name)
+    check_input(args)
+
+
+    ## Devide jobs based on memory
+    smiles_per_G = args.capacity  #Allow 40 SMILES per GB memory
+    if args.use_gpu:
+        gpu_idx = int(args.gpu_idx)
+        t = int(math.ceil(torch.cuda.get_device_properties(gpu_idx).total_memory/(1024**3)))
+    else:
+        t = psutil.virtual_memory().total/(1024**3)
+    chunk_size = t * smiles_per_G
+
+    #Get indexes for each chunk
+    df = pd.read_csv(path, sep=" ", header=None)
+    data_size = len(df)
+    num_chunks = int(data_size // chunk_size + 1)
+    print(f"There are {len(df)} SMILES, available memory is {t} GB.")
+    print(f"The task will be divided into {num_chunks} jobs.")
+    chunk_idxes = [[] for _ in range(num_chunks)]
+    for i in range(num_chunks):
+        idx = i
+        while idx < data_size:
+            chunk_idxes[i].append(idx)
+            idx += num_chunks
+
+    #Save each chunk as smi
+    chunk_info = []
+    basename = os.path.basename(path).split(".")[0].strip()
+    for i in range(num_chunks):
+        dir = os.path.join(job_name, f"job{i+1}")
+        os.mkdir(dir)
+        new_basename = basename + "_" + str(i+1) + ".smi"
+        new_name = os.path.join(dir, new_basename)
+        df_i = df.iloc[chunk_idxes[i], :]
+        df_i.to_csv(new_name, header=None, index=None, sep=" ")
+        path = new_name
+
+        print(f"Job{i+1}, number of inputs: {len(df_i)}")
+        chunk_info.append((path, dir))
+
+    p1 = mp.Process(target=isomer_wraper, args=(chunk_info, args, chunk_line))
+    p2 = mp.Process(target=optim_rank_wrapper, args=(args, chunk_line,))
+    p1.start()
+    p2.start()
+    p1.join()
+    p2.join()
+
+    #Combine jobs into a single sdf
+    data = []
+    paths = os.path.join(job_name, "job*/*.sdf")
+    files = glob.glob(paths)
+    for file in files:
+        with open(file, "r") as f:
+            data_i = f.readlines()
+        data += data_i
+    combined_basename = basename + "_out.sdf"
+    path_combined = os.path.join(job_name, combined_basename)
+    with open(path_combined, "w+") as f:
+        for line in data:
+            f.write(line)
+
+    # Program ends
+    end = time.time()
+    print("Energy unit: Hartree if implicit.")
+    print(f'Program running time: {end - start} seconds')
diff --git a/src/Auto3D/auto3D.yml b/src/Auto3D/auto3D.yml
new file mode 100644
index 0000000..3c32734
--- /dev/null
+++ b/src/Auto3D/auto3D.yml
@@ -0,0 +1,15 @@
+name: ani_rep5
+channels:
+  - rmg
+  - conda-forge
+  - rdkit
+  - pytorch
+  - defaults
+dependencies:
+  - python=3.7
+  - numpy
+  - pytorch
+  - rdkit[version='>=2021.09.02']
+  - openbabel[version='>=3.1.1']
+  - tqdm
+prefix: /home/jack/miniconda3/envs/ani_rep5
diff --git a/src/Auto3D/auto3d_logo.png b/src/Auto3D/auto3d_logo.png
new file mode 100644
index 0000000..88447f1
Binary files /dev/null and b/src/Auto3D/auto3d_logo.png differ
diff --git a/src/Auto3D/batch_opt/ANI2xt energies and repulsion.ipynb b/src/Auto3D/batch_opt/ANI2xt energies and repulsion.ipynb
new file mode 100644
index 0000000..786123a
--- /dev/null
+++ b/src/Auto3D/batch_opt/ANI2xt energies and repulsion.ipynb
@@ -0,0 +1,656 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "9ed8a52a",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/liucmu/miniconda3/envs/py37/lib/python3.7/site-packages/torchani-2.3.dev568+geffb1bb-py3.7.egg/torchani/infer.py:16: UserWarning: mnp not installed\n",
+      "  warnings.warn(\"mnp not installed\")\n",
+      "/home/liucmu/miniconda3/envs/py37/lib/python3.7/site-packages/torchani-2.3.dev568+geffb1bb-py3.7.egg/torchani/aev/aev_computer.py:24: UserWarning: cuaev not installed\n",
+      "  warnings.warn(\"cuaev not installed\")\n"
+     ]
+    }
+   ],
+   "source": [
+    "from openbabel import pybel\n",
+    "import numpy as np\n",
+    "import torch\n",
+    "import torch.nn as nn\n",
+    "import torchani\n",
+    "from torchani.repulsion import RepulsionCalculator, StandaloneRepulsionCalculator\n",
+    "from torchani.units import hartree2kcalmol \n",
+    "\n",
+    "device = torch.device(\"cuda:0\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "fa98015c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def mols2arrays(sdf_fil):\n",
+    "    species_order = (\"H\", 'C', 'N', 'O', 'S', 'F', 'Cl')\n",
+    "    ani2xt_index = {1:0, 6:1, 7:2, 8:3, 16:4, 9:5, 17:6}\n",
+    "    coord = np.array([[a.coords for a in mol.atoms] for mol in mols])\n",
+    "    numbers = np.array([[ani2xt_index[a.atomicnum] for a in mol.atoms] for mol in mols])\n",
+    "    return coord, numbers\n",
+    "\n",
+    "class ANI2xt(nn.Module):\n",
+    "    def __init__(self, device, state_dict='models/ani2xt.pt'):\n",
+    "        super().__init__()\n",
+    "        self.device = device\n",
+    "        self.state_dict = state_dict\n",
+    "        # setup constants and construct an AEV computer\n",
+    "        Rcr = 5.2000e+00\n",
+    "        Rca = 3.5000e+00\n",
+    "        EtaR = torch.tensor([1.6000000e+01], device=device)\n",
+    "        ShfR = torch.tensor([9.0000000e-01, 1.1687500e+00, 1.4375000e+00, 1.7062500e+00, 1.9750000e+00, 2.2437500e+00, 2.5125000e+00, 2.7812500e+00, 3.0500000e+00, 3.3187500e+00, 3.5875000e+00, 3.8562500e+00, 4.1250000e+00, 4.3937500e+00, 4.6625000e+00, 4.9312500e+00], device=device)\n",
+    "        Zeta = torch.tensor([3.2000000e+01], device=device)\n",
+    "        ShfZ = torch.tensor([1.9634954e-01, 5.8904862e-01, 9.8174770e-01, 1.3744468e+00, 1.7671459e+00, 2.1598449e+00, 2.5525440e+00, 2.9452431e+00], device=device)\n",
+    "        EtaA = torch.tensor([8.0000000e+00], device=device)\n",
+    "        ShfA = torch.tensor([9.0000000e-01, 1.5500000e+00, 2.2000000e+00, 2.8500000e+00], device=device)\n",
+    "#         species_order = [b'H', b'C', b'N', b'O', b'S', b'F', b'Cl']\n",
+    "        species_order = [\"H\", 'C', 'N', 'O', 'S', 'F', 'Cl']\n",
+    "\n",
+    "        num_species = len(species_order)\n",
+    "        aev_computer = torchani.AEVComputer(Rcr, Rca, EtaR, ShfR, EtaA, Zeta, ShfA, ShfZ, num_species)\n",
+    "        aev_dim = aev_computer.aev_length\n",
+    "\n",
+    "        H_network = torch.nn.Sequential(\n",
+    "            torch.nn.Linear(aev_dim, 256),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(256, 192),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(192, 160),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(160, 1)\n",
+    "        )\n",
+    "\n",
+    "        C_network = torch.nn.Sequential(\n",
+    "            torch.nn.Linear(aev_dim, 224),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(224, 192),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(192, 160),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(160, 1)\n",
+    "        )\n",
+    "\n",
+    "        N_network = torch.nn.Sequential(\n",
+    "            torch.nn.Linear(aev_dim, 192),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(192, 160),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(160, 128),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(128, 1)\n",
+    "        )\n",
+    "\n",
+    "        O_network = torch.nn.Sequential(\n",
+    "            torch.nn.Linear(aev_dim, 192),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(192, 160),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(160, 128),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(128, 1)\n",
+    "        )\n",
+    "\n",
+    "        S_network = torch.nn.Sequential(\n",
+    "            torch.nn.Linear(aev_dim, 160),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(160, 128),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(128, 96),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(96, 1)\n",
+    "        )\n",
+    "\n",
+    "        F_network = torch.nn.Sequential(\n",
+    "            torch.nn.Linear(aev_dim, 160),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(160, 128),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(128, 96),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(96, 1)\n",
+    "        )\n",
+    "\n",
+    "        Cl_network = torch.nn.Sequential(\n",
+    "            torch.nn.Linear(aev_dim, 160),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(160, 128),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(128, 96),\n",
+    "            torch.nn.CELU(0.1),\n",
+    "            torch.nn.Linear(96, 1)\n",
+    "        )\n",
+    "\n",
+    "        nn = torchani.ANIModel([H_network, C_network, N_network, O_network, S_network, F_network, Cl_network])\n",
+    "        checkpoint = torch.load(self.state_dict, map_location=self.device)\n",
+    "        nn.load_state_dict(checkpoint)\n",
+    "        self.shifter = torchani.utils.EnergyShifter(torch.tensor([ -0.6002,  -38.1224,  -54.7241,  -75.2086, \n",
+    "                                                                    -398.1251,  -99.8006, -460.1375], \n",
+    "                                                                    device=self.device, dtype=torch.float64))\n",
+    "        self.model = torchani.nn.Sequential(aev_computer, nn).to(device)\n",
+    "        \n",
+    "#         self.species_order = (\"H\", 'C', 'N', 'O', 'S', 'F', 'Cl')\n",
+    "        self.rep = StandaloneRepulsionCalculator(elements=species_order, periodic_table_index=False).to(device)\n",
+    "#         self.model.train(False)\n",
+    "\n",
+    "    \n",
+    "    def forward(self, species, coords):\n",
+    "        #Electronic part\n",
+    "        energy = self.shifter(self.model((species, coords))).energies\n",
+    "        gradient = torch.autograd.grad([energy.sum()], [coords])[0]\n",
+    "        force = -gradient\n",
+    "    \n",
+    "        #Repulsion part\n",
+    "        rep_energy = self.rep((species, coords)).energies\n",
+    "        derivative = torch.autograd.grad(rep_energy.sum(), coords)[0]\n",
+    "        rep_force = -derivative\n",
+    "        \n",
+    "        final_energy = energy + rep_energy\n",
+    "        final_force = force + rep_force\n",
+    "        return (final_energy, final_force)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "be934878",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "(2,) (2,)\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/liucmu/miniconda3/envs/py37/lib/python3.7/site-packages/ipykernel_launcher.py:4: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
+      "  after removing the cwd from sys.path.\n",
+      "/home/liucmu/miniconda3/envs/py37/lib/python3.7/site-packages/ipykernel_launcher.py:5: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
+      "  \"\"\"\n"
+     ]
+    },
+    {
+     "ename": "TypeError",
+     "evalue": "can't convert np.ndarray of type numpy.object_. The only supported types are: float64, float32, float16, complex64, complex128, int64, int32, int16, int8, uint8, and bool.",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
+      "\u001b[0;32m/tmp/ipykernel_2711/458750939.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcoord\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnumbers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0mcoord\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcoord\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrequires_grad\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdevice\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdevice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0mspecies\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnumbers\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlong\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdevice\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdevice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mTypeError\u001b[0m: can't convert np.ndarray of type numpy.object_. The only supported types are: float64, float32, float16, complex64, complex128, int64, int32, int16, int8, uint8, and bool."
+     ]
+    }
+   ],
+   "source": [
+    "# Repulsion energy and forces\n",
+    "# species_order = (\"H\", 'C', 'N', 'O', 'S', 'F', 'Cl')\n",
+    "# rep = StandaloneRepulsionCalculator(elements=species_order, periodic_table_index=False).to(device)\n",
+    "\n",
+    "f = \"/home/liucmu/auto3D/workflow/batch_opt/test100.sdf\"\n",
+    "f2 = \"/home/liucmu/auto3D/workflow/batch_opt/test2.sdf\"\n",
+    "mols = list(pybel.readfile('sdf', f2))\n",
+    "coord, numbers = mols2arrays(mols)\n",
+    "print(coord.shape, numbers.shape)\n",
+    "\n",
+    "coord = torch.tensor(coord, dtype=torch.float64, requires_grad=True, device=device).float()\n",
+    "species = torch.tensor(numbers, dtype=torch.long, device=device)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "8cddfdb0",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/liucmu/miniconda3/envs/py37/lib/python3.7/site-packages/torchani-2.3.dev568+geffb1bb-py3.7.egg/torchani/utils.py:294: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
+      "  self_energies = torch.tensor(self_energies, dtype=torch.double)\n"
+     ]
+    }
+   ],
+   "source": [
+    "# ANI2xt energy and forces\n",
+    "ani2xt = ANI2xt(device, state_dict='/home/liucmu/auto3D/workflow/models/ani2xt.pt')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "8d5d0e12",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "ANI2xt(\n",
+      "  (shifter): EnergyShifter()\n",
+      "  (model): Sequential(\n",
+      "    (0): AEVComputer(\n",
+      "      (angular_terms): StandardAngular(\n",
+      "        (cutoff_fn): CutoffCosine()\n",
+      "      )\n",
+      "      (radial_terms): StandardRadial(\n",
+      "        (cutoff_fn): CutoffCosine()\n",
+      "      )\n",
+      "      (neighborlist): FullPairwise()\n",
+      "    )\n",
+      "    (1): ANIModel(\n",
+      "      (0): Sequential(\n",
+      "        (0): Linear(in_features=1008, out_features=256, bias=True)\n",
+      "        (1): CELU(alpha=0.1)\n",
+      "        (2): Linear(in_features=256, out_features=192, bias=True)\n",
+      "        (3): CELU(alpha=0.1)\n",
+      "        (4): Linear(in_features=192, out_features=160, bias=True)\n",
+      "        (5): CELU(alpha=0.1)\n",
+      "        (6): Linear(in_features=160, out_features=1, bias=True)\n",
+      "      )\n",
+      "      (1): Sequential(\n",
+      "        (0): Linear(in_features=1008, out_features=224, bias=True)\n",
+      "        (1): CELU(alpha=0.1)\n",
+      "        (2): Linear(in_features=224, out_features=192, bias=True)\n",
+      "        (3): CELU(alpha=0.1)\n",
+      "        (4): Linear(in_features=192, out_features=160, bias=True)\n",
+      "        (5): CELU(alpha=0.1)\n",
+      "        (6): Linear(in_features=160, out_features=1, bias=True)\n",
+      "      )\n",
+      "      (2): Sequential(\n",
+      "        (0): Linear(in_features=1008, out_features=192, bias=True)\n",
+      "        (1): CELU(alpha=0.1)\n",
+      "        (2): Linear(in_features=192, out_features=160, bias=True)\n",
+      "        (3): CELU(alpha=0.1)\n",
+      "        (4): Linear(in_features=160, out_features=128, bias=True)\n",
+      "        (5): CELU(alpha=0.1)\n",
+      "        (6): Linear(in_features=128, out_features=1, bias=True)\n",
+      "      )\n",
+      "      (3): Sequential(\n",
+      "        (0): Linear(in_features=1008, out_features=192, bias=True)\n",
+      "        (1): CELU(alpha=0.1)\n",
+      "        (2): Linear(in_features=192, out_features=160, bias=True)\n",
+      "        (3): CELU(alpha=0.1)\n",
+      "        (4): Linear(in_features=160, out_features=128, bias=True)\n",
+      "        (5): CELU(alpha=0.1)\n",
+      "        (6): Linear(in_features=128, out_features=1, bias=True)\n",
+      "      )\n",
+      "      (4): Sequential(\n",
+      "        (0): Linear(in_features=1008, out_features=160, bias=True)\n",
+      "        (1): CELU(alpha=0.1)\n",
+      "        (2): Linear(in_features=160, out_features=128, bias=True)\n",
+      "        (3): CELU(alpha=0.1)\n",
+      "        (4): Linear(in_features=128, out_features=96, bias=True)\n",
+      "        (5): CELU(alpha=0.1)\n",
+      "        (6): Linear(in_features=96, out_features=1, bias=True)\n",
+      "      )\n",
+      "      (5): Sequential(\n",
+      "        (0): Linear(in_features=1008, out_features=160, bias=True)\n",
+      "        (1): CELU(alpha=0.1)\n",
+      "        (2): Linear(in_features=160, out_features=128, bias=True)\n",
+      "        (3): CELU(alpha=0.1)\n",
+      "        (4): Linear(in_features=128, out_features=96, bias=True)\n",
+      "        (5): CELU(alpha=0.1)\n",
+      "        (6): Linear(in_features=96, out_features=1, bias=True)\n",
+      "      )\n",
+      "      (6): Sequential(\n",
+      "        (0): Linear(in_features=1008, out_features=160, bias=True)\n",
+      "        (1): CELU(alpha=0.1)\n",
+      "        (2): Linear(in_features=160, out_features=128, bias=True)\n",
+      "        (3): CELU(alpha=0.1)\n",
+      "        (4): Linear(in_features=128, out_features=96, bias=True)\n",
+      "        (5): CELU(alpha=0.1)\n",
+      "        (6): Linear(in_features=96, out_features=1, bias=True)\n",
+      "      )\n",
+      "    )\n",
+      "  )\n",
+      "  (rep): StandaloneRepulsionCalculator(\n",
+      "    (cutoff_function): CutoffSmooth()\n",
+      "    (neighborlist): FullPairwise()\n",
+      "    (species_converter): SpeciesConverter()\n",
+      "  )\n",
+      ")\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(ani2xt)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "afa0f71e",
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "torch.Size([100, 23]) torch.Size([100, 23, 3])\n",
+      "torch.Size([100]) torch.Size([100, 23, 3])\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(species.shape, coord.shape)\n",
+    "e, f = ani2xt(species, coord)\n",
+    "print(e.size(), f.size())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "f1fc3391",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "tensor([-612.4149, -612.4152, -612.4145, -612.4130, -612.4192, -612.4322,\n",
+      "        -612.4315, -612.4190, -612.4199, -612.4145, -612.4167, -612.4210,\n",
+      "        -612.4157, -612.4252, -612.4227, -612.4123, -612.4212, -612.4182,\n",
+      "        -612.4103, -612.4162, -612.4165, -612.4093, -612.4110, -612.4005,\n",
+      "        -612.4043, -612.4265, -612.4108, -612.4109, -612.4178, -612.4153,\n",
+      "        -612.4110, -612.4025, -612.4054, -612.4193, -612.4091, -612.4187,\n",
+      "        -612.4244, -612.4065, -612.4345, -612.4116, -612.4200, -612.4200,\n",
+      "        -612.4094, -612.4120, -612.3883, -612.4274, -612.4202, -612.4151,\n",
+      "        -612.4162, -612.4085, -612.4070, -612.4196, -612.4221, -612.3979,\n",
+      "        -612.4320, -612.4278, -612.4259, -612.4216, -612.4101, -612.4235,\n",
+      "        -612.4080, -612.4258, -612.4057, -612.4034, -612.4145, -612.4074,\n",
+      "        -612.4117, -612.4316, -612.4200, -612.4236, -612.4276, -612.4061,\n",
+      "        -612.4100, -612.4095, -612.4084, -612.4228, -612.4117, -612.4105,\n",
+      "        -612.4178, -612.4219, -612.4192, -612.4172, -612.4164, -612.4264,\n",
+      "        -612.3984, -612.4131, -612.4219, -612.4150, -612.4082, -612.4163,\n",
+      "        -612.4381, -612.4219, -612.4095, -612.4171, -612.4276, -612.4048,\n",
+      "        -612.4083, -612.4205, -612.4185, -612.4095], device='cuda:0',\n",
+      "       dtype=torch.float64, grad_fn=<AddBackward0>)\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(e)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "774a93da",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "tensor([[[ 0.0874,  0.0021,  0.0089],\n",
+      "         [-0.0821, -0.1736,  0.0256],\n",
+      "         [-0.0326, -0.0429,  0.0072],\n",
+      "         ...,\n",
+      "         [ 0.0127,  0.0706, -0.0476],\n",
+      "         [-0.0626, -0.0512, -0.0202],\n",
+      "         [-0.0546,  0.0329,  0.0982]],\n",
+      "\n",
+      "        [[-0.0811,  0.0759, -0.0304],\n",
+      "         [-0.0949, -0.1333,  0.0583],\n",
+      "         [ 0.1338, -0.0816,  0.0456],\n",
+      "         ...,\n",
+      "         [-0.0902, -0.0747, -0.0263],\n",
+      "         [ 0.0138,  0.0724, -0.0881],\n",
+      "         [-0.0367,  0.0611,  0.0664]],\n",
+      "\n",
+      "        [[ 0.1318,  0.0148,  0.0071],\n",
+      "         [-0.1010, -0.1883,  0.0903],\n",
+      "         [-0.0427, -0.0106,  0.0014],\n",
+      "         ...,\n",
+      "         [-0.0477,  0.0214,  0.0535],\n",
+      "         [ 0.0128,  0.0891, -0.0776],\n",
+      "         [-0.0583, -0.0734, -0.0442]],\n",
+      "\n",
+      "        ...,\n",
+      "\n",
+      "        [[-0.0802, -0.0682, -0.0255],\n",
+      "         [-0.1031,  0.1351,  0.0745],\n",
+      "         [ 0.1780,  0.0580,  0.0304],\n",
+      "         ...,\n",
+      "         [-0.0174, -0.0627, -0.0996],\n",
+      "         [-0.0277, -0.0399,  0.0647],\n",
+      "         [-0.0729,  0.0808, -0.0146]],\n",
+      "\n",
+      "        [[ 0.1036, -0.0283,  0.0526],\n",
+      "         [ 0.0402,  0.1487, -0.0794],\n",
+      "         [-0.1046,  0.0499, -0.0746],\n",
+      "         ...,\n",
+      "         [-0.0700,  0.0340,  0.0516],\n",
+      "         [-0.0015, -0.1039,  0.0370],\n",
+      "         [-0.0326,  0.0115, -0.0937]],\n",
+      "\n",
+      "        [[ 0.1105, -0.0161, -0.0162],\n",
+      "         [-0.1313,  0.1647,  0.1625],\n",
+      "         [-0.0278,  0.0069,  0.0030],\n",
+      "         ...,\n",
+      "         [-0.0841,  0.0571,  0.0061],\n",
+      "         [-0.0460, -0.0909,  0.0634],\n",
+      "         [ 0.0212, -0.0355, -0.1094]]], device='cuda:0')\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(f)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f256f1d0",
+   "metadata": {},
+   "source": [
+    "# Ani2xt on different molecules"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "6015fd4a",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/liucmu/miniconda3/envs/py37/lib/python3.7/site-packages/torchani-2.3.dev568+geffb1bb-py3.7.egg/torchani/infer.py:16: UserWarning: mnp not installed\n",
+      "  warnings.warn(\"mnp not installed\")\n",
+      "/home/liucmu/miniconda3/envs/py37/lib/python3.7/site-packages/torchani-2.3.dev568+geffb1bb-py3.7.egg/torchani/aev/aev_computer.py:24: UserWarning: cuaev not installed\n",
+      "  warnings.warn(\"cuaev not installed\")\n"
+     ]
+    }
+   ],
+   "source": [
+    "from openbabel import pybel\n",
+    "import numpy as np\n",
+    "import torch\n",
+    "import torch.nn as nn\n",
+    "import torchani\n",
+    "from torchani.repulsion import RepulsionCalculator, StandaloneRepulsionCalculator\n",
+    "from torchani.units import hartree2kcalmol\n",
+    "from utils import *\n",
+    "\n",
+    "device = torch.device(\"cuda:0\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "9a8bcca2",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[(3.5291, 0.1498, 0.2911), (2.752, -0.9304, 0.6543), (1.3866, -0.9348, 0.4238), (0.7977, 0.1598, -0.1788), (-0.6492, 0.1758, -0.4363), (-1.159, 1.1824, -0.979), (-1.5111, -0.9655, -0.0709), (-2.8731, -0.6811, -0.4321), (-3.4636, 0.3714, 0.3598), (1.5829, 1.2328, -0.5372), (2.9427, 1.2552, -0.3158), (4.589, 0.0881, 0.4994), (3.196, -1.8143, 1.1359), (0.7444, -1.7585, 0.6924), (-1.1934, -1.8247, -0.7339), (-1.4173, -1.2787, 0.9786), (-2.951, -0.5103, -1.457), (-2.804, 0.8004, 1.1135), (-3.8451, 1.1851, -0.2935), (-4.336, -0.0946, 0.8937), (1.1195, 2.0844, -1.0067), (3.5627, 2.1078, -0.6015)]] [[1, 1, 1, 1, 1, 3, 1, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Repulsion energy and forces\n",
+    "# species_order = (\"H\", 'C', 'N', 'O', 'S', 'F', 'Cl')\n",
+    "# rep = StandaloneRepulsionCalculator(elements=species_order, periodic_table_index=False).to(device)\n",
+    "\n",
+    "f = \"/home/liucmu/auto3D/workflow/batch_opt/test100.sdf\"\n",
+    "f2 = \"/home/liucmu/auto3D/workflow/batch_opt/test2.sdf\"\n",
+    "mols = list(pybel.readfile('sdf', f2))\n",
+    "coord, numbers = mols2lists(mols)\n",
+    "\n",
+    "coord_padded = padding_coords(coord[1:])\n",
+    "numbers_padded = padding_species(numbers[1:])\n",
+    "print(coord_padded, numbers_padded)\n",
+    "\n",
+    "# coord = torch.tensor(coord, dtype=torch.float64, requires_grad=True, device=device).float()\n",
+    "# species = torch.tensor(numbers, dtype=torch.long, device=device)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "e06bcbb8",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "torch.Size([1, 22]) torch.Size([1, 22, 3])\n",
+      "torch.Size([1]) torch.Size([1, 22, 3])\n",
+      "tensor([-479.2707], device='cuda:0', dtype=torch.float64,\n",
+      "       grad_fn=<AddBackward0>) tensor([[[ 0.0936,  0.0158,  0.0086],\n",
+      "         [ 0.0319, -0.1198,  0.0603],\n",
+      "         [-0.0293, -0.0727,  0.0289],\n",
+      "         [-0.1148, -0.0251, -0.0072],\n",
+      "         [ 0.0396,  0.0324,  0.0035],\n",
+      "         [-0.0553,  0.1102, -0.0629],\n",
+      "         [-0.0458, -0.0487, -0.0639],\n",
+      "         [-0.0360, -0.0847,  0.0275],\n",
+      "         [-0.0630, -0.0320, -0.0316],\n",
+      "         [-0.0255,  0.1054, -0.0522],\n",
+      "         [ 0.0796,  0.0623, -0.0150],\n",
+      "         [ 0.1309,  0.0017,  0.0215],\n",
+      "         [ 0.0437, -0.0735,  0.0412],\n",
+      "         [-0.0687, -0.1123,  0.0389],\n",
+      "         [ 0.0232, -0.0503, -0.0206],\n",
+      "         [ 0.0106, -0.0283,  0.1032],\n",
+      "         [-0.0035,  0.0213, -0.0652],\n",
+      "         [ 0.0779,  0.0474,  0.0817],\n",
+      "         [-0.0279,  0.0646, -0.0426],\n",
+      "         [-0.0541, -0.0123,  0.0387],\n",
+      "         [-0.0639,  0.1116, -0.0625],\n",
+      "         [ 0.0568,  0.0871, -0.0303]]], device='cuda:0')\n"
+     ]
+    }
+   ],
+   "source": [
+    "coords = torch.tensor(coord_padded, dtype=torch.float64, requires_grad=True, device=device).float()\n",
+    "species = torch.tensor(numbers_padded, dtype=torch.long, device=device)\n",
+    "\n",
+    "ani2xt = ANI2xt(device, state_dict='/home/liucmu/auto3D/workflow/models/ani2xt.pt')\n",
+    "\n",
+    "print(species.shape, coords.shape)\n",
+    "e, f = ani2xt(species, coords)\n",
+    "print(e.size(), f.size())\n",
+    "print(e, f)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "9a70215d",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[(3.5291, 0.1498, 0.2911),\n",
+       " (2.752, -0.9304, 0.6543),\n",
+       " (1.3866, -0.9348, 0.4238),\n",
+       " (0.7977, 0.1598, -0.1788),\n",
+       " (-0.6492, 0.1758, -0.4363),\n",
+       " (-1.159, 1.1824, -0.979),\n",
+       " (-1.5111, -0.9655, -0.0709),\n",
+       " (-2.8731, -0.6811, -0.4321),\n",
+       " (-3.4636, 0.3714, 0.3598),\n",
+       " (1.5829, 1.2328, -0.5372),\n",
+       " (2.9427, 1.2552, -0.3158),\n",
+       " (4.589, 0.0881, 0.4994),\n",
+       " (3.196, -1.8143, 1.1359),\n",
+       " (0.7444, -1.7585, 0.6924),\n",
+       " (-1.1934, -1.8247, -0.7339),\n",
+       " (-1.4173, -1.2787, 0.9786),\n",
+       " (-2.951, -0.5103, -1.457),\n",
+       " (-2.804, 0.8004, 1.1135),\n",
+       " (-3.8451, 1.1851, -0.2935),\n",
+       " (-4.336, -0.0946, 0.8937),\n",
+       " (1.1195, 2.0844, -1.0067),\n",
+       " (3.5627, 2.1078, -0.6015)]"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "coord[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3099feb5",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.7.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/src/Auto3D/batch_opt/ANI2xt.py b/src/Auto3D/batch_opt/ANI2xt.py
new file mode 100644
index 0000000..ca769d8
--- /dev/null
+++ b/src/Auto3D/batch_opt/ANI2xt.py
@@ -0,0 +1,124 @@
+#!/bin/usr/env python
+import torch
+import torch.nn as nn
+import torchani
+from torchani.repulsion import StandaloneRepulsionCalculator
+
+
+hartree2eV = 27.211385
+class ANI2xt(nn.Module):
+    def __init__(self, device, state_dict='models/ani2xt.pt'):
+        super().__init__()
+        self.device = device
+        self.state_dict = state_dict
+        # setup constants and construct an AEV computer
+        Rcr = 5.2000e+00
+        Rca = 3.5000e+00
+        EtaR = torch.tensor([1.6000000e+01], device=device)
+        ShfR = torch.tensor([9.0000000e-01, 1.1687500e+00, 1.4375000e+00, 1.7062500e+00, 1.9750000e+00, 2.2437500e+00, 2.5125000e+00, 2.7812500e+00, 3.0500000e+00, 3.3187500e+00, 3.5875000e+00, 3.8562500e+00, 4.1250000e+00, 4.3937500e+00, 4.6625000e+00, 4.9312500e+00], device=device)
+        Zeta = torch.tensor([3.2000000e+01], device=device)
+        ShfZ = torch.tensor([1.9634954e-01, 5.8904862e-01, 9.8174770e-01, 1.3744468e+00, 1.7671459e+00, 2.1598449e+00, 2.5525440e+00, 2.9452431e+00], device=device)
+        EtaA = torch.tensor([8.0000000e+00], device=device)
+        ShfA = torch.tensor([9.0000000e-01, 1.5500000e+00, 2.2000000e+00, 2.8500000e+00], device=device)
+#         species_order = [b'H', b'C', b'N', b'O', b'S', b'F', b'Cl']
+        # species_order = ["H", 'C', 'N', 'O', 'S', 'F', 'Cl']
+        species_order = ["H", 'C', 'N', 'O', 'F', 'S', 'Cl']
+
+        num_species = len(species_order)
+        aev_computer = torchani.AEVComputer(Rcr, Rca, EtaR, ShfR, EtaA, Zeta, ShfA, ShfZ, num_species)
+        aev_dim = aev_computer.aev_length
+
+        H_network = torch.nn.Sequential(
+            torch.nn.Linear(aev_dim, 256),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(256, 192),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(192, 160),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(160, 1)
+        )
+
+        C_network = torch.nn.Sequential(
+            torch.nn.Linear(aev_dim, 224),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(224, 192),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(192, 160),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(160, 1)
+        )
+
+        N_network = torch.nn.Sequential(
+            torch.nn.Linear(aev_dim, 192),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(192, 160),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(160, 128),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(128, 1)
+        )
+
+        O_network = torch.nn.Sequential(
+            torch.nn.Linear(aev_dim, 192),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(192, 160),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(160, 128),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(128, 1)
+        )
+
+        S_network = torch.nn.Sequential(
+            torch.nn.Linear(aev_dim, 160),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(160, 128),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(128, 96),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(96, 1)
+        )
+
+        F_network = torch.nn.Sequential(
+            torch.nn.Linear(aev_dim, 160),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(160, 128),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(128, 96),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(96, 1)
+        )
+
+        Cl_network = torch.nn.Sequential(
+            torch.nn.Linear(aev_dim, 160),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(160, 128),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(128, 96),
+            torch.nn.CELU(0.1),
+            torch.nn.Linear(96, 1)
+        )
+
+        nn = torchani.ANIModel([H_network, C_network, N_network, O_network, F_network, S_network, Cl_network])
+        checkpoint = torch.load(self.state_dict, map_location=self.device)
+        nn.load_state_dict(checkpoint)
+        self.shifter = torchani.utils.EnergyShifter(torch.tensor([ -0.6002,  -38.1224,  -54.7241,  -75.2086, 
+                                                                   -99.8006, -398.1251,  -460.1375], 
+                                                                    device=self.device, dtype=torch.float64))
+        self.model = torchani.nn.Sequential(aev_computer, nn).to(device)
+        self.rep = StandaloneRepulsionCalculator(elements=species_order, periodic_table_index=False).to(device)
+
+    
+    def forward(self, species, coords):
+        #Electronic part
+        energy = self.shifter(self.model((species, coords))).energies * hartree2eV
+        gradient = torch.autograd.grad([energy.sum()], [coords])[0]
+        force = -gradient
+    
+        #Repulsion part
+        rep_energy = self.rep((species, coords)).energies * hartree2eV
+        derivative = torch.autograd.grad(rep_energy.sum(), coords)[0]
+        rep_force = -derivative
+        
+        final_energy = energy + rep_energy
+        final_force = force + rep_force
+        return (final_energy, final_force)
diff --git a/src/Auto3D/batch_opt/__init__.py b/src/Auto3D/batch_opt/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/Auto3D/batch_opt/batchopt.py b/src/Auto3D/batch_opt/batchopt.py
new file mode 100644
index 0000000..613731b
--- /dev/null
+++ b/src/Auto3D/batch_opt/batchopt.py
@@ -0,0 +1,362 @@
+# Original source: /labspace/models/aimnet/batch_opt_script/
+
+import torch
+import warnings
+import os
+try:
+    import torchani
+except:
+    pass
+from collections import defaultdict
+from openbabel import pybel
+try:
+    from .ANI2xt import ANI2xt
+except:
+    pass
+from tqdm import tqdm
+
+
+hartree2eV = 27.211385
+@torch.jit.script
+class FIRE():
+    """a general optimization program """
+    # For a list of documentation for different optimization programs: https://wiki.fysik.dtu.dk/ase/ase/optimize.html
+    def __init__(self, coord):
+        ## default parameters
+        self.dt_max = 1.0
+        self.Nmin = 5
+        self.maxstep = 0.2
+        self.finc = 1.1
+        self.fdec = 0.5
+        self.astart = 0.1
+        self.fa = 0.99
+        self.v = torch.zeros_like(coord)
+        self.Nsteps = torch.zeros(coord.shape[0], dtype=torch.long, device=coord.device)
+        self.dt = torch.full(coord.shape[:1], 0.1, device=coord.device)
+        self.a = torch.full(coord.shape[:1], 0.1, device=coord.device)
+
+    def __call__(self, coord, forces):
+        """Moving atoms based on forces
+        
+        Arguments:
+            coord: coordinates of atoms. Size (Batch, N, 3), where Batch is
+                   the number of structures, N is the number of atom in each structure.
+            forces: forces on each atom. Size (Batch, N, 3).
+            
+        Return:
+            new coordinates that are moved based on input forces. Size (Batch, N, 3)"""
+        vf = (forces * self.v).flatten(-2, -1).sum(-1)
+        w_vf = vf > 0.0
+        if w_vf.all():
+            a = self.a.unsqueeze(-1).unsqueeze(-1)
+            v = self.v
+            f = forces
+            self.v = (1.0 - a) * v + a * v.flatten(-2, -1).norm(p=2, dim=-1).unsqueeze(-1).unsqueeze(-1) * f / f.flatten(-2, -1).norm(p=2, dim=-1).unsqueeze(-1).unsqueeze(-1)
+            self.Nsteps += 1
+        elif w_vf.any():
+            a = self.a[w_vf].unsqueeze(-1).unsqueeze(-1)
+            v = self.v[w_vf]
+            f = forces[w_vf]
+            self.v[w_vf] = (1.0 - a) * v + a * v.flatten(-2, -1).norm(p=2, dim=-1).unsqueeze(-1).unsqueeze(-1) * f / f.flatten(-2, -1).norm(p=2, dim=-1).unsqueeze(-1).unsqueeze(-1)
+
+            w_N = self.Nsteps > self.Nmin
+            w_vfN = w_vf & w_N
+            self.dt[w_vfN] = (self.dt[w_vfN] * self.finc).clamp(max=self.dt_max)
+            self.a[w_vfN] *= self.fa
+            self.Nsteps[w_vfN] += 1
+
+        w_vf = ~w_vf
+        if w_vf.all():
+            self.v[:] = 0.0
+            self.a[:] = torch.tensor(self.astart, device=self.a.device)
+            self.dt[:] *= self.fdec
+            self.Nsteps[:] = 0
+        elif w_vf.any():
+            self.v[w_vf] = torch.tensor(0.0, device=self.v.device)
+            self.a[w_vf] = torch.tensor(self.astart, device=self.a.device)
+            self.dt[w_vf] *= self.fdec
+            self.Nsteps[w_vf] = torch.tensor(0, device=self.v.device)
+
+        dt = self.dt.unsqueeze(-1).unsqueeze(-1)
+        self.v += dt * forces
+        dr = dt * self.v
+        normdr = dr.flatten(-2, -1).norm(p=2, dim=-1).unsqueeze(-1).unsqueeze(-1)
+        dr *= (self.maxstep / normdr).clamp(max=1.0)
+        return coord + dr
+
+    def clean(self, mask):
+        # types: (Tensor) -> bool
+        self.v = self.v[mask]
+        self.Nsteps = self.Nsteps[mask]
+        self.dt = self.dt[mask]
+        self.a = self.a[mask]
+        return True
+
+
+class EnForce_ANI(torch.nn.Module):
+    """Takes in an torch model, then defines two forward functions for it.
+    The input model should be able to calculate energy and disp_energy given
+    coordiantes, species and charges of a molecule. 
+
+    Arguments:
+        ani: model
+        batchsize_atoms: the maximum nmber atoms that can be handled in one batch.
+
+    Returns:
+        the energies and forces for the input molecules. One time calculation.
+    """
+    def __init__(self, ani, name, batchsize_atoms=1024*16):
+        super().__init__()
+        self.add_module('ani', ani)
+        self.name = name
+        self.batchsize_atoms = batchsize_atoms
+
+    def forward(self, coord, numbers, charges):
+        """Calculate the energies and forces for input molecules. Called by self.forward_batched
+        
+        Arguments:
+            coord: coordinates for all input structures. size (B, N, 3), where
+                  B is the number of structures in coord, N is the number of
+                  atoms in each structure, 3 represents xyz dimensions.
+            numbers: the periodic numbers for all atoms.
+            
+        Returns:
+            energies
+            forces
+        """
+        charge = torch.zeros_like(numbers[:, 0])
+
+        if self.name == "AIMNET":
+            d = self.ani(dict(coord=coord, numbers=numbers, charge=charges))  #Output from the model
+            e = (d['energy'] + d['disp_energy']).to(torch.double)
+            g = torch.autograd.grad([e.sum()], [coord])[0]  # size(100, 23, 3)
+            assert g is not None
+            f = -g
+        elif self.name == "ANI2xt":
+            e, f = self.ani(numbers, coord)
+        elif self.name == "ANI2x":
+            e = self.ani((numbers, coord)).energies
+            e = e * hartree2eV  #ANI2x (torch.models.ANI2x()) output energy unit is Hatree;
+                                     #ANI ASE interface unit is eV
+            g = torch.autograd.grad([e.sum()], [coord])[0]
+            f = -g
+        return e, f
+
+#    @torch.jit.script_method
+    def forward_batched(self, coord, numbers, charges):
+        """Calculate the energies and forces for input molecules.
+        
+        Arguments:
+            coord: coordinates for all input structures. size (B, N, 3), where
+                  B is the number of structures in coord, N is the number of
+                  atoms in each structure, 3 represents xyz dimensions.
+            numbers: the periodic numbers for all atoms. size (N, )
+            
+        Returns:
+            energies
+            forces
+        """
+        B, N = coord.shape[:2]
+        e = []
+        f = []
+        idx = torch.arange(B, device=coord.device)
+        # for batch in idx.split(self.batchsize_atoms // N):  #How was the batchsize_atoms decided?
+        for batch in idx.split(self.batchsize_atoms // N):  #How was the batchsize_atoms decided?
+            _e, _f = self(coord[batch], numbers[batch], charges[batch])
+            e.append(_e)
+            f.append(_f)
+        return torch.cat(e, dim=0), torch.cat(f, dim=0)
+
+def print_stats(state):
+    """Print the optimization status"""
+    numbers = state['numbers']
+    num_total = numbers.size()[0]
+    num_converged = torch.sum(state['converged_mask']).to('cpu')
+    num_active = num_total - num_converged
+    print("Total 3D structures: %i  Converged: %i   Active: %i" % (num_total, num_converged, num_active))
+
+def n_steps(state, n, opttol):
+    """Doing n steps optimization for each input. Only converged structures are 
+    modified at each step.
+    n_steps does not change input conformer order
+    
+    Argument:
+        state: an dictionary containing all information about this optimization step
+        n: optimization step"""
+    # t0 = perf_counter()
+    numbers = state['numbers']
+    charges = state['charges']
+    num_total = numbers.size()[0]
+    coord = state['coord']
+    optimizer = FIRE(coord)
+    for istep in tqdm(range(1, (n+1), 1)):
+        not_converged = ~ state['converged_mask']  #Size fixed
+        
+        # stop optimization if all structures converged.
+        if not not_converged.any():
+            break
+
+        coord = state['coord'][not_converged]  #Subset coordinates, size=not_converged.
+        numbers = state['numbers'][not_converged]
+        charges = state['charges'][not_converged]
+
+        coord.requires_grad_(True)
+        e, f = state['nn'].forward_batched(coord, numbers, charges)  #Key step to calculate all coordinates and forces.
+        coord.requires_grad_(False)
+
+        coord = optimizer(coord, f)
+        fmax = f.norm(dim=-1).max(dim=-1)[0]  #Size(100), Norm is the length of each vector. Here it returns the maximum force length for ecah conformer. Size (100)
+        not_converged_post = fmax > opttol
+
+        optimizer.clean(not_converged_post)  #Subset v, a in FIRE for next optimization
+        
+        state['converged_mask'][not_converged] = ~ not_converged_post  #Update converged_mask, so that converged structures will not be updated in future steps.
+        state['fmax'][not_converged] = fmax  #Update fmax for conformers that are optimized in this iteration
+        state['energy'][not_converged] = e.detach()  #Update energy for conformers that are optimized in this iteration
+        state['coord'][not_converged] = coord  #Update coordinates for conformers that are optimized in this iteration
+    
+        if (istep % (n//10)) == 0:
+            print_stats(state)
+    if istep == (n):
+        print("Reaching maximum optimization step:   ", end="")
+    else:
+        print(f"All structures converged at step {istep}:   ", end="")
+    print_stats(state)
+
+
+def ensemble_opt(net, coord, numbers, charges, param, model, device):
+    """Optimizing a group of molecules
+    
+    Arguments:
+    net: a EnForce_ANI object
+    coord: coordinates of input molecules (N, m, 3). N is the number of structures
+           m is the number of atoms in each structure.
+    numbers: periodict number of atoms in the molecule (include H). size is m.
+    param: a dictionary containing parameters
+    
+    """
+    coord = torch.tensor(coord, dtype=torch.float, device=device)
+    numbers = torch.tensor(numbers, dtype=torch.long, device=device)
+    charges = torch.tensor(charges, dtype=torch.long, device=device)
+    converged_mask = torch.zeros(coord.shape[0], dtype=torch.bool, device=device)
+    fmax = torch.full(coord.shape[:1], 999.0, device=coord.device)  # size=N, a tensored filled with 999.0, representing the current maximum forces at each conformer.
+    energy = torch.full(coord.shape[:1], 999.0, dtype=torch.double, device=coord.device)
+    ids = torch.arange(coord.shape[0], device=coord.device)  #Returns a 1D tensor
+    optimizer = FIRE(coord)
+
+    state = dict(
+        ids=ids,
+        coord=coord, numbers=numbers, converged_mask=converged_mask,
+        optimizer=optimizer, nn=net, fmax=fmax, energy=energy,
+        timing=defaultdict(float), charges=charges,
+        he=list(), close=list()  #!!! he and close?
+    )
+
+    n_steps(state, param['opt_steps'], param['opttol'])
+
+    return dict(
+        coord=state['coord'].tolist(),
+        ids=state['ids'].tolist(),
+        energy=state['energy'].tolist(),
+        fmax=state['fmax'].tolist(),
+        he=state['he'],
+        close=state['close'],
+        timing=dict(state['timing']),
+        numbers=state['numbers'].tolist()
+    )
+
+
+def padding_coords(lists, pad_value=0.0):
+    lengths = [len(lst) for lst in lists]
+    pad_length = [max(lengths) - len(lst) for lst in lists]
+    assert (len(pad_length) == len(lists))
+    
+    lists_padded = []
+    for i in range(len(pad_length)):
+        lst_i = lists[i]
+        pad_i = [(pad_value, pad_value, pad_value) for _ in range(pad_length[i])]
+        lst_i_padded = lst_i + pad_i
+        lists_padded.append(lst_i_padded)
+    return lists_padded
+
+def padding_species(lists, pad_value=-1):
+    lengths = [len(lst) for lst in lists]
+    pad_length = [max(lengths) - len(lst) for lst in lists]
+    assert (len(pad_length) == len(lists))
+    
+    lists_padded = []
+    for i in range(len(pad_length)):
+        lst_i = lists[i]
+        pad_i = [pad_value for _ in range(pad_length[i])]
+        lst_i_padded = lst_i + pad_i
+        lists_padded.append(lst_i_padded)
+    return lists_padded
+
+def mols2lists(mols, model):
+    species_order = ("H", 'C', 'N', 'O', 'S', 'F', 'Cl')
+    # ani2xt_index = {1:0, 6:1, 7:2, 8:3, 16:4, 9:5, 17:6}
+    ani2xt_index = {1:0, 6:1, 7:2, 8:3, 9:4, 16:5, 17:6}
+    coord = [[a.coords for a in mol.atoms] for mol in mols]
+    charges = [mol.charge for mol in mols]
+    if model == "ANI2xt":
+        numbers = [[ani2xt_index[a.atomicnum] for a in mol.atoms] for mol in mols]
+    else:
+        numbers = [[a.atomicnum for a in mol.atoms] for mol in mols]
+    return coord, numbers, charges
+
+
+class optimizing(object):
+    def __init__(self, in_f, out_f, model, device, config):
+        self.in_f = in_f
+        self.out_f = out_f
+        # self.energies = energies
+        self.model = model
+        self.device = device
+        self.config = config
+        root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+        if model == "AIMNET":
+            self.ani = torch.jit.load(os.path.join(root, "models/aimnet2nqed_pc14iall_b97m_sae.jpt"), map_location=device) # Return a ScriptModule object
+        elif model == "ANI2xt":
+            self.ani = ANI2xt(device, state_dict=os.path.join(root, "models/ani2xt_seed0.pt"))
+        elif model == "ANI2x":
+            self.ani = torchani.models.ANI2x(periodic_table_index=True).to(device)
+        else:
+            raise ValueError("Model has to be ANI2x, ANI2xt or AIMNET.")
+        
+    def run(self):
+        print("Preparing for parallel optimizing... (Max optimization steps: %i)" % self.config["opt_steps"])
+        mols = list(pybel.readfile('sdf', self.in_f))
+        print(f"Total 3D conformers: {len(mols)}")
+        coord, numbers, charges = mols2lists(mols, self.model)
+        if self.model == "AIMNET":
+            coord_padded = padding_coords(coord, 0)
+            numbers_padded = padding_species(numbers, 0)
+        else:
+            coord_padded = padding_coords(coord, 0)
+            numbers_padded = padding_species(numbers, -1)
+
+        for p in self.ani.parameters():  
+            p.requires_grad_(False)
+
+        ani = EnForce_ANI(self.ani, self.model)  # Interesting, EnForce_ANI inherites nn.module, bu can still accept a ScriptModule object as the input
+
+        with torch.jit.optimized_execution(False):
+            optdict = ensemble_opt(ani, coord_padded, numbers_padded, charges,
+                                   self.config, self.model, self.device)  #Magic step
+
+        energies = optdict['energy']
+        fmax = optdict['fmax']
+        convergence_mask = list(map(lambda x: (x <= self.config['opttol']), fmax))
+        with open(self.out_f, 'w') as f:
+            for i in range(len(mols)):
+                mol = mols[i]
+                idx = mol.title
+                fmax_i = fmax[i]
+                mol.data['E_tot'] = energies[i]
+                mol.data['fmax'] = fmax_i
+                mol.data['Converged'] = str(convergence_mask[i])
+                mol.data['ID'] = idx
+                coord = optdict['coord'][i]
+                for atom, c in zip(mol.atoms, coord):
+                    atom.OBAtom.SetVector(*c)
+                f.write(mol.write('sdf'))
diff --git a/src/Auto3D/input/case1.smi b/src/Auto3D/input/case1.smi
new file mode 100644
index 0000000..1a7aa64
--- /dev/null
+++ b/src/Auto3D/input/case1.smi
@@ -0,0 +1,2 @@
+N=C1OC(CN2CC(C)OC(C)C2)CN1 case1
+CCC1CN(C2COCC2)CC(CC)O1 case2
\ No newline at end of file
diff --git a/src/Auto3D/input/case1_2.smi b/src/Auto3D/input/case1_2.smi
new file mode 100644
index 0000000..9cd5009
--- /dev/null
+++ b/src/Auto3D/input/case1_2.smi
@@ -0,0 +1,3 @@
+N=C1OC(CN2CC(CC)OC(C)C2)CN1 smi1
+CC12CCC(C1)C(C)(C)C2O CHEMBL1907996
+NC1C2C3C4C1C1C2C3C41 CHEMBL2283501
\ No newline at end of file
diff --git a/src/Auto3D/input/case1_3.smi b/src/Auto3D/input/case1_3.smi
new file mode 100644
index 0000000..ee580b9
--- /dev/null
+++ b/src/Auto3D/input/case1_3.smi
@@ -0,0 +1,1 @@
+[H]/N=C1OC(CN2CC(C)OCC2C)CN/1 smi1
\ No newline at end of file
diff --git a/src/Auto3D/input/cyclo_stereo.smi b/src/Auto3D/input/cyclo_stereo.smi
new file mode 100644
index 0000000..52a0353
--- /dev/null
+++ b/src/Auto3D/input/cyclo_stereo.smi
@@ -0,0 +1,4 @@
+CC1CCC(CO1)C smi2
+CC1CCCC(CC)C1 smi3
+C1CCCCCC=C1 cyclooctene
+C1CCCC=CCCCCCCCCOCCCCC1 smi4
\ No newline at end of file
diff --git a/src/Auto3D/input/examples.smi b/src/Auto3D/input/examples.smi
new file mode 100644
index 0000000..8a036aa
--- /dev/null
+++ b/src/Auto3D/input/examples.smi
@@ -0,0 +1,9 @@
+OC1(C)CCCCO1 rs0
+C1=CC(=CC=C1)C2[C@H](N[C@H](O2)C3=CC=C(C=C3)Cl)C rs1
+CC1(C)O[C@@]2([H])[C@@]([C@H](N(O)C)OC2)([H])O1 rs2
+CCC(C=C(CC(C)C)O)=S ez0
+C1=CC(=CC=C1)C(O)[C@H](C)N=CC2=CC=C(C=C2)F ez1
+CC/C=C(OC)\O ez2
+C(C=C(C(C)CC)O)(C)=S both0
+C1=CC(=CC=C1)C(O)[C@@H](C)N=CC2=CC=C(C=C2)C both1
+OC(C)(C)C([H])[C@@H](C)/C(N)=C/C(OC)=O both2
diff --git a/src/Auto3D/input/r_s.smi b/src/Auto3D/input/r_s.smi
new file mode 100644
index 0000000..a5d2223
--- /dev/null
+++ b/src/Auto3D/input/r_s.smi
@@ -0,0 +1,5 @@
+CC(O)/C=C/C1=CC=CC=C1 smi1
+CCCCCC(O)C=C smi2
+C/C(O)=C/CCC(C)CC smi3
+CCCCCC(O)C=CCC smi4
+BrC=CC1OC(C2)(F)C2(Cl)C1 smi5
\ No newline at end of file
diff --git a/src/Auto3D/input/smi4.smi b/src/Auto3D/input/smi4.smi
new file mode 100644
index 0000000..c650ad6
--- /dev/null
+++ b/src/Auto3D/input/smi4.smi
@@ -0,0 +1,1 @@
+[C@H]1(C2=CC=CC=C2)CC/C=C\CC[C@H]1C3=CC=CC=C3 keystone
diff --git a/src/Auto3D/input/smiles.smi b/src/Auto3D/input/smiles.smi
new file mode 100644
index 0000000..eba6749
--- /dev/null
+++ b/src/Auto3D/input/smiles.smi
@@ -0,0 +1,4 @@
+CNCC(C1=CC=CC=C1)=O smi1
+CC(CC)=O smi2
+COC(/C=C/C)=O smi3
+CC=CCC(CC)=O smi4
\ No newline at end of file
diff --git a/src/Auto3D/input/smiles2.smi b/src/Auto3D/input/smiles2.smi
new file mode 100644
index 0000000..e1e65fe
--- /dev/null
+++ b/src/Auto3D/input/smiles2.smi
@@ -0,0 +1,2 @@
+CNC(C)C(C1=CC=CC=C1)=O smi1
+N=C1OC(CN2CC(C)OC(C)C2)CN1 case1
\ No newline at end of file
diff --git a/src/Auto3D/input/smiles3.smi b/src/Auto3D/input/smiles3.smi
new file mode 100644
index 0000000..65370db
--- /dev/null
+++ b/src/Auto3D/input/smiles3.smi
@@ -0,0 +1,1 @@
+[C@H]1(C2=CC=CC=C2)CC/C=C\CC[C@H]1C3=CC=CC=C3 keystone
\ No newline at end of file
diff --git a/src/Auto3D/input/smiles4.smi b/src/Auto3D/input/smiles4.smi
new file mode 100644
index 0000000..98d099d
--- /dev/null
+++ b/src/Auto3D/input/smiles4.smi
@@ -0,0 +1,12 @@
+CC12CCC(C1)C(C)(C)C2O CHEMBL1907996
+C1C2C3CC4C1C1C2NC3C41 CHEMBL353418
+CC1(C)C2CCC(CN)C1C2 CHEMBL1886231
+CC(O)C(C)NCCNC(C)C(C)O CHEMBL34447#
+CC(N)COCC(C)OCC(C)OCC(C)N CHEMBL273942#
+CC(C)C1C2CCC(C2)C1N CHEMBL430112
+N#CCC1C2CC3C(C2=O)C13 CHEMBL1982767
+CCC(O)OCC(C)OC(O)CC CHEMBL1900555
+OC1COC(O)C(O)C1O CHEMBL14343
+OC(C1CO1)C(O)C1CO1 CHEMBL1987719
+OC1CNC(O)C(O)C1O CHEMBL1187927
+CC1CC(O)C(C)CN1 CHEMBL1344552
\ No newline at end of file
diff --git a/src/Auto3D/input/stereo.smi b/src/Auto3D/input/stereo.smi
new file mode 100644
index 0000000..94512c9
--- /dev/null
+++ b/src/Auto3D/input/stereo.smi
@@ -0,0 +1,3 @@
+CC(CC)CCC(CC)OC smi1
+N=C1OC(CN2CC(C)OC(C)C2)CN1 smi2
+CC(CC)CC(O)C(F)CC smi3
\ No newline at end of file
diff --git a/src/Auto3D/input/sugars.smi b/src/Auto3D/input/sugars.smi
new file mode 100644
index 0000000..26fc1d2
--- /dev/null
+++ b/src/Auto3D/input/sugars.smi
@@ -0,0 +1,2 @@
+OC(C1O)C(OC1(CO)O)CO honey
+OCC1OC(OC2OC(CO)C(O)C2O)(CO)C(O)C(O)C1O sugar
\ No newline at end of file
diff --git a/src/Auto3D/input/taut.smi b/src/Auto3D/input/taut.smi
new file mode 100644
index 0000000..93367c2
--- /dev/null
+++ b/src/Auto3D/input/taut.smi
@@ -0,0 +1,3 @@
+O=C(CC(C)=O)CC smi2
+O=C(C)CC smi1
+N/C=C/C smi3
diff --git a/src/Auto3D/isomer_engine.py b/src/Auto3D/isomer_engine.py
new file mode 100644
index 0000000..db681d2
--- /dev/null
+++ b/src/Auto3D/isomer_engine.py
@@ -0,0 +1,305 @@
+#!/usr/bin/env python
+"""
+Enumerating stereoisomers for each SMILES representation with RDKit.
+"""
+import warnings
+import os
+import glob
+import collections
+from openbabel import pybel
+from rdkit import Chem
+from rdkit.Chem import AllChem
+from rdkit.Chem.EnumerateStereoisomers import EnumerateStereoisomers
+from rdkit.Chem.EnumerateStereoisomers import StereoEnumerationOptions
+from rdkit.Chem.MolStandardize import rdMolStandardize
+from .utils import hash_enumerated_smi_IDs, combine_smi, amend_configuration_w
+from .utils import remove_enantiomers
+try:
+    from openeye import oechem
+    from openeye import oequacpac
+    from openeye import oeomega
+except:
+    pass
+from tqdm import tqdm
+
+
+class tautomer_engine(object):
+    """Enemerate possible tautomers for the input
+    
+    Arguments:
+        mode: rdkit or oechem
+        input: smi file
+        output: smi file
+        
+    """
+    def __init__(self, mode, input, out):
+        self.mode = mode
+        self.input = input
+        self.output = out
+
+    def oe_taut(self):
+        """OEChem enumerating tautomers, modified from
+        https://docs.eyesopen.com/toolkits/python/quacpactk/examples_summary_getreasonabletautomers.html"""
+        ifs = oechem.oemolistream()
+        ifs.open(self.input)
+
+        ofs = oechem.oemolostream()
+        ofs.open(self.output)
+
+        tautomerOptions = oequacpac.OETautomerOptions()
+        pKaNorm = False
+
+        for mol in ifs.GetOEGraphMols():
+            for tautomer in oequacpac.OEGetReasonableTautomers(mol, tautomerOptions, pKaNorm):
+                oechem.OEWriteMolecule(ofs, tautomer)
+        
+        # Appending input smiles into output
+        combine_smi([self.input, self.output], self.output)
+
+    def rd_taut(self):
+        """RDKit enumerating tautomers"""
+        enumerator = rdMolStandardize.TautomerEnumerator()
+        smiles = []
+        with open(self.input, 'r') as f:
+            data = f.readlines()
+            for line in data:
+                line = line.strip().split()
+                smi, idx = line[0], line[1]
+                smiles.append((smi, idx))
+        tautomers = []
+        for smi_idx in smiles:
+            smi, idx = smi_idx
+            mol = Chem.MolFromSmiles(smi)
+            tauts = enumerator.Enumerate(mol)
+            for taut in tauts:
+                tautomers.append((Chem.MolToSmiles(taut), idx))
+        with open(self.output, 'w+') as f:
+            for smi_idx in tautomers:
+                smi, idx = smi_idx
+                line = smi.strip() + ' ' + str(idx.strip()) + '\n'
+                f.write(line)
+
+    def run(self):
+        if self.mode == 'oechem':
+            self.oe_taut()
+        elif self.mode == 'rdkit':
+            self.rd_taut()
+        else:
+            raise ValueError(f'{self.mode} must be one of "oechem" or "rdkit".')
+
+class rd_isomer(object):
+    """
+    Enumerating stereoisomers for each SMILES representation with RDKit.
+
+    Arguments:
+        smi: A smi file containing SMILES and IDs.
+        smiles_enumerated: A smi containing cis/trans isomers for the smi file.
+        smiles_hashed: For smiles_enumerated, each ID is hashed.
+        enumerated_sdf: for smiles_hashed, generating possible 3D structures.
+        job_name: as the name suggests.
+        max_confs: maximum number of conformers for each smi.
+        threshold: Maximum RMSD to be considered as duplicates.
+    """
+    def __init__(self, smi, smiles_enumerated, smiles_enumerated_reduced,
+                 smiles_hashed, enumerated_sdf, job_name, max_confs, threshold, np,
+                 flipper=True):
+        self.input = smi
+        self.n_conformers = max_confs
+        self.enumerate = {}
+        self.enumerated_smi_path = smiles_enumerated
+        self.enumerated_smi_path_reduced = smiles_enumerated_reduced
+        self.enumerated_smi_hashed_path = smiles_hashed
+        self.enumerated_sdf = enumerated_sdf
+        self.num2sym = {1: 'H', 6: 'C', 8: 'O', 7: 'N',
+                        9: 'F', 16: 'S', 17: 'Cl'}
+        self.rdk_tmp = os.path.join(job_name, 'rdk_tmp')
+        os.mkdir(self.rdk_tmp)
+        self.threshold = threshold
+        self.np = np
+        self.flipper = flipper
+
+    @staticmethod
+    def read(input):
+        outputs = {}
+        with open(input, 'r') as f:
+            data = f.readlines()
+        for line in data:
+            smiles, name = tuple(line.strip().split())
+            outputs[name.strip()] = smiles.strip()
+        return outputs
+
+    @staticmethod
+    def enumerate_func(mol):
+        """Enumerate the R/S and cis/trans isomers
+        
+        Argument:
+            mol: rd mol object
+            
+        Return:
+            isomers: a list of SMILES"""
+        opts = StereoEnumerationOptions(unique=True)
+        isomers = tuple(EnumerateStereoisomers(mol, options=opts))
+        isomers = sorted(Chem.MolToSmiles(x, isomericSmiles=True, doRandom=False) for x in isomers)
+        return isomers
+
+    def write_enumerated_smi(self):
+        with open(self.enumerated_smi_path, 'w+') as f:
+            for name, smi in self.enumerate.items():
+                for i, isomer in enumerate(smi):
+                    new_name = str(name).strip() + '_' + str(i)
+                    line = isomer.strip() + '\t' + new_name + '\n'
+                    f.write(line)
+
+    def conformer_func(self, smi_name):
+        """smi_name is a tuple (smiles, name)"""
+        smi, name = smi_name
+        mol = Chem.AddHs(Chem.MolFromSmiles(smi))
+        AllChem.EmbedMultipleConfs(mol, numConfs=self.n_conformers,
+                                   randomSeed=42, numThreads=self.np,
+                                   pruneRmsThresh=self.threshold)
+        files = []
+        for i in range(mol.GetNumConformers()):
+            basename = name.strip() + f"_{i}.sdf"
+            mol_ID = basename.split(".")[0].strip()
+            mol.SetProp('ID', mol_ID)
+            file_path = os.path.join(self.rdk_tmp, basename)
+            
+            writer = Chem.SDWriter(file_path)
+            writer.write(mol, confId=i)
+            writer.close()
+            files.append(file_path)
+
+    def combine_SDF(self, SDFs, out):
+        """Combine and sort SDF files in folder into a single file"""
+        paths = os.path.join(SDFs, '*.sdf')
+        files = glob.glob(paths)
+        dict0 = {}
+        for file in files:
+            mols = pybel.readfile('sdf', file)
+            for mol in mols:
+                idx = mol.data['ID']
+                mol.title = idx
+                dict0[idx] = mol
+            dict0 = collections.OrderedDict(sorted(dict0.items()))
+
+        f = pybel.Outputfile('sdf', out)
+        for idx, mol in sorted(dict0.items()):
+            f.write(mol)
+        f.close()  
+
+    def run(self):
+        """
+        When called, enumerate 3 dimensional structures for the input file and
+        writes all structures in 'job_name/smiles_enumerated.sdf'
+        """
+        if self.flipper:
+            print("Enumerating cis/tran isomers for unspecified double bonds...")
+            print("Enumerating R/S isomers for unspecified atomic centers...")
+            smiles_og = self.read(self.input)
+            for name, smiles in smiles_og.items():
+                # mol = Chem.AddHs(Chem.MolFromSmiles(smiles))
+                mol = Chem.MolFromSmiles(smiles)
+                isomers = self.enumerate_func(mol)
+                self.enumerate[name] = isomers
+            self.write_enumerated_smi()
+            print("Removing enantiomers...")
+            amend_configuration_w(self.enumerated_smi_path)
+            remove_enantiomers(self.enumerated_smi_path, self.enumerated_smi_path_reduced)
+            hash_enumerated_smi_IDs(self.enumerated_smi_path_reduced,
+                                    self.enumerated_smi_hashed_path)
+        else:
+            hash_enumerated_smi_IDs(self.input,
+                                    self.enumerated_smi_hashed_path)
+
+        print("Enumerating conformers/rotamers, removing duplicates...")
+        smiles2 = self.read(self.enumerated_smi_hashed_path)
+
+        smi_name_tuples = [(smi, name) for name, smi in smiles2.items()]
+        # if __name__ == "__main__":
+        #     with Pool(self.np) as p:
+        #         p.map(self.conformer_func, smi_name_tuples)
+        for smi_name in smi_name_tuples:
+            self.conformer_func(smi_name)
+
+        self.combine_SDF(self.rdk_tmp, self.enumerated_sdf)
+
+
+def oe_flipper(input, out):
+    """helper function for oe_isomer"""
+    ifs = oechem.oemolistream()
+    ifs.open(input)
+    ofs = oechem.oemolostream()
+    ofs.open(out)
+
+    flipperOpts = oeomega.OEFlipperOptions()
+    flipperOpts.SetWarts(True)
+    flipperOpts.SetMaxCenters(12)
+    flipperOpts.SetEnumNitrogen(True)
+    flipperOpts.SetEnumBridgehead(True)
+    flipperOpts.SetEnumEZ(False)
+    flipperOpts.SetEnumRS(False)
+    for mol in ifs.GetOEMols():
+        # oechem.OEThrow.Info("Title: %s" % mol.GetTitle())
+        for enantiomer in oeomega.OEFlipper(mol.GetActive(), flipperOpts):
+            enantiomer = oechem.OEMol(enantiomer)
+            oechem.OEWriteMolecule(ofs, enantiomer)
+
+def oe_isomer(mode, input, smiles_enumerated, smiles_reduced, smiles_hashed, output, max_confs, threshold, flipper=True):
+    """Generating R/S, cis/trans and conformers using omega
+    Arguments:
+        mode: 'classic' or 'macrocycle'
+        input: input smi file
+        output: output SDF file
+        flipper: optional R/S and cis/trans enumeration"""
+    
+    if mode == "classic":
+        omegaOpts = oeomega.OEOmegaOptions()
+    elif mode == "macrocycle":
+        omegaOpts = oeomega.OEMacrocycleOmegaOptions()
+    else:
+        raise ValueError(f"mode has to be 'classic' or 'macrocycle', but received {mode}.")
+    omegaOpts.SetParameterVisibility(oechem.OEParamVisibility_Hidden) 
+    omegaOpts.SetParameterVisibility("-rms", oechem.OEParamVisibility_Simple)
+    omegaOpts.SetParameterVisibility("-ewindow", oechem.OEParamVisibility_Simple)
+    omegaOpts.SetParameterVisibility("-maxconfs", oechem.OEParamVisibility_Simple)
+    
+    omegaOpts.SetMaxConfs(max_confs)
+    omegaOpts.SetEnergyWindow(999)
+    if mode == "classic":
+        omegaOpts.SetFixRMS(threshold)  #macrocycle mode does not have the attribute 'SetFixRMS'
+        omegaOpts.SetStrictStereo(False)
+        omegaOpts.SetWarts(True)
+    else:
+        omegaOpts.SetIterCycleSize(1000)
+        omegaOpts.SetMaxIter(2000)
+    opts = oechem.OESimpleAppOptions(omegaOpts, "Omega", oechem.OEFileStringType_Mol, oechem.OEFileStringType_Mol3D)
+
+    omegaOpts.UpdateValues(opts)
+    if mode == "classic":
+        omega = oeomega.OEOmega(omegaOpts)
+    else:
+        omega = oeomega.OEMacrocycleOmega(omegaOpts)
+
+    if flipper:
+        print("Enumerating stereoisomers.")
+        oe_flipper(input, smiles_enumerated)
+        amend_configuration_w(smiles_enumerated)
+        remove_enantiomers(smiles_enumerated, smiles_reduced)
+        ifs = oechem.oemolistream()
+        ifs.open(smiles_reduced)
+    else:
+        ifs = oechem.oemolistream()
+        ifs.open(input)
+    ofs = oechem.oemolostream()
+    ofs.open(output)
+
+    print("Enumerating conformers.")
+    for mol in tqdm(ifs.GetOEMols()):
+        # oechem.OEThrow.Info("Title: %s" % mol.GetTitle())
+        ret_code = omega.Build(mol)
+        if ret_code == oeomega.OEOmegaReturnCode_Success:
+            oechem.OEWriteMolecule(ofs, mol)
+        else:
+            oechem.OEThrow.Warning("%s: %s" % (mol.GetTitle(), oeomega.OEGetOmegaError(ret_code)))
+
+    return 0
diff --git a/src/Auto3D/models/__init__.py b/src/Auto3D/models/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/Auto3D/models/aimnet2nqed_pc14iall_b97m.jpt b/src/Auto3D/models/aimnet2nqed_pc14iall_b97m.jpt
new file mode 100644
index 0000000..563190e
Binary files /dev/null and b/src/Auto3D/models/aimnet2nqed_pc14iall_b97m.jpt differ
diff --git a/src/Auto3D/models/aimnet2nqed_pc14iall_b97m_sae.jpt b/src/Auto3D/models/aimnet2nqed_pc14iall_b97m_sae.jpt
new file mode 100644
index 0000000..7ff475d
Binary files /dev/null and b/src/Auto3D/models/aimnet2nqed_pc14iall_b97m_sae.jpt differ
diff --git a/src/Auto3D/models/ani2xt_seed0.pt b/src/Auto3D/models/ani2xt_seed0.pt
new file mode 100644
index 0000000..f95482c
Binary files /dev/null and b/src/Auto3D/models/ani2xt_seed0.pt differ
diff --git a/src/Auto3D/ranking.py b/src/Auto3D/ranking.py
new file mode 100644
index 0000000..d8705d5
--- /dev/null
+++ b/src/Auto3D/ranking.py
@@ -0,0 +1,205 @@
+#!/usr/bin/env python
+'''
+Finding 3D structures that satisfy the input requirement.
+'''
+import os
+import pandas as pd
+from openbabel import pybel
+from .utils import guess_file_type, filter_unique
+
+eV2kcalmol = 23.06055  #1 eV = 23.06055 kcal/mol
+hartree2eV = 27.211385
+class ranking(object):
+    '''
+    Finding 3D structures that satisfy the user-defined requirements.
+
+    Arguments:
+        input_path: An SDF file that contains all isomers.
+        energies: A txt file that contains the IDs and energies.
+        out_path: An SDF file that stores the optimical structures.
+        k: Outputs the top-k structures for each SMILES.
+        window: Outputs the structures whose energies are within
+                window (kcal/mol) from the lowest energy
+    Returns:
+        None
+    '''
+    def __init__(self, input_path,
+                 out_path, threshold, k=False, window=False):
+        self.input_path = input_path
+        # self.energies = energies
+        self.out_path = out_path
+        self.threshold = threshold
+        self.atomic_number2symbol = {1: 'H', 6: 'C', 8: 'O',
+                                     7: 'N', 9: 'F', 16: 'S', 17: 'Cl'}
+        self.k = k
+        self.window = window
+
+    @staticmethod
+    def similar(name, names):
+        name2 = name.strip().split("_")[0]
+        names2 = names[0].strip().split('_')[0]
+        return (name2 == names2)
+
+    @staticmethod
+    def add_relative_e(list0):
+        """Adding relative energies compared with lowest-energy structure
+        
+        Argument:
+            list: a list of tuple (idx, name, energy)
+            
+        Return:
+            list of tuple (idx, name, energy, relative_energy)
+        """
+        list0_ = []
+        _, _, e_m = list0[0]
+        for idx_name_e in list0:
+            idx_i, name_i, e_i = idx_name_e
+            e_relative = e_i - e_m
+            list0_.append((idx_i, name_i, e_i, e_relative))
+        return list0_
+
+    def top_k(self, energies, names, mols, k=1):
+        '''
+        Given a group of energy_name_idxes,
+        return the top-k lowest name-energies pairs with idxes as keys.
+        '''
+        assert(len(energies) == len(names))
+        assert(len(energies) == len(mols))
+        assert(len(set(names)) == 1)
+
+        df = pd.DataFrame({"names": names, "energies": energies, "mols": mols})
+        df2 = df.sort_values(by=['energies'])
+        
+        out_mols_ = filter_unique(list(df2["mols"]), self.threshold)
+        if k < len(out_mols_):
+            out_mols = out_mols_[:k]
+        else:
+            out_mols = out_mols_
+
+        if len(out_mols) == 0:
+            name = names[0].split("_")[0].strip()
+            print(f"No structure converged for {name}, try a larger convergence threshold.")
+        else:
+            #Adding relative energies
+            ref_energy = out_mols[0].data['E_tot']
+            for mol in out_mols:
+                my_energy = mol.data['E_tot']
+                rel_energy = float(my_energy) - float(ref_energy)
+                mol.data['E_rel(eV)'] = rel_energy
+        return out_mols
+
+    def top_window(self, energies, names, mols, window=1):
+        '''
+        Given a group of energy_name_idxes,
+        return all (idx, name, e) tuples whose energies are within
+        window (Hatree) from the lowest energy. Unit table is based on: 
+        http://wild.life.nctu.edu.tw/class/common/energy-unit-conv-table.html
+        '''
+        window = (window/eV2kcalmol)  # convert energy window into eV unit
+        assert(len(energies) == len(names))
+        assert(len(energies) == len(mols))
+        assert(window >= 0)
+        assert(len(set(names)) == 1)
+
+
+        df = pd.DataFrame({"names": names, "energies": energies, "mols": mols})
+        df2 = df.sort_values(by=['energies'])
+
+        out_mols_ = filter_unique(list(df2['mols']), self.threshold)
+        out_mols = []
+
+        if len(out_mols_) == 0:
+            name = names[0].split("_")[0].strip()
+            print(f"No structure converged for {name}, try a larger convergence threshold.")
+        else:
+            ref_energy = float(out_mols_[0].data['E_tot'])
+            for mol in out_mols_:
+                my_energy = float(mol.data['E_tot'])
+                rel_energy = my_energy - ref_energy
+                if rel_energy <= window:
+                    mol.data['E_rel(eV)'] = rel_energy
+                    out_mols.append(mol)
+                else:
+                    break
+        return out_mols
+
+    def run(self):
+        """
+        When runs, lowest-energy structure will be stored in out_path folder.
+        """
+        print("Beggin to slelect structures that satisfy the requirements...")
+        results = []
+
+        file_type = guess_file_type(self.input_path)
+        data2 = pybel.readfile(file_type, self.input_path)
+
+        mols = [mol for mol in data2]
+        names = [mol.title.strip() for mol in mols]
+        energies = [float(mol.data['E_tot']) for mol in mols]
+        assert(len(mols) == len(names))
+        assert(len(mols) == len(energies))
+
+        #Grouping, ranking
+        names2 = map(lambda x: x.strip().split("_")[0].strip(), names)
+        df = pd.DataFrame({"names": names2, "energies": energies, "mols": mols})
+        df2 = df.groupby("names")
+        for group_name in df2.indices:
+            group = df2.get_group(group_name)
+            # Filter similar structures
+            g_mols = list(group.loc[:, 'mols'])
+            g_names = list(group.loc[:, 'names'])
+            g_energies = list(group.loc[:, 'energies'])
+
+            if self.k:
+                top_results = self.top_k(g_energies, g_names,
+                                            g_mols, self.k)
+            elif self.window:
+                top_results = self.top_window(g_energies, g_names,
+                                                g_mols, self.window)
+            else:
+                raise ValueError(('Parameter k or window needs to be '
+                                    'specified. Append "--k=1" if you'
+                                    'only want one structure per SMILES'))
+            results += top_results
+
+        f = pybel.Outputfile('sdf', self.out_path)
+        for mol in results:
+            # Change the energy unit from eV back to Hartree
+            mol.data['E_tot'] = (float(mol.data['E_tot'])/hartree2eV)
+            mol.data['E_rel(kcal/mol)'] = (float(mol.data['E_rel(eV)']) * eV2kcalmol)
+            del mol.data['E_rel(eV)']
+            #Remove _ in the molecule title
+            t = mol.title
+            t_simplified = t.split("_")[0].strip()
+            mol.title = t_simplified
+            f.write(mol)
+        f.close()
+
+
+if __name__ == "__main__":
+    """Utils for select top conformers after running Auto3D"""
+    import argparse
+    import os
+
+
+    parser = argparse.ArgumentParser(description="Utils for select top conformers after running Auto3D")
+    parser.add_argument("path", type=str, help="Input sdf/xyz file")
+    parser.add_argument('--k', type=int, default=False,
+                        help='Outputs the top-k structures for each SMILES.')
+    parser.add_argument('--window', type=float, default=False,
+                        help=('Outputs the structures whose energies are within '
+                            'window (kcal/mol) from the lowest energy'))
+    parser.add_argument("--threshold", type=float, default=0.1, help="RMSD threshold for two conformers to be considered as duplicates")
+
+    args = parser.parse_args()
+    path = args.path
+    k = args.k
+    window = args.window
+    threshold = args.threshold
+
+    dir = os.path.dirname(os.path.abspath(path))
+    basename = os.path.basename(path).split(".")[0] + "_out." +  os.path.basename(path).split(".")[1]
+    out_path = os.path.join(dir, basename)
+
+    rank_engine = ranking(path, out_path, threshold, k=k, window=window)
+    rank_engine.run()
diff --git a/src/Auto3D/utils.py b/src/Auto3D/utils.py
new file mode 100644
index 0000000..3ffbfa2
--- /dev/null
+++ b/src/Auto3D/utils.py
@@ -0,0 +1,708 @@
+#!/usr/bin/env python
+"""
+Providing utilities for the workflow package.
+"""
+import math
+import warnings
+import os
+import sys
+import re
+import glob
+import torch
+import collections
+from collections import defaultdict, OrderedDict
+import shutil
+from openbabel import pybel
+from tqdm import tqdm
+from io import StringIO
+from rdkit import Chem
+from rdkit.Chem.rdMolDescriptors import CalcNumAtomStereoCenters
+from rdkit.Chem.rdMolDescriptors import CalcNumUnspecifiedAtomStereoCenters
+
+
+def guess_file_type(filename):
+    """Returns the extension for the filename"""
+    assert '.' in filename
+    return os.path.splitext(filename)[1][1:]
+
+def check_input(args):
+    """
+    Check the input file and give recommendations.
+
+    Arguments:
+        args: Arguments to auto3d.
+
+    Returns:
+        This function checks the format of the input file, the properties for
+        each SMILES in the input file.
+    """
+    print("Checking input file...")
+    ANI_elements = {1, 6, 7, 8, 9, 16, 17}
+    ANI = True
+    # Check --use_gpu
+    gpu_flag = args.use_gpu
+    if gpu_flag:
+        if torch.cuda.is_available() == False:
+            sys.exit("No cuda device was detected. Please set --use_gpu=False.")
+    # Check the availability of omega
+    if "OE_LICENSE" not in os.environ:
+        warnings.warn("OpenEye software license is not detected. Please use RDKit for your program where applicable.")
+        isomer_engine = args.isomer_engine
+        if isomer_engine == "omega":
+            sys.exit("Omega is used as isomer engine, but OE_LICENSE is not detected.")
+
+    # Check the installation for open toolkits, torchani
+    if args.isomer_engine == "omega":
+        try:
+            from openeye import oechem
+        except:
+            sys.exit("Omega is used as isomer engine, but openeye toolkits are not installed.")
+    
+    if args.optimizing_engine == "ANI2x":
+        try:
+            import torchani
+        except:
+            sys.exit("ANI2x is used as optimizing engine, but TorchANI is not installed.")
+
+    if args.optimizing_engine == "ANI2xt":
+        try:
+            from torchani.repulsion import StandaloneRepulsionCalculator
+        except:
+            sys.exit("ANI2xt is used as optimizing engine, but TorchANI with repulsion calculator is not installed.")
+
+    if int(args.opt_steps) < 10:
+        sys.exit(f"Number of optimization steps cannot be smaller than 10, but received {args.opt_steps}")
+
+    # Check the input format
+    smiles_all = []
+    with open(args.path, 'r') as f:
+        data = f.readlines()
+    for line in data:
+        smiles, id = tuple(line.strip().split())
+        assert len(smiles) > 0, \
+               "Empty SMILES string"
+        assert len(id) > 0, \
+               "Empty ID"
+        smiles_all.append(smiles)
+    print(f"\tThere are {len(data)} SMILES in the input file {args.path}. ")
+    print("\tAll SMILES and IDs are valid.")
+
+
+    # Check the properties of molecules
+    only_aimnet_smiles = []
+    for smiles in smiles_all:
+        mol = Chem.MolFromSmiles(smiles)
+        charge = Chem.rdmolops.GetFormalCharge(mol)
+        elements = set([a.GetAtomicNum() for a in mol.GetAtoms()])
+        if ((elements.issubset(ANI_elements) is False) or (charge != 0)):
+            ANI = False
+            only_aimnet_smiles.append(smiles)
+
+    print("Suggestions for choosing isomer_engine and optimizing_engine: ")
+    if ANI:
+        print("\tIsomer engine options: RDKit and Omega.\n"
+              "\tOptimizing engine options: ANI2x, ANI_2xt and AIMNET.")
+    else:
+        print("\tIsomer engine options: RDKit and Omega.\n"
+              "\tOptimizing engine options: AIMNET.")
+        optimizing_engine = args.optimizing_engine
+        if optimizing_engine != "AIMNET":
+            sys.exit(f"Only AIMNET can handle: {only_aimnet_smiles}, but {optimizing_engine} was parsed to Auto3D.")
+
+
+class NullIO(StringIO):
+    """
+    Place holder for a clean terminal
+    """
+    def write(self, txt):
+        pass
+
+def countXYZ(xyz):
+    """Counting the number of structures in XYZ file"""
+    c = 0
+    for _ in pybel.readfile('xyz', xyz):
+        c += 1
+    return c
+
+def countSDF(sdf):
+    """Counting the number of structures in SDF file"""
+    mols = pybel.readfile('sdf', sdf)
+    mols2 = [mol for mol in mols]
+    c = len(mols2)
+    return c
+
+def hash_enumerated_smi_IDs(smi, out):
+    '''
+    Writes all SMILES with hashed IDs into smiles_enumerated_hashed.smi
+
+    Arguments:
+        smi: a .smi File path
+        out: the path for the new .smi file where original IDs are hashed.
+    Returns:
+        writes all SMILES with hashed IDs into smiles_enumerated_hashed.smi
+    '''
+    with open(smi, 'r') as f:
+        data = f.readlines()
+
+    dict0 = {}
+    for line in data:
+        smiles, id = line.strip().split()
+        while (id in dict0.keys()):
+            id += '_0'
+        dict0[id] = smiles
+
+    dict0 = collections.OrderedDict(sorted(dict0.items()))
+
+    new_smi = out
+    with open(new_smi, 'w+') as f:
+        for id, smiles in dict0.items():
+            molecule = smiles.strip() + ' ' + id.strip() + '\n'
+            f.write(molecule)
+
+def hash_taut_smi(smi, out):
+    '''
+    Writes all SMILES with hashed IDs for tautomers
+
+    Arguments:
+        smi: a .smi File path
+        out: the path for the new .smi file where original IDs are hashed.
+    '''
+    with open(smi, 'r') as f:
+        data = f.readlines()
+
+    dict0 = {}
+    for line in data:
+        smiles, id = line.strip().split()
+        c = 1
+        id_ = id
+        while (id_ in dict0.keys()):
+            id_ = id + f"taut{c}"
+            c += 1
+        dict0[id_] = smiles
+
+    dict0 = collections.OrderedDict(sorted(dict0.items()))
+
+    with open(out, 'w+') as f:
+        for id, smiles in dict0.items():
+            molecule = smiles.strip() + ' ' + id.strip() + '\n'
+            f.write(molecule)
+
+
+def sort_helper(mol):
+    '''
+    Converting the pybel mol object into an array containing
+    atomic number and coordinates.
+    '''
+    num2symbol = {1: 'H', 6: 'C', 8: 'O', 7: 'N', 9: 'F', 16: 'S', 17: 'Cl'}
+    unit = []
+    for atom in mol.atoms:
+        num = atom.atomicnum
+        symbol = num2symbol[num]
+        x, y, z = atom.coords
+        symbol_coord = [symbol, x, y, z]
+        unit.append(symbol_coord)
+    return unit
+
+
+def sort_enumerated_xyz(xyz, out):
+    '''
+    Sort an xyz file based on the IDs of SMILES in the input file.
+
+    Arguments:
+        xyz: the input xyz file.
+        out: the path for the sorted xyz file
+    Returns:
+        writes the sorted molecules into out.
+    '''
+    dict0 = {}
+    mols = pybel.readfile('xyz', xyz)
+    for mol in tqdm(mols):
+        id = str(mol).strip().split('\t')[1].strip()
+        unit = sort_helper(mol)
+        dict0[id] = unit
+    dict0 = collections.OrderedDict(sorted(dict0.items()))
+
+    # new_xyz = out[:-4] + '_ordered.xyz'
+    with open(out, 'w+') as f:
+        for id, unit in dict0.items():
+            length = str(len(unit)) + '\n'
+            id = id.strip() + '\n'
+            f.write(length)
+            f.write(id)
+            for line in unit:
+                s, x, y, z = line
+                atom = (str(s).strip() + '\t' + str(x) + '\t' +
+                        str(y) + '\t' + str(z) + '\n')
+                f.write(atom)
+
+def sort_enumerated_sdf(sdf, out):
+    """
+    Sort an SDF file based on the IDs of SMILES in the input file.
+
+    Arguments:
+        sdf: the input SDF file.
+        out: the path for the sorted SDF file
+    Returns:
+        writes the sorted molecules into out.
+    """
+    dict0 = {}
+    
+    mols = pybel.readfile('sdf', sdf)
+    for mol in mols:
+        idx = str(mol).split('\t')[1].strip().split(' ')[0].strip()
+        mol.data['ID'] = idx
+        dict0[idx] = mol
+    dict0 = collections.OrderedDict(sorted(dict0.items()))
+    
+    f = pybel.Outputfile('sdf', out)
+    for idx, mol in dict0.items():
+        f.write(mol)
+    f.close()
+
+def combine_xyz(in_folder, out_path):
+    """
+    Combining all xyz files in the in_folder into a single xyz file (out_path).
+
+    Arguemnts:
+        in_folder: a folder contains all xyz files.
+        out_path: a path of xyz file to store every structure in the in_folder
+
+    Returns:
+        Combining all xyz files in the in_folder into out_path.
+    """
+    file_paths = os.path.join(f"{in_folder}/*.xyz")
+    files = glob.glob(file_paths)
+    # print(f'There are {len(files)} single xyz files...')
+
+    results = []
+    for file in files:
+        with open(file, 'r') as f:
+            data = f.readlines()
+        assert(len(data) == (int(data[0]) + 2))
+        results += data
+
+    with open(out_path, 'w+') as f:
+        for line in results:
+            f.write(line)
+    # print(f'Combined in a singl file {out_path}!')
+
+def combine_smi(smies, out):
+    """Combine smi files into a single file"""
+    data = []
+    for smi in smies:
+        with open(smi, 'r') as f:
+            datai = f.readlines()
+        data += datai
+    data = list(set(data))
+    with open(out, 'w+') as f2:
+        for line in data:
+            if not line.isspace():
+                f2.write((line.strip() + '\n'))
+
+
+def housekeeping_helper(folder, file):
+    basename = os.path.basename(file)
+    new_name = os.path.join(folder, basename)
+    shutil.move(file, new_name)
+
+
+def housekeeping(job_name, folder, optimized_structures):
+    """
+    Moving all meta data into a folder
+
+    Arguments:
+        folder: a folder name to contain all meta data
+        out: the resulting SDF output
+    Returns:
+        whe the function is called, it moves all meta data into a folder.
+    """
+    paths = os.path.join(job_name, '*')
+    files = glob.glob(paths)
+    for file in files:
+        if file != optimized_structures:
+            shutil.move(file, folder)
+
+    try:
+        paths1 = os.path.join('', 'oeomega_*')
+        files1 = glob.glob(paths1)
+        paths2 = os.path.join('', 'flipper_*')
+        files2 = glob.glob(paths2)
+        files = files1 + files2
+        for file in files:
+            shutil.move(file, folder)
+    except:
+        pass
+
+def enantiomer(l1,  l2):
+    """Check if two lists of stereo centers are enantiomers"""
+    indicator = True
+    assert (len(l1) == len(l2))
+    for i in range(len(l1)):
+        tp1 = l1[i]
+        tp2 = l2[i]
+        idx1, stereo1 = tp1
+        idx2, stereo2 = tp2
+        assert(idx1 == idx2)
+        if (stereo1 == stereo2):
+            indicator = False
+            return indicator
+    return indicator
+            
+
+def enantiomer_helper(smiles):
+    """get non-enantiomer SMILES from given smiles"""
+    mols = [Chem.MolFromSmiles(smi) for smi in smiles]
+    stereo_centers = [Chem.FindMolChiralCenters(mol, useLegacyImplementation=False) for mol in mols]
+    non_enantiomers = []
+    non_centers = []
+    for i in range(len(stereo_centers)):
+        smi = smiles[i]
+        stereo = stereo_centers[i]
+        indicator = True
+        for j in range(len(non_centers)):
+            stereo_j = non_centers[j]
+            if enantiomer(stereo_j, stereo):
+                indicator = False
+        if indicator:
+            non_centers.append(stereo)
+            non_enantiomers.append(smi)
+    return non_enantiomers
+            
+
+def remove_enantiomers(inpath, out):
+    """Removing enantiomers for the input file
+    Arguments:
+        inpath: input smi
+        output: output smi
+    """
+    with open(inpath, 'r') as f:
+        data = f.readlines()
+    
+    smiles = defaultdict(lambda: [])
+    for line in data:
+        vals = line.split()
+        smi, name = vals[0].strip(), vals[1].strip().split("_")[0].strip()
+        smiles[name].append(smi)
+
+    for key, values in smiles.items():
+        try:
+            new_values = enantiomer_helper(values)
+        except:
+            new_values = values
+            print(f"Enantiomers not removed for {key}")
+        smiles[key] = new_values
+        
+    with open(out, 'w+') as f:
+        for key, val in smiles.items():
+            for i in range(len(val)):
+                new_key = key + "_" + str(i)
+                line = val[i].strip() + ' ' + new_key + '\n'
+                f.write(line)
+    return smiles
+
+
+def filter_unique(mols, crit=0.3):
+    """Remove structures that are very similar.
+       Remove unconverged structures.
+    
+    Arguments:
+        mols: pybel mol objects
+    Returns:
+        unique_mols: unique molecules
+    """
+
+    #Remove unconverged structures
+    mols_ = []
+    for mol in mols:
+        convergence_flag = str(mol.data['Converged']).lower() == "true"
+        if convergence_flag:
+            mols_.append(mol)
+    mols = mols_
+
+    #Remove similar structures
+    unique_mols = []
+    aligner = pybel.ob.OBAlign()
+    for mol_i in mols:
+        aligner.SetRefMol(mol_i.OBMol)
+        unique = True
+        for mol_j in unique_mols:
+            aligner.SetTargetMol(mol_j.OBMol)
+            aligner.Align()
+            rmsd = aligner.GetRMSD()
+            if rmsd < crit:
+                unique = False
+                break
+        if unique:
+            unique_mols.append(mol_i)
+    return unique_mols
+
+
+def unique_conformers(files, crit=0.5):
+    """Removing conformers whose RMSD is within crit
+    
+    Arguments:
+        files: sdf files
+    """
+    unique_files = []
+    duplicate_files = []
+    aligner = pybel.ob.OBAlign()
+    for file in files:
+        mol = next(pybel.readfile("sdf", file))
+        aligner.SetRefMol(mol.OBMol)
+        unique = True
+        for f in unique_files:
+            mol_j = next(pybel.readfile("sdf", f))
+            aligner.SetTargetMol(mol_j.OBMol)
+            aligner.Align()
+            rmsd = aligner.GetRMSD()
+            if rmsd < crit:
+                unique = False
+                break
+        if unique:
+            unique_files.append(file)
+        else:
+            duplicate_files.append(file)
+    c = len(unique_files) + len(duplicate_files)
+    assert(c == len(files))
+    for file in duplicate_files:
+        os.remove(file)
+
+def sort_output(mols, out):
+    """Sort molecules based on energies
+    
+    Arguments:
+        mols: a list of pybel mol objects
+        out: a SDF file to store the molecules
+    """
+    l = []
+    for mol in mols:
+        name = mol.data['ID'].split('_')[0].strip()
+        e = mol.data['E_tot']
+        e_rel = mol.data['E_relative']
+
+        l.append((name, e_rel, mol))
+
+    l = sorted(l)
+
+    f = pybel.Outputfile('sdf', out)
+    for n_er_m in l:
+        name, e_relative, mol = n_er_m
+        f.write(mol)
+    f.close()
+
+def no_enantiomer_helper(info1, info2):
+    """Return true if info1 and info2 are enantiomers"""
+    assert (len(info1) == len(info2))
+    for i in range(len(info1)):
+        if info1[i].strip() == info2[i].strip():
+            return False
+    return True
+
+def get_stereo_info(smi):
+    "Return a dictionary of @@ or @  in smi"
+    dct = {}
+    regex1 = re.compile("[^@]@[^@]")
+    regex2 = re.compile("@@")
+
+    # match @
+    for m in regex1.finditer(smi):
+        dct[m.start()+1] = '@'
+
+    #match  @@
+    for m in regex2.finditer(smi):
+        dct[m.start()] = "@@"
+
+    dct2 = OrderedDict(sorted(dct.items()))
+    return dct2
+
+
+def no_enantiomer(smi, smiles):
+    """Return True if there is no enantiomer for smi in smiles"""
+
+    stereo_infoi = list(get_stereo_info(smi).values())
+    for i in range(len(smiles)):
+        tar = smiles[i]
+        if tar != smi:
+            stereo_infoj = list(get_stereo_info(tar).values())
+            if no_enantiomer_helper(stereo_infoi, stereo_infoj):
+                return False
+    return True
+
+def create_enantiomer(smi):
+    """Create an enantiomer SMILES for input smi"""
+    stereo_info = get_stereo_info(smi)
+    new_smi = ""
+    # for key in stereo_info.keys():
+    #     val = stereo_info[key]
+    #     if val == '@':
+    keys = list(stereo_info.keys())
+    if len(keys) == 1:
+        key = keys[0]
+        val = stereo_info[key]
+        if val == "@":
+            new_smi += smi[:key]
+            new_smi += "@@"
+            new_smi += smi[(key+1):]
+        elif val == "@@":
+            new_smi += smi[:key]
+            new_smi += "@"
+            new_smi += smi[(key+2):]
+        else:
+            raise ValueError("Invalid %s" % smi)
+        return new_smi
+
+    for i in range(len(keys)):
+        if i == 0:
+            key = keys[i]
+            new_smi += smi[:key]
+        else:
+            key1 = keys[i-1]
+            key2 = keys[i]
+            val1 = stereo_info[key1]
+            if val1 == "@":
+                new_smi += "@@"
+                new_smi += smi[int(key1+1): key2]
+            elif val1 == "@@":
+                new_smi += "@"
+                new_smi += smi[int(key1+2): key2]
+    val2 = stereo_info[key2]
+    if val2 == "@":
+        new_smi += "@@"
+        new_smi += smi[int(key2+1):]
+    elif val2 == "@@":
+        new_smi += "@"
+        new_smi += smi[int(key2+2):]
+    return new_smi
+
+def check_value(n):
+    """Return True if n is a power of 2. 2^-2 and 2^-1 should be accessible,
+    because not all stereo centers can be enumerated. For example: CC12CCC(C1)C(C)(C)C2O"""
+    
+    power = math.log(n, 2)
+    decimal, integer = math.modf(power)
+    i = abs(power - integer)
+    if (i < 0.0001):
+        return True
+    return False
+
+def amend_configuration(smis):
+    """Adding the missing configurations. 
+    Example: N=C1OC(CN2CC(C)OC(C)C2)CN1"""
+
+    with open(smis, 'r') as f:
+        data = f.readlines()
+    dct = defaultdict(lambda: [])
+    for line in data:
+        smi, idx = tuple(line.strip().split())
+        idx = idx.split("_")[0].strip()
+        dct[idx].append(smi)
+    
+    for key in dct.keys():
+        value = dct[key]
+        smi = value[0]
+        mol = Chem.MolFromSmiles(smi)
+        num_centers = CalcNumAtomStereoCenters(mol)
+        num_unspecified_centers = CalcNumUnspecifiedAtomStereoCenters(mol)
+        # num_configurations = 2 ** num_unspecified_centers
+        num_configurations = 2 ** num_centers
+        num = len(value)/num_configurations
+
+        if ((check_value(num) == False) and ("@" in smi)):  # Missing configurations
+            try:
+                new_value = []
+                for val in value:
+                    if no_enantiomer(val, value):
+                        new_val = create_enantiomer(val)
+                        new_value.append(new_val)
+                value += new_value
+                
+                # assert 
+                new_num = len(value)/num_configurations
+                assert (check_value(new_num) == True)
+                dct[key] = value
+            except:
+                print(f"Stereo centers for {key} are not fully enumerated.")
+    return dct
+
+def amend_configuration_w(smi):
+    """Write the output from dictionary"""
+    dct = amend_configuration(smi)
+    with open(smi, 'w+') as f:
+        for key in dct.keys():
+            val = dct[key]
+            for i, smi in enumerate(val):
+                idx = str(key).strip() + "_" + str(i+1)
+                line = smi + ' ' + idx + '\n'
+                f.write(line)
+
+
+def is_macrocycle(smi, size=10):
+    """Check is a SMIELS is contains a macrocycle part (a 10-membered or larger 
+    ring regardless of their aromaticity and hetero atoms content)"""
+    mol = Chem.MolFromSmiles(smi)
+    ring = mol.GetRingInfo()
+    ring_bonds = ring.BondRings()
+    for bonds in ring_bonds:
+        if len(bonds) >= 10:
+            return True
+    return False
+
+
+def split_smi(smi):
+    """Split an input .smi file into two files:
+    one contains small SMILES, the other contain macrocycle smiles
+    """
+    #Prepare out file path
+    dir = os.path.dirname(os.path.realpath(smi))
+    basename = os.path.basename(smi)
+    normal_name = basename.split(".")[0].strip() + "_normal.smi"
+    macrocycle_name = basename.split(".")[0].strip() + "_macrocycle.smi"
+    normal_path = os.path.join(dir, normal_name)
+    macrocycle_path = os.path.join(dir, macrocycle_name)
+
+    normal = []
+    macrocycle = []
+    with open(smi, "r") as f:
+        data = f.readlines()
+    for line in tqdm(data):
+        smi_idx = line.strip().split()
+        smi = smi_idx[0]
+        if is_macrocycle(smi):
+            macrocycle.append(line)
+        else:
+            normal.append(line)
+    
+    l_ = len(normal) + len(macrocycle)
+    l = len(data)
+    assert(l == l_)
+
+    with open(normal_path, "w+") as f:
+        for line in normal:
+            f.write(line)
+
+    with open(macrocycle_path, "w+") as f:
+        for line in macrocycle:
+            f.write(line)
+
+
+class my_name_space(dict):
+    """A modified dictionary whose keys can be accessed via dot. This dict is
+    similar to a NameSpace object from parser.
+    """
+    def __getattr__(self, key):
+        try:
+            return self[key]
+        except KeyError as k:
+            raise AttributeError(k)
+
+    def __setattr__(self, key, value):
+        self[key] = value
+
+    def __delattr__(self, key):
+        try:
+            del self[key]
+        except KeyError as k:
+            raise AttributeError(k)
+
+    def __repr__(self):
+        return '<my_name_space ' + dict.__repr__(self) + '>'


Test patch:
diff --git a/src/Auto3D/input/test.smi b/src/Auto3D/input/test.smi
new file mode 100644
index 0000000..bd439c7
--- /dev/null
+++ b/src/Auto3D/input/test.smi
@@ -0,0 +1,6 @@
+CC1(C2=CC=CC=C2)C=C1 smi1
+O=C(C)OCC1(C2=CC=CC=C2)C=C1 smi2
+C1=C\C=C/C=C\C=C/1 smi3
+C1(C2=CC=CC=C2)CC/C=C\CCC1C3=CC=CC=C3 smi4
+C1CCOC1 smi5
+CC(=Cc1ccc(cc1)O)C(=O)[O-] smi6
\ No newline at end of file
diff --git a/tests/files/b973c_ANI2xt_E.sdf b/tests/files/b973c_ANI2xt_E.sdf
deleted file mode 100644
index 61c6121..0000000
--- a/tests/files/b973c_ANI2xt_E.sdf
+++ /dev/null
@@ -1,143 +0,0 @@
-/Users/liu5/Documents/tautomer/dG/dG3/b973c/analyze/817-2-473.xyz
- OpenBabel04272223063D
-
- 20 19  0  0  0  0  0  0  0  0999 V2000
-    1.0833   -0.6247   -0.2391 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.8271    0.6201   -0.6549 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.5840   -2.1116    0.9509 C   0  0  0  0  0  0  0  0  0  0  0  0
-    2.5883    3.8140   -0.6118 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.1483   -1.7836   -0.3495 C   0  0  0  0  0  0  0  0  0  0  0  0
-    1.4111    2.9365   -0.9446 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.3132    1.0368   -1.2776 O   0  0  0  0  0  0  0  0  0  0  0  0
-    1.7275    1.6120   -0.4878 O   0  0  0  0  0  0  0  0  0  0  0  0
-    2.0335   -0.7789    0.2514 H   0  0  0  0  0  0  0  0  0  0  0  0
-    0.1266   -2.3245    1.7484 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.2324   -2.9820    0.8384 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.1946   -1.2703    1.2775 H   0  0  0  0  0  0  0  0  0  0  0  0
-    2.7717    3.8316    0.4607 H   0  0  0  0  0  0  0  0  0  0  0  0
-    3.4906    3.4612   -1.1071 H   0  0  0  0  0  0  0  0  0  0  0  0
-    2.3924    4.8333   -0.9425 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.5930   -1.6280   -1.1438 H   0  0  0  0  0  0  0  0  0  0  0  0
-    0.7059   -2.6657   -0.6738 H   0  0  0  0  0  0  0  0  0  0  0  0
-    1.2224    2.9220   -2.0190 H   0  0  0  0  0  0  0  0  0  0  0  0
-    0.5040    3.2908   -0.4525 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.9011    0.2769   -1.3594 H   0  0  0  0  0  0  0  0  0  0  0  0
-  1  9  1  0  0  0  0
-  2  8  1  0  0  0  0
-  2  1  2  0  0  0  0
-  3 12  1  0  0  0  0
-  3 10  1  0  0  0  0
-  4 13  1  0  0  0  0
-  5  1  1  0  0  0  0
-  5  3  1  0  0  0  0
-  6  4  1  0  0  0  0
-  6  8  1  0  0  0  0
-  6 19  1  0  0  0  0
-  7  2  1  0  0  0  0
- 11  3  1  0  0  0  0
- 14  4  1  0  0  0  0
- 15  4  1  0  0  0  0
- 16  5  1  0  0  0  0
- 17  5  1  0  0  0  0
- 18  6  1  0  0  0  0
- 20  7  1  0  0  0  0
-M  END
->  <SPE(Hatree)>
--386.110819707842
-
->  <Enthalpy(Hatree)>
--385.91602271
-
->  <G(Hatree)>
--386.00075446
-
->  <Runtime(m)>
-3
-
->  <ID>
-817-2-473
-
->  <E_hatree>
--386.11389232197996
-
-$$$$
-/Users/liu5/Documents/tautomer/dG/dG3/b973c/analyze/510-2-443.xyz
- OpenBabel04272223063D
-
- 27 29  0  0  0  0  0  0  0  0999 V2000
-   -0.9627    3.7501   -0.3367 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.7986    4.2297    1.2298 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.1480    5.1066   -0.5303 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.6062    5.5874    1.0364 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.0175    3.2909    0.5489 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.3679    6.0228    0.1565 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.8853   -1.4005   -1.1159 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.0504   -1.5246    0.0822 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.2291    0.8631    0.0391 C   0  0  0  0  0  0  0  0  0  0  0  0
-    1.1949    1.4613    1.7786 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.9970   -0.0521    1.8271 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.3163   -0.2868   -1.6135 N   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.9781    0.9027   -1.0199 N   0  0  0  0  0  0  0  0  0  0  0  0
-    0.2174    1.9246    0.7926 N   0  0  0  0  0  0  0  0  0  0  0  0
-    0.2544   -0.2980    0.6017 N   0  0  0  0  0  0  0  0  0  0  0  0
-    0.3235   -2.5944    0.5568 O   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.2189   -2.5581   -1.7045 O   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.6120    7.7420   -0.0926 Cl  0  0  0  0  0  0  0  0  0  0  0  0
-   -1.5605    3.0426   -0.8852 H   0  0  0  0  0  0  0  0  0  0  0  0
-    1.5721    3.9207    1.9156 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.9061    5.4494   -1.2187 H   0  0  0  0  0  0  0  0  0  0  0  0
-    1.2182    6.3006    1.5681 H   0  0  0  0  0  0  0  0  0  0  0  0
-    1.0099    1.9220    2.7451 H   0  0  0  0  0  0  0  0  0  0  0  0
-    2.2075    1.7203    1.4614 H   0  0  0  0  0  0  0  0  0  0  0  0
-    1.9271   -0.6111    1.8276 H   0  0  0  0  0  0  0  0  0  0  0  0
-    0.4038   -0.3724    2.6821 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.8060   -3.2572   -1.1667 H   0  0  0  0  0  0  0  0  0  0  0  0
-  1  5  1  0  0  0  0
-  2 20  1  0  0  0  0
-  3  1  2  0  0  0  0
-  3  6  1  0  0  0  0
-  4  2  1  0  0  0  0
-  4 22  1  0  0  0  0
-  5 14  1  0  0  0  0
-  5  2  2  0  0  0  0
-  6  4  2  0  0  0  0
-  7  8  1  0  0  0  0
-  8 16  2  0  0  0  0
-  8 15  1  0  0  0  0
-  9 15  1  0  0  0  0
-  9 14  1  0  0  0  0
- 10 11  1  0  0  0  0
- 10 23  1  0  0  0  0
- 11 25  1  0  0  0  0
- 11 26  1  0  0  0  0
- 12  7  2  0  0  0  0
- 12 13  1  0  0  0  0
- 13  9  2  0  0  0  0
- 14 10  1  0  0  0  0
- 15 11  1  0  0  0  0
- 17 27  1  0  0  0  0
- 17  7  1  0  0  0  0
- 18  6  1  0  0  0  0
- 19  1  1  0  0  0  0
- 21  3  1  0  0  0  0
- 24 10  1  0  0  0  0
-M  END
->  <SPE(Hatree)>
--1253.811921710955
-
->  <Enthalpy(Hatree)>
--1253.58519689
-
->  <G(Hatree)>
--1253.68784878
-
->  <Runtime(m)>
-18
-
->  <ID>
-510-2-443
-
->  <E_hatree>
--1253.8134542668342
-
-$$$$
diff --git a/tests/files/wb97x_dz_ANI2x_E.sdf b/tests/files/wb97x_dz_ANI2x_E.sdf
deleted file mode 100644
index 461ebe7..0000000
--- a/tests/files/wb97x_dz_ANI2x_E.sdf
+++ /dev/null
@@ -1,143 +0,0 @@
-/Users/liu5/Documents/tautomer/dG/dG3/DFT_output_analyze2/817-2-473.xyz
- OpenBabel04272223063D
-
- 20 19  0  0  0  0  0  0  0  0999 V2000
-    1.0662   -0.6046   -0.2205 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.8149    0.6422   -0.6377 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.5667   -2.1485    0.9594 C   0  0  0  0  0  0  0  0  0  0  0  0
-    2.5912    3.8169   -0.6111 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.1288   -1.7733   -0.3538 C   0  0  0  0  0  0  0  0  0  0  0  0
-    1.4154    2.9308   -0.9652 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.3072    1.0636   -1.2787 O   0  0  0  0  0  0  0  0  0  0  0  0
-    1.7077    1.6333   -0.4471 O   0  0  0  0  0  0  0  0  0  0  0  0
-    2.0141   -0.7590    0.2850 H   0  0  0  0  0  0  0  0  0  0  0  0
-    0.1738   -2.3808    1.7325 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.2148   -3.0241    0.8369 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.1755   -1.3160    1.3291 H   0  0  0  0  0  0  0  0  0  0  0  0
-    2.7247    3.8626    0.4738 H   0  0  0  0  0  0  0  0  0  0  0  0
-    3.5124    3.4286   -1.0555 H   0  0  0  0  0  0  0  0  0  0  0  0
-    2.4234    4.8319   -0.9855 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.6372   -1.5887   -1.1242 H   0  0  0  0  0  0  0  0  0  0  0  0
-    0.6855   -2.6435   -0.7254 H   0  0  0  0  0  0  0  0  0  0  0  0
-    1.2725    2.8769   -2.0506 H   0  0  0  0  0  0  0  0  0  0  0  0
-    0.4864    3.3094   -0.5239 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.9014    0.3083   -1.3735 H   0  0  0  0  0  0  0  0  0  0  0  0
-  1  9  1  0  0  0  0
-  2  8  1  0  0  0  0
-  2  1  2  0  0  0  0
-  3 12  1  0  0  0  0
-  3 10  1  0  0  0  0
-  4 13  1  0  0  0  0
-  5  1  1  0  0  0  0
-  5  3  1  0  0  0  0
-  6  4  1  0  0  0  0
-  6 19  1  0  0  0  0
-  6  8  1  0  0  0  0
-  7  2  1  0  0  0  0
- 11  3  1  0  0  0  0
- 14  4  1  0  0  0  0
- 15  4  1  0  0  0  0
- 16  5  1  0  0  0  0
- 17  5  1  0  0  0  0
- 18  6  1  0  0  0  0
- 20  7  1  0  0  0  0
-M  END
->  <SPE(Hatree)>
--386.177798664469
-
->  <Enthalpy(Hatree)>
--385.97801882
-
->  <G(Hatree)>
--386.0620377
-
->  <Runtime(m)>
-11
-
->  <ID>
-817-2-473
-
->  <E_hatree>
--386.1793964525865
-
-$$$$
-/Users/liu5/Documents/tautomer/dG/dG3/DFT_output_analyze2/510-2-443.xyz
- OpenBabel04272223063D
-
- 27 29  0  0  0  0  0  0  0  0999 V2000
-   -0.9493    3.7513   -0.3628 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.7709    4.2300    1.2607 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.1236    5.1135   -0.5636 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.5899    5.5930    1.0576 C   0  0  0  0  0  0  0  0  0  0  0  0
-    0.0053    3.2937    0.5554 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.3577    6.0311    0.1450 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.9174   -1.3978   -1.1115 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.0446   -1.5210    0.0770 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.2534    0.8603    0.0514 C   0  0  0  0  0  0  0  0  0  0  0  0
-    1.1911    1.4641    1.7851 C   0  0  0  0  0  0  0  0  0  0  0  0
-    1.0350   -0.0643    1.8043 C   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.3661   -0.2931   -1.5857 N   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.0242    0.9021   -0.9817 N   0  0  0  0  0  0  0  0  0  0  0  0
-    0.2030    1.9202    0.8080 N   0  0  0  0  0  0  0  0  0  0  0  0
-    0.2601   -0.3010    0.5986 N   0  0  0  0  0  0  0  0  0  0  0  0
-    0.3452   -2.5900    0.5234 O   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.2385   -2.5639   -1.6833 O   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.5857    7.7373   -0.1142 Cl  0  0  0  0  0  0  0  0  0  0  0  0
-   -1.5384    3.0387   -0.9254 H   0  0  0  0  0  0  0  0  0  0  0  0
-    1.5231    3.9149    1.9760 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -1.8627    5.4624   -1.2782 H   0  0  0  0  0  0  0  0  0  0  0  0
-    1.1893    6.3110    1.6084 H   0  0  0  0  0  0  0  0  0  0  0  0
-    0.9887    1.9010    2.7671 H   0  0  0  0  0  0  0  0  0  0  0  0
-    2.2020    1.7606    1.4760 H   0  0  0  0  0  0  0  0  0  0  0  0
-    1.9863   -0.5976    1.7549 H   0  0  0  0  0  0  0  0  0  0  0  0
-    0.4783   -0.4220    2.6757 H   0  0  0  0  0  0  0  0  0  0  0  0
-   -0.7994   -3.2531   -1.1526 H   0  0  0  0  0  0  0  0  0  0  0  0
-  1  5  1  0  0  0  0
-  2 20  1  0  0  0  0
-  3  1  2  0  0  0  0
-  3  6  1  0  0  0  0
-  4  2  1  0  0  0  0
-  4 22  1  0  0  0  0
-  5 14  1  0  0  0  0
-  5  2  2  0  0  0  0
-  6  4  2  0  0  0  0
-  7  8  1  0  0  0  0
-  8 16  2  0  0  0  0
-  8 15  1  0  0  0  0
-  9 15  1  0  0  0  0
-  9 14  1  0  0  0  0
- 10 11  1  0  0  0  0
- 10 23  1  0  0  0  0
- 11 26  1  0  0  0  0
- 12  7  2  0  0  0  0
- 12 13  1  0  0  0  0
- 13  9  2  0  0  0  0
- 14 10  1  0  0  0  0
- 15 11  1  0  0  0  0
- 17 27  1  0  0  0  0
- 17  7  1  0  0  0  0
- 18  6  1  0  0  0  0
- 19  1  1  0  0  0  0
- 21  3  1  0  0  0  0
- 24 10  1  0  0  0  0
- 25 11  1  0  0  0  0
-M  END
->  <SPE(Hatree)>
--1254.007233026973
-
->  <Enthalpy(Hatree)>
--1253.7738878
-
->  <G(Hatree)>
--1253.87479791
-
->  <Runtime(m)>
-111
-
->  <ID>
-510-2-443
-
->  <E_hatree>
--1254.017388286501
-
-$$$$