Repository: ifilot/pyqint (language: python)

Issue:
Molecule class should safeguard for vector input
Encapsulate variables in `float` on https://github.com/ifilot/pyqint/blob/master/pyqint/molecule.py#L24


Discussion / hints:

Fix patch:
diff --git a/MANIFEST.in b/MANIFEST.in
index 5696f9c..f67a7de 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,1 +1,2 @@
-include pyqint/basissets/*.json
\ No newline at end of file
+include pyqint/basissets/*.json
+include pyqint/molecules/*.xyz
\ No newline at end of file
diff --git a/README.md b/README.md
index ca88fd9..8020d3a 100644
--- a/README.md
+++ b/README.md
@@ -6,31 +6,60 @@
 [![PyPI](https://img.shields.io/pypi/v/pyqint?style=flat-square)](https://pypi.org/project/pyqint/)
 [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
 
-__Table of Contents__
-
-* [Purpose](#purpose)
-* [Documentation](#documentation)
-* [Examples](#documentation)
-    - [Contour plots of the molecular orbitals of CO](#contour-plots-of-the-molecular-orbitals-of-co)
-
 ## Purpose
 
 PyQInt is a Python package for calculating one- and two-electron integrals as
 encountered in electronic structure calculations. Since integral evaluation can
-be quite computationally intensive, they are programmed in C++ and connected to
-Python using Cython. PyQInt offers additional features such as performed restricted **Hartree-Fock
-calculations**, calculation of **Crystal Orbital Hamilton Population** coefficients
-and constructing localized orbitals using the **Boys-Foster** method.
+be quite computationally intensive, the evaluation is programmed in C++ and
+connected to Python using Cython.
+
+PyQInt mainly serves as an educational package to teach students how to perform
+(simple) electronic structure calculations wherein the most difficult task,
+i.e. the integral evaluation, is already encapsulated in a handy set of
+routines. With PyQInt, the student can for example build their own Hartree-Fock
+routine. Some common electronic structure routine, most notably the
+Hartree-Fock algorithm, is also readily available.
+
+> **Note**
+> Although PyQInt connects to a C++ backend, it is certainly not optimized for
+> speed and might be (too) slow for anything outside of the calculation of the
+> electronic structure of simple molecules.
 
 ## Documentation
 
 PyQInt comes with detailed documentation and examples, which can be found
 at https://pyqint.imc-tue.nl.
 
-## Examples
+## Features
+
+The following molecular integrals are supported by PyQInt
+
+- [x] Overlap integral
+- [x] Kinetic integral
+- [x] Dipole integral
+- [x] Nuclear integral
+- [x] Two-electron repulsion integral
+
+as well as the following geometric derivatives
+
+- [x] Overlap integral
+- [x] Kinetic integral
+- [x] Nuclear integral
+- [x] Two-electron repulsion integral
+
+PyQInt offers additional features such as
+* Performing [restricted Hartree-Fock](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method)
+  calculations using [DIIS](https://en.wikipedia.org/wiki/DIIS)
+* Calculation of [Crystal Orbital Hamilton Population](http://www.cohp.de/) coefficients
+* Construction of localized orbitals using the [Boys-Foster method](https://en.wikipedia.org/wiki/Localized_molecular_orbitals#Foster-Boys)
+* Visualization of molecular orbitals
 
-### Contour plots of the molecular orbitals of CO
+All routines are (automatically) tested and verified against several open-source
+as well as commercial programs that use cartesian Gaussian orbitals. Nevertheless,
+if you spot any mistake, please kindly open an [issue](https://github.com/ifilot/pyqint/issues)
+in this Github repository.
 
-Below, the molecular orbitals for CO using a `sto-3g` basis set are shown.
+In the image below, the (canonical) molecular orbitals as found using a restricted
+Hartree-Fock calculation for the CO molecule are shown.
 
 ![Molecular orbitals of CO](img/co.jpg)
diff --git a/examples/.gitignore b/examples/.gitignore
index 6708eb6..ef28ff6 100644
--- a/examples/.gitignore
+++ b/examples/.gitignore
@@ -1,2 +1,3 @@
 *.npy
 *.jpg
+*.ply
diff --git a/examples/cohp_lih.py b/examples/cohp_lih.py
new file mode 100644
index 0000000..b5d1264
--- /dev/null
+++ b/examples/cohp_lih.py
@@ -0,0 +1,110 @@
+from pyqint import Molecule, HF
+import numpy as np
+import scipy.optimize
+import matplotlib.pyplot as plt
+from matplotlib.patches import Rectangle
+
+#
+# Perform a (simplified) geometry optimization on the CO molecule and
+# calculate the COHP coefficients for the molecule.
+#
+
+def main():
+    # find optimum for CO molecule
+    res = scipy.optimize.minimize(optimize_lih, [1.14], tol=1e-4)
+    print('Optimal distance found at d = %f' % res.x)
+    
+    # calculate sto-3g coefficients for h2o
+    result = calculate_lih(res.x)
+    
+    energies = result['orbe']
+    coeff = result['orbc']
+    H = result['fock']
+    
+    cohp = np.zeros(len(coeff))
+    for k in range(0, len(coeff)): # loop over molecular orbitals
+        for i in range(0,len(coeff)-1): # loop over orbitals on Li
+            for j in range(len(coeff)-1,len(coeff)): # loop over orbitals on H
+                cohp[k] += 2.0 * H[i,j] * coeff[i,k] * coeff[j,k]
+                
+    for i,(chi,e) in enumerate(zip(cohp, energies)):
+        print('%02i %+8.4f %+8.4f' % (i+1, e, chi))
+
+    # labels = [
+    #     'Li,1s',
+    #     'Li,2s',
+    #     'Li,2p$_{x}$',
+    #     'Li,2p$_{y}$',
+    #     'Li,2p$_{z}$',
+    #     'H,1s'
+    # ]
+
+    # # fig, ax = plt.subplots(1,2, dpi=300)
+    # # plot_matrix(ax[0], H, xlabels=labels, ylabels=labels,
+    # #             title='Hamiltonian matrix')
+    # # plot_matrix(ax[1], coeff, ylabels=labels, title='Coefficient matrix',
+    # #             orbe=energies)
+
+def optimize_lih(d):
+    """
+    Optimization function for scipy.optimize.minimize
+    """
+    mol = Molecule()
+    mol.add_atom('Li', 0.0, 0.0, -d[0]/2, unit='angstrom')
+    mol.add_atom('H', 0.0, 0.0,  d[0]/2, unit='angstrom')
+    
+    result = HF().rhf(mol, 'sto3g')
+    
+    return result['energy']
+
+def calculate_lih(d):
+    """
+    Full function for evaluation
+    """
+    mol = Molecule()
+    mol.add_atom('Li', 0.0, 0.0, -d[0]/2, unit='angstrom')
+    mol.add_atom('H', 0.0, 0.0,  d[0]/2, unit='angstrom')
+    
+    result = HF().rhf(mol, 'sto3g')
+    
+    return result
+
+def plot_matrix(ax, mat, xlabels=None, ylabels=None, title = None, orbe=None):
+    """
+    Produce plot of matrix
+    """
+    ax.imshow(mat, vmin=-1, vmax=1, cmap='PiYG')
+    for i in range(mat.shape[0]):
+        for j in range(mat.shape[1]):
+            ax.text(i, j, '%.3f' % mat[j,i], ha='center', va='center',
+                    fontsize=7, zorder=5)
+    ax.set_xticks([])
+    ax.set_yticks([])
+    ax.hlines(np.arange(1, mat.shape[0])-0.5, -0.5, mat.shape[0] - 0.5,
+              color='black', linestyle='--', linewidth=1)
+    ax.vlines(np.arange(1, mat.shape[0])-0.5, -0.5, mat.shape[0] - 0.5,
+              color='black', linestyle='--', linewidth=1)
+    
+    # add xlabels if available
+    if xlabels:
+        ax.set_xticks(np.arange(0, mat.shape[0]))
+        ax.set_xticklabels(ylabels)
+        ax.tick_params(axis='both', which='major', labelsize=7)
+    
+    # add basis functions as axes labels
+    if ylabels:
+        ax.set_yticks(np.arange(0, mat.shape[0]))
+        ax.set_yticklabels(ylabels)
+        ax.tick_params(axis='both', which='major', labelsize=7)
+    
+    if orbe is not None:
+        orbe = ['(%i) / %.4f' % ((i+1),e) for i,e in enumerate(orbe)]
+        ax.set_xticks(np.arange(0, mat.shape[0]))
+        ax.set_xticklabels(orbe, rotation=90)
+    
+    # add title if supplied
+    if title:
+        ax.set_title(title)
+
+if __name__ == '__main__':
+    main()
diff --git a/examples/foster_boys.py b/examples/foster_boys.py
new file mode 100644
index 0000000..b0d2849
--- /dev/null
+++ b/examples/foster_boys.py
@@ -0,0 +1,50 @@
+from pyqint import Molecule, HF, PyQInt, FosterBoys
+import pyqint
+import numpy as np
+from pytessel import PyTessel
+
+#
+# Plot the isosurfaces for the CO molecule
+#
+
+def main():   
+    res = calculate_co(1.145414)
+    resfb = FosterBoys(res).run()
+
+    for i in range(len(res['cgfs'])):
+        build_isosurface('MO_%03i' % (i+1), 
+                         res['cgfs'],
+                         resfb['orbc'][:,i],
+                         0.1)
+
+def calculate_co(d):
+    """
+    Full function for evaluation
+    """
+    mol = Molecule()
+    mol.add_atom('C', 0.0, 0.0, -d/2, unit='angstrom')
+    mol.add_atom('O', 0.0, 0.0,  d/2, unit='angstrom')
+
+    result = HF().rhf(mol, 'sto3g')
+
+    return result
+
+def build_isosurface(filename, cgfs, coeff, isovalue, sz=5, npts=100):
+    # generate some data
+    isovalue = np.abs(isovalue)
+    integrator = PyQInt()
+    grid = integrator.build_rectgrid3d(-sz, sz, npts)
+    scalarfield = np.reshape(integrator.plot_wavefunction(grid, coeff, cgfs), (npts, npts, npts))
+    unitcell = np.diag(np.ones(3) * 2 * sz)
+
+    pytessel = PyTessel()
+    vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), isovalue)
+    fname = filename + '_pos.ply'
+    pytessel.write_ply(fname, vertices, normals, indices)
+    
+    vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), -isovalue)
+    fname = filename + '_neg.ply'
+    pytessel.write_ply(fname, vertices, normals, indices)
+
+if __name__ == '__main__':
+    main()
diff --git a/examples/molecular_orbitals_co.py b/examples/molecular_orbitals_co.py
index 6b98b0d..f9bea40 100644
--- a/examples/molecular_orbitals_co.py
+++ b/examples/molecular_orbitals_co.py
@@ -14,7 +14,7 @@ def main():
         
     fig, ax = plt.subplots(2, 5, dpi=144, figsize=(12,5))           
     for i,(chi,e) in enumerate(zip(coeff.transpose(), energies)):
-        res, x, y = build_isosurface(result['cgfs'], chi, sz=3, plane='xz')
+        res, x, y = build_contourplot(result['cgfs'], chi, sz=3, plane='xz')
         vmax = np.max(np.abs(res))
         vmin = -vmax
         ax[i//5,i%5].contourf(x, y, res, levels=15, cmap='PiYG',
@@ -53,7 +53,7 @@ def calculate_co(d):
     
     return result
 
-def build_isosurface(cgfs, coeff, sz=2, npts=50, plane='xy'):
+def build_contourplot(cgfs, coeff, sz=2, npts=50, plane='xy'):
     integrator = PyQInt()
     
     # build grid
diff --git a/meta.yaml b/meta.yaml
index 22450f8..0b02a3d 100644
--- a/meta.yaml
+++ b/meta.yaml
@@ -1,6 +1,6 @@
 package:
   name: "pyqint"
-  version: "0.11.0.0"
+  version: "0.12.0.0"
 
 source:
   path: .
@@ -33,6 +33,9 @@ test:
   source_files:
     - tests/*.py
     - tests/results/*.npy
+    - tests/results/*.txt
+  commands:
+    - nosetests tests
 
 about:
   home: https://github.com/ifilot/pyqint
diff --git a/pyqint/__init__.py b/pyqint/__init__.py
index 21d34de..c201978 100644
--- a/pyqint/__init__.py
+++ b/pyqint/__init__.py
@@ -4,5 +4,8 @@
 from .gto import gto
 from .hf import HF
 from .foster_boys import FosterBoys
+from .cohp import COHP
+from .molecule_builder import MoleculeBuilder
+from .geometry_optimization import GeometryOptimization
 
 from ._version import __version__
diff --git a/pyqint/_version.py b/pyqint/_version.py
index 9a3ad98..5e5199d 100644
--- a/pyqint/_version.py
+++ b/pyqint/_version.py
@@ -1,2 +1,2 @@
-__version__ = "0.11.0.0"
+__version__ = "0.12.0.0"
 
diff --git a/pyqint/cgf.cpp b/pyqint/cgf.cpp
index ca89dee..fd41564 100644
--- a/pyqint/cgf.cpp
+++ b/pyqint/cgf.cpp
@@ -119,7 +119,7 @@ vec3 GTO::get_grad(const vec3& r) const {
  * @return void
  */
 void GTO::calculate_normalization_constant() {
-    static const double pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286;
+    static const double pi = boost::math::constants::pi<double>();
 
     double nom =   std::pow(2.0, 2.0 * (l + m + n) + 3.0 / 2.0) *
                    std::pow(alpha, (l + m + n) + 3.0 / 2.0);
diff --git a/pyqint/cohp.py b/pyqint/cohp.py
new file mode 100644
index 0000000..a6a1275
--- /dev/null
+++ b/pyqint/cohp.py
@@ -0,0 +1,50 @@
+# -*- coding: utf-8 -*-
+
+from .pyqint import PyQInt
+import numpy as np
+
+class COHP:
+    """
+    Class for performing Crystal Orbital Hamilton Population Analysis
+    """
+    def __init__(self, res):
+        # copy objects from Hartree-Fock result dictionary
+        self.orbc_canonical = res['orbc']
+        self.orbe_canonical = res['orbe']
+        self.nelec = res['nelec']
+        self.nuclei = res['nuclei']
+        self.H = res['fock']
+        self.cgfs = res['cgfs']
+        self.maxiter = 1000
+        self.occ = [1 if i < self.nelec//2 else 0 for i in range(0, len(self.cgfs))]
+
+    def run(self, C, n1, n2):
+        """
+        Run the COHP algorithm, take a coefficient matrix as input
+
+        C: coefficient matrix
+        n1: index of nucleus 1
+        n2: index of nucleus 2
+        """
+        if n1 == n2:
+            raise Exception('Cannot perform COHP for the same atoms')
+
+        # figure out which basis functions below to which nucleus
+        cgfs1 = []
+        cgfs2 = []
+        nuc1 = self.nuclei[n1][0]
+        nuc2 = self.nuclei[n2][0]
+        for i,cgf in enumerate(self.cgfs):
+            if np.linalg.norm(cgf.p - nuc1) < 1e-3:
+                cgfs1.append(i)
+            if np.linalg.norm(cgf.p - nuc2) < 1e-3:
+                cgfs2.append(i)
+
+        N = len(C)
+        cohp = np.zeros(N)
+        for k in range(0, N): # loop over molecular orbitals
+            for i in cgfs1: # loop over atom A
+                for j in cgfs2: # loop over atom B
+                    cohp[k] += 2.0 * self.H[i,j] * C[i,k] * C[j,k]
+
+        return cohp
diff --git a/pyqint/foster_boys.py b/pyqint/foster_boys.py
index 319e2d2..0e22db8 100644
--- a/pyqint/foster_boys.py
+++ b/pyqint/foster_boys.py
@@ -4,12 +4,13 @@
 import numpy as np
 import random
 import scipy.optimize
+from scipy.stats import ortho_group
 
 class FosterBoys:
     """
     Routine for constructing localized orbitals using Foster-Boys procedure
     """
-    def __init__(self, res):
+    def __init__(self, res, seed=None):
         # copy objects from Hartree-Fock result dictionary
         self.orbc_canonical = res['orbc']
         self.orbe_canonical = res['orbe']
@@ -18,16 +19,36 @@ def __init__(self, res):
         self.cgfs = res['cgfs']
         self.maxiter = 1000
         self.occ = [1 if i < self.nelec//2 else 0 for i in range(0, len(self.cgfs))]
+        self.rngseed = seed
+        self.rng = np.random.default_rng(seed=self.rngseed)
         
         # construct dipole tensor
         self.dipol = self.__build_dipole_tensor(self.cgfs)
     
-    def run(self):
+    def run(self, nr_runners=1):
+        """
+        Perform the Foster-Boys localization procedure
+
+        Because the Foster-Boys procedure uses random initialization, it is possible that
+        the algorithm ends up in a local minima. To circumvent this, the user can use
+        a number of runners and automatically use the best result found.
+        """
+        bestres = None
+        bestr2 = 0.0
+
+        for i in range(0, nr_runners):
+            res = self.__single_runner()
+            if res['r2final'] > bestr2:
+                bestres = res
+
+        return res
+
+    def __single_runner(self):
         """
         Perform the Foster-Boys optimization routine
         """
         # start with a random unitary transformation
-        C = self.__construct_random_orthogonal_transformation(self.orbc_canonical, 
+        C = self.__construct_random_orthogonal_transformation(self.orbc_canonical,
                                                               self.nelec)
 
         old_r2 = 0.0 # start assuming orbital centroids lie at origin
@@ -66,9 +87,12 @@ def __perform_foster_boys_mixing_scipy(self, C):
         for i in range(N):
             for j in range(i+1, N):
                 
+                # optimize local pair of MOs on the interval -(pi,pi)
                 res = scipy.optimize.minimize(self.__evaluate_orbital_centroids, 
-                                              np.pi/4, 
-                                              args=(C, i, j))
+                                              0.0,
+                                              args=(C, i, j),
+                                              bounds=[(-np.pi, np.pi)],
+                                              tol=1e-12)
                 
                 alpha = res.x[0]
                 Cnew = C.copy()
@@ -116,10 +140,12 @@ def __construct_random_orthogonal_transformation(self, C, nelec, nops=100):
         """
         for i in range(0,nops):
             Cnew = C.copy()
+
             # randomly pick two numbers among the occupied orbitals
-            n = random.sample(range(nelec//2),2)
+            n = self.rng.choice(range(nelec//2), size=2, replace=False)
             
-            gamma = random.random() * 2.0 * np.pi
+            # and mix them by a random angle
+            gamma = self.rng.uniform() * 2.0 * np.pi
             
             for j in range(0,len(C)):
                 Cnew[j,n[0]] = np.cos(gamma) * C[j,n[0]] + np.sin(gamma) * C[j,n[1]]
@@ -127,7 +153,22 @@ def __construct_random_orthogonal_transformation(self, C, nelec, nops=100):
 
             C = Cnew.copy()
 
-        return C
+        return Cnew
+
+    def __construct_random_orthogonal_transformation_from_orthogroup(self, C, nelec, nops=100):
+        """
+        Create orthogonal transformation usign the scipy ortho_group function
+
+        This routine should in theory act as a substitute for the
+        __construct_random_orthogonal_transformation() routine, but it is not working
+        appropriately. Leaving it here for potential future development.
+        """
+        N = nelec//2
+        orthomat = ortho_group.rvs(N)
+        identitymat = np.identity(len(C) - N)
+        T = scipy.linalg.block_diag(orthomat, identitymat)
+
+        return T @ C @ T.transpose()
     
     def __build_dipole_tensor(self, cgfs):
         """
diff --git a/pyqint/geometry_optimization.py b/pyqint/geometry_optimization.py
new file mode 100644
index 0000000..d227f8d
--- /dev/null
+++ b/pyqint/geometry_optimization.py
@@ -0,0 +1,171 @@
+# -*- coding: utf-8 -*-
+
+from . import HF, Molecule
+import numpy as np
+import scipy.optimize
+
+class GeometryOptimization:
+    """
+    Class to perform geometry optimizaton using the Conjugate Gradient algorithm
+    as implemented in the scipy library
+    """
+    def __init__(self, verbose=False):
+        self.cinit = None
+        self.P = None
+        self.orbe = None
+        self.verbose = verbose
+        self.iter = 0
+        self.energies_history = []
+        self.forces_history = []
+        self.coordinates_history = []
+
+    def run(self, mol, basis):
+        """
+        Perform geometry optimization
+        """
+        x0 = self.__unpack(mol) # unpack coordinates from molecule class
+        self.iter = 0
+
+        # reset arrays
+        self.energies_history = []
+        self.forces_history = []
+        self.coordinates_history = []
+
+        if self.verbose:
+            self.__print_break(newline=False)
+            print("START GEOMETRY OPTIMIZATION")
+            print("USING CONJUGATE GRADIENT PROCEDURE")
+            self.__print_break(newline=True)
+
+        res_opt = scipy.optimize.minimize(self.energy, x0, args=(mol, basis), method='CG',
+                                          jac=self.jacobian)
+
+        res = {
+            'opt': res_opt,
+            'energies': self.energies_history,
+            'forces': self.forces_history,
+            'coordinates': self.coordinates_history
+        }
+
+        return res
+
+    def energy(self, x, mol, basis):
+        self.iter += 1
+
+        mol = self.__pack(mol, x) # pack positions into new molecule class
+
+        if self.verbose:
+            self.__print_break(newline=False)
+            print('  START GEOMETRY OPTIMIZATION STEP %03i' % self.iter)
+            self.__print_break()
+            self.__print_positions(mol)
+            print() # newline
+
+        res = HF().rhf(mol, basis, orbc_init=self.cinit, calc_forces=True)
+
+        if self.verbose:
+            self.__print_energies(res)
+            print() # print newline
+
+        # cache matrices
+        self.cinit = res['orbc']
+        self.P = res['density']
+        self.orbe = res['orbe']
+        self.forces = res['forces']
+        self.coord = x
+
+        # append history
+        self.forces_history.append(self.forces)
+        self.coordinates_history.append(self.coord.reshape(len(mol.atoms),3))
+        self.energies_history.append(res['energies'][-1])
+
+        if self.verbose:
+            self.__print_forces(mol, self.forces)
+            print() # print newline
+            self.__print_break(newline=False)
+            print('  END GEOMETRY OPTIMIZATION STEP %03i' % self.iter)
+            self.__print_break(newline=True)
+
+        return res['energies'][-1]
+
+    def jacobian(self, x, mol, basis):
+        """
+        Calculate the forces. Note that the forces are typically already
+        calculated in the energy step so here only a simple check is done
+        to see if the coordinates match. If so, the forces are simply returned.
+        """
+        # check if forces are already calculated
+        if np.max(x - self.coord) < 1e-5:
+            return self.forces.flatten()
+        else:
+            res = HF().rhf(mol, basis, orbc_init=self.cinit, calc_forces=True)
+            return res['forces'].flatten()
+
+    def __unpack(self, mol):
+        """
+        Unpack coordinates from molecule class
+        """
+        coords = []
+        for atom in mol.atoms:
+            coords.append(atom[1])
+
+        return np.array(coords).flatten()
+
+    def __pack(self, mol, coords):
+        """
+        Pack coordinates into new molecule class
+        """
+
+        newmol = Molecule(mol.name)
+        coords = coords.reshape((len(mol.atoms), 3))
+        for i,atom in enumerate(mol.atoms):
+            newmol.add_atom(mol.atoms[i][0], coords[i][0], coords[i][1], coords[i][2])
+
+        return newmol
+
+    def __print_positions(self, mol):
+        """
+        Print atomic positions in nice formatting
+        """
+        print('-------------')
+        print('  POSITIONS  ')
+        print('-------------')
+        for atom in mol.atoms:
+            print('  %2s %12.8f %12.8f %12.8f' % (atom[0],
+                                                  atom[1][0],
+                                                  atom[1][1],
+                                                  atom[1][2]))
+
+    def __print_forces(self, mol, forces):
+        """
+        Print forces using nice formatting
+        """
+        print('----------')
+        print('  FORCES  ')
+        print('----------')
+
+        for atom,force in zip(mol.atoms, forces):
+            print('  %2s %12.4e %12.4e %12.4e' % (atom[0],
+                                                  force[0],
+                                                  force[1],
+                                                  force[2]))
+
+    def __print_energies(self, res):
+        """
+        Print the energy terms in each iteration
+        """
+        print('------------')
+        print('  ENERGIES  ')
+        print('------------')
+
+        print('  Kinetic:                     %12.8f' % res['ekin'])
+        print('  Nuclear:                     %12.8f' % res['enuc'])
+        print('  Electron-electron repulsion: %12.8f' % res['erep'])
+        print('  Exchange:                    %12.8f' % res['ex'])
+        print('  Nuclear repulsion:           %12.8f' % res['enucrep'])
+        print('  TOTAL:                       %12.8f' % res['energies'][-1])
+
+    def __print_break(self, newline=True):
+        print('=' * 80)
+        if newline:
+            print()
diff --git a/pyqint/hf.py b/pyqint/hf.py
index cfc1c8f..92bf50f 100644
--- a/pyqint/hf.py
+++ b/pyqint/hf.py
@@ -2,7 +2,6 @@
 
 import json
 import os
-from .cgf import cgf
 import numpy as np
 from . import PyQInt
 import time
@@ -16,7 +15,8 @@ class HF:
     Routines to perform a restricted Hartree-Fock calculations
     """
     def rhf(self, mol, basis, calc_forces=False, itermax=100,
-            use_diis=True, verbose=False, tolerance=1e-7):
+            use_diis=True, verbose=False, tolerance=1e-7,
+            orbc_init=None):
         """
         Performs a Hartree-Fock type calculation
 
@@ -32,14 +32,23 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100,
         # build cgfs, nuclei and calculate nr of electrons
         cgfs, nuclei = mol.build_basis(basis)
         nelec = int(np.sum([at[1] for at in nuclei]))
+        N = len(cgfs)
+        occ = [2 if i < nelec//2 else 0 for i in range(N)]
 
         # build integrals
         integrator = PyQInt()
         start = time.time()
-        S, T, V, teint = integrator.build_integrals_openmp(cgfs, nuclei)
+        S, T, V, teints = integrator.build_integrals_openmp(cgfs, nuclei)
         end = time.time()
         time_stats['integral_evaluation'] = end - start
 
+        # calculate nuclear repulsion
+        nuc_rep = 0.0
+        for i in range(0, len(nuclei)):
+            for j in range(i+1, len(nuclei)):
+                r = np.linalg.norm(np.array(nuclei[i][0]) - np.array(nuclei[j][0]))
+                nuc_rep += nuclei[i][1] * nuclei[j][1] / r
+
         # diagonalize S
         s, U = np.linalg.eigh(S)
 
@@ -47,7 +56,10 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100,
         X = U.dot(np.diag(1.0/np.sqrt(s)))
 
         # create empty P matrix as initial guess
-        P = np.zeros(S.shape)
+        if orbc_init is None:
+            P = np.zeros(S.shape)
+        else:
+            P = np.einsum('ik,jk,k->ij', orbc_init, orbc_init, occ)
 
         # keep track of time
         start = time.time()
@@ -65,27 +77,28 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100,
             iterstart = time.time()
 
             if niter > SUBSPACE_START and use_diis:
-                diis_coeff = self.calculate_diis_coefficients(evs_diis)
+                try:
+                    diis_coeff = self.calculate_diis_coefficients(evs_diis)
+                except np.linalg.LinAlgError:
+                    # stop diis procedure and revert to linear stepping
+                    use_diis = False
+                    continue
 
                 F = self.extrapolate_fock_from_diis_coefficients(fmats_diis, diis_coeff)
                 Fprime = X.transpose().dot(F).dot(X)
                 e, Cprime = np.linalg.eigh(Fprime)
                 C = X.dot(Cprime)
-                P = np.zeros(S.shape)
-                for i in range(S.shape[0]):
-                    for j in range(S.shape[0]):
-                        for k in range(0,int(nelec/2)):
-                            P[i,j] += 2.0 * C[i,k] * C[j,k]
+                P = np.einsum('ik,jk,k->ij', C, C, occ)
 
             # calculate G
-            G = np.zeros(S.shape)
-            for i in range(S.shape[0]):
-                for j in range(S.shape[0]):
-                    for k in range(S.shape[0]):
-                        for l in range(S.shape[0]):
+            G = np.zeros_like(S)
+            for i in range(N):
+                for j in range(N):
+                    for k in range(N):
+                        for l in range(N):
                             idx_rep = integrator.teindex(i,j,l,k)
                             idx_exc = integrator.teindex(i,k,l,j)
-                            G[i,j] += P[k,l] * (teint[idx_rep] - 0.5 * teint[idx_exc])
+                            G[i,j] += P[k,l] * (teints[idx_rep] - 0.5 * teints[idx_exc])
 
             # build Fock matrix
             F = T + V + G
@@ -106,11 +119,8 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100,
                 for j in range(S.shape[0]):
                     energy += 0.5 * P[j,i] * M[i,j]
 
-            # calculate repulsion of the nuclei
-            for i in range(0, len(nuclei)):
-                for j in range(i+1, len(nuclei)):
-                    r = np.linalg.norm(np.array(nuclei[i][0]) - np.array(nuclei[j][0]))
-                    energy += nuclei[i][1] * nuclei[j][1] / r
+            # add nuclear repulsion
+            energy += nuc_rep
 
             # store energy for next iteration
             energies.append(energy)
@@ -169,6 +179,15 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100,
         end = time.time()
         time_stats['self_convergence'] = end - start
 
+        # build two-electron integral tensor
+        tetens = np.zeros((N,N,N,N))
+        for i in range(N):
+            for j in range(N):
+                for k in range(N):
+                    for l in range(N):
+                        idx = integrator.teindex(i,j,k,l)
+                        tetens[i,j,k,l] = teints[idx]
+
         # build solution dictionary
         sol = {
             "energy": energies[-1],
@@ -184,9 +203,15 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100,
             "kinetic": T,
             "nuclear": V,
             'hcore': T+V,
+            'tetens': tetens,
             "time_stats" : time_stats,
             "ecore": np.sum(P * (T + V)),
-            "teint": teint,
+            "ekin": np.einsum('ij,ji', T, P),
+            "enuc": np.einsum('ij,ji', V, P),
+            "erep": np.einsum('ijlk,ij,kl', tetens, P, P),
+            "ex": np.einsum('ijlk,ij,kl', tetens, P, P),
+            "enucrep": nuc_rep,
+            "teint": teints,
             "nelec": nelec,
             "forces": self.rhf_forces(mol, basis, C, P, orbe) if calc_forces else None
         }
@@ -242,70 +267,46 @@ def rhf_forces(self, mol, basis, C, P, e):
         # intialization
         integrator = PyQInt()
         cgfs, nuclei = mol.build_basis(basis)
-        n_elec = np.sum([nucleus[1] for nucleus in nuclei])
-        forces = np.zeros(shape = [len(mol.atoms), 3])
-        n = len(cgfs)
+        nelec = np.sum([nucleus[1] for nucleus in nuclei])
+        forces = np.zeros((len(nuclei),3))
+        N = len(cgfs)
+        occ = [2 if i < nelec//2 else 0 for i in range(N)]
     
-        # calculate energy weighted density matrix:
-        # It might be preferrable to recalculate P as well
-        ew_density = np.zeros(shape = [n, n])
-        for i in range(n):
-            for j in range(n):
-                for k in range(0, n_elec//2):
-                    ew_density[i,j] += 2.0 * e[k] * C[i,k] * C[j,k]
+        # calculate energy weighted density matrix
+        ew_density = np.einsum('ik,jk,k,k->ij', C, C, occ, e)
+
+        # collect derivatives
+        S, T, V, teints = integrator.build_geometric_derivatives_openmp(cgfs, nuclei)
 
         # Loop over all cartesian direction for every nucleus
         # This could be made more efficient by incorporating symmetry
-        for n_nuc, deriv_nucleus in enumerate(nuclei):            
-            for deriv_direction in range(3): 
-                
-                overlap = np.zeros(shape=[n, n])
-                kinetic = np.zeros(shape=[n, n])
-                nuclear = np.zeros(shape=[n, n])
-                repulsion = np.zeros(shape=[n, n, n, n])
-
-                # Loop over every permutation of cgfs in the basis
-                for i, cgf_1 in enumerate(cgfs):
-                    for j, cgf_2 in enumerate(cgfs):
-
-                        # derivative of overlap matrix
-                        overlap[i, j] += integrator.overlap_deriv(cgf_1, cgf_2, deriv_nucleus[0], deriv_direction)
-
-                        # derivative of kinetic matrix
-                        kinetic[i, j] += integrator.kinetic_deriv(cgf_1, cgf_2, deriv_nucleus[0], deriv_direction)
-                        
-                        # derivative nuclear electron attraction
-                        for nucleus in mol.nuclei:
-                            
-                            # nuclear_deriv_op returns wrong values if both cgfs and nucleus is on deriv_nucleus
-                            # the following if statment eliminates that term using translational symmetry 
-                            if not (np.linalg.norm(cgf_1.p - deriv_nucleus[0]) < 0.0001 and  np.linalg.norm(cgf_2.p - deriv_nucleus[0]) < 0.0001 and np.linalg.norm(nucleus[0] - deriv_nucleus[0]) < 0.0001):
-                                nuclear[i, j] += integrator.nuclear_deriv(cgf_1, cgf_2, nucleus[0], nucleus[1], deriv_nucleus[0], deriv_direction)
-
-                        # derivative of electron-electron repulsions
-                        for k, cgf_3 in enumerate(cgfs):
-                            for l, cgf_4 in enumerate(cgfs):
-                                repulsion[i, j, k, l] += integrator.repulsion_deriv(cgf_1, cgf_2, cgf_3, cgf_4, deriv_nucleus[0], deriv_direction)
-                        
+        for n, deriv_nucleus in enumerate(nuclei):
+            for d in range(3):
+
                 # derivate nucleus-nucleus repulsion
                 term_nn = 0
                 for nucleus in nuclei:
                     if np.linalg.norm(nucleus[0] - deriv_nucleus[0]) > 0.0001:
-                        term_nn += deriv_nucleus[1] * nucleus[1] * (nucleus[0][deriv_direction] - deriv_nucleus[0][deriv_direction]) / (np.linalg.norm(nucleus[0]- deriv_nucleus[0])**3)
+                        term_nn += deriv_nucleus[1] * nucleus[1] * \
+                                   (nucleus[0][d] - deriv_nucleus[0][d]) / \
+                                   (np.linalg.norm(nucleus[0]- deriv_nucleus[0])**3)
 
-                hcore = kinetic + nuclear
+                hcore = T[n,d] + V[n,d]
                 term_hcore = np.sum(np.multiply(P, hcore))
 
                 term_repulsion = 0
-                for i in range(n):
-                    for j in range(n):
-                        for k in range(n):
-                            for l in range(n):
-                                term_repulsion += 1/2 * P[i, j] * P[k, l]*(repulsion[i, j, k, l] - 1/2*repulsion[i, k, j, l]) 
+                for i in range(N):
+                    for j in range(N):
+                        for k in range(N):
+                            for l in range(N):
+                                idx1 = integrator.teindex(i,j,k,l)
+                                idx2 = integrator.teindex(i,k,j,l)
+
+                                term_repulsion += 0.5 * P[i,j] * P[k,l] * (teints[n,d,idx1] - 0.5 * teints[n,d,idx2])
                 
-                term_overlap = - np.sum(np.multiply(ew_density, overlap))
+                term_overlap = - np.sum(np.multiply(ew_density, S[n,d]))
            
                 # F = - d/dR E
-                forces[n_nuc, deriv_direction] -= term_hcore + term_repulsion + term_overlap + term_nn 
+                forces[n,d] = term_hcore + term_repulsion + term_overlap + term_nn
 
         return forces
diff --git a/pyqint/integrals.cpp b/pyqint/integrals.cpp
index 6b502cc..e727dae 100644
--- a/pyqint/integrals.cpp
+++ b/pyqint/integrals.cpp
@@ -56,7 +56,7 @@ std::vector<double> Integrator::evaluate_cgfs(const std::vector<CGF>& cgfs,
             S(i,j) = S(j,i) = this->overlap(cgfs[i], cgfs[j]);
             T(i,j) = T(j,i) = this->kinetic(cgfs[i], cgfs[j]);
             for(unsigned int k=0; k<charges.size(); k++) {
-                // there is a race condition here!!
+                // this is to avoid the race condition
                 Vn[k](i,j) = Vn[k](j,i) = this->nuclear(cgfs[i], cgfs[j], vec3(px[k], py[k], pz[k]), charges[k]);
             }
         }
@@ -127,6 +127,124 @@ std::vector<double> Integrator::evaluate_cgfs(const std::vector<CGF>& cgfs,
     return results;
 }
 
+/**
+ * @brief      Evaluate the geometric derivates for all cgfs in buffer
+ */
+std::vector<double> Integrator::evaluate_geometric_derivatives(const std::vector<CGF>& cgfs,
+                                                               const std::vector<int>& charges,
+                                                               const std::vector<double>& px,
+                                                               const std::vector<double>& py,
+                                                               const std::vector<double>& pz) const {
+    size_t sz = cgfs.size();
+
+    // Construct 2x2 matrices to hold values for the overlap,
+    // kinetic and two nuclear integral values, respectively.
+    auto S = std::vector<Eigen::MatrixXd>(charges.size() * 3, Eigen::MatrixXd::Zero(sz, sz));
+    auto T = std::vector<Eigen::MatrixXd>(charges.size() * 3, Eigen::MatrixXd::Zero(sz, sz));
+    auto V = std::vector<Eigen::MatrixXd>(charges.size() * 3, Eigen::MatrixXd::Zero(sz, sz));
+
+    // calculate the integral values using the integrator class
+    for(unsigned int n=0; n<charges.size(); n++) { // loop over nuclei
+        for(unsigned int k=0; k<3; k++) { // loop over directions
+
+            // build container for summation over nuclei
+            std::vector<Eigen::MatrixXd> Vn(charges.size(), Eigen::MatrixXd::Zero(sz, sz));
+
+            #pragma omp parallel for schedule(dynamic)
+            for(int i=0; i<(int)sz; i++) {  // have to use signed int for MSVC OpenMP here
+                for(int j=0; j<(int)sz; j++) {
+                    S[n*3+k](i,j) = this->overlap_deriv(cgfs[i], cgfs[j], vec3(px[n], py[n], pz[n]), k);
+                    T[n*3+k](i,j) = this->kinetic_deriv(cgfs[i], cgfs[j], vec3(px[n], py[n], pz[n]), k);
+
+
+                    for(unsigned int l=0; l<charges.size(); l++) {
+                        Vn[l](i,j) = this->nuclear_deriv(cgfs[i], cgfs[j],
+                                                         vec3(px[l], py[l], pz[l]),
+                                                         charges[l],
+                                                         vec3(px[n], py[n], pz[n]),
+                                                         k);
+                    }
+                }
+            }
+
+            // combine all nuclear attraction integrals
+            for(unsigned int i=0; i<sz; i++) {
+                for(unsigned int j=0; j<sz; j++) {
+                    for(unsigned int l=0; l<charges.size(); l++) {
+                        V[n*3+k](i,j) += Vn[l](i,j);
+                    }
+                }
+            }
+        }
+    }
+
+    // // calculate all two-electron integrals
+    const unsigned int max_teints = this->teindex(sz-1,sz-1,sz-1,sz-1) + 1;
+    std::vector<double> tedouble(max_teints * charges.size() * 3, -1.0);
+
+    // it is more efficient to first 'unroll' the fourfold nested loop
+    // into a single vector of jobs to execute
+    std::vector<std::array<size_t, 7>> jobs;
+    for(size_t n=0; n<charges.size(); n++) {        // nuclear derivatives
+        for(size_t d=0; d<3; d++) {                 // directions
+            for(size_t i=0; i<sz; i++) {
+                for(size_t j=0; j<sz; j++) {
+                    size_t ij = i*(i+1)/2 + j;
+                    for(size_t k=0; k<sz; k++) {
+                        for(size_t l=0; l<sz; l++) {
+                            size_t kl = k * (k+1)/2 + l;
+                            if(ij <= kl) {
+                                size_t idx = this->teindex(i,j,k,l) + (n*3 + d) * max_teints;
+
+                                if(idx >= tedouble.size()) {
+                                    throw std::runtime_error("Process tried to access illegal array position");
+                                }
+
+                                if(tedouble[idx] < 0.0) {
+                                    tedouble[idx] = 1.0;
+                                    jobs.push_back({idx, i, j, k, l, n, d});
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // // evaluate jobs
+    #pragma omp parallel for schedule(dynamic)
+    for(int s=0; s<(int)jobs.size(); s++) {  // have to use signed int for MSVC OpenMP here
+        const size_t idx = jobs[s][0];
+        const size_t i = jobs[s][1];
+        const size_t j = jobs[s][2];
+        const size_t k = jobs[s][3];
+        const size_t l = jobs[s][4];
+        const size_t n = jobs[s][5];
+        const size_t d = jobs[s][6];
+        tedouble[idx] = this->repulsion_deriv(cgfs[i], cgfs[j], cgfs[k], cgfs[l], vec3(px[n], py[n], pz[n]), d);
+    }
+
+    // package everything into results vector, will be unpacked in
+    // connected Python class
+    std::vector<double> results(charges.size() * 3 * sz * sz * 3 + tedouble.size());
+    for(size_t n=0; n<charges.size(); n++) { // nuclear derivatives
+        for(size_t d=0; d<3; d++) { // directions
+            size_t spos = (n * 3 + d) * sz * sz;
+            std::memcpy(&results[spos], S[n*3+d].data(), sz * sz * sizeof(double));
+
+            size_t tpos = spos + charges.size() * 3 * sz * sz;
+            std::memcpy(&results[tpos], T[n*3+d].data(), sz * sz * sizeof(double));
+
+            size_t vpos = tpos + charges.size() * 3 * sz * sz;
+            std::memcpy(&results[vpos], V[n*3+d].data(), sz * sz * sizeof(double));
+        }
+    }
+    std::memcpy(&results[charges.size() * 3 * sz * sz * 3], tedouble.data(), tedouble.size() * sizeof(double));
+
+    return results;
+}
+
 /**************************************************************************
  * OVERLAP INTEGRALS
  **************************************************************************/
@@ -549,84 +667,60 @@ double Integrator::nuclear_gto(const GTO& gto1,
  *
  * Calculates the value of < cgf1 | V | cgf2 >
  *
+ * In contrast to the other derivatives, this will be done using a finite
+ * difference procedure.
+ *
  * @return double value of the nuclear integral
  */
 double Integrator::nuclear_deriv(const CGF& cgf1, const CGF& cgf2, const vec3 &nucleus, unsigned int charge,
                                  const vec3& nucderiv, unsigned int coord) const {
-    double sum = 0.0;
 
     // check if cgf originates from nucleus
     bool n1 = (cgf1.get_r() - nucderiv).squaredNorm() < 0.0001;
     bool n2 = (cgf2.get_r() - nucderiv).squaredNorm() < 0.0001;
     bool n3 = (nucleus - nucderiv).squaredNorm() < 0.0001;
 
-    for(unsigned int k = 0; k < cgf1.size(); k++) {
-        for(unsigned int l = 0; l < cgf2.size(); l++) {
-
-            // take the derivative towards the basis functions
-            double t1 = n1 ? this->nuclear_deriv_bf(cgf1.get_gto(k), cgf2.get_gto(l), nucleus, coord) : 0.0;
-            double t2 = n2 ? this->nuclear_deriv_bf(cgf2.get_gto(l), cgf1.get_gto(k), nucleus, coord) : 0.0;
-
-            // take the derivative of the operator towards the coordinate
-            double t3 = n3 ? this->nuclear_deriv_op(cgf1.get_gto(k), cgf2.get_gto(l), nucleus, coord) : 0.0;
-
-            sum += cgf1.get_norm_gto(k) *
-                   cgf2.get_norm_gto(l) *
-                   cgf1.get_coefficient_gto(k) *
-                   cgf2.get_coefficient_gto(l) *
-                   (t1 + t2 + t3);
-        }
+    if(n1 == n2 && n2 == n3) {
+        return 0.0;
     }
 
-    return sum * (double)charge;
-}
+    static const double delta = 1e-5;
+
+    // create copies of objects
+    CGF cgf1m = cgf1;
+    CGF cgf1p = cgf1;
+    CGF cgf2m = cgf2;
+    CGF cgf2p = cgf2;
+    vec3 nucm = nucleus;
+    vec3 nucp = nucleus;
+
+    if(n1) {
+        vec3 rm = cgf1.get_r();
+        rm[coord] -= delta;
+        cgf1m.set_position(rm);
+        vec3 rp = cgf1.get_r();
+        rp[coord] += delta;
+        cgf1p.set_position(rp);
+    }
 
-/**
- * @brief Calculates nuclear integral of two CGF
- *
- * @param const GTO& gto1       Contracted Gaussian Function
- * @param const GTO& gto2       Contracted Gaussian Function
- * @param unsigned int charge   charge of the nucleus in a.u.
- *
- * Calculates the value of < d/dc * gto1 | V | gto2 >
- *
- * @return double value of the nuclear integral
- */
-double Integrator::nuclear_deriv_bf(const GTO& gto1, const GTO& gto2, const vec3 &nucleus, unsigned int coord) const {
-    std::array<unsigned int, 3> gto_ang = {gto1.get_l(), gto1.get_m(), gto1.get_n()};
-    if(gto_ang[coord] != 0) {
-        gto_ang[coord] += 1; // calculate l+1 term
-        double term_plus = this->nuclear(gto1.get_position(), gto_ang[0], gto_ang[1], gto_ang[2], gto1.get_alpha(),
-                                         gto2.get_position(), gto2.get_l(), gto2.get_m(), gto2.get_n(), gto2.get_alpha(), nucleus);
-        gto_ang[coord] -= 2; // calculate l-1 term
-        double term_min = this->nuclear(gto1.get_position(), gto_ang[0], gto_ang[1], gto_ang[2], gto1.get_alpha(),
-                                        gto2.get_position(), gto2.get_l(), gto2.get_m(), gto2.get_n(), gto2.get_alpha(), nucleus);
-        gto_ang[coord] += 1; // recover l
+    if(n2) {
+        vec3 rm = cgf2.get_r();
+        rm[coord] -= delta;
+        cgf2m.set_position(rm);
+        vec3 rp = cgf2.get_r();
+        rp[coord] += delta;
+        cgf2p.set_position(rp);
+    }
 
-        return 2.0 * gto1.get_alpha() * term_plus - gto_ang[coord] * term_min;
-    } else { // s-type GTO
-        gto_ang[coord] += 1;
-        double term1 = this->nuclear(gto1.get_position(), gto_ang[0], gto_ang[1], gto_ang[2], gto1.get_alpha(),
-                                     gto2.get_position(), gto2.get_l(), gto2.get_m(), gto2.get_n(), gto2.get_alpha(), nucleus);
-        return 2.0 * gto1.get_alpha() * term1;
+    if(n3) {
+        nucm[coord] -= delta;
+        nucp[coord] += delta;
     }
-}
 
-/**
- * @brief Calculates nuclear integral of two CGF
- *
- * @param const GTO& gto1       Contracted Gaussian Function
- * @param const GTO& gto2       Contracted Gaussian Function
- * @param unsigned int charge   charge of the nucleus in a.u.
- *
- * Calculates the value of < gto1 | V | gto2 >
- *
- * @return double value of the nuclear integral
- */
-double Integrator::nuclear_deriv_op(const GTO& gto1, const GTO& gto2, const vec3 &nucleus, unsigned int coord) const {
-    return nuclear_deriv_op(gto1.get_position(), gto1.get_l(), gto1.get_m(), gto1.get_n(), gto1.get_alpha(),
-                            gto2.get_position(), gto2.get_l(), gto2.get_m(), gto2.get_n(), gto2.get_alpha(),
-                            nucleus, coord);
+    double vm = this->nuclear(cgf1m, cgf2m, nucm, charge);
+    double vp = this->nuclear(cgf1p, cgf2p, nucp, charge);
+
+    return (vp - vm) / (2.0 * delta);
 }
 
 /**************************************************************************
@@ -658,7 +752,7 @@ double Integrator::repulsion(const CGF &cgf1,const CGF &cgf2,const CGF &cgf3,con
                     const double n4 = cgf4.get_norm_gto(l);
 
                     double pre = cgf1.get_coefficient_gto(i) * cgf2.get_coefficient_gto(j) * cgf3.get_coefficient_gto(k) * cgf4.get_coefficient_gto(l);
-                    sum += n1 * n2 * n3 * n4 * pre * repulsion(cgf1.get_gto(i), cgf2.get_gto(j), cgf3.get_gto(k), cgf4.get_gto(l));
+                    sum += n1 * n2 * n3 * n4 * pre * this->repulsion(cgf1.get_gto(i), cgf2.get_gto(j), cgf3.get_gto(k), cgf4.get_gto(l));
                 }
             }
         }
@@ -964,64 +1058,6 @@ double Integrator::nuclear(const vec3& a, int l1, int m1, int n1, double alpha1,
     return -2.0 * pi / gamma * std::exp(-alpha1*alpha2*rab2/gamma) * sum;
 }
 
-double Integrator::nuclear_deriv_op(const vec3& a, int l1, int m1, int n1, double alpha1,
-                                    const vec3& b, int l2, int m2, int n2,
-                                    double alpha2, const vec3& c, unsigned int coord) const {
-
-    static const double pi = boost::math::constants::pi<double>();
-
-    double gamma = alpha1 + alpha2;
-
-    vec3 p = gaussian_product_center(alpha1, a, alpha2, b);
-    double rab2 = (a-b).squaredNorm();
-    double rcp2 = (c-p).squaredNorm();
-    double rcpcoord = (c-p)[coord];
-
-    std::vector<double> ax = A_array(l1, l2, p[0]-a[0], p[0]-b[0], p[0]-c[0], gamma);
-    std::vector<double> ay = A_array(m1, m2, p[1]-a[1], p[1]-b[1], p[1]-c[1], gamma);
-    std::vector<double> az = A_array(n1, n2, p[2]-a[2], p[2]-b[2], p[2]-c[2], gamma);
-
-    // build array of doubles for the derivatives towards C[coord] appearing
-    // in each of the terms
-    std::vector<double> ad;
-    switch(coord) {
-        case 0: // x coordinate
-            ad = A_array_deriv(l1, l2, p[0]-a[0], p[0]-b[0], p[0]-c[0], gamma);
-        break;
-        case 1: // y coordinate
-            ad = A_array_deriv(m1, m2, p[1]-a[1], p[1]-b[1], p[1]-c[1], gamma);
-        break;
-        case 2: // z coordinate
-            ad = A_array_deriv(n1, n2, p[2]-a[2], p[2]-b[2], p[2]-c[2], gamma);
-        break;
-    }
-
-    double sum = 0.0;
-
-    // build arrays of all values; based on which coordinate derivative is requested,
-    // different indices are generated to ensure by which all different cases
-    // can be handled in one set of nested loops
-    const std::array<int, 3> itmax = {l1 + l2, m1 + m2, n1 + n2};
-    const std::array<std::vector<double>, 3> v = {ax, ay, az};
-    const unsigned int v0 = coord;
-    const unsigned int v1 = (coord+1)%3;
-    const unsigned int v2 = (coord+2)%3;
-
-    for(int i=0; i<=itmax[v0];i++) {
-        for(int j=0; j<=itmax[v1];j++) {
-            for(int k=0; k<=itmax[v2];k++) {
-
-                // apply product rule as both the terms as well as the incomplete gamma function
-                // have terms in C[coord]
-                sum += (v[v0][i] * -2.0 * gamma * rcpcoord * this->gamma_inc.Fgamma(i+j+k+1,rcp2*gamma)
-                        + ad[i] * this->gamma_inc.Fgamma(i+j+k,rcp2*gamma)) * v[v1][j] * v[v2][k];
-            }
-        }
-    }
-
-    return -2.0 * pi / gamma * std::exp(-alpha1*alpha2*rab2/gamma) * sum;
-}
-
 std::vector<double> Integrator::A_array(const int l1, const int l2, const double pa, const double pb, const double cp, const double g) const {
     int imax = l1 + l2 + 1;
     std::vector<double> arrA(imax, 0);
@@ -1038,29 +1074,6 @@ std::vector<double> Integrator::A_array(const int l1, const int l2, const double
     return arrA;
 }
 
-std::vector<double> Integrator::A_array_deriv(const int l1, const int l2, const double pa, const double pb, const double cp, const double g) const {
-    int imax = l1 + l2 + 1;
-    std::vector<double> arrA(imax, 0);
-
-    for(int i=0; i<imax; i++) {
-        for(int r=0; r<=i/2; r++) {
-            for(int u=0; u<=(i-2*r)/2; u++) {
-                int iI = i - 2 * r - u;
-                int cppow = i-2*r-2*u; // power in the coefficient cp
-
-                double term = A_term(i, r, u, l1, l2, pa, pb, cp, g);
-
-                if(cppow != 0 && cp != 0.0) {
-                    term *= -1.0 * cppow / cp;
-                    arrA[iI] += term;
-                }
-            }
-        }
-    }
-
-    return arrA;
-}
-
 double Integrator::A_term(const int i, const int r, const int u, const int l1, const int l2, const double pax, const double pbx, const double cpx, const double gamma) const {
     return  std::pow(-1,i) * this->binomial_prefactor(i,l1,l2,pax,pbx)*
             std::pow(-1,u) * boost::math::factorial<double>(i)*std::pow(cpx,i-2*r-2*u)*
@@ -1116,7 +1129,7 @@ std::vector<double> Integrator::B_array(const int l1,const int l2,const int l3,c
                 for(int r2=0; r2 < i2/2+1; r2++) {
                     for(int u=0; u<(i1+i2)/2-r1-r2+1; u++) {
                         int i = i1+i2-2*(r1+r2)-u;
-                        arrB[i] += B_term(i1,i2,r1,r2,u,l1,l2,l3,l4,
+                        arrB[i] += this->B_term(i1,i2,r1,r2,u,l1,l2,l3,l4,
                                                             p,a,b,q,c,d,g1,g2,delta);
                     }
                 }
diff --git a/pyqint/integrals.h b/pyqint/integrals.h
index f7b73b6..838cf60 100644
--- a/pyqint/integrals.h
+++ b/pyqint/integrals.h
@@ -23,8 +23,8 @@
 
 #include <boost/math/special_functions/binomial.hpp>
 #include <boost/math/special_functions/gamma.hpp>
-#include <boost/lexical_cast.hpp>
 #include <boost/math/constants/constants.hpp>
+#include <boost/lexical_cast.hpp>
 #include <string>
 #include <unordered_map>
 
@@ -115,6 +115,15 @@ class Integrator {
                                       const std::vector<double>& py,
                                       const std::vector<double>& pz) const;
 
+    /**
+     * @brief      Evaluate all integrals for cgfs in buffer
+     */
+    std::vector<double> evaluate_geometric_derivatives(const std::vector<CGF>& cgfs,
+                                                       const std::vector<int>& charges,
+                                                       const std::vector<double>& px,
+                                                       const std::vector<double>& py,
+                                                       const std::vector<double>& pz) const;
+
     /**************************************************************************
      * OVERLAP INTEGRALS
      **************************************************************************/
@@ -383,73 +392,6 @@ class Integrator {
         return this->nuclear_deriv(cgf1, cgf2, vec3(cx, cy, cz), charge, vec3(dx, dy, dz), coord);
     }
 
-    /**
-     * @brief Calculates derivative towards one of the basis functions in the
-     * nuclear integral of two GTOs
-     *
-     * @param const GTO& gto1       Contracted Gaussian Function
-     * @param const GTO& gto2       Contracted Gaussian Function
-     * @param const vec3 nucleus    Position of the nucleus
-     * @param unsigned int charge   charge of the nucleus in a.u.
-     * @param unsigned int coord    Cartesian direction to take derivative towards
-     *
-     * Calculates the value of < d/dx1 gto1 | V | gto2 >
-     *
-     * @return double value of the nuclear integral
-     */
-    double nuclear_deriv_bf(const GTO &gto1, const GTO &gto2, const vec3& nucleus, unsigned int coord) const;
-
-    /**
-     * @brief Calculates derivative of the nuclear attraction operator
-     * in a nuclear attraction integral of two CGF
-     *
-     * @param const GTO& gto1       Contracted Gaussian Function
-     * @param const GTO& gto2       Contracted Gaussian Function
-     * @param const vec3 nucleus    Position of the nucleus
-     * @param unsigned int charge   charge of the nucleus in a.u.
-     * @param unsigned int coord    Cartesian direction to take derivative towards
-     *
-     * Calculates the value of < gto1 | d/dcx V | gto2 >
-     *
-     * @return double value of the nuclear integral
-     */
-    inline double nuclear_deriv_bf(const GTO &gto1, const GTO &gto2, double rx, double ry, double rz, unsigned int coord) const {
-        return this->nuclear_deriv_bf(gto1, gto2, vec3(rx, ry, rz), coord);
-    }
-
-    /**
-     * @brief Calculates nuclear integral of two CGF
-     *
-     * @param const GTO& gto1       Contracted Gaussian Function
-     * @param const GTO& gto2       Contracted Gaussian Function
-     * @param const vec3 nucleus    Position of the nucleus
-     * @param unsigned int charge   charge of the nucleus in a.u.
-     * @param unsigned int coord    Cartesian direction to take derivative towards
-     *
-     * Calculates the value of < gto1 | V | gto2 >
-     *
-     * @return double value of the nuclear integral
-     */
-    inline double nuclear_deriv_op(const GTO &gto1, const GTO &gto2, const vec3& nucleus, unsigned int coord) const;
-
-    /**
-     * @brief Calculates derivative of the nuclear attraction operator
-     * in a nuclear attraction integral of two CGF
-     *
-     * @param const GTO& gto1       Contracted Gaussian Function
-     * @param const GTO& gto2       Contracted Gaussian Function
-     * @param const vec3 nucleus    Position of the nucleus
-     * @param unsigned int charge   charge of the nucleus in a.u.
-     * @param unsigned int coord    Cartesian direction to take derivative towards
-     *
-     * Calculates the value of < gto1 | d/dcx V | gto2 >
-     *
-     * @return double value of the nuclear integral
-     */
-    inline double nuclear_deriv_op(const GTO &gto1, const GTO &gto2, double rx, double ry, double rz, unsigned int coord) const {
-        return this->nuclear_deriv_op(gto1, gto2, vec3(rx, ry, rz), coord);
-    }
-
     /**
      * @brief Calculates nuclear integral of two CGF
      *
diff --git a/pyqint/molecule.py b/pyqint/molecule.py
index a7700dc..d3caccb 100644
--- a/pyqint/molecule.py
+++ b/pyqint/molecule.py
@@ -22,7 +22,11 @@ def __str__(self):
         return res
 
     def add_atom(self, atom, x, y, z, unit='bohr'):
-        ang2bohr = 1.88973
+        ang2bohr = 1.8897259886
+
+        x = float(x)
+        y = float(y)
+        z = float(z)
 
         if unit == "bohr":
             self.atoms.append([atom, np.array([x, y, z])])
diff --git a/pyqint/molecule_builder.py b/pyqint/molecule_builder.py
new file mode 100644
index 0000000..f1f0d3d
--- /dev/null
+++ b/pyqint/molecule_builder.py
@@ -0,0 +1,36 @@
+import os
+from .molecule import Molecule 
+
+class MoleculeBuilder:
+    """
+    Class that builds molecules from templates or from point group descriptions
+    """
+    def __init__(self):
+        pass
+
+    def from_name(self, molname):
+        """
+        Build molecule from file and return it
+        """
+        fname = os.path.join(os.path.dirname(__file__), 'molecules', molname.lower() + '.xyz')
+        with open(fname, 'r') as f:
+            lines = f.readlines()
+            
+            nratoms = int(lines[0].strip())
+
+            mol = Molecule(molname)
+            for line in lines[2:2+nratoms]:
+                pieces = line.split()
+                mol.add_atom(pieces[0], float(pieces[1]), float(pieces[2]), float(pieces[3]), unit='angstrom')
+
+            return mol
+
+    def build_complex_td(self, r, at1, at2, unit='bohr'):
+        mol = Molecule()
+        mol.add_atom(at1,  0,  0,  0, unit=unit)
+        mol.add_atom(at2,  1,  1,  1, unit=unit)
+        mol.add_atom(at2,  1, -1, -1, unit=unit)
+        mol.add_atom(at2, -1, -1,  1, unit=unit)
+        mol.add_atom(at2, -1,  1,  1, unit=unit)
+
+        return mol
diff --git a/pyqint/molecules/benzene.xyz b/pyqint/molecules/benzene.xyz
new file mode 100644
index 0000000..ecc96c8
--- /dev/null
+++ b/pyqint/molecules/benzene.xyz
@@ -0,0 +1,14 @@
+12
+
+  C        0.0000000        1.4024230        0.0000000
+  C        1.2145340        0.7012110        0.0000000
+  C        1.2145340       -0.7012110        0.0000000
+  C        0.0000000       -1.4024230        0.0000000
+  C       -1.2145340       -0.7012110        0.0000000
+  C       -1.2145340        0.7012110        0.0000000
+  H        0.0000000        2.5016720        0.0000000
+  H        2.1665120        1.2508360        0.0000000
+  H        2.1665120       -1.2508360        0.0000000
+  H        0.0000000       -2.5016720        0.0000000
+  H       -2.1665120       -1.2508360        0.0000000
+  H       -2.1665120        1.2508360        0.0000000
diff --git a/pyqint/molecules/bf3.xyz b/pyqint/molecules/bf3.xyz
new file mode 100644
index 0000000..9d5dbf3
--- /dev/null
+++ b/pyqint/molecules/bf3.xyz
@@ -0,0 +1,6 @@
+4
+
+  B        0.0000000        0.0000000        0.0000000
+  F        0.0000000        1.3304260        0.0000000
+  F        1.1521830       -0.6652130        0.0000000
+  F       -1.1521830       -0.6652130        0.0000000
diff --git a/pyqint/molecules/ch4.xyz b/pyqint/molecules/ch4.xyz
new file mode 100644
index 0000000..07a9cc7
--- /dev/null
+++ b/pyqint/molecules/ch4.xyz
@@ -0,0 +1,7 @@
+5
+
+  C        0.0000000        0.0000000        0.0000000
+  H        0.6327670        0.6327670        0.6327670
+  H       -0.6327670       -0.6327670        0.6327670
+  H       -0.6327670        0.6327670       -0.6327670
+  H        0.6327670       -0.6327670       -0.6327670
diff --git a/pyqint/molecules/co.xyz b/pyqint/molecules/co.xyz
new file mode 100644
index 0000000..f3257ab
--- /dev/null
+++ b/pyqint/molecules/co.xyz
@@ -0,0 +1,4 @@
+2
+
+  O        0.0000000        0.0000000        0.5120840
+  C        0.0000000        0.0000000       -0.6827790
diff --git a/pyqint/molecules/co2.xyz b/pyqint/molecules/co2.xyz
new file mode 100644
index 0000000..af66fa6
--- /dev/null
+++ b/pyqint/molecules/co2.xyz
@@ -0,0 +1,5 @@
+3
+
+  C        0.0000000        0.0000000        0.0000000
+  O        0.0000000        0.0000000        1.2292730
+  O        0.0000000        0.0000000       -1.2292730
diff --git a/pyqint/molecules/ethylene.xyz b/pyqint/molecules/ethylene.xyz
new file mode 100644
index 0000000..9bf64aa
--- /dev/null
+++ b/pyqint/molecules/ethylene.xyz
@@ -0,0 +1,8 @@
+6
+
+  C        0.0000000        0.0000000        0.6655810
+  H        0.0000000        0.9321350        1.2456070
+  H        0.0000000       -0.9321350        1.2456070
+  C        0.0000000        0.0000000       -0.6655810
+  H        0.0000000       -0.9321350       -1.2456070
+  H        0.0000000        0.9321350       -1.2456070
diff --git a/pyqint/molecules/h2.xyz b/pyqint/molecules/h2.xyz
new file mode 100644
index 0000000..95de1fa
--- /dev/null
+++ b/pyqint/molecules/h2.xyz
@@ -0,0 +1,4 @@
+2
+
+  H        0.0000000        0.0000000        0.3679020
+  H        0.0000000        0.0000000       -0.3679020
diff --git a/pyqint/molecules/h2o.xyz b/pyqint/molecules/h2o.xyz
new file mode 100644
index 0000000..e540033
--- /dev/null
+++ b/pyqint/molecules/h2o.xyz
@@ -0,0 +1,5 @@
+3
+
+  O        0.0000000        0.0000000        0.1368220
+  H        0.0000000        0.7697660       -0.5472890
+  H        0.0000000       -0.7697660       -0.5472890
diff --git a/pyqint/molecules/he.xyz b/pyqint/molecules/he.xyz
new file mode 100644
index 0000000..71256b0
--- /dev/null
+++ b/pyqint/molecules/he.xyz
@@ -0,0 +1,3 @@
+1
+
+He 0.0 0.0 0.0
diff --git a/pyqint/molecules/lih.xyz b/pyqint/molecules/lih.xyz
new file mode 100644
index 0000000..a56821b
--- /dev/null
+++ b/pyqint/molecules/lih.xyz
@@ -0,0 +1,4 @@
+2
+
+  Li        0.0000000        0.0000000        0.3834270
+  H        0.0000000        0.0000000       -1.1502800
diff --git a/pyqint/molecules/nh3.xyz b/pyqint/molecules/nh3.xyz
new file mode 100644
index 0000000..9f35c83
--- /dev/null
+++ b/pyqint/molecules/nh3.xyz
@@ -0,0 +1,7 @@
+5
+
+  H        0.7034100        0.7034100        0.7034100
+  H       -0.7034100       -0.7034100        0.7034100
+  H       -0.7034100        0.7034100       -0.7034100
+  H        0.7034100       -0.7034100       -0.7034100
+  N        0.0000000        0.0000000        0.0000000
diff --git a/pyqint/pyqint.pxd b/pyqint/pyqint.pxd
index c1a2496..9caf48f 100644
--- a/pyqint/pyqint.pxd
+++ b/pyqint/pyqint.pxd
@@ -86,8 +86,12 @@ cdef extern from "integrals.h":
         # two-electron indexing
         int teindex(int, int, int, int) except +
 
+        # openmp routine for integral evaluation
         vector[double] evaluate_cgfs(vector[CGF], vector[int], vector[double], vector[double], vector[double]) except +
 
+        # openmp routine for geometric derivatives evaluation
+        vector[double] evaluate_geometric_derivatives(vector[CGF], vector[int], vector[double], vector[double], vector[double]) except +
+
         string get_compiler_version() except +
         string get_openmp_version() except +
         string get_compile_date() except +
diff --git a/pyqint/pyqint.pyx b/pyqint/pyqint.pyx
index c8e645a..1bc77e7 100644
--- a/pyqint/pyqint.pyx
+++ b/pyqint/pyqint.pyx
@@ -195,26 +195,6 @@ cdef class PyQInt:
 
         return self.integrator.nuclear_gto(c_gto1, c_gto2, rc[0], rc[1], rc[2])
 
-    def nuclear_gto_deriv_bf(self, gto1, gto2, rc, coord):
-
-        cdef GTO c_gto1
-        cdef GTO c_gto2
-
-        c_gto1 = GTO(gto1.c, gto1.p[0], gto1.p[1], gto1.p[2], gto1.alpha, gto1.l, gto1.m, gto1.n)
-        c_gto2 = GTO(gto2.c, gto2.p[0], gto2.p[1], gto2.p[2], gto2.alpha, gto2.l, gto2.m, gto2.n)
-
-        return self.integrator.nuclear_deriv_bf(c_gto1, c_gto2, rc[0], rc[1], rc[2], coord)
-
-    def nuclear_gto_deriv_op(self, gto1, gto2, rc, coord):
-
-        cdef GTO c_gto1
-        cdef GTO c_gto2
-
-        c_gto1 = GTO(gto1.c, gto1.p[0], gto1.p[1], gto1.p[2], gto1.alpha, gto1.l, gto1.m, gto1.n)
-        c_gto2 = GTO(gto2.c, gto2.p[0], gto2.p[1], gto2.p[2], gto2.alpha, gto2.l, gto2.m, gto2.n)
-
-        return self.integrator.nuclear_deriv_op(c_gto1, c_gto2, rc[0], rc[1], rc[2], coord)
-
     def nuclear(self, cgf1, cgf2, rc, zc):
 
         cdef CGF c_cgf1
@@ -409,6 +389,40 @@ cdef class PyQInt:
 
         return S, T, V, teint
 
+    def build_geometric_derivatives_openmp(self, cgfs, nuclei):
+        cdef vector[CGF] c_cgfs
+        cdef vector[int] charges
+        cdef vector[double] px
+        cdef vector[double] py
+        cdef vector[double] pz
+
+        # build CGFS objects
+        for cgf in cgfs:
+            c_cgfs.push_back(CGF(cgf.p[0], cgf.p[1], cgf.p[2]))
+            for gto in cgf.gtos:
+                c_cgfs.back().add_gto(gto.c, gto.alpha, gto.l, gto.m, gto.n)
+
+        # add nuclei to buffer
+        for nucleus in nuclei:
+            charges.push_back(nucleus[1])
+            px.push_back(nucleus[0][0])
+            py.push_back(nucleus[0][1])
+            pz.push_back(nucleus[0][2])
+
+        results = np.array(self.integrator.evaluate_geometric_derivatives(c_cgfs, charges, px, py, pz))
+
+        sz = len(cgfs)
+        ntei = self.teindex(sz-1,sz-1,sz-1,sz-1)+1 # calculate number of two-electron integrals
+
+        tsize = len(nuclei) * 3 * sz * sz
+        S = results[0:tsize].reshape((len(nuclei), 3, sz,sz))
+        T = results[tsize:2*tsize].reshape((len(nuclei), 3, sz,sz))
+        V = results[2*tsize:3*tsize].reshape((len(nuclei), 3, sz,sz))
+
+        teints = results[3*tsize:].reshape(len(nuclei), 3, ntei)
+
+        return S, T, V, teints
+
     def build_rectgrid3d(self, xmin, xmax, sz):
         """
         Build rectangular grid with z the slowest moving index
diff --git a/setup.py b/setup.py
index 5382ab6..d00e100 100644
--- a/setup.py
+++ b/setup.py
@@ -124,7 +124,7 @@ def find_windows_versions():
     ext_modules=cythonize(ext_modules[0],
                           language_level = "3",
                           build_dir="build"),
-    packages=['pyqint', 'pyqint.basissets'],
+    packages=['pyqint', 'pyqint.basissets', 'pyqint.molecules'],
     include_package_data=True,
     classifiers=[
         "Programming Language :: Python :: 3",


Test patch:
diff --git a/tests/results/nuclear_deriv_h2o.txt b/tests/results/nuclear_deriv_h2o.txt
new file mode 100644
index 0000000..53aea80
--- /dev/null
+++ b/tests/results/nuclear_deriv_h2o.txt
@@ -0,0 +1,441 @@
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+6.890920502529956714e-12
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.197937775332391652e+00
+-9.359368253303300689e-01
+0.000000000000000000e+00
+1.197937775332391652e+00
+9.359368253969434504e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+1.197937775332391652e+00
+-9.359368253303300689e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.197937775332391652e+00
+9.359368253969434504e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+2.601951087285533972e-16
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-6.270579197834706584e-12
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+2.601951087285533972e-16
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.918855460925072753e+00
+-1.499182616537808377e+00
+0.000000000000000000e+00
+1.918855460969481674e+00
+1.499182616537808377e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+1.918855460925072753e+00
+-1.499182616537808377e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.918855460969481674e+00
+1.499182616537808377e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+2.601951087285535451e-16
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-5.829522794043029421e-01
+-1.156442681460312372e+00
+0.000000000000000000e+00
+5.829522794043029421e-01
+1.156442681438107911e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-5.829522794043029421e-01
+1.156442681460312372e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+5.829522794043029421e-01
+-1.156442681438107911e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+6.890920502529956714e-12
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-6.270579197834704161e-12
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.156442681438107911e+00
+-6.298587673647658945e-03
+0.000000000000000000e+00
+1.156442681438107911e+00
+6.298587651443198453e-03
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+1.156442681438107911e+00
+-6.298587673647658945e-03
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.156442681438107911e+00
+6.298587651443198453e-03
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+2.601951087285535451e-16
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+8.972185350599605203e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+-8.972185350599607423e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+8.972185350599605203e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-8.972185350599607423e-01
+-1.197937775310187192e+00
+-9.359368253303300689e-01
+0.000000000000000000e+00
+1.197937775332391652e+00
+9.359368253747389899e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.918855460925072753e+00
+-1.499182616493399456e+00
+0.000000000000000000e+00
+1.918855460969481674e+00
+1.499182616448990535e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-5.829522793820984816e-01
+-1.156442681482516832e+00
+0.000000000000000000e+00
+5.829522794043029421e-01
+1.156442681482516832e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.156442681460312372e+00
+-6.298587673647658945e-03
+0.000000000000000000e+00
+1.156442681460312372e+00
+6.298587629238737960e-03
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+8.972185350599604092e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+-8.972185350599607423e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.297366463592908303e+00
+-1.013619467027382370e+00
+0.000000000000000000e+00
+1.297366463548499382e+00
+1.013619467027382370e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-2.941662268107769762e-01
+0.000000000000000000e+00
+7.465604617840070345e-01
+1.470831134220418335e-01
+0.000000000000000000e+00
+-7.465604617840070345e-01
+1.470831134109396032e-01
+0.000000000000000000e+00
+1.197937775310187192e+00
+-9.359368253303300689e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.197937775332391652e+00
+9.359368253747389899e-01
+0.000000000000000000e+00
+1.918855460925072753e+00
+-1.499182616493399456e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.918855460969481674e+00
+1.499182616448990535e+00
+0.000000000000000000e+00
+-5.829522793820984816e-01
+1.156442681482516832e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+5.829522794043029421e-01
+-1.156442681482516832e+00
+0.000000000000000000e+00
+1.156442681460312372e+00
+-6.298587673647658945e-03
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.156442681460312372e+00
+6.298587629238737960e-03
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+8.972185350599604092e-01
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-8.972185350599607423e-01
+0.000000000000000000e+00
+-2.941662268107769762e-01
+0.000000000000000000e+00
+7.465604617840070345e-01
+1.470831134109396032e-01
+0.000000000000000000e+00
+-7.465604617840070345e-01
+1.470831134220418335e-01
+0.000000000000000000e+00
+1.297366463592908303e+00
+-1.013619467027382370e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+0.000000000000000000e+00
+-1.297366463548499382e+00
+1.013619467027382370e+00
+0.000000000000000000e+00
diff --git a/tests/test_cohp.py b/tests/test_cohp.py
new file mode 100644
index 0000000..502c06f
--- /dev/null
+++ b/tests/test_cohp.py
@@ -0,0 +1,40 @@
+import unittest
+from pyqint import Molecule, HF, COHP
+import numpy as np
+
+class TestCOHP(unittest.TestCase):
+
+    def testCO(self):
+        """
+        Test construction of localized orbitals using Foster-Boys procedure
+        for the CO molecule
+        """
+        d = 1.145414
+        mol = Molecule()
+        mol.add_atom('C', 0.0, 0.0, -d/2, unit='angstrom')
+        mol.add_atom('O', 0.0, 0.0,  d/2, unit='angstrom')
+
+        res = HF().rhf(mol, 'sto3g')
+        cohp = COHP(res).run(res['orbc'], 0, 1)
+
+        cohp_ref = np.array([
+            0.0399,
+            0.0104,
+           -0.4365,
+            0.2051,
+           -0.2918,
+           -0.2918,
+            0.1098,
+            0.5029,
+            0.5029,
+            6.4827
+        ])
+
+        # note that Foster-Boys optimization is somewhat random and thus
+        # we use relatively loose testing criteria
+        np.testing.assert_almost_equal(cohp,
+                                       cohp_ref,
+                                       decimal=4)
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/test_derivatives_openmp.py b/tests/test_derivatives_openmp.py
new file mode 100644
index 0000000..042a087
--- /dev/null
+++ b/tests/test_derivatives_openmp.py
@@ -0,0 +1,61 @@
+import unittest
+from pyqint import PyQInt, Molecule
+import numpy as np
+
+class TestDerivativesOpenMP(unittest.TestCase):
+
+    def test_hartree_fock_forces_h2o(self):
+        """
+        Test Hartree-Fock calculation on water using STO-3G basis set
+        """
+        mol = Molecule()
+        mol.add_atom('O', 0.0, 0.1, 0.0)
+        mol.add_atom('H', 0.7570, 0.5860, 0.0)
+        mol.add_atom('H', -0.7570, 0.5860, 0.0)
+
+        # calculate forces using analytical derivatives
+        cgfs, nuclei = mol.build_basis('sto3g')
+        N = len(cgfs)
+
+        # build containers
+        S = np.zeros((3,3,N,N))
+        T = np.zeros_like(S)
+        V = np.zeros_like(S)
+
+        # build integrator object
+        integrator = PyQInt()
+
+        ntei = integrator.teindex(N-1,N-1,N-1,N-1)+1 # calculate number of two-electron integrals
+        teints = np.zeros((3,3,ntei))
+
+        for n,deriv_nucleus in enumerate(nuclei):
+            for d in range(3):
+                for i,cgf1 in enumerate(cgfs):
+                    for j,cgf2 in enumerate(cgfs):
+
+                        # derivative of overlap matrix
+                        S[n,d,i,j] += integrator.overlap_deriv(cgf1, cgf2, deriv_nucleus[0], d)
+
+                        # derivative of kinetic matrix
+                        T[n,d,i,j] += integrator.kinetic_deriv(cgf1, cgf2, deriv_nucleus[0], d)
+
+                        # derivative nuclear electron attraction
+                        for nucleus in nuclei:
+                            V[n,d,i,j] += integrator.nuclear_deriv(cgf1, cgf2, nucleus[0], nucleus[1], deriv_nucleus[0], d)
+
+                        ij = i*(i+1)//2+j;
+
+                        for k,cgf3 in enumerate(cgfs):
+                            for l,cgf4 in enumerate(cgfs):
+                                kl = k * (k+1)//2+l;
+                                if ij <= kl:
+                                    idx = integrator.teindex(i,j,k,l)
+                                    teints[n,d,idx] = integrator.repulsion_deriv(cgf1, cgf2, cgf3, cgf4, deriv_nucleus[0], d)
+
+
+        S2, T2, V2, teints2 = integrator.build_geometric_derivatives_openmp(cgfs, nuclei)
+
+        np.testing.assert_almost_equal(S, S2)
+        np.testing.assert_almost_equal(T, T2)
+        np.testing.assert_almost_equal(V, V2)
+        np.testing.assert_almost_equal(teints, teints2)
diff --git a/tests/test_dipole.py b/tests/test_dipole.py
index a79b80f..41b3585 100644
--- a/tests/test_dipole.py
+++ b/tests/test_dipole.py
@@ -35,7 +35,7 @@ def test_cgf_dipole(self):
             exact = np.load(os.path.join(os.path.dirname(__file__),
                                          'results',
                                          'h2o_dipole_%s.npy' % cart))
-            np.testing.assert_almost_equal(D[:,:,i], exact, decimal=8)
+            np.testing.assert_almost_equal(D[:,:,i], exact, decimal=4)
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_foster_boys.py b/tests/test_foster_boys.py
index baba529..0aa13ee 100644
--- a/tests/test_foster_boys.py
+++ b/tests/test_foster_boys.py
@@ -15,7 +15,9 @@ def testCO(self):
         mol.add_atom('O', 0.0, 0.0,  d/2, unit='angstrom')
 
         res = HF().rhf(mol, 'sto3g')
-        res_fb = FosterBoys(res).run()
+
+        # note that a seed is given here for reproducibility purposes
+        res_fb = FosterBoys(res, seed=0).run(nr_runners=5)
 
         orbe_ref = np.array([
             -20.30750217,
@@ -30,9 +32,11 @@ def testCO(self):
             res['orbe'][9]
         ])
 
+        # note that Foster-Boys optimization is somewhat random and thus
+        # we use relatively loose testing criteria
         np.testing.assert_almost_equal(res_fb['orbe'],
                                        orbe_ref,
-                                       decimal=3)
+                                       decimal=2)
 
     def testCH4(self):
         """
@@ -48,7 +52,9 @@ def testCH4(self):
         mol.add_atom('H', dist, -dist, -dist, unit='angstrom')
 
         res = HF().rhf(mol, 'sto3g')
-        res_fb = FosterBoys(res).run()
+
+        # note that a seed is given here for reproducibility purposes
+        res_fb = FosterBoys(res, seed=0).run(nr_runners=5)
 
         orbe_ref = np.array([
             -11.050113,
@@ -63,14 +69,16 @@ def testCH4(self):
         ])
 
         # assert orbital energies
-        np.testing.assert_almost_equal(res_fb['orbe'], orbe_ref, decimal=3)
+        np.testing.assert_almost_equal(res_fb['orbe'], orbe_ref, decimal=2)
 
         # specifically test for quadruple degenerate orbital
         for i in range(0,4):
             for j in range(i+1,4):
+                # note that Foster-Boys optimization is somewhat random and thus
+                # we use relatively loose testing criteria
                 np.testing.assert_almost_equal(res_fb['orbe'][i+1],
                                                res_fb['orbe'][j+1],
-                                               decimal=4)
+                                               decimal=2)
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_geometry_optimization.py b/tests/test_geometry_optimization.py
new file mode 100644
index 0000000..865ce2d
--- /dev/null
+++ b/tests/test_geometry_optimization.py
@@ -0,0 +1,106 @@
+import unittest
+from pyqint import Molecule,GeometryOptimization
+import numpy as np
+
+class TestGeometryOptimization(unittest.TestCase):
+    """
+    Test the calculation of analytical gradients by comparing these
+    with the finite difference results
+    """
+
+    def test_optimization_h2_sto3g(self):
+        """
+        Optimize dihydrogen molecule and assess that the energy corresponds to
+        an optimized structure
+        """
+        # create new H2 molecule with perturbed geometry
+        mol = Molecule()
+        mol.add_atom('H', 0.9, 0.0, 0.0)
+        mol.add_atom('H', -0.9, 0.0, 0.0)
+
+        res = GeometryOptimization(verbose=False).run(mol, 'sto3g')
+        np.testing.assert_almost_equal(res['opt'].fun, -1.117506)
+
+        self.assertEqual(len(res['energies']), len(res['forces']))
+        self.assertEqual(len(res['energies']), len(res['coordinates']))
+        self.assertEqual(res['coordinates'][0].shape, (len(mol.atoms),3))
+
+    def test_optimization_h2_p321(self):
+        """
+        Optimize dihydrogen molecule and assess that the energy corresponds to
+        an optimized structure
+        """
+        # create new H2 molecule with perturbed geometry
+        mol = Molecule()
+        mol.add_atom('H', 0.9, 0.0, 0.0)
+        mol.add_atom('H', -0.9, 0.0, 0.0)
+
+        res = GeometryOptimization(verbose=False).run(mol, 'p321')
+        np.testing.assert_almost_equal(res['opt'].fun, -1.1229598312004625, decimal=4)
+
+        self.assertEqual(len(res['energies']), len(res['forces']))
+        self.assertEqual(len(res['energies']), len(res['coordinates']))
+        self.assertEqual(res['coordinates'][0].shape, (len(mol.atoms),3))
+
+    def test_optimization_ch4(self):
+        """
+        Optimize methane molecule and assess that the energy corresponds to
+        an optimized structure
+        """
+        # create new CH4 molecule with perturbed geometry
+        mol = Molecule()
+        dist = 2.0/2
+        mol.add_atom('C', 0.1, 0.0, 0.1, unit='angstrom')
+        mol.add_atom('H', dist, dist, dist, unit='angstrom')
+        mol.add_atom('H', -dist, -dist, dist, unit='angstrom')
+        mol.add_atom('H', -dist, dist, -dist, unit='angstrom')
+        mol.add_atom('H', dist, -dist, -dist, unit='angstrom')
+
+        res = GeometryOptimization(verbose=False).run(mol, 'sto3g')
+        np.testing.assert_almost_equal(res['opt'].fun, -39.7268635)
+
+        self.assertEqual(len(res['energies']), len(res['forces']))
+        self.assertEqual(len(res['energies']), len(res['coordinates']))
+        self.assertEqual(res['coordinates'][0].shape, (len(mol.atoms),3))
+
+    def test_optimization_h2o(self):
+        """
+        Optimize the water molecule and assess that the energy corresponds to
+        an optimized structure
+        """
+        # create new H2O molecule with perturbed geometry
+        mol = Molecule()
+        mol.add_atom('O', 0.0,  0.0, -0.2271310707, unit='angstrom')
+        mol.add_atom('H', 0.0,  -0.8580158822, 0.5085242828, unit='angstrom')
+        mol.add_atom('H', 0.0,  0.8580158822, 0.5085242828, unit='angstrom')
+
+        res = GeometryOptimization(verbose=False).run(mol, 'sto3g')
+        np.testing.assert_almost_equal(res['opt'].fun, -74.96590013836217)
+
+        self.assertEqual(len(res['energies']), len(res['forces']))
+        self.assertEqual(len(res['energies']), len(res['coordinates']))
+        self.assertEqual(res['coordinates'][0].shape, (len(mol.atoms),3))
+
+    def test_optimization_c2h4(self):
+        """
+        Optimize ethylene molecule and assess that the energy corresponds to
+        an optimized structure
+        """
+        # create new C2H4 molecule with perturbed geometry
+        mol = Molecule()
+        mol.add_atom('C', -0.5530176758,  0.0000000000 ,0.0000000000, unit='angstrom')
+        mol.add_atom('C',  0.4530176758,  0.0000000000 ,0.0000000000, unit='angstrom')
+        mol.add_atom('H', -1.3288875372, -0.9556191261 ,0.1000000000, unit='angstrom')
+        mol.add_atom('H', -1.1288875372,  0.9356191261 ,0.0000000000, unit='angstrom')
+        mol.add_atom('H',  1.3288875372,  0.9556191261 ,0.0000000000, unit='angstrom')
+        mol.add_atom('H',  1.1288875372, -0.9156191261 ,0.1000000000, unit='angstrom')
+
+        res = GeometryOptimization(verbose=False).run(mol, 'sto3g')
+        np.testing.assert_almost_equal(res['opt'].fun, -77.0739544)
+
+        self.assertEqual(len(res['energies']), len(res['forces']))
+        self.assertEqual(len(res['energies']), len(res['coordinates']))
+        self.assertEqual(res['coordinates'][0].shape, (len(mol.atoms),3))
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/test_hf.py b/tests/test_hf.py
index 43b72e5..e4e457c 100644
--- a/tests/test_hf.py
+++ b/tests/test_hf.py
@@ -7,7 +7,7 @@
 
 class TestHF(unittest.TestCase):
 
-    def testHartreeFockWater(self):
+    def test_hartree_fock_h2o(self):
         """
         Test Hartree-Fock calculation on water using STO-3G basis set
         """
@@ -16,7 +16,7 @@ def testHartreeFockWater(self):
         mol.add_atom('H', 0.7570, 0.5860, 0.0)
         mol.add_atom('H', -0.7570, 0.5860, 0.0)
 
-        results = perform_hf(mol)
+        results = HF().rhf(mol, 'sto3g')
 
         # check that energy matches
         np.testing.assert_almost_equal(results['energy'], -73.21447132, 4)
@@ -25,7 +25,7 @@ def testHartreeFockWater(self):
         #self.assertTrue(results['time_stats']['integral_evaluation'] > 0)
         #self.assertTrue(results['time_stats']['self_convergence'] > 0)
 
-    def testHartreeFockCH4Properties(self):
+    def test_hartree_fock_ch4(self):
         """
         Test Hartree-Fock calculation on water using STO-3G basis set
         """
@@ -37,7 +37,7 @@ def testHartreeFockCH4Properties(self):
         mol.add_atom('H', -dist, dist, -dist, unit='angstrom')
         mol.add_atom('H', dist, -dist, -dist, unit='angstrom')
 
-        results = perform_hf(mol)
+        results = HF().rhf(mol, 'sto3g')
 
         # check that orbital energies are correctly approximated
         ans = np.array([-11.0707,
@@ -54,9 +54,55 @@ def testHartreeFockCH4Properties(self):
         en = -39.35007843284954
         np.testing.assert_almost_equal(results['energies'][-1], en, 4)
 
-def perform_hf(mol):
-    results = HF().rhf(mol, 'sto3g')
-    return results
+    def test_hartree_fock_restart(self):
+        """
+        Test Hartree-Fock calculation on water using STO-3G basis set
+        """
+        mol = Molecule()
+        dist = 1.78/2
+        mol.add_atom('C', 0.0, 0.0, 0.0, unit='angstrom')
+        mol.add_atom('H', dist, dist, dist, unit='angstrom')
+        mol.add_atom('H', -dist, -dist, dist, unit='angstrom')
+        mol.add_atom('H', -dist, dist, -dist, unit='angstrom')
+        mol.add_atom('H', dist, -dist, -dist, unit='angstrom')
+        results1 = HF().rhf(mol, 'sto3g')
+
+        # check that orbital energies are correctly approximated
+        ans = np.array([-11.0707,
+                         -0.7392,
+                         -0.3752,
+                         -0.3752,
+                         -0.3752,
+                          0.2865,
+                          0.4092,
+                          0.4092,
+                          0.4092])
+        np.testing.assert_almost_equal(results1['orbe'], ans, 4)
+
+        en = -39.35007843284954
+        np.testing.assert_almost_equal(results1['energies'][-1], en, 4)
+
+        # create new CH4 molecule with slight adjustment in geometry and
+        # seed the calculation with the previous converged result
+        mol = Molecule()
+        dist = 1.78/2
+        mol.add_atom('C', 0.0, 0.0, 0.1, unit='angstrom')
+        mol.add_atom('H', dist, dist, dist, unit='angstrom')
+        mol.add_atom('H', -dist, -dist, dist, unit='angstrom')
+        mol.add_atom('H', -dist, dist, -dist, unit='angstrom')
+        mol.add_atom('H', dist, -dist, -dist, unit='angstrom')
+
+        # perform HF calculation with coefficient matrix from previous
+        # calculation to speed up the convergence
+        results2 = HF().rhf(mol, 'sto3g', orbc_init=results1['orbc'])
+
+        # assess that the energy of the perturbed result is different
+        # (and also higher)
+        en = -39.34538546003782
+        np.testing.assert_almost_equal(results2['energies'][-1], en, 3)
+
+        # check that the convergence is quicker
+        self.assertTrue(len(results1['energies']) > len(results2['energies']))
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_hf_deriv.py b/tests/test_hf_deriv.py
index e9e1499..0dec4dd 100644
--- a/tests/test_hf_deriv.py
+++ b/tests/test_hf_deriv.py
@@ -8,58 +8,96 @@
 
 class TestHFDeriv(unittest.TestCase):
 
-    @nottest
-    def testHartreeFockForcesCore(self):
+    def test_hartree_fock_forces_h2(self):
+        """
+        Test Hartree-Fock calculation on H2 using STO-3G basis set
+        """
+        for i in range(0,5):
+            mol = Molecule()
+            mol.add_atom('H', 0.0000, 0.0000, 0.30 + i * 0.50, unit='angstrom')
+            mol.add_atom('H', 0.0000, 0.0000, -0.30 - i * 0.50, unit='angstrom')
+
+            # calculate forces using analytical derivatives
+            solver = HF()
+            res = solver.rhf(mol, 'sto3g', calc_forces=True)
+
+            # calculate forces using finite difference
+            forces = calculate_forces_finite_difference(mol)
+
+            np.testing.assert_almost_equal(res['forces'], forces, decimal=4)
+
+    def test_hartree_fock_forces_h2o(self):
         """
         Test Hartree-Fock calculation on water using STO-3G basis set
         """
         mol = Molecule()
-        mol.add_atom('O', 0.0, 0.0, 0.0)
+        mol.add_atom('O', 0.0, 0.1, 0.0)
         mol.add_atom('H', 0.7570, 0.5860, 0.0)
         mol.add_atom('H', -0.7570, 0.5860, 0.0)
-        basis = 'sto3g'
 
         # calculate forces using analytical derivatives
         solver = HF()
-        res = solver.rhf(mol, basis, calc_forces=False)
+        res = solver.rhf(mol, 'sto3g', calc_forces=True)
+
+        # calculate forces using finite difference
+        forces = calculate_forces_finite_difference(mol)
 
-        # calculate forces via routines
-        forces = np.zeros((len(mol.atoms), 3))
-        for i in range(0, len(mol.atoms)): # loop over nuclei
-            for j in range(0, 3): # loop over directions
-                forces[i,j] = np.sum(res['density'] * solver.rhf_force_core(mol, basis, i, j))
+        np.testing.assert_almost_equal(res['forces'], forces, decimal=4)
+
+    def test_hartree_fock_forces_ch4(self):
+        """
+        Test Hartree-Fock calculation on CH4 using STO-3G basis set
+
+        Note that the CH4 molecule is slightly perturbed
+        """
+        mol = Molecule()
+        mol.add_atom('C', 0.0,  0.1, 0.0, unit='angstrom')
+        mol.add_atom('H', 0.1,  0.15, 1.0830098121, unit='angstrom')
+        mol.add_atom('H', 0.0,  -1.0210714424, -0.3610032723, unit='angstrom')
+        mol.add_atom('H', -0.8842738057,  0.5105357237, -0.3610032748, unit='angstrom')
+        mol.add_atom('H', 0.8842738117,  0.5105357203, -0.3610032651, unit='angstrom')
+
+        # calculate forces using analytical derivatives
+        solver = HF()
+        res = solver.rhf(mol, 'sto3g', calc_forces=True)
 
         # calculate forces using finite difference
-        forces_ans = calculate_forces_core_finite_difference(mol)
+        forces = calculate_forces_finite_difference(mol)
 
-        np.testing.assert_almost_equal(forces, forces_ans)
+        np.testing.assert_almost_equal(res['forces'], forces, decimal=4)
 
-    # def testHartreeFockForces(self):
-    #     """
-    #     Test Hartree-Fock calculation on water using STO-3G basis set
-    #     """
-    #     mol = Molecule()
-    #     mol.add_atom('O', 0.0, 0.0, 0.0)
-    #     mol.add_atom('H', 0.7570, 0.5860, 0.0)
-    #     mol.add_atom('H', -0.7570, 0.5860, 0.0)
+    def test_hartree_fock_forces_co2(self):
+        """
+        Test Hartree-Fock calculation on CH4 using STO-3G basis set
 
-    #     # calculate forces using analytical derivatives
-    #     solver = HF()
-    #     res = solver.rhf(mol, 'sto3g', calc_forces=True)
+        Note that the CO2 molecule is slightly perturbed
+        """
+        mol = Molecule()
+        mol.add_atom('C', 0.0,  0.0, 0.0, unit='angstrom')
+        mol.add_atom('O', 0.0,  0.0, 1.2879700928, unit='angstrom')
+        mol.add_atom('O', 0.0,  0.0, -1.2879700928, unit='angstrom')
 
-    #     # calculate forces using finite difference
-    #     forces = calculate_forces_finite_difference(mol)
+        # calculate forces using analytical derivatives
+        solver = HF()
+        res = solver.rhf(mol, 'sto3g', calc_forces=True)
 
-    #     np.testing.assert_almost_equal(res['forces'], forces)
+        # calculate forces using finite difference
+        forces = calculate_forces_finite_difference(mol)
+
+        np.testing.assert_almost_equal(res['forces'], forces, decimal=3)
 
 def perform_hf(mol):
     sol = HF().rhf(mol, 'sto3g')
     return sol
 
 def calculate_forces_finite_difference(mol):
-    forces = np.zeros((3,3))
+    """
+    Calculates the forces on each of the atoms using a finite difference
+    approach.
+    """
+    forces = np.zeros((len(mol.atoms),3))
 
-    sz = 0.0001
+    sz = 1e-4
 
     for i in range(0, len(mol.atoms)): # loop over nuclei
         for j in range(0, 3): # loop over directions
@@ -75,24 +113,5 @@ def calculate_forces_finite_difference(mol):
 
     return forces
 
-def calculate_forces_core_finite_difference(mol):
-    forces = np.zeros((3,3))
-
-    sz = 0.0001
-
-    for i in range(0, len(mol.atoms)): # loop over nuclei
-        for j in range(0, 3): # loop over directions
-            mol1 = deepcopy(mol)
-            mol1.atoms[i][1][j] -= sz / 2
-            mol2 = deepcopy(mol)
-            mol2.atoms[i][1][j] += sz / 2
-
-            energy1 = perform_hf(mol1)['ecore']
-            energy2 = perform_hf(mol2)['ecore']
-
-            forces[i,j] = (energy2 - energy1) / sz
-
-    return forces
-
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_molecule_builder.py b/tests/test_molecule_builder.py
new file mode 100644
index 0000000..1e9a917
--- /dev/null
+++ b/tests/test_molecule_builder.py
@@ -0,0 +1,26 @@
+import unittest
+from pyqint import MoleculeBuilder
+import numpy as np
+
+class TestMoleculeBuilder(unittest.TestCase):
+    
+    def test_construction(self):
+        """
+        Build a molecule from point group
+        """
+        mol = MoleculeBuilder().build_complex_td(1.0, 'C', 'H')
+
+    def test_load_molecule_from_file(self):
+        """
+        Build a molecule from file
+        """
+        mol = MoleculeBuilder().from_name('ch4')
+
+        np.testing.assert_almost_equal(mol.atoms[0][1], 
+                                       np.array([0.0,0.0,0.0], dtype=np.float64))
+        np.testing.assert_almost_equal(mol.atoms[1][1], 
+                                       np.array([0.6327670,0.6327670,0.6327670]) * 1.8897259886)
+        np.testing.assert_equal(mol.atoms[0][0], 'C')
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/test_nuclear_deriv.py b/tests/test_nuclear_deriv.py
index baf683c..07c74dc 100644
--- a/tests/test_nuclear_deriv.py
+++ b/tests/test_nuclear_deriv.py
@@ -7,350 +7,6 @@
 
 class TestNuclearDeriv(unittest.TestCase):
 
-    def testDerivBF_s1(self):
-        """
-        Test Derivative of the s-type basis functions in x-direction
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = [-0.5, 0.0, 0.0]
-        gto1 = gto(0.154329, [-0.50, 0.0, 0.0], 3.425251, 0, 0, 0)
-        gto2 = gto(0.154329, [ 0.50, 0.0, 0.0], 3.425251, 0, 0, 0)
-
-        # derivative towards bf coordinates
-        t1 = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 0)
-
-        # test the first derivative of the second hydrogen atom in the
-        # Cartesian direction
-        t2 = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 0)
-
-        # also calculate this integral using finite difference
-        fd_01 = calculate_fx_bf_finite_difference(gto1, nucleus)
-        fd_02 = calculate_fx_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(-2 * t1, fd_01, 4)
-        #np.testing.assert_almost_equal(-2 * t2, fd_02, 4)
-
-    def testDerivBF_s2(self):
-        """
-        Test Derivative of the s-type basis functions in y-direction
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = [-0.5, 0.0, 0.0]
-        gto1 = gto(0.154329, [-0.50, 0.0, 0.0], 3.425251, 0, 0, 0)
-        gto2 = gto(0.154329, [ 0.50, 0.0, 0.0], 3.425251, 0, 0, 0)
-
-        # derivative towards bf coordinates
-        t1 = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 1)
-
-        # test the first derivative of the second hydrogen atom in the
-        # Cartesian direction
-        t2 = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 1)
-
-        # also calculate this integral using finite difference
-        fd_01 = calculate_fy_bf_finite_difference(gto1, nucleus)
-        fd_02 = calculate_fy_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(-2 * t1, fd_01, 4)
-        np.testing.assert_almost_equal(-2 * t2, fd_02, 4)
-
-    def testDerivBF_px(self):
-        """
-        Test Derivative of the px-type basis functions in x-direction
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = [-0.5, 0.0, 0.0]
-        gto1 = gto(0.154329, [-0.50, 0.0, 0.0], 3.425251, 1, 0, 0)
-        gto2 = gto(0.154329, [ 0.50, 0.0, 0.0], 3.425251, 1, 0, 0)
-
-        # derivative towards bf coordinates
-        t1 = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 0)
-
-        # test the first derivative of the second hydrogen atom in the
-        # Cartesian direction
-        t2 = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 0)
-
-        # also calculate this integral using finite difference
-        fx_fd_01 = calculate_fx_bf_finite_difference(gto1, nucleus)
-        fx_fd_02 = calculate_fx_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(-2 * t1, fx_fd_01, 4)
-        #np.testing.assert_almost_equal(-2 * t2, fx_fd_02, 4)
-
-    def testDerivBF_pxy(self):
-        """
-        Test Derivative of the px-type basis functions in y-direction
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = [-0.5, 0.0, 0.0]
-        gto1 = gto(0.154329, [-0.50, 0.0, 0.0], 3.425251, 1, 0, 0)
-        gto2 = gto(0.154329, [ 0.50, 0.0, 0.0], 3.425251, 1, 0, 0)
-
-        # derivative towards bf coordinates
-        t1 = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 1)
-
-        # test the first derivative of the second hydrogen atom in the
-        # Cartesian direction
-        t2 = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 1)
-
-        # also calculate this integral using finite difference
-        fx_fd_01 = calculate_fy_bf_finite_difference(gto1, nucleus)
-        fx_fd_02 = calculate_fy_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(-2 * t1, fx_fd_01, 4)
-        np.testing.assert_almost_equal(-2 * t2, fx_fd_02, 4)
-
-    def testDerivBF_pyy(self):
-        """
-        Test Derivative of the py-type basis functions in y-direction
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = [-0.5, 0.0, 0.0]
-        gto1 = gto(0.154329, [-0.50, 0.0, 0.0], 3.425251, 0, 1, 0)
-        gto2 = gto(0.154329, [ 0.50, 0.0, 0.0], 3.425251, 0, 1, 0)
-
-        # derivative towards bf coordinates
-        t1 = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 1)
-
-        # test the first derivative of the second hydrogen atom in the
-        # Cartesian direction
-        t2 = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 1)
-
-        # also calculate this integral using finite difference
-        fx_fd_01 = calculate_fy_bf_finite_difference(gto1, nucleus)
-        fx_fd_02 = calculate_fy_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(-2 * t1, fx_fd_01, 4)
-        np.testing.assert_almost_equal(-2 * t2, fx_fd_02, 4)
-
-    def testDerivBF_d(self):
-        """
-        Test Derivative of the dx2-type basis functions in x-direction
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = [-0.5, 0.0, 0.0]
-        gto1 = gto(0.154329, [-0.50, 0.0, 0.0], 3.425251, 2, 0, 0)
-        gto2 = gto(0.154329, [ 0.50, 0.0, 0.0], 3.425251, 2, 0, 0)
-
-        # derivative towards bf coordinates
-        t1 = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 0)
-
-        # test the first derivative of the second hydrogen atom in the
-        # Cartesian direction
-        t2 = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 0)
-
-        # also calculate this integral using finite difference
-        fx_fd_01 = calculate_fx_bf_finite_difference(gto1, nucleus)
-        fx_fd_02 = calculate_fx_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(-2 * t1, fx_fd_01, 4)
-        #np.testing.assert_almost_equal(-2 * t2, fx_fd_02, 4)
-
-    def testDerivOpt_s1(self):
-        """
-        Test Derivative of the nuclear attraction operator in x-direction
-        for s-type orbitals
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = np.array([-0.5, 0.0, 0.0])
-        gto1 = gto(0.154329,  nucleus, 3.425251, 0, 0, 0)
-        gto2 = gto(0.154329, -nucleus, 3.425251, 0, 0, 0)
-
-        t1a = integrator.nuclear_gto_deriv_op(gto1, gto1, nucleus, 0)
-        t1b = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 0)
-
-        t2a = integrator.nuclear_gto_deriv_op(gto2, gto2, nucleus, 0)
-        t2b = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 0)
-
-        # also calculate this integral using finite difference
-        fd_01 = calculate_fx_op_finite_difference(gto1, nucleus)
-        fd_02a = calculate_fx_op_finite_difference(gto2, nucleus)
-        fd_02b = calculate_fx_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(t1a, 0.0, 4)
-        np.testing.assert_almost_equal(t1b, 0.0, 4)
-
-        #np.testing.assert_almost_equal(-2.0 * t2b, fd_02b, 4)
-        np.testing.assert_almost_equal(t2a, fd_02a, 4)
-
-    def testDerivOpt_s2(self):
-        """
-        Test Derivative of the nuclear attraction operator in y-direction
-        for s-type orbitals
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = np.array([-0.5, 0.0, 0.0])
-        gto1 = gto(0.154329,  nucleus, 3.425251, 0, 0, 0)
-        gto2 = gto(0.154329, -nucleus, 3.425251, 0, 0, 0)
-
-        t1a = integrator.nuclear_gto_deriv_op(gto1, gto1, nucleus, 1)
-        t1b = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 1)
-
-        t2a = integrator.nuclear_gto_deriv_op(gto2, gto2, nucleus, 1)
-        t2b = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 1)
-
-        # also calculate this integral using finite difference
-        fd_01 = calculate_fy_op_finite_difference(gto1, nucleus)
-        fd_02a = calculate_fy_op_finite_difference(gto2, nucleus)
-        fd_02b = calculate_fy_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(t1a, 0.0, 4)
-        np.testing.assert_almost_equal(t1b, 0.0, 4)
-
-        np.testing.assert_almost_equal(-2.0 * t2b, fd_02b, 4)
-        np.testing.assert_almost_equal(t2a, fd_02a, 4)
-
-    def testDerivOpt_p1(self):
-        """
-        Test Derivative of the nuclear attraction operator in x-direction
-        for p-type orbitals
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = np.array([-0.5, 0.0, 0.0])
-        gto1 = gto(0.154329,  nucleus, 3.425251, 1, 0, 0)
-        gto2 = gto(0.154329, -nucleus, 3.425251, 1, 0, 0)
-
-        t1a = integrator.nuclear_gto_deriv_op(gto1, gto1, nucleus, 0)
-        t1b = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 0)
-
-        t2a = integrator.nuclear_gto_deriv_op(gto2, gto2, nucleus, 0)
-        t2b = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 0)
-
-        # also calculate this integral using finite difference
-        fd_01 = calculate_fx_op_finite_difference(gto1, nucleus)
-        fd_02a = calculate_fx_op_finite_difference(gto2, nucleus)
-        fd_02b = calculate_fx_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(t1a, 0.0, 4)
-        np.testing.assert_almost_equal(t1b, 0.0, 4)
-
-        #np.testing.assert_almost_equal(-2.0 * t2b, fd_02b, 4)
-        np.testing.assert_almost_equal(t2a, fd_02a, 4)
-
-    def testDerivOpt_d1(self):
-        """
-        Test Derivative of the nuclear attraction operator in x-direction
-        for d-type orbitals
-        """
-
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        nucleus = np.array([-0.5, 0.0, 0.0])
-        gto1 = gto(0.154329,  nucleus, 3.425251, 2, 0, 0)
-        gto2 = gto(0.154329, -nucleus, 3.425251, 2, 0, 0)
-
-        t1a = integrator.nuclear_gto_deriv_op(gto1, gto1, nucleus, 0)
-        t1b = integrator.nuclear_gto_deriv_bf(gto1, gto1, nucleus, 0)
-
-        t2a = integrator.nuclear_gto_deriv_op(gto2, gto2, nucleus, 0)
-        t2b = integrator.nuclear_gto_deriv_bf(gto2, gto2, nucleus, 0)
-
-        # also calculate this integral using finite difference
-        fd_01 = calculate_fx_op_finite_difference(gto1, nucleus)
-        fd_02a = calculate_fx_op_finite_difference(gto2, nucleus)
-        fd_02b = calculate_fx_bf_finite_difference(gto2, nucleus)
-
-        # testing
-        np.testing.assert_almost_equal(t1a, 0.0, 4)
-        np.testing.assert_almost_equal(t1b, 0.0, 4)
-
-        #np.testing.assert_almost_equal(-2.0 * t2b, fd_02b, 4)
-        np.testing.assert_almost_equal(t2a, fd_02a, 4)
-
-    def testDerivH2(self):
-        """
-        Test Derivatives of dihydrogen
-        """
-
-        # build integrator object
-        integrator = PyQInt()
-
-        # build hydrogen molecule
-        mol = Molecule('H2')
-        mol.add_atom('H', -0.5, 0.0, 0.0)
-        mol.add_atom('H',  0.5, 0.0, 0.0)
-        cgfs, nuclei = mol.build_basis('sto3g')
-
-        fx1 = integrator.nuclear_deriv(cgfs[0], cgfs[0], nuclei[0][0], nuclei[0][1], nuclei[0][0], 0)
-        fx2 = integrator.nuclear_deriv(cgfs[1], cgfs[1], nuclei[0][0], nuclei[0][1], nuclei[0][0], 0)
-
-        ans1 = calculate_force_finite_difference(cgfs[0], cgfs[0], nuclei[0][0], nuclei[0][1], 0)
-        ans2 = calculate_force_finite_difference(cgfs[1], cgfs[1], nuclei[0][0], nuclei[0][1], 0)
-
-        np.testing.assert_almost_equal(fx1, ans1, 4)
-        np.testing.assert_almost_equal(fx2, ans2, 4)
-
-        # assert that the nuclear gradient in the x direction is zero because the
-        # basis functions spawn from this atom
-        self.assertTrue(np.abs(fx1) < 1e-8)
-
-        # assert that the nuclear gradient has a meaningful number
-        self.assertTrue(np.abs(fx2) > 1e-1)
-
-        fx3 = integrator.nuclear_deriv(cgfs[0], cgfs[0], nuclei[1][0], nuclei[1][1], nuclei[1][0], 0)
-        fx4 = integrator.nuclear_deriv(cgfs[1], cgfs[1], nuclei[1][0], nuclei[1][1], nuclei[1][0], 0)
-
-        ans3 = calculate_force_finite_difference(cgfs[0], cgfs[0], nuclei[1][0], nuclei[1][1], 0)
-        ans4 = calculate_force_finite_difference(cgfs[1], cgfs[1], nuclei[1][0], nuclei[1][1], 0)
-
-        np.testing.assert_almost_equal(fx3, ans3, 4)
-        np.testing.assert_almost_equal(fx4, ans4, 4)
-
-        fy1 = integrator.nuclear_deriv(cgfs[0], cgfs[0], nuclei[0][0], nuclei[0][1], nuclei[0][0], 1)
-        fy2 = integrator.nuclear_deriv(cgfs[1], cgfs[1], nuclei[0][0], nuclei[0][1], nuclei[0][0], 1)
-
-        ans5 = calculate_force_finite_difference(cgfs[0], cgfs[0], nuclei[0][0], nuclei[0][1], 1)
-        ans6 = calculate_force_finite_difference(cgfs[1], cgfs[1], nuclei[0][0], nuclei[0][1], 1)
-
-        np.testing.assert_almost_equal(fy1, ans5, 4)
-        np.testing.assert_almost_equal(fy2, ans6, 4)
-
     def testNuclearRepulsionDerivatives(self):
         """
         Test Hartree-Fock calculation on water using STO-3G basis set
@@ -371,86 +27,62 @@ def testNuclearRepulsionDerivatives(self):
 
         np.testing.assert_almost_equal(forces, forces_fd, 4)
 
-def calculate_force_finite_difference(cgf1, cgf2, nucleus, charge, coord):
-    # build integrator object
-    integrator = PyQInt()
-
-    # distance
-    diff = 0.00001
-
-    # build hydrogen molecule
-    nucleus[coord] -= diff / 2.0
-    left = integrator.nuclear(cgf1, cgf2, nucleus, charge)
-    nucleus[coord] += diff
-    right = integrator.nuclear(cgf1, cgf2, nucleus, charge)
-
-    return (right - left) / diff
-
-def calculate_fx_bf_finite_difference(_gto, nucleus):
-    """
-    Calculate the basis function derivative in the x-direction using
-    a finite difference method
-    """
-    # build integrator object
-    integrator = PyQInt()
-
-    diff = 0.01
-    gto1 = gto(0.154329, [ _gto.p[0] - diff / 2, _gto.p[1], _gto.p[2]], 3.425251, _gto.l, _gto.m, _gto.n)
-    gto2 = gto(0.154329, [ _gto.p[0] + diff / 2, _gto.p[1], _gto.p[2]], 3.425251, _gto.l, _gto.m, _gto.n)
-
-    left = integrator.nuclear_gto(gto1, gto1, nucleus)
-    right = integrator.nuclear_gto(gto2, gto2, nucleus)
+    def test_derivatives_h2o_fulltest(self):
+        """
+        Test Derivatives of water
+        """
+        # build integrator object
+        integrator = PyQInt()
 
-    return (right - left) / diff
+        # build hydrogen molecule
+        mol = Molecule('H2O')
+        mol.add_atom('O', 0.00000, -0.07579, 0.00000, unit='angstrom')
+        mol.add_atom('H', 0.86681, 0.60144, 0.00000, unit='angstrom')
+        mol.add_atom('H',  -0.86681, 0.60144, 0.00000, unit='angstrom')
+        cgfs, nuclei = mol.build_basis('sto3g')
 
-def calculate_fy_bf_finite_difference(_gto, nucleus):
-    """
-    Calculate the basis function derivative in the y-direction using
-    a finite difference method
-    """
+        # get position and charge
+        probe_atom = 0
+        O = nuclei[probe_atom][0]
+        Ochg = nuclei[probe_atom][1]
+
+        # load results from file
+        fname = os.path.join(os.path.dirname(__file__), 'results', 'nuclear_deriv_h2o.txt')
+        vals = np.loadtxt(fname).reshape((len(cgfs), len(cgfs), 3, 3))
+        for i in range(0, len(cgfs)): # loop over cgfs
+            for j in range(0, len(cgfs)): # loop over cgfs
+                for k in range(0,3):  # loop over nuclei
+                    for l in range(0,3):  # loop over directions
+                        force = integrator.nuclear_deriv(cgfs[i], cgfs[j], O, Ochg, nuclei[k][0], l)
+                        np.testing.assert_almost_equal(force, vals[i,j,k,l], 4)
+
+def calculate_force_finite_difference(cgf_id1, cgf_id2, nuc_id, coord, probe_atom):
     # build integrator object
     integrator = PyQInt()
 
-    diff = 0.01
-    gto1 = gto(0.154329, [ _gto.p[0], _gto.p[1] - diff / 2, _gto.p[2]], 3.425251, _gto.l, _gto.m, _gto.n)
-    gto2 = gto(0.154329, [ _gto.p[0], _gto.p[1] + diff / 2, _gto.p[2]], 3.425251, _gto.l, _gto.m, _gto.n)
-
-    left = integrator.nuclear_gto(gto1, gto1, nucleus)
-    right = integrator.nuclear_gto(gto2, gto2, nucleus)
-
-    return (right - left) / diff
-
-def calculate_fx_op_finite_difference(_gto, nucleus):
-    """
-    Calculate the operator derivative in the x-direction using
-    a finite difference method
-    """
-    # build integrator object
-    integrator = PyQInt()
+    # distance
+    diff = 0.00001
 
-    diff = 0.000001
-    nucleus[0] -= diff / 2.0
-    left = integrator.nuclear_gto(_gto, _gto, nucleus)
-    nucleus[0] += diff
-    right = integrator.nuclear_gto(_gto, _gto, nucleus)
+    vals = np.zeros(2)
+    for i,v in enumerate([-1,1]):
+        # build hydrogen molecule
+        mol = Molecule('H2O')
+        mol.add_atom('O', 0.00000, -0.07579, 0.00000, unit='angstrom')
+        mol.add_atom('H', 0.86681, 0.60144, 0.00000, unit='angstrom')
+        mol.add_atom('H',  -0.86681, 0.60144, 0.00000, unit='angstrom')
 
-    return (right - left) / diff
+        # adjust molecule
+        mol.atoms[nuc_id][1][coord] += v * diff / 2
 
-def calculate_fy_op_finite_difference(_gto, nucleus):
-    """
-    Calculate the operator derivative in the y-direction using
-    a finite difference method
-    """
-    # build integrator object
-    integrator = PyQInt()
+        # build basis
+        cgfs, nuclei = mol.build_basis('sto3g')
 
-    diff = 0.000001
-    nucleus[1] -= diff / 2.0
-    left = integrator.nuclear_gto(_gto, _gto, nucleus)
-    nucleus[1] += diff
-    right = integrator.nuclear_gto(_gto, _gto, nucleus)
+        # calculate values
+        O = nuclei[probe_atom][0]
+        Ochg = nuclei[probe_atom][1]
+        vals[i] = integrator.nuclear(cgfs[cgf_id1], cgfs[cgf_id2], O, Ochg)
 
-    return (right - left) / diff
+    return (vals[1] - vals[0]) / diff
 
 def energy_nuclear_repulsion(nuclei):
     energy = 0.0
diff --git a/tests/test_repulsion_deriv.py b/tests/test_repulsion_deriv.py
index 65cad52..4f903f5 100644
--- a/tests/test_repulsion_deriv.py
+++ b/tests/test_repulsion_deriv.py
@@ -39,6 +39,8 @@ def testDerivH2O(self):
         fx5 = integrator.repulsion_deriv(cgfs[5], cgfs[3], cgfs[5], cgfs[3], nuclei[0][0], 0)
         fx6 = integrator.repulsion_deriv(cgfs[3], cgfs[5], cgfs[3], cgfs[5], nuclei[1][0], 0)
 
+        print(fx3)
+
         # mol | nuc_id | cgf_id1 | cgf_id2 | cgf_id3 | cgf_id4 | coord
         ans3 = calculate_force_finite_difference(mol, 0, 3, 3, 5, 5, 0)
         ans4 = calculate_force_finite_difference(mol, 1, 3, 3, 5, 5, 0)
@@ -122,25 +124,6 @@ def testDerivH2(self):
         np.testing.assert_almost_equal(fx5, ans5, 4)
         np.testing.assert_almost_equal(fx6, ans6, 4)
 
-    def testDerivRepulsionGTO(self):
-        # construct integrator object
-        integrator = PyQInt()
-
-        # build gtos
-        gto1 = gto(0.154329, [-0.50, 0.0, 0.0], 3.425251, 0, 0, 0)
-        gto2 = gto(0.154329, [-0.50, 0.0, 0.0], 3.425251, 0, 0, 0)
-        gto3 = gto(0.154329, [ 0.50, 0.0, 0.0], 3.425251, 0, 0, 0)
-        gto4 = gto(0.154329, [ 0.50, 0.0, 0.0], 3.425251, 0, 0, 0)
-
-        # perform integration
-        rep1 = integrator.repulsion_gto_deriv(gto1, gto2, gto3, gto4, 0)
-        rep2 = integrator.repulsion_gto_deriv(gto2, gto1, gto3, gto4, 0)
-
-        # perform finite difference
-        ans1 = calculate_deriv_gto(gto1, gto2, gto3, gto4, 0)
-
-        np.testing.assert_almost_equal(rep1, ans1, 4)
-
 def calculate_force_finite_difference(mol, nuc_id, cgf_id1, cgf_id2, cgf_id3, cgf_id4, coord):
     # build integrator object
     integrator = PyQInt()