diff --git a/geometric/data/hcn_minimized.tcin b/geometric/data/hcn_minimized.tcin
new file mode 100644
index 00000000..bb451b6b
--- /dev/null
+++ b/geometric/data/hcn_minimized.tcin
@@ -0,0 +1,6 @@
+run gradient
+coordinates hcn_minimized.xyz
+method rhf
+basis 3-21g
+charge 0
+spinmult 1
diff --git a/geometric/data/hcn_minimized.xyz b/geometric/data/hcn_minimized.xyz
new file mode 100644
index 00000000..705457a5
--- /dev/null
+++ b/geometric/data/hcn_minimized.xyz
@@ -0,0 +1,5 @@
+3    
+Iteration 0 Energy -92.35408417
+C        0.0000000000    0.0000000000    0.0000000000
+N        0.0000000000    0.0000000000    1.1371177000
+H        0.0000000000    0.0000000000   -1.0502281000
diff --git a/geometric/data/hcn_minimized_ts.qcin b/geometric/data/hcn_minimized_ts.qcin
new file mode 100644
index 00000000..5476e2a5
--- /dev/null
+++ b/geometric/data/hcn_minimized_ts.qcin
@@ -0,0 +1,12 @@
+$molecule
+0 1
+C       -0.1088783634   -0.6365101639    0.0043221742
+N       -0.6393457902    0.4205365638    0.0052498438
+H        0.7532976101    0.2173493463   -0.0090384631
+$end
+
+$rem
+jobtype force
+exchange hf
+basis 3-21g
+$end
diff --git a/geometric/data/hcn_neb_input.qcin b/geometric/data/hcn_neb_input.qcin
new file mode 100644
index 00000000..53539c75
--- /dev/null
+++ b/geometric/data/hcn_neb_input.qcin
@@ -0,0 +1,13 @@
+$molecule
+0 1
+C       -0.1088783634   -0.6365101639    0.0043221742
+N       -0.6393457902    0.4205365638    0.0052498438
+H        0.7532976101    0.2173493463   -0.0090384631
+$end
+
+$rem
+jobtype force
+exchange hf
+basis 3-21g
+$end
+
diff --git a/geometric/data/hcn_neb_input.tcin b/geometric/data/hcn_neb_input.tcin
new file mode 100644
index 00000000..dce588bb
--- /dev/null
+++ b/geometric/data/hcn_neb_input.tcin
@@ -0,0 +1,6 @@
+run gradient
+coordinates hcn_neb_input.xyz
+method rhf
+basis 3-21g
+charge 0
+spinmult 1
diff --git a/geometric/engine.py b/geometric/engine.py
index db009205..69a1b87e 100644
--- a/geometric/engine.py
+++ b/geometric/engine.py
@@ -383,6 +383,11 @@ def __init__(self, molecule, tcin, dirname=None, pdb=None):
             self.scr = self.tcin['scrdir']
         else:
             self.scr = 'scr'
+
+        # Getting method and basis
+        self.method = self.tcin['method']
+        self.basis = self.tcin['basis']
+
         # Always specify the TeraChem scratch folder name
         self.tcin['scrdir'] = self.scr
         # A few notes about the electronic structure method
@@ -1069,7 +1074,8 @@ def load_gaussian_input(self, gaussian_input):
 
                 elif line.startswith('#'):
                     self.route_line = line
-
+                    method_basis = line.split(' ')[1]
+                    self.method, self.basis = method_basis.split('/')
                     found_force_nostep = "force=nostep" in line.lower()
                     found_force = "force=" in line.lower()
                     found_symmetry_none = "symmetry=none" in line.lower()
@@ -1388,7 +1394,10 @@ def load_psi4_input(self, psi4in):
             else:
                 psi4_temp.append(line)
             if "gradient(" in line:
+                self.method = re.findall(r'\((.*?)\)', line)[0].replace("'","").replace('"','')
                 found_gradient = True
+            if "set basis" in line:
+                self.basis = line.split(" ")[-1].replace('\n', '')
         if found_gradient == False:
             raise RuntimeError("Psi4 inputfile %s should have gradient() command." % psi4in)
         self.M = Molecule()
@@ -1413,7 +1422,7 @@ def calc_new(self, coords, dirname):
                     outfile.write(line)
         try:
             # Run Psi4
-            subprocess.check_call('psi4%s input.dat' % self.nt(), cwd=dirname, shell=True)
+            subprocess.check_call('psi4%s input.dat run.out' % self.nt(), cwd=dirname, shell=True)
             # Read energy and gradients from Psi4 output
             result = self.read_result(dirname)
         except (OSError, IOError, RuntimeError, subprocess.CalledProcessError):
@@ -1439,17 +1448,22 @@ def calc_wq_new(self, coords, dirname):
         # self.M.edit_qcrems({'jobtype':'force'})
         # self.M[0].write(os.path.join(dirname, 'run.in'))
         in_files = [('%s/input.dat' % dirname, 'input.dat')]
-        out_files = [('%s/output.dat' % dirname, 'output.dat'), ('%s/run.log' % dirname, 'run.log')]
+        out_files = [('%s/run.out' % dirname, 'run.out'), ('%s/run.log' % dirname, 'run.log')]
         # We will assume that the number of threads on the worker is 1, as this maximizes efficiency
         # in the limit of large numbers of jobs, although it may be controlled via environment variables.
-        queue_up_src_dest(wq, 'psi4 input.dat > run.log 2>&1', in_files, out_files, verbose=False)
-    
+        queue_up_src_dest(wq, 'psi4 input.dat run.out 2>&1', in_files, out_files, verbose=False)
+
+    def number_output(self, dirname, calcNum):
+        if not os.path.exists(os.path.join(dirname, 'run.out')):
+            raise RuntimeError('run.out does not exist')
+        shutil.copy2(os.path.join(dirname,'run.out'), os.path.join(dirname,'run_%03i.out' % calcNum))
+
     def read_result(self, dirname, check_coord=None):
         """ Read Psi4 calculation output. """
         if check_coord is not None:
             raise CheckCoordError("Coordinate checking not implemented")
         energy, gradient = None, None
-        psi4out = os.path.join(dirname, 'output.dat')
+        psi4out = os.path.join(dirname, 'run.out')
         with open(psi4out) as outfile:
             found_grad = False
             found_num_grad = False
@@ -1518,6 +1532,11 @@ def __init__(self, molecule, dirname=None, qcdir=None, threads=None):
         self.threads = threads
         self.prep_temp_folder(dirname, qcdir)
 
+        self.method = self.M.qcrems[0].get('method')
+        if self.method is None:
+            self.method = self.M.qcrems[0].get('exchange')
+        self.basis = self.M.qcrems[0].get('basis')
+
     def prep_temp_folder(self, dirname, qcdir):
         # Provide an initial qchem scratch folder (e.g. supplied initial guess)
         if not qcdir:
@@ -1855,13 +1874,21 @@ def calc_new(self, coords, dirname):
         new_schema["molecule"]["geometry"] = coords.tolist()
         new_schema.pop("program", None)
         ret = qcengine.compute(new_schema, self.program, return_dict=True)
-        # store the schema_traj for run_json to pick up
-        self.schema_traj.append(ret)
+
         if ret["success"] is False:
             raise QCEngineAPIEngineError("QCEngineAPI computation did not execute correctly. Message: " + ret["error"]["error_message"])
         # Unpack the energy and gradient
+        if "return_energy" not in ret["properties"]:
+            # SinglepointRecord dictionary from terachem.py in QCEngine won't have 'return_energy' key.
+            # Setting 'scf_total_energy' as 'return_energy' here, otherwise energies in OptimizationResult will be None.
+            ret["properties"]["return_energy"] = ret["properties"]["scf_total_energy"]
+
         energy = ret["properties"]["return_energy"]
         gradient = np.array(ret["return_result"])
+
+        # store the schema_traj for run_json to pick up
+        self.schema_traj.append(ret)
+
         return {'energy':energy, 'gradient':gradient}
 
     def calc(self, coords, dirname, **kwargs):
diff --git a/geometric/neb.py b/geometric/neb.py
index 16c2eebf..f547e118 100644
--- a/geometric/neb.py
+++ b/geometric/neb.py
@@ -262,6 +262,7 @@ def __init__(self, molecule, engine, tmpdir, params, coords=None):
         self.tmpdir = tmpdir
         # HP 5/2/2023: coordtype is set to 'cart' here.
         self.coordtype='cart'
+
         # coords is a 2D array with dimensions (N_image, N_atomx3) in atomic units
         if coords is not None:
             # Reshape the array if we passed in a flat array
@@ -292,19 +293,49 @@ def ComputeEnergyGradient(self, cyc=None, result=None):
         """Compute energies and gradients for each structure."""
         # This is the parallel point.
         wq = getWorkQueue()
-        if wq is None:
-            for i in range(len(self)):
-                if result:
-                    self.Structures[i].ComputeEnergyGradient(result=result[i])
-                else:
-                    self.Structures[i].ComputeEnergyGradient()
-        else:  # If work queue is available, handle jobs with the work queue.
+        if wq:  # If work queue is available, handle jobs with the work queue.
             for i in range(len(self)):
                 self.Structures[i].QueueEnergyGradient()
             wq_wait(wq, print_time=600)
             for i in range(len(self)):
                 self.Structures[i].GetEnergyGradient()
-        if cyc is not None:
+        elif self.params.bigchem:
+            # If BigChem is available, it will be used to carry the single-point calculations.
+            from qcio import Molecule as qcio_Molecule, ProgramInput
+            from bigchem import compute, group
+            elems = self.Structures[0].M.elem
+            molecules = []
+            
+            # Creating a list with qcio Molecule objects and submitting calculations.
+            for Structure in self.Structures:
+                molecules.append(qcio_Molecule(symbols=elems, geometry=Structure.cartesians.reshape(-1,3)))
+
+            outputs = group(compute.s(self.engine.__class__.__name__.lower(),
+                            ProgramInput(molecule=qcio_M, calctype="gradient",
+                                         model={"method":self.engine.method, "basis": self.engine.basis},
+                                         extras={"order":i}),
+                            ) for i, qcio_M in enumerate(molecules)).apply_async()
+
+            # Getting the records
+            records = outputs.get()
+            
+            # Passing the results to chain.ComputeEnergyGradient()
+            for i in range(len(self)):
+                assert records[i].input_data.extras["order"] == i
+                result = {"energy": records[i].results.energy, "gradient": np.array(records[i].results.gradient).ravel()}
+                self.Structures[i].ComputeEnergyGradient(result=result)
+
+            # Deleting the records
+            outputs.forget()
+
+        else:
+            for i in range(len(self)):
+                if result:
+                    self.Structures[i].ComputeEnergyGradient(result=result[i])
+                else:
+                    self.Structures[i].ComputeEnergyGradient()
+        # When BigChem is used, it does not write output files to disk.
+        if cyc is not None and not self.params.bigchem:
             for i in range(len(self)):
                 self.Structures[i].engine.number_output(self.Structures[i].tmpdir, cyc)
         self.haveCalcs = True
@@ -1682,8 +1713,11 @@ def main():
 
     M, engine = get_molecule_engine(**args)
 
-    if params.port != 0:
-        createWorkQueue(params.port, debug=params.verbose)
+    if params.bigchem:
+        logger.info("BigChem will be used to carry singlepoint calculations. \n")
+
+    if args.get('wqport', 0):
+        createWorkQueue(args.get('wqport'), debug=params.verbose)
 
     if params.prefix is None:
         tmpdir = os.path.splitext(args["input"])[0] + ".tmp"
diff --git a/geometric/normal_modes.py b/geometric/normal_modes.py
index 7fed9e26..8fec7be0 100644
--- a/geometric/normal_modes.py
+++ b/geometric/normal_modes.py
@@ -42,7 +42,7 @@
 from .molecule import Molecule, PeriodicTable
 from .nifty import logger, kb, kb_si, hbar, au2kj, au2kcal, ang2bohr, bohr2ang, c_lightspeed, avogadro, cm2au, amu2au, ambervel2au, wq_wait, getWorkQueue, commadash, bak
 
-def calc_cartesian_hessian(coords, molecule, engine, dirname, read_data=True, verbose=0):
+def calc_cartesian_hessian(coords, molecule, engine, dirname, read_data=True, bigchem=False, verbose=0):
     """ 
     Calculate the Cartesian Hessian using finite difference, and/or read data from disk. 
     Data is stored in a folder <prefix>.tmp/hessian, with gradient calculations found in
@@ -51,7 +51,7 @@ def calc_cartesian_hessian(coords, molecule, engine, dirname, read_data=True, ve
     Parameters
     ----------
     coords : np.ndarray
-        Nx3 array of Cartesian coordinates in atomic units
+        1-dimensional array of shape (3*N_atoms) containing atomic coordinates in Bohr
     molecule : Molecule
         Molecule object
     engine : Engine
@@ -134,6 +134,58 @@ def calc_cartesian_hessian(coords, molecule, engine, dirname, read_data=True, ve
             gbak = engine.read_wq(coords, dirname_d)['gradient']
             coords[i] += h
             Hx[i] = (gfwd-gbak)/(2*h)
+
+    elif bigchem:
+        # If BigChem is ready, it will be used to parallelize the Hessian calculation.
+        logger.info("BigChem will be used to calculate the Hessian. \n")
+        from qcio import Molecule as qcio_Molecule, ProgramInput
+        from bigchem import compute, group
+
+        elems = molecule[0].elem
+
+        # Uncommenting the following 6 lines and commenting out the rest of the lines will use BigChem's parallel_hessian function.
+        #from bigchem.algos import parallel_hessian
+        #qcio_M = qcio_Molecule(symbols=elems, geometry=coords.reshape(-1,3))
+        #input = ProgramInput(molecule=qcio_M, model={"method":engine.method, "basis": engine.basis}, calctype='hessian')
+        #output = parallel_hessian("psi4", input).delay()
+        #rec = output.get()
+        #Hx = rec.results.hessian
+
+        # Creating a list containing qcio Molecule obejcts with different geometries.
+        molecules = []
+        for i in range(nc):
+            coords[i] += h
+            molecules.append(qcio_Molecule(symbols=elems, geometry=coords.reshape(-1,3)))
+            coords[i] -= 2*h
+            molecules.append(qcio_Molecule(symbols=elems, geometry=coords.reshape(-1,3)))
+            coords[i] += h
+
+        # Submitting calculations
+        outputs = group(compute.s(engine.__class__.__name__.lower(),
+                        ProgramInput(molecule=qcio_M, calctype='gradient',
+                                     model={"method":engine.method, "basis": engine.basis},
+                                     extras={"order":i}),
+                        ) for i, qcio_M in enumerate(molecules)).apply_async()
+
+        # Getting the records
+        records = outputs.get(disable_sync_subtasks=False)
+        assert len(records) == nc*2
+
+        # Grouping the recrods
+        grouped_records = list(zip(records[::2], records[1::2]))
+
+        # Iterating through the grouped records to calculate the Hessian
+        for i, (fwd_rec, bak_rec) in enumerate(grouped_records):
+            # Double checking the order
+            assert fwd_rec.input_data.extras["order"] == i*2
+            assert bak_rec.input_data.extras["order"] == fwd_rec.input_data.extras["order"] + 1
+            gfwd = fwd_rec.results.gradient.ravel()
+            gbak = bak_rec.results.gradient.ravel()
+            Hx[i] = (gfwd-gbak)/(2*h)
+
+        # Deleting the records in the backend
+        outputs.forget()
+
     else:
         # First calculate a gradient at the central point, for linking scratch files.
         engine.calc(coords, dirname, read_data=read_data)
diff --git a/geometric/optimize.py b/geometric/optimize.py
index 323384e7..4e0a3993 100644
--- a/geometric/optimize.py
+++ b/geometric/optimize.py
@@ -5,7 +5,7 @@
 
 Authors: Lee-Ping Wang, Chenchen Song
 
-Contributors: Yudong Qiu, Daniel G. A. Smith, Alberto Gobbi, Josh Horton
+Contributors: Heejune Park, Yudong Qiu, Daniel G. A. Smith, Alberto Gobbi, Josh Horton
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
@@ -67,7 +67,7 @@ def __init__(self, coords, molecule, IC, engine, dirname, params, print_info=Tru
         Parameters
         ----------
         coords : np.ndarray
-            Nx3 array of Cartesian coordinates in atomic units
+            1-dimensional array of shape (3*N_atoms) containing atomic coordinates in Bohr
         molecule : Molecule
             Molecule object (Units Angstrom)
         IC : InternalCoordinates
@@ -333,7 +333,7 @@ def calcEnergyForce(self):
         if self.params.hessian == 'each' or self.recalcHess:
             # Hx is assumed to be the Cartesian Hessian at the current step.
             # Otherwise we use the variable name Hx0 to avoid almost certain confusion.
-            self.Hx = calc_cartesian_hessian(self.X, self.molecule, self.engine, self.dirname, read_data=True, verbose=self.params.verbose)
+            self.Hx = calc_cartesian_hessian(self.X, self.molecule, self.engine, self.dirname, read_data=True, bigchem=self.params.bigchem, verbose=self.params.verbose)
             if self.params.frequency:
                 self.frequency_analysis(self.Hx, 'iter%03i' % self.Iteration, False)
             if self.recalcHess:
@@ -345,7 +345,7 @@ def calcEnergyForce(self):
                 self.recalcHess = False
         elif self.Iteration == 0:
             if self.params.hessian in ['first', 'stop', 'first+last'] and not hasattr(self.params, 'hess_data'):
-                self.Hx0 = calc_cartesian_hessian(self.X, self.molecule, self.engine, self.dirname, read_data=True, verbose=self.params.verbose)
+                self.Hx0 = calc_cartesian_hessian(self.X, self.molecule, self.engine, self.dirname, read_data=True, bigchem=self.params.bigchem, verbose=self.params.verbose)
                 logger.info(">> Initial Cartesian Hessian Eigenvalues\n")
                 self.SortedEigenvalues(self.Hx0)
                 if self.params.frequency:
@@ -467,7 +467,7 @@ def IRC_step(self):
             logger.info('First, following the imaginary mode vector\n')
             if self.TSWavenum[1] < 0:
                 raise IRCError("There are more than one imaginary vibrational mode. Please optimize the structure and try again.\n")
-            elif self.TSWavenum[0] > 0:
+            elif self.TSWavenum[0] >= 0:
                 raise IRCError("No imaginary mode detected. Please optimize the structure and try again.\n")
 
             self.IRC_adjfactor = np.linalg.norm(self.TSNormal_modes_x[0] * np.sqrt(self.IC.mass))
@@ -1050,7 +1050,7 @@ def optimizeGeometry(self):
             Hx = self.IC.calcHessCart(self.X, self.G, self.H)
             np.savetxt(self.params.write_cart_hess, Hx, fmt='% 14.10f')
         if self.params.hessian in ['last', 'first+last', 'each']:
-            Hx = calc_cartesian_hessian(self.X, self.molecule, self.engine, self.dirname, read_data=False, verbose=self.params.verbose)
+            Hx = calc_cartesian_hessian(self.X, self.molecule, self.engine, self.dirname, read_data=False, bigchem=self.params.bigchem, verbose=self.params.verbose)
             if self.params.frequency:
                 self.frequency_analysis(Hx, 'last', True)
         return self.progress
@@ -1095,7 +1095,7 @@ def Optimize(coords, molecule, IC, engine, dirname, params, print_info=True):
     Parameters
     ----------
     coords : np.ndarray
-        Nx3 array of Cartesian coordinates in atomic units
+        1-dimensional array of shape (3*N_atoms) containing atomic coordinates in Bohr
     molecule : Molecule
         Molecule object
     IC : InternalCoordinates
@@ -1187,9 +1187,9 @@ def run_optimizer(**kwargs):
     M, engine = get_molecule_engine(**kwargs)
     
     # Create Work Queue object
-    if kwargs.get('port', 0):
+    if kwargs.get('wqport', 0):
         logger.info("Creating Work Queue object for distributed Hessian calculation\n")
-        createWorkQueue(kwargs['port'], debug=verbose>1)
+        createWorkQueue(kwargs['wqport'], debug=verbose>1)
 
     # Get initial coordinates in bohr
     coords = M.xyzs[0].flatten() * ang2bohr
@@ -1337,7 +1337,7 @@ def run_optimizer(**kwargs):
             if params.qdata is not None: Mfinal.write('qdata-final.txt')
     print_citation(logger)
     logger.info("Time elapsed since start of run_optimizer: %.3f seconds\n" % (time.time()-t0))
-    if kwargs.get('port', 0):
+    if kwargs.get('wqport', 0):
         destroyWorkQueue()
     return progress
 
diff --git a/geometric/params.py b/geometric/params.py
index 11397a62..465343c9 100644
--- a/geometric/params.py
+++ b/geometric/params.py
@@ -5,7 +5,7 @@
 
 Authors: Lee-Ping Wang, Chenchen Song
 
-Contributors: Yudong Qiu, Daniel G. A. Smith, Alberto Gobbi, Josh Horton
+Contributors: Heejune Park, Yudong Qiu, Daniel G. A. Smith, Alberto Gobbi, Josh Horton
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
@@ -107,6 +107,9 @@ def __init__(self, **kwargs):
         self.bothre = kwargs.get('bothre', 0.6 if self.transition else 0.0)
         # Whether to calculate or read a Hessian matrix.
         self.hessian = kwargs.get('hessian', None)
+        # Whether to use BigChem to carry the Hessian and NEB calculations.
+        self.bigchem = kwargs.get('bigchem', False)
+
         if self.hessian is None:
             # Default is to calculate Hessian in the first step if searching for a transition state.
             # Otherwise the default is to never calculate the Hessian.
@@ -226,13 +229,11 @@ def printInfo(self):
         elif self.hessian.startswith('file+last:'):
             logger.info(' Hessian data will be read from file: %s, then computed for the last step.\n' % self.hessian[5:])
 
-
 class NEBParams(object):
     """
     Container for optimization parameters.
     """
     def __init__(self, **kwargs):
-        self.port = kwargs.get('port', 0)
         self.prefix = kwargs.get('prefix', None)
         self.images = kwargs.get('images', 11)
         self.plain = kwargs.get('plain', 0)
@@ -253,6 +254,8 @@ def __init__(self, **kwargs):
         self.tmax = kwargs.get('tmax', 0.3)
         self.tmin = kwargs.get('tmin', 1.2e-3)
         self.skip = kwargs.get('skip', False)
+        self.bigchem = kwargs.get('bigchem', False)
+
         # Sanity checks on trust radius
         if self.tmax < self.tmin:
             raise ParamError("Max trust radius must be larger than min")
@@ -355,8 +358,10 @@ def parse_optimizer_args(*args):
                              '"stop" : Calculate Hessian for initial structure, then exit.\n'
                              'file:<path> : Read initial Hessian data in NumPy format from path, e.g. file:folder/hessian.txt\n'
                              'file+last:<path> : Read initial Hessian data as above, then compute for the last structure.\n ')
-    grp_hessian.add_argument('--port', type=int, help='Work Queue port used to distribute Hessian calculations. Workers must be started separately.')
-    grp_hessian.add_argument('--frequency', type=str2bool, help='Perform frequency analysis whenever Hessian is calculated, default is yes/true.\n ')
+    grp_hessian.add_argument('--wqport', type=int, help='Work Queue port used to distribute Hessian calculations. Workers must be started separately. \n ')
+    grp_hessian.add_argument('--bigchem', type=str2bool, help='Provide "Yes" to use BigChem for performing the Hessian calculation in parallel. \n'
+                                                              'Please ensure that BigChem is running with workers properly. \n ')
+    grp_hessian.add_argument('--frequency', type=str2bool, help='Perform frequency analysis whenever Hessian is calculated, default is yes/true. \n')
     grp_hessian.add_argument('--thermo', type=float, nargs=2, help='Temperature (K) and pressure (bar) for harmonic free energy\n'
                              'following frequency analysis, default is 300 K and 1.0 bar.\n ')
     grp_hessian.add_argument('--wigner', type=int, help='Number of desired samples from Wigner distribution after frequency analysis.\n'
@@ -495,8 +500,9 @@ def parse_neb_args(*args):
     grp_nebparam.add_argument('--tmin', type=float, help='Minimum trust radius, do not reject steps trust radius is below this threshold (method-dependent).\n ')
     grp_nebparam.add_argument('--skip', action='store_true', help='Skip Hessian updates that would introduce negative eigenvalues.\n ')
     grp_nebparam.add_argument('--epsilon', type=float, help='Small eigenvalue threshold for resetting Hessian, default 1e-5.\n ')
-    grp_nebparam.add_argument('--port', type=int, help='Work Queue port used to distribute singlepoint calculations. Workers must be started separately.\n')
-
+    grp_nebparam.add_argument('--wqport', type=int, help='Work Queue port used to distribute singlepoint calculations. Workers must be started separately.\n ')
+    grp_nebparam.add_argument('--bigchem', type=str2bool, help='Provide "Yes" to use BigChem for performing the NEB calculation in parallel. \n'
+                                                               'Please ensure that BigChem is running with workers properly. \n')
     grp_output = parser.add_argument_group('output', 'Control the format and amount of the output')
     grp_output.add_argument('--prefix', type=str, help='Specify a prefix for log file and temporary directory.\n'
                             'Defaults to the input file path (incl. file name with extension removed).\n ')
diff --git a/geometric/tests/addons.py b/geometric/tests/addons.py
index c230f2bb..a4cac00c 100644
--- a/geometric/tests/addons.py
+++ b/geometric/tests/addons.py
@@ -46,6 +46,20 @@ def get_gaussian_version():
 def workqueue_found():
     return (_plugin_import("work_queue") is True) and (geometric.nifty.which('work_queue_worker'))
 
+def bigchem_available(max_retries: int = 2):
+    """ 
+    A function written by Colton to check if BigChem is working properly.
+    """
+    try:
+        from bigchem.app import bigchem as bigchem_app
+        # Use the connection() method to get a broker connection object
+        with bigchem_app.connection() as conn:
+            # Open a connection to the broker
+            conn.ensure_connection(max_retries=max_retries)
+            return True
+    except Exception as e:
+        return False
+
 # Modify paths for testing
 os.environ["DQM_CONFIG_PATH"] = os.path.dirname(os.path.abspath(__file__))
 os.environ["TMPDIR"] = "/tmp/"
@@ -63,6 +77,8 @@ def workqueue_found():
     _plugin_import("openmm") is False and _plugin_import("simtk.openmm") is False, reason="could not find openmm. please install the package to enable tests")
 using_workqueue = pytest.mark.skipif(
     (not workqueue_found()), reason="could not find work_queue module or work_queue_worker executable. please install the package to enable tests")
+using_bigchem = pytest.mark.skipif(
+    (not bigchem_available()), reason="BigChem is not working. please install the package to enable tests")
 using_terachem = pytest.mark.skipif(
     not geometric.nifty.which("terachem"), reason="could not find terachem. please make sure TeraChem is installed for these tests")
 using_qchem = pytest.mark.skipif(
diff --git a/geometric/tests/conftest.py b/geometric/tests/conftest.py
index 2f8d275d..f85b126a 100644
--- a/geometric/tests/conftest.py
+++ b/geometric/tests/conftest.py
@@ -1,5 +1,9 @@
 import logging
 import _pytest.logging
+import pytest
+import geometric
+import os, tempfile
+from . import addons 
 
 # 2022-09-22: This modification allows live and captured logs to be 
 # formatted in the expected way when running tests using pytest. 
@@ -9,3 +13,33 @@
 # It took a long time to figure out and I'm glad it only took one line to fix it!
 
 logging.StreamHandler.terminator = ""
+
+datad = addons.datad
+
+@pytest.fixture
+def molecule_engine_hcn_neb():
+    """Return the Molecule and Engine for an NEB Calculation."""
+    input_ext = {'psi4': 'psi4in', 'qchem': 'qcin', 'tera': 'tcin'}
+    def get_molecule_engine(engine: str, images: int):
+        
+        return geometric.prepare.get_molecule_engine(
+            input=os.path.join(datad, "hcn_neb_input.%s" %input_ext.get(engine)),
+            chain_coords=os.path.join(datad, "hcn_neb_input.xyz"),
+            images=images,
+            neb=True,
+            engine=engine,
+    )
+
+    return get_molecule_engine
+
+@pytest.fixture
+def bigchem_frequency():
+    """Return frequency calculation results carried by BigChem"""
+    def get_freqs(engine: str, input: str):
+        molecule, engine = geometric.prepare.get_molecule_engine(engine=engine, input=os.path.join(datad, input))
+        coords = molecule.xyzs[0].flatten()*geometric.nifty.ang2bohr
+        hessian = geometric.normal_modes.calc_cartesian_hessian(coords, molecule, engine, tempfile.mkdtemp(), bigchem=True)
+        freqs, modes, G = geometric.normal_modes.frequency_analysis(coords, hessian, elem=molecule.elem, verbose=True)
+        return freqs
+
+    return get_freqs
diff --git a/geometric/tests/test_hessian.py b/geometric/tests/test_hessian.py
index 7d459d2f..f2a444e5 100644
--- a/geometric/tests/test_hessian.py
+++ b/geometric/tests/test_hessian.py
@@ -7,12 +7,14 @@
 import subprocess
 import os, shutil
 import numpy as np
+import tempfile
 from . import addons
 from geometric.internal import *
 
 datad = addons.datad
 localizer = addons.in_folder
 
+
 def test_hessian_assort():
     M = geometric.molecule.Molecule(os.path.join(datad, 'assort.xyz'))
     coords = M.xyzs[0].flatten() * ang2bohr
@@ -42,6 +44,37 @@ def test_hessian_assort():
     assert np.allclose(fgrad_dlc, agrad_dlc, atol=1.e-6)
     assert np.allclose(fhess_dlc, ahess_dlc, atol=1.e-6)
 
+@addons.using_psi4
+@addons.using_bigchem
+def test_psi4_bigchem_hessian(bigchem_frequency):
+    freqs = bigchem_frequency("psi4", "hcn_minimized.psi4in")
+
+    np.testing.assert_almost_equal(freqs, [989.5974, 989.5992, 2394.0352, 3690.5745], decimal=0)
+    assert len(freqs) == 4
+
+
+@addons.using_qchem
+@addons.using_bigchem
+def test_qchem_bigchem_hessian(bigchem_frequency):
+    # Optimized TS structure of HCN -> CNH reaction. 
+    freqs = bigchem_frequency("qchem", "hcn_minimized_ts.qcin")
+
+    np.testing.assert_almost_equal(freqs, [-1215.8587, 2126.6551, 2451.8599], decimal=0)
+    assert freqs[0] < 0
+    assert len(freqs) == 3
+
+
+@addons.using_terachem
+@addons.using_bigchem
+def test_terachem_bigchem_hessian(bigchem_frequency):
+    shutil.copy2(os.path.join(datad,"hcn_minimized.xyz"), os.path.join(os.getcwd(), "hcn_minimized.xyz"))
+    freqs = bigchem_frequency("tera", "hcn_minimized.tcin")
+    os.remove("hcn_minimized.xyz")
+
+    np.testing.assert_almost_equal(freqs, [989.5482, 989.7072, 2394.4201, 3687.4741], decimal=0)
+    assert len(freqs) == 4
+
+
 class TestPsi4WorkQueueHessian:
 
     """ Tests are put into class so that the fixture can terminate the worker process. """
@@ -74,6 +107,7 @@ def test_psi4_work_queue_hessian(self, localizer):
         assert len(freqs) == 4
         geometric.nifty.destroyWorkQueue()
 
+
 class TestASEWorkQueueHessian:
 
     """ Tests are put into class so that the fixture can terminate the worker process. """
diff --git a/geometric/tests/test_irc.py b/geometric/tests/test_irc.py
index d574bf1c..f48037ad 100644
--- a/geometric/tests/test_irc.py
+++ b/geometric/tests/test_irc.py
@@ -9,12 +9,12 @@
 
 localizer = addons.in_folder
 datad = addons.datad
-exampled = addons.exampled
+
 
 @addons.using_psi4
 def test_hcn_irc_psi4(localizer):
     """
-    IRC test
+    IRC test without BigChem
     """
     shutil.copy2(os.path.join(datad, 'hcn_irc.psi4in'), os.path.join(os.getcwd(), 'hcn_irc.psi4in'))
     progress = geometric.optimize.run_optimizer(engine='psi4', input='hcn_irc.psi4in', converge=['set', 'GAU_LOOSE'],
@@ -35,3 +35,31 @@ def test_hcn_irc_psi4(localizer):
 
     # Check that the IRC converged in less than 100 iterations
     assert len(progress) < 100
+
+
+@addons.using_psi4
+@addons.using_bigchem
+def test_hcn_irc_psi4_bigchem(localizer):
+    """
+    IRC test with BigChem
+    """
+    shutil.copy2(os.path.join(datad, 'hcn_irc.psi4in'), os.path.join(os.getcwd(), 'hcn_irc.psi4in'))
+    progress = geometric.optimize.run_optimizer(engine='psi4', input='hcn_irc.psi4in', converge=['set', 'GAU_LOOSE'],
+                                                nt=4, reset=False, trust=0.05, irc=True, maxiter=50, bigchem=True)
+    e_ref1 = -92.35408411
+    e_ref2 = -92.33971205
+    max_e = np.max(progress.qm_energies)
+    reac_e = progress.qm_energies[0]
+    prod_e = progress.qm_energies[-1]
+
+    # Check the mex_e is not from the end-points
+    assert reac_e < max_e
+    assert prod_e < max_e
+
+    # Check the end-point energies
+    assert np.isclose(min(reac_e, prod_e), e_ref1)
+    assert np.isclose(max(reac_e, prod_e), e_ref2)
+
+    # Check that the IRC converged in less than 100 iterations
+    assert len(progress) < 100
+
diff --git a/geometric/tests/test_neb.py b/geometric/tests/test_neb.py
index c4265245..137e3154 100644
--- a/geometric/tests/test_neb.py
+++ b/geometric/tests/test_neb.py
@@ -14,47 +14,28 @@
 datad = addons.datad
 exampled = addons.exampled
 
-
-def test_hcn_neb_input(localizer):
+def test_hcn_neb_input(localizer, molecule_engine_hcn_neb):
     """
     Test lengths of input chains
     """
     chain_M = geometric.molecule.Molecule(os.path.join(datad, "hcn_neb_input.xyz"))
-
     nimg = 7
-    M1, engine = geometric.prepare.get_molecule_engine(
-        input=os.path.join(datad, "hcn_neb_input.psi4in"),
-        chain_coords=os.path.join(datad, "hcn_neb_input.xyz"),
-        images=nimg,
-        neb=True,
-        engine="psi4",
-    )
+
+    M1, engine = molecule_engine_hcn_neb('psi4', nimg)
 
     # The number of images can't exceed the maximum number of images in the input chain
-    M2, engine = geometric.prepare.get_molecule_engine(
-        input=os.path.join(datad, "hcn_neb_input.psi4in"),
-        chain_coords=os.path.join(datad, "hcn_neb_input.xyz"),
-        images=9999,
-        neb=True,
-        engine="psi4",
-    )
+    M2, engine = molecule_engine_hcn_neb('psi4', 9999)
 
     assert nimg == len(M1)
     assert len(M2) == len(chain_M)
 
 
 @addons.using_psi4
-def test_hcn_neb_optimize_1(localizer):
+def test_hcn_neb_optimize_1(localizer, molecule_engine_hcn_neb):
     """
     Optimize a HCN chain without alignment
     """
-    M, engine = geometric.prepare.get_molecule_engine(
-        input=os.path.join(datad, "hcn_neb_input.psi4in"),
-        chain_coords=os.path.join(datad, "hcn_neb_input.xyz"),
-        images=11,
-        neb=True,
-        engine="psi4",
-    )
+    M, engine = molecule_engine_hcn_neb('psi4', 11)
 
     params = geometric.params.NEBParams(**{"optep": True, "align": False, "verbose": 1})
     chain = geometric.neb.ElasticBand(
@@ -65,22 +46,17 @@ def test_hcn_neb_optimize_1(localizer):
 
     final_chain, optCycle = geometric.neb.OptimizeChain(chain, engine, params)
 
-    assert optCycle < 10
+    assert optCycle <= 10
     assert final_chain.maxg < params.maxg
     assert final_chain.avgg < params.avgg
 
+
 @addons.using_psi4
-def test_hcn_neb_optimize_2(localizer):
+def test_hcn_neb_optimize_2(localizer, molecule_engine_hcn_neb):
     """
     Optimize a HCN chain with alignment
     """
-    M, engine = geometric.prepare.get_molecule_engine(
-        input=os.path.join(datad, "hcn_neb_input.psi4in"),
-        chain_coords=os.path.join(datad, "hcn_neb_input.xyz"),
-        images=7,
-        neb=True,
-        engine="psi4",
-    )
+    M, engine = molecule_engine_hcn_neb('psi4', 11)
 
     # maxg and avgg are increased here to make them converge faster after the alignment
     params = geometric.params.NEBParams(**{"verbose": 1, "maxg": 3.0, "avgg": 2.0})
@@ -92,10 +68,83 @@ def test_hcn_neb_optimize_2(localizer):
 
     final_chain, optCycle = geometric.neb.OptimizeChain(chain, engine, params)
 
-    assert optCycle < 10
+    assert optCycle <= 10
     assert final_chain.maxg < params.maxg
     assert final_chain.avgg < params.avgg
 
+
+@addons.using_psi4
+@addons.using_bigchem
+def test_psi4_bigchem(localizer, molecule_engine_hcn_neb):
+    """
+    Optimize a HCN chain without alignment with BigChem and Psi4
+    """
+    M, engine = molecule_engine_hcn_neb('psi4', 11)
+
+    params = geometric.params.NEBParams(**{"align": False, "verbose": 1, "bigchem": True})
+    chain = geometric.neb.ElasticBand(
+        M, engine=engine, tmpdir=tempfile.mkdtemp(), params=params, plain=0
+    )
+
+    assert chain.coordtype == "cart"
+
+    final_chain, optCycle = geometric.neb.OptimizeChain(chain, engine, params)
+
+    assert optCycle <= 10
+    assert final_chain.maxg < params.maxg
+    assert final_chain.avgg < params.avgg
+
+
+
+@addons.using_qchem
+@addons.using_bigchem
+def test_qchem_bigchem(localizer, molecule_engine_hcn_neb):
+    """
+    Optimize a HCN chain without alignment with BigChem and QChem
+    """
+    M, engine = molecule_engine_hcn_neb('qchem', 11)
+
+    params = geometric.params.NEBParams(**{"align": False, "verbose": 1, "bigchem": True})
+    chain = geometric.neb.ElasticBand(
+        M, engine=engine, tmpdir=tempfile.mkdtemp(), params=params, plain=0
+    )
+
+    assert chain.coordtype == "cart"
+
+    final_chain, optCycle = geometric.neb.OptimizeChain(chain, engine, params)
+
+    assert optCycle <= 10
+    assert final_chain.maxg < params.maxg
+    assert final_chain.avgg < params.avgg
+
+
+@addons.using_terachem
+@addons.using_bigchem
+def test_terachem_bigchem(localizer, molecule_engine_hcn_neb):
+    """
+    Optimize a HCN chain without alignment with BigChem and TeraChem
+    """
+    shutil.copy2(os.path.join(datad, 'hcn_neb_input.xyz'), os.path.join(os.getcwd(), 'hcn_neb_input.xyz'))
+
+    M, engine = molecule_engine_hcn_neb('tera', 11)
+
+    params = geometric.params.NEBParams(**{"align": False, "verbose": 1, "bigchem": True})
+    chain = geometric.neb.ElasticBand(
+        M, engine=engine, tmpdir=tempfile.mkdtemp(), params=params, plain=0
+    )
+
+    os.remove('hcn_neb_input.xyz')
+
+    assert chain.coordtype == "cart"
+
+    final_chain, optCycle = geometric.neb.OptimizeChain(chain, engine, params)
+
+    assert optCycle <= 10
+    assert final_chain.maxg < params.maxg
+    assert final_chain.avgg < params.avgg
+
+
+
 class TestPsi4WorkQueueNEB:
 
     """ Tests are put into class so that the fixture can terminate the worker process. """
@@ -110,14 +159,8 @@ def work_queue_cleanup(self):
 
     @addons.using_psi4
     @addons.using_workqueue
-    def test_psi4_work_queue_neb(self, localizer):
-        M, engine = geometric.prepare.get_molecule_engine(
-            input=os.path.join(datad, "hcn_neb_input.psi4in"),
-            chain_coords=os.path.join(datad, "hcn_neb_input.xyz"),
-            images=11,
-            neb=True,
-            engine="psi4",
-        )
+    def test_psi4_work_queue_neb(self, localizer, molecule_engine_hcn_neb):
+        M, engine = molecule_engine_hcn_neb('psi4', 11)
         params = geometric.params.NEBParams(**{"optep": True, "align": False, "verbose": 1})
         chain = geometric.neb.ElasticBand(
             M, engine=engine, tmpdir=tempfile.mkdtemp(), params=params, plain=0
@@ -137,6 +180,7 @@ def test_psi4_work_queue_neb(self, localizer):
         assert final_chain.avgg < params.avgg
         geometric.nifty.destroyWorkQueue()
 
+
 @addons.using_qcelemental
 def test_hcn_neb_service_arrange(localizer):
     """
diff --git a/setup.py b/setup.py
index 426bc310..b4f43526 100644
--- a/setup.py
+++ b/setup.py
@@ -31,6 +31,9 @@
         'pytest',
         'pytest-cov',
     ],
+    extra_require={
+        'bigchem': ['bigchem']        
+    },
     version=versioneer.get_version(),
     cmdclass=versioneer.get_cmdclass(),
     classifiers=[