Repository: qubekit/QUBEKit (language: python)

Issue:
General QM support through QCEngine
We currently use QCEngine for optimization and hessian calculations, but only for very specific settings (PSI4 and Geometric) it would be good to allow more QC programs via qcengine for optimizations, hessians and torsiondrives. The single point calculations used to derive charges and LJ parameters are more specialized since they involve the use of an implicit solvent so we can be more restrictive on the use of program there. This could be achieved by changing the check done [here](https://github.com/qubekit/QUBEKit/blob/267f5fbb3ed0f5bdec913d8ffe71bbbf18f361fb/QUBEKit/run.py#L754) to check the program is not `gaussian` as this is the only engine we support that is not currently available through QCEngine. A similar check could also be done at the hessian and torsion scan stages. 
Change MM optimise to initial optimise
It might be a good idea to change the MM optimise stage to be a general initial optimisation stage which aims to help speed up the full QM optimisation. If we allow the use of QCEngine at this stage we get access to QM/MM/ML methods and semi-empirical methods including xtb this way we expand the range of usable methods while keeping the current options and use a single interface which should keep development cost low. The downside of this would be that the openforcefield toolkit is needed to run gaff/openff optimisations via the qcengine interface, however, as these packages are now all on conda we can just add these to the recipe. 


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml
index 1dd0e5e5..b78b190b 100644
--- a/.github/workflows/CI.yaml
+++ b/.github/workflows/CI.yaml
@@ -12,20 +12,36 @@ on:
 
 jobs:
   test:
-    name: Test on ${{ matrix.os }}, 🐍=${{ matrix.python-version }}, Program=${{ matrix.conda-env }}
-    runs-on: ${{ matrix.os }}
+    name: Test on ${{ matrix.cfg.os }}, 🐍=${{ matrix.cfg.python-version }}, Program=${{ matrix.cfg.conda-env }}
+    runs-on: ${{ matrix.cfg.os }}
     strategy:
       fail-fast: false
       matrix:
-        python-version: [3.6, 3.7, 3.8]
-        os: [ubuntu-latest, macOS-latest]
-        conda-env: [basic]
+        cfg:
+          - conda-env: basic
+            python-version: 3.7
+            os: ubuntu-latest
 
+          - conda-env: basic
+            python-version: 3.7
+            os: macOS-latest
+
+          - conda-env: everything
+            python-version: 3.7
+            os: ubuntu-latest
+
+          - conda-env: basic
+            python-version: 3.8
+            os: ubuntu-latest
+
+          - conda-env: basic
+            python-version: 3.8
+            os: macOS-latest
 
     env:
-      CI_OS: ${{ matrix.os }}
-      ENV: ${{ matrix.conda-env }}
-      PYVER: ${{ matrix.python-version }}
+      CI_OS: ${{ matrix.cfg.os }}
+      ENV: ${{ matrix.cfg.conda-env }}
+      PYVER: ${{ matrix.cfg.python-version }}
 
     steps:
       - uses: actions/checkout@v2
@@ -41,8 +57,8 @@ jobs:
         uses: conda-incubator/setup-miniconda@v2
         with:
           activate-environment: test
-          environment-file: devtools/conda-envs/${{ matrix.conda-env }}.yaml
-          python-version: ${{ matrix.python-version }}
+          environment-file: devtools/conda-envs/${{ matrix.cfg.conda-env }}.yaml
+          python-version: ${{ matrix.cfg.python-version }}
           auto-activate-base: false
 
       - name: Environment Information
diff --git a/QUBEKit/GUI/gui.py b/QUBEKit/GUI/gui.py
index 525edfbc..4b8b8109 100644
--- a/QUBEKit/GUI/gui.py
+++ b/QUBEKit/GUI/gui.py
@@ -6,7 +6,7 @@
 from PyQt5.QtCore import QUrl
 from PyQt5.QtWebEngineWidgets import QWebEnginePage, QWebEngineView
 
-from QUBEKit.ligand import Ligand
+from QUBEKit.molecules.ligand import Ligand
 
 
 class MainWindow(QtWidgets.QMainWindow):
diff --git a/QUBEKit/data/capped_leu.pdb b/QUBEKit/data/capped_leu.pdb
index 365ac6fb..56490951 100644
--- a/QUBEKit/data/capped_leu.pdb
+++ b/QUBEKit/data/capped_leu.pdb
@@ -7,31 +7,31 @@ ATOM      3  HY2 ACE P   1      -0.575   3.023  -2.673  1.00  0.00           H
 ATOM      4  HY3 ACE P   1      -2.186   2.880  -1.953  1.00  0.00           H  
 ATOM      5  CY  ACE P   1      -0.657   1.876  -0.935  1.00  0.00           C  
 ATOM      6  OY  ACE P   1       0.440   2.278  -0.627  1.00  0.00           O  
-ATOM      7  N   LEU P   1      -1.424   1.034  -0.270  1.00  0.00           N  
-ATOM      8  HN  LEU P   1      -2.276   0.814  -0.764  1.00  0.00           H  
-ATOM      9  CA  LEU P   1      -1.219   0.343   1.054  1.00  0.00           C  
-ATOM     10  HA  LEU P   1      -2.149   0.453   1.600  1.00  0.00           H  
-ATOM     11  CB  LEU P   1      -1.074  -1.163   0.752  1.00  0.00           C  
-ATOM     12  HB1 LEU P   1      -0.879  -1.684   1.698  1.00  0.00           H  
-ATOM     13  HB2 LEU P   1      -2.054  -1.529   0.429  1.00  0.00           H  
-ATOM     14  CG  LEU P   1      -0.019  -1.546  -0.276  1.00  0.00           C  
-ATOM     15  HG  LEU P   1      -0.260  -1.078  -1.222  1.00  0.00           H  
-ATOM     16  CD1 LEU P   1      -0.085  -3.047  -0.537  1.00  0.00           C  
-ATOM     17 HD11 LEU P   1       0.638  -3.327  -1.308  1.00  0.00           H  
-ATOM     18 HD12 LEU P   1       0.163  -3.615   0.380  1.00  0.00           H  
-ATOM     19 HD13 LEU P   1      -1.093  -3.350  -0.881  1.00  0.00           H  
-ATOM     20  CD2 LEU P   1       1.373  -1.129   0.147  1.00  0.00           C  
-ATOM     21 HD21 LEU P   1       2.104  -1.566  -0.552  1.00  0.00           H  
-ATOM     22 HD22 LEU P   1       1.468  -0.057   0.099  1.00  0.00           H  
-ATOM     23 HD23 LEU P   1       1.597  -1.493   1.157  1.00  0.00           H  
-ATOM     24  C   LEU P   1      -0.071   0.913   2.074  1.00  0.00           C  
-ATOM     25  O   LEU P   1       0.040   0.459   3.189  1.00  0.00           O  
-ATOM     26  NT  NME P   1       0.702   1.942   1.724  1.00  0.00           N  
-ATOM     27  HNT NME P   1       0.624   2.190   0.705  1.00  0.00           H  
-ATOM     28  CAT NME P   1       1.656   2.626   2.587  1.00  0.00           C  
-ATOM     29  HT1 NME P   1       2.536   2.007   2.725  1.00  0.00           H  
-ATOM     30  HT2 NME P   1       1.949   3.565   2.144  1.00  0.00           H  
-ATOM     31  HT3 NME P   1       1.215   2.817   3.556  1.00  0.00           H  
+ATOM      7  N   LEU P   2      -1.424   1.034  -0.270  1.00  0.00           N
+ATOM      8  HN  LEU P   2      -2.276   0.814  -0.764  1.00  0.00           H
+ATOM      9  CA  LEU P   2      -1.219   0.343   1.054  1.00  0.00           C
+ATOM     10  HA  LEU P   2      -2.149   0.453   1.600  1.00  0.00           H
+ATOM     11  CB  LEU P   2      -1.074  -1.163   0.752  1.00  0.00           C
+ATOM     12  HB1 LEU P   2      -0.879  -1.684   1.698  1.00  0.00           H
+ATOM     13  HB2 LEU P   2      -2.054  -1.529   0.429  1.00  0.00           H
+ATOM     14  CG  LEU P   2      -0.019  -1.546  -0.276  1.00  0.00           C
+ATOM     15  HG  LEU P   2      -0.260  -1.078  -1.222  1.00  0.00           H
+ATOM     16  CD1 LEU P   2      -0.085  -3.047  -0.537  1.00  0.00           C
+ATOM     17 HD11 LEU P   2       0.638  -3.327  -1.308  1.00  0.00           H
+ATOM     18 HD12 LEU P   2       0.163  -3.615   0.380  1.00  0.00           H
+ATOM     19 HD13 LEU P   2      -1.093  -3.350  -0.881  1.00  0.00           H
+ATOM     20  CD2 LEU P   2       1.373  -1.129   0.147  1.00  0.00           C
+ATOM     21 HD21 LEU P   2       2.104  -1.566  -0.552  1.00  0.00           H
+ATOM     22 HD22 LEU P   2       1.468  -0.057   0.099  1.00  0.00           H
+ATOM     23 HD23 LEU P   2       1.597  -1.493   1.157  1.00  0.00           H
+ATOM     24  C   LEU P   2      -0.071   0.913   2.074  1.00  0.00           C
+ATOM     25  O   LEU P   2       0.040   0.459   3.189  1.00  0.00           O
+ATOM     26  NT  NME P   3       0.702   1.942   1.724  1.00  0.00           N
+ATOM     27  HNT NME P   3       0.624   2.190   0.705  1.00  0.00           H
+ATOM     28  CAT NME P   3       1.656   2.626   2.587  1.00  0.00           C
+ATOM     29  HT1 NME P   3       2.536   2.007   2.725  1.00  0.00           H
+ATOM     30  HT2 NME P   3       1.949   3.565   2.144  1.00  0.00           H
+ATOM     31  HT3 NME P   3       1.215   2.817   3.556  1.00  0.00           H
 CONECT    1    2                                                      
 CONECT    2    1    3    4    5                                       
 CONECT    2                                                           
diff --git a/QUBEKit/data/gaussian.fchk b/QUBEKit/data/gaussian.fchk
new file mode 100644
index 00000000..610dc46d
--- /dev/null
+++ b/QUBEKit/data/gaussian.fchk
@@ -0,0 +1,129 @@
+Virial Ratio                               R      2.004558755529153E+00
+SCF Energy                                 R     -1.931393770857046E+02
+Total Energy                               R     -1.931393770857046E+02
+RMS Density                                R      7.648940185449268E-09
+External E-field                           R   N=          35
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+IOpCl                                      I                0
+IROHF                                      I                0
+Atom Modified MM Charges                   R   N=          10
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+Link Distances                             R   N=          40
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
+Cartesian Gradient                         R   N=          30
+ -1.01049212E-02  1.66760481E-02 -6.67413060E-03 -1.81952308E-02 -4.23809696E-02
+  1.30554873E-02  1.52557150E-02  1.73802654E-02 -7.35784961E-03  2.03435545E-02
+ -1.02864200E-03  1.32130300E-02  3.41955166E-03  9.08805029E-03  1.71492181E-02
+  9.54261690E-04  5.75134689E-03 -1.39221813E-02 -1.01400161E-02 -1.34540923E-02
+  6.96887563E-05 -6.59511817E-03  1.80285533E-03 -2.02985913E-02 -9.77845533E-04
+  7.88438894E-03  1.20844530E-02  6.04004889E-03 -1.71925106E-03 -7.31912443E-03
+Cartesian Force Constants                  R   N=         465
+  4.87352968E-01  9.84610254E-03  5.48875488E-01 -2.85051990E-02 -2.31530236E-02
+  4.94849324E-01 -2.04922090E-01  6.27711490E-02  5.49889686E-03  5.86994326E-01
+  5.09202709E-02 -1.04303881E-01 -2.40146168E-03  2.29437759E-02  8.56110638E-01
+  1.95674624E-03 -1.52525982E-03 -8.17702107E-02 -4.87195918E-02 -1.67537969E-01
+  3.03371301E-01  6.19459027E-03  4.37709374E-02 -8.42863795E-03 -1.15403499E-01
+ -3.12262203E-02  1.34973686E-02  1.07878579E-01  2.41827858E-02 -5.13107887E-02
+  1.86356868E-02 -2.73550701E-02 -6.29274481E-01  1.52829244E-01  2.46643824E-02
+  7.23114521E-01 -3.51273854E-03  1.63674891E-02  1.55094776E-02  1.34670001E-02
+  1.52979374E-01 -1.24078057E-01 -1.65480115E-02 -1.90729660E-01  8.74757145E-02
+ -1.99888385E-02 -8.86345266E-03  9.08486133E-03 -2.11678387E-01 -4.18458752E-02
+  3.10990007E-02  4.84752633E-03 -2.25666819E-02  5.28597603E-03  4.89334596E-01
+  1.09287903E-02  8.97990929E-03  5.70239948E-04 -5.88986717E-02 -1.05931569E-01
+  1.26981512E-02 -4.01383790E-02 -3.13137047E-02  1.97574907E-02 -1.66598335E-02
+  5.85240656E-01  2.32418493E-03  3.56049034E-03  2.57476793E-02  2.97504581E-02
+  8.26920529E-03 -8.50936264E-02  1.05590000E-02  1.61697229E-02  1.01874758E-02
+ -2.07020701E-02  5.73834300E-03  4.83052898E-01 -9.42579190E-02 -3.20395204E-02
+ -8.67966929E-02 -1.66068911E-02 -8.78216111E-03 -2.75666884E-02  1.48675369E-03
+  1.41175934E-04  3.52025839E-03 -1.00141418E-03  1.64982244E-03  4.76720203E-03
+  1.04020542E-01 -3.04747974E-02 -6.82266184E-02 -5.29072249E-02  6.58636262E-03
+  3.91803980E-03  1.53394390E-02  1.03270113E-03 -2.45579727E-03 -1.86566379E-03
+ -4.22300489E-05 -4.46312971E-04 -1.09220532E-03  3.26411121E-02  7.12468966E-02
+ -8.12123592E-02 -5.35591974E-02 -1.90799825E-01 -2.16499686E-03  1.28819367E-05
+ -2.78408425E-03  1.46314079E-04  3.77940047E-04  2.12964993E-03 -2.63437448E-05
+  2.63878103E-04 -1.29729232E-03  9.41278621E-02  5.92095339E-02  2.06781108E-01
+ -1.07558006E-01 -2.93021620E-02  1.02600982E-01 -2.16768136E-02 -6.98455989E-03
+  2.68779951E-02  3.10548640E-03 -2.95768529E-03 -2.33531744E-03 -3.01508068E-03
+  1.99029924E-03 -3.23876708E-03  9.49110606E-03  4.70352751E-03 -1.33171666E-02
+  1.23168210E-01 -3.21174525E-02 -6.78532575E-02  5.15574122E-02  7.20604121E-03
+  3.56760648E-03 -1.46335748E-02  1.45248622E-03 -9.78211430E-04  1.20307501E-03
+  6.28906766E-04 -8.73011433E-04  1.07748600E-03  6.03097590E-03  1.46233262E-03
+ -7.43208643E-03  3.15653010E-02  6.74340418E-02  1.03923085E-01  4.86769818E-02
+ -2.07161928E-01 -1.43161939E-03  2.05508432E-03 -2.69503914E-04 -1.60537519E-03
+ -4.21595251E-04  1.65349153E-03 -5.63654370E-04 -5.10458629E-04 -9.43581578E-05
+  1.38730388E-02  5.94164876E-03 -1.66040553E-02 -1.11962525E-01 -5.51223505E-02
+  2.20892120E-01 -6.14689668E-02 -4.23513784E-02  5.09171628E-03  1.08401429E-02
+  2.61892160E-02  3.04409974E-05 -6.20844369E-03  7.28334735E-03 -2.62837724E-03
+ -7.46434725E-04  1.07075714E-04 -3.36252962E-04 -2.33924539E-03 -1.44740005E-02
+  2.43726188E-03 -2.35399099E-03 -1.49398417E-02 -1.94803095E-03  6.41456669E-02
+ -3.39363174E-02 -2.65036110E-01  8.19539972E-03 -1.07014976E-02 -1.63049302E-02
+  1.03998651E-03  7.01818509E-03 -3.77573267E-03  8.28799139E-04 -8.12437298E-04
+ -9.40466499E-05  6.33488997E-04  1.08202956E-04 -6.01225377E-03  9.92186570E-04
+  5.07416715E-04 -3.66437405E-03 -6.36226174E-04  3.81508702E-02  2.92290695E-01
+  4.46922736E-03  9.33691053E-03 -5.41299064E-02 -2.69598073E-04  1.85917440E-03
+ -3.21272523E-03 -1.28426647E-03  1.01263744E-03  2.65547954E-03  3.50500673E-04
+ -1.74486332E-04 -2.78239022E-04 -1.59764100E-03 -2.44229435E-02  2.50060545E-03
+  8.71290128E-04  2.33583846E-02  1.58803129E-03 -2.91656471E-03 -1.08320671E-02
+  5.03797290E-02 -2.01915688E-03 -1.47329290E-03  4.82191992E-03 -2.48421996E-02
+  6.88310222E-03 -2.40024784E-02  3.56725247E-03  1.47718705E-03  2.06662877E-03
+ -1.26953119E-01  2.36853875E-02 -9.94984923E-02  4.79550203E-04 -3.15987788E-04
+ -3.07389247E-04 -1.32221560E-03  2.86729575E-04 -8.83787126E-05  1.23962243E-04
+ -1.73216630E-04 -6.48305529E-05  1.44356511E-01 -2.82371785E-04 -4.71398865E-04
+  8.75583350E-04 -8.44808275E-03  3.11219641E-03 -1.34989138E-02 -1.29863034E-04
+ -6.69370300E-04  1.89841401E-03  2.43410718E-02 -6.09469339E-02  3.02525051E-02
+  2.66639119E-04 -1.86938644E-04  2.55919748E-05 -1.68429677E-04  4.70263523E-04
+  3.18863393E-05 -1.28432415E-04  4.13947310E-04  1.36636837E-04 -2.71218440E-02
+  6.07492980E-02  7.89798438E-04 -3.59440507E-04 -1.40060035E-03  4.12007917E-03
+ -2.09319235E-04  3.86880460E-03  7.57227917E-04 -3.13387340E-04  9.64823272E-04
+ -9.77575061E-02  3.00273846E-02 -1.65153629E-01 -2.19603255E-04  5.62977133E-05
+  1.18066654E-05  3.15697548E-04 -1.17290258E-04  4.23786207E-05  1.33386586E-04
+ -1.01955186E-04 -7.35448439E-05  1.08841668E-01 -3.51065733E-02  1.78450073E-01
+ -2.44486150E-03 -1.97619247E-03 -3.56840182E-03 -1.15117655E-02  1.25170253E-02
+  2.63569518E-02  1.25772583E-03 -1.77114452E-04 -2.89675898E-03 -7.52482123E-02
+  3.21239838E-02  7.81939839E-02 -1.34912326E-03  2.70236061E-04  1.81869705E-04
+ -3.25276791E-05 -2.36357911E-04  1.30139062E-05  2.38661089E-04 -2.80277016E-04
+  6.87968655E-05  9.10008959E-03 -8.21068346E-03 -1.68799105E-02  8.24665822E-02
+ -3.95778962E-04 -6.56979197E-04 -1.04798003E-03 -2.38045984E-03  6.06516578E-03
+  1.47347834E-02 -8.89362757E-04 -2.18903346E-03 -1.82719135E-03  3.03365134E-02
+ -8.35681778E-02 -7.47243926E-02 -2.34493112E-04  4.11667098E-04  3.51873600E-05
+  2.62864806E-04 -3.96584499E-05  5.35492577E-06 -4.78276985E-05  3.91742669E-04
+ -1.72084530E-04 -3.31350666E-03  1.50351388E-03  5.75441200E-03 -3.37573693E-02
+  8.39087390E-02 -1.83089485E-04  6.23046554E-04 -6.81811945E-04  2.03335970E-03
+ -2.54605522E-03 -6.85418262E-03  1.67571056E-04  8.80453520E-04  1.97085564E-03
+  7.42661007E-02 -7.62099443E-02 -2.17993700E-01 -2.93212444E-05 -3.07584559E-05
+  1.16025117E-04  1.96516129E-04  8.62546727E-06 -7.20323810E-06 -1.96678542E-04
+ -2.00555651E-05  4.05521116E-06  9.81213461E-03 -7.73822830E-03 -1.58021691E-02
+ -8.31414329E-02  8.30836787E-02  2.34790818E-01 -8.87719158E-04 -3.82190133E-04
+  2.00555592E-04  8.80717650E-03 -3.06145738E-02  4.70255100E-04 -6.72597196E-03
+ -4.69232681E-03  3.58134037E-03 -5.55506353E-02  4.52115252E-02 -1.81924649E-03
+  7.66413838E-05  7.30763390E-05  1.34947857E-04  1.93832050E-04  1.23211394E-04
+ -2.09553869E-04 -2.23135152E-03  1.19070972E-04  3.73085766E-04 -2.49067449E-03
+  1.98819962E-02 -1.00837292E-04 -2.47656846E-03  1.04194201E-02 -2.92515998E-03
+  6.12852710E-02  1.32876848E-03  3.63621785E-06 -3.24631911E-04  8.27645334E-03
+ -1.69587852E-02  5.54112556E-04 -5.55486706E-03 -1.14740154E-03  1.38787364E-03
+  3.54840185E-02 -3.11046808E-01  1.01153563E-02  2.18246209E-04  2.88984923E-04
+  7.40838892E-05  3.83427513E-04  4.74268446E-04 -2.03255484E-05  2.10971471E-04
+  1.79106237E-03 -1.02162275E-04  6.54415594E-05 -3.97457741E-03  3.69871496E-04
+ -2.73250608E-04 -5.82697954E-03  1.94923759E-03 -4.01392095E-02  3.36396600E-01
+ -4.96559112E-05  3.20030353E-05 -1.62199731E-04 -2.28398782E-03  7.51908506E-03
+ -3.17771534E-03  2.73880952E-03  1.55895783E-03  1.53108967E-03 -1.03686512E-03
+  7.83940167E-03 -4.90772088E-02 -7.84146428E-05 -2.28123459E-04 -5.39388041E-05
+ -8.70411269E-06  1.00318645E-04 -3.89737362E-05  3.33098671E-04 -9.95568674E-05
+  5.66515005E-04 -1.58078159E-03  2.31230978E-02 -9.07943603E-04  1.67188808E-03
+ -2.58417679E-02  4.45731319E-03  2.94612941E-04 -1.40034158E-02  4.68630622E-02
+Dipole Moment                              R   N=           3
+ -1.02402519E-01 -1.17188635E+00  3.46976141E-01
\ No newline at end of file
diff --git a/QUBEKit/data/gaussian.log b/QUBEKit/data/gaussian.log
new file mode 100644
index 00000000..072f7441
--- /dev/null
+++ b/QUBEKit/data/gaussian.log
@@ -0,0 +1,45 @@
+
+ ******************************************
+ Gaussian 09:  ES64L-G09RevD.01 24-Apr-2013
+                25-Feb-2021
+ ******************************************
+ %Mem=10GB
+ %NProcShared=10
+ Will use up to   10 processors via shared memory.
+ %Chk=lig
+ ----------------------------------------------------------------------
+ # wB97XD/6-311++G(d,p) SCF=(XQC,MaxConventionalCycles=350) nosymm Forc
+ e=Nostep geom=check
+ ----------------------------------------------------------------------
+ 1/10=10,29=2,30=1,38=1/1;
+ 2/12=2,15=1,40=1/2;
+ 3/5=4,6=6,7=1111,11=2,14=-4,16=1,25=1,30=1,71=1,74=-58/1,2,3;
+ 4//1;
+ 5/5=2,8=3,13=1,38=5,85=350/2,8;
+ 6/7=2,8=2,9=2,10=2,28=1/1;
+ 7/30=1/1,2,3,16;
+ 99//99;
+ ---
+ MOL
+ ---
+ Structure from the checkpoint file:  "lig.chk"
+ Charge =  0 Multiplicity = 1
+ Redundant internal coordinates found in file.
+ C,0,0.8954058983,0.0094381372,0.0017746914
+ C,0,0.4563749929,0.9961724141,-1.0864832887
+ C,0,0.3790349334,-1.3928398142,-0.3414115381
+ C,0,0.3203546134,0.4499237371,1.3532702958
+ C,0,2.4263469757,-0.0160429706,0.08207858
+ H,0,0.809267297,2.0083703738,-0.864472636
+ H,0,-0.6346566891,1.0311989956,-1.1689764192
+ H,0,0.8582062449,0.7075442888,-2.0630052294
+ H,0,0.6745306299,-2.1181453729,0.4233277044
+ H,0,0.7802629978,-1.7350981568,-1.3007940612
+ H,0,-0.7133710068,-1.4033610104,-0.4112307365
+ H,0,0.6175037941,-0.240431324,2.1491119719
+ H,0,-0.7735780286,0.4788473522,1.3237691689
+ H,0,0.6760014052,1.4491890905,1.6239143311
+ H,0,2.767393633,-0.7158305746,0.8516966543
+ H,0,2.8243884356,0.9738027398,0.3272241876
+
+ Normal termination of Gaussian
\ No newline at end of file
diff --git a/QUBEKit/data/templates/gaussian.com b/QUBEKit/data/templates/gaussian.com
new file mode 100644
index 00000000..4f04e74d
--- /dev/null
+++ b/QUBEKit/data/templates/gaussian.com
@@ -0,0 +1,13 @@
+%Mem={{ memory }}GB
+%NProcShared={{ threads }}
+%Chk=lig
+# {{ theory }}/{{ basis }}  SCF=(XQC, MaxConventionalCycles={{ scf_maxiter }}) nosymm {{ driver }}
+
+{{ title }}
+
+{{ charge }} {{ multiplicity }}
+{%- for element, coords in data %}
+{{ element }}  {{ '{: .10f}'.format(coords[0]) }} {{ '{: .10f}'.format(coords[1]) }} {{ '{: .10f}'.format(coords[2]) }}
+{%- endfor %}
+
+
diff --git a/QUBEKit/dihedrals.py b/QUBEKit/dihedrals.py
index 248e95e2..73d74c4b 100644
--- a/QUBEKit/dihedrals.py
+++ b/QUBEKit/dihedrals.py
@@ -14,8 +14,8 @@
 from scipy.optimize import differential_evolution, minimize
 from scipy.stats import linregress
 
-from QUBEKit.engines import PSI4, Gaussian, OpenMM, RDKit
-from QUBEKit.ligand import Ligand
+from QUBEKit.engines import Gaussian, OpenMM
+from QUBEKit.molecules.ligand import Ligand
 from QUBEKit.utils import constants
 from QUBEKit.utils.exceptions import TorsionDriveFailed
 from QUBEKit.utils.file_handling import make_and_change_into
@@ -48,9 +48,9 @@ def __init__(self, molecule: Ligand, constraints_made=None):
         self.molecule.convergence = "GAU"
         self.constraints_made = constraints_made
 
-        self.qm_engine = {"psi4": PSI4, "g09": Gaussian, "g16": Gaussian}.get(
-            molecule.bonds_engine
-        )(molecule)
+        self.qm_engine = {"g09": Gaussian, "g16": Gaussian}.get(molecule.bonds_engine)(
+            molecule
+        )
         self.native_opt = True
 
         # Ensure geometric can only be used with psi4 so far
@@ -1066,6 +1066,7 @@ def _create_rdkit_molecules(self, coordinates):
         :param coordinates: A list of numpy arrays used to generate the conformers
         :return: a list of rdkit molecules each corresponding to a point on the torsionscan
         """
+        from QUBEKit.molecules.utils import RDKit
 
         self.rmsd_atoms = []
         for coord in coordinates:
@@ -1082,6 +1083,8 @@ def scan_rmsd(self, coordinates):
         :param coordinates:  A list of numpy coordinate arrays
         :return: a list containing the rmsd values for each pair of coordinates
         """
+        from QUBEKit.molecules.utils import RDKit
+
         # Make sure the number of coordinates we pass is the same as the number of reference positions that we have
         if len(coordinates) != len(self.rmsd_atoms):
             print(
diff --git a/QUBEKit/engines/__init__.py b/QUBEKit/engines/__init__.py
index 487a4e36..2b7d090d 100644
--- a/QUBEKit/engines/__init__.py
+++ b/QUBEKit/engines/__init__.py
@@ -1,14 +1,10 @@
-#!/usr/bin/env python3
 """
 Ensures that errors with particular engines are only raised if said engine is being used.
 No point raising "cannot execute psi4 command" if user is using g09. See try_load docstring for info.
 """
-from QUBEKit.utils.helpers import try_load
-
-from .chargemol import Chargemol
-from .gaussian import Gaussian
-from .rdkit import RDKit
-
-OpenMM = try_load("OpenMM", "QUBEKit.engines.openmm")
-PSI4 = try_load("PSI4", "QUBEKit.engines.psi4")
-QCEngine = try_load("QCEngine", "QUBEKit.engines.qcengine")
+from QUBEKit.engines.chargemol import Chargemol
+from QUBEKit.engines.gaussian import Gaussian
+from QUBEKit.engines.gaussian_harness import GaussianHarness
+from QUBEKit.engines.geometry_optimiser import GeometryOptimiser
+from QUBEKit.engines.openmm import OpenMM
+from QUBEKit.engines.qcengine import QCEngine
diff --git a/QUBEKit/engines/base_engine.py b/QUBEKit/engines/base_engine.py
index 2dd9912d..de7fbf1c 100644
--- a/QUBEKit/engines/base_engine.py
+++ b/QUBEKit/engines/base_engine.py
@@ -1,4 +1,70 @@
-#!/usr/bin/env python3
+from typing import Dict, Optional
+
+import qcelemental as qcel
+import qcengine as qcng
+from pydantic import BaseModel, PositiveInt, validator
+
+from QUBEKit.utils.exceptions import SpecificationError
+
+
+class BaseEngine(BaseModel):
+    """
+    A common base model for engines to run computations via external programs.
+    """
+
+    program: str = "psi4"
+    basis: Optional[str] = "6-311++G(d,p)"
+    method: str = "wB97X-D"
+    cores: PositiveInt = 4
+    memory: PositiveInt = 4
+
+    class Config:
+        validate_assignment = True
+        arbitrary_types_allowed = True
+        fields = {
+            "program": {
+                "description": "The name of the program which should be used to carry out the computation, such as psi4"
+            },
+            "basis": {
+                "description": "The basis that should be used in the computation."
+            },
+            "method": {
+                "description": "The method that should be used for the computation."
+            },
+            "cores": {
+                "description": "The number of cores to be allocated to the computation."
+            },
+            "memory": {
+                "description": "The amount of memory that should be allocated to the computation in GB."
+            },
+        }
+
+    @validator("program")
+    def validate_program(cls, program: str) -> str:
+        """
+        Validate the choice of program against those supported by QCEngine and QUBEKit.
+        """
+        programs = qcng.list_available_programs()
+        programs.discard("dftd3")
+
+        if program.lower() not in programs:
+            raise SpecificationError(
+                f"The program {program} is not available, available programs are {programs}"
+            )
+        return program.lower()
+
+    @property
+    def local_options(self) -> Dict[str, int]:
+        """return the local options."""
+        return {"memory": self.memory, "ncores": self.cores}
+
+    @property
+    def qc_model(self) -> qcel.models.common_models.Model:
+        """
+        Build the QC model for the computation.
+        """
+        model = qcel.models.common_models.Model(method=self.method, basis=self.basis)
+        return model
 
 
 class Engines:
diff --git a/QUBEKit/engines/gaussian_harness.py b/QUBEKit/engines/gaussian_harness.py
new file mode 100644
index 00000000..5e0abfb5
--- /dev/null
+++ b/QUBEKit/engines/gaussian_harness.py
@@ -0,0 +1,302 @@
+"""
+A custom gaussian harness for QCEngine which should be registered with qcengine.
+"""
+from typing import Any, Dict, List, Optional, Tuple
+
+from qcelemental.models import AtomicInput, AtomicResult
+from qcelemental.util import which
+from qcengine import register_program
+from qcengine.exceptions import UnknownError
+from qcengine.programs import ProgramHarness
+from qcengine.util import execute
+
+from QUBEKit.utils import constants
+from QUBEKit.utils.file_handling import get_data
+
+
+class GaussianHarness(ProgramHarness):
+    """
+    A very minimal Gaussian wrapper for optimisations and torsiondrives.
+    Note a lot of key information is not extracted need for QCArchive.
+    """
+
+    name = "gaussian"
+    scratch = True
+    thread_safe = True
+    thread_parallel = True
+    node_parallel = False
+    managed_memory = True
+
+    @staticmethod
+    def found(raise_error: bool = False) -> bool:
+        """
+        Check if gaussian can be found via the command line, check for g09 and g16.
+        """
+        gaussian09 = which("g09", return_bool=True)
+        gaussian16 = which("g16", return_bool=True)
+        return gaussian09 or gaussian16
+
+    def get_version(self) -> str:
+        """
+        Work out which version of gaussian we should be running.
+        """
+        if which("g09", return_bool=True, raise_error=False):
+            return "g09"
+        else:
+            which(
+                "g16",
+                return_bool=True,
+                raise_error=True,
+                raise_msg="Gaussian 09/16 can not be found make sure they are available.",
+            )
+            return "g16"
+
+    def compute(self, input_data: AtomicInput, config: "TaskConfig") -> AtomicResult:
+        """
+        Run the compute job via the gaussian CLI.
+        """
+        # we always use an internal template
+        job_inputs = self.build_input(input_model=input_data, config=config)
+        exe_success, proc = self.execute(job_inputs)
+        if exe_success:
+            result = self.parse_output(proc["outfiles"], input_data)
+            return result
+        else:
+            raise UnknownError(proc["stderr"])
+
+    def execute(
+        self,
+        inputs: Dict[str, Any],
+        extra_outfiles: Optional[List[str]] = None,
+        extra_commands: Optional[List[str]] = None,
+        scratch_name: Optional[str] = None,
+        timeout: Optional[int] = None,
+    ) -> Tuple[bool, Dict[str, Any]]:
+        """
+        Run the gaussian single point job and fchk conversion
+        """
+        # collect together inputs
+        infiles = inputs["infiles"]
+        outfiles = ["gaussian.log", "lig.chk"]
+        gaussian_version = self.get_version()
+        commands = [gaussian_version, "gaussian.com"]
+        scratch_directory = inputs["scratch_directory"]
+
+        # remember before formatting lig.chk is binary, run calculation
+        exe_success, proc = execute(
+            command=commands,
+            infiles=infiles,
+            outfiles=outfiles,
+            scratch_directory=scratch_directory,
+            as_binary=["lig.chk"],
+        )
+        if exe_success:
+            # now we need to run the conversion of the chk file
+            commands = ["formchk", "lig.chk", "lig.fchk"]
+            infiles = {"lig.chk": proc["outfiles"]["lig.chk"]}
+            chk_success, proc_chk = execute(
+                command=commands,
+                infiles=infiles,
+                outfiles=["lig.fchk"],
+                scratch_directory=scratch_directory,
+                as_binary=["lig.chk"],
+            )
+            # now we want to put this out file into the main proc
+            proc["outfiles"]["lig.fchk"] = proc_chk["outfiles"]["lig.fchk"]
+            # we can dump the chk file
+            del proc["outfiles"]["lig.chk"]
+
+        return exe_success, proc
+
+    @staticmethod
+    def functional_converter(method: str) -> str:
+        """
+        Convert the given function to the correct format for gaussian.
+        """
+        functionals = {"pbe": "PBEPBE", "wb97x-d": "wB97XD"}
+
+        # check for dispersion
+        if "-d3" in method:
+            theory = method.split("-d3")[0]
+            dispersion = "GD3"
+            if "bj" in method:
+                dispersion += "BJ"
+        else:
+            theory = method
+            dispersion = None
+        # convert the theory
+        theory = functionals.get(theory, theory)
+        if dispersion is not None:
+            return f"EmpiricalDispersion={dispersion.upper()} {theory}"
+        return theory
+
+    @staticmethod
+    def driver_conversion(driver: str) -> str:
+        """
+        Convert the qcengine driver enum to the specific keyword for gaussian.
+        """
+        drivers = {"energy": "SP", "gradient": "Force=NoStep", "hessian": "FREQ"}
+        return drivers[driver.lower()]
+
+    def build_input(
+        self,
+        input_model: AtomicInput,
+        config: "TaskConfig",
+        template: Optional[str] = None,
+    ) -> Dict[str, Any]:
+        """
+        Use the template files stored in QUBEKit to build a gaussian input file for the given driver.
+        """
+        import os
+
+        from jinja2 import Template
+
+        template_file = get_data(os.path.join("templates", "gaussian.com"))
+        with open(template_file) as file:
+            template = Template(file.read())
+
+        template_data = {
+            "memory": int(config.memory),
+            "threads": config.ncores,
+            "driver": self.driver_conversion(driver=input_model.driver),
+            "title": "gaussian job",
+        }
+        molecule = input_model.molecule
+        spec = input_model.model
+        theory = self.functional_converter(method=spec.method)
+        template_data["theory"] = theory
+        template_data["basis"] = spec.basis
+        template_data["charge"] = int(molecule.molecular_charge)
+        template_data["multiplicity"] = molecule.molecular_multiplicity
+        template_data["scf_maxiter"] = input_model.extras.get("maxiter", 300)
+        # now we need to build the coords data
+        data = []
+        for i, symbol in enumerate(molecule.symbols):
+            # we must convert the atomic input back to angstroms
+            data.append((symbol, molecule.geometry[i] * constants.BOHR_TO_ANGS))
+        template_data["data"] = data
+
+        rendered_template = template.render(**template_data)
+        return {
+            "infiles": {"gaussian.com": rendered_template},
+            "scratch_directory": config.scratch_directory,
+            "input_result": input_model.copy(deep=True),
+        }
+
+    @staticmethod
+    def check_convergence(logfile: str) -> None:
+        """
+        Check the gaussian log file to make sure we have normal termination.
+        """
+        if "Normal termination of Gaussian" not in logfile:
+            # raise an error with the log file as output
+            raise UnknownError(message=logfile)
+
+    def parse_output(
+        self, outfiles: Dict[str, str], input_model: AtomicInput
+    ) -> AtomicResult:
+        """
+        For the set of output files parse them to extract as much info as possible and return the atomic result.
+        From the fchk file we get the energy and hessian, the gradient is taken from the log file.
+        """
+        properties = {}
+        # make sure we got vaild exit status
+        self.check_convergence(logfile=outfiles["gaussian.log"])
+        version = self.parse_version(logfile=outfiles["gaussian.log"])
+        # build the main data dict
+        output_data = input_model.dict()
+        provenance = {"version": version, "creator": "gaussian", "routine": "CLI"}
+        # collect the total energy from the fchk file
+        logfile = outfiles["lig.fchk"]
+        for line in logfile.split("\n"):
+            if "Total Energy" in line:
+                energy = float(line.split()[3])
+                properties["return_energy"] = energy
+                properties["scf_total_energy"] = energy
+                if input_model.driver == "energy":
+                    output_data["return_result"] = energy
+        if input_model.driver == "gradient":
+            # now we need to parse out the forces
+            gradient = self.parse_gradient(fchfile=outfiles["lig.fchk"])
+            output_data["return_result"] = gradient
+        elif input_model.driver == "hessian":
+            hessian = self.parse_hessian(fchkfile=outfiles["lig.fchk"])
+            output_data["return_result"] = hessian
+
+        output_data["properties"] = properties
+        output_data["schema_name"] = "qcschema_output"
+        output_data["stdout"] = outfiles["gaussian.log"]
+        output_data["success"] = True
+        output_data["provenance"] = provenance
+        return AtomicResult(**output_data)
+
+    @staticmethod
+    def parse_version(logfile: str) -> str:
+        """
+        Try and parse the gaussian version from the logfile.
+        """
+        # the version is printed after the first ***** line
+        lines = logfile.split("\n")
+        star_line = 0
+        for i, line in enumerate(lines):
+            if "******************************************" in line:
+                star_line = i
+                break
+        # now extract the line
+        version = lines[star_line + 1].strip()
+        return version
+
+    @staticmethod
+    def parse_gradient(fchfile: str) -> List[float]:
+        """
+        Parse the cartesian gradient from a gaussian fchk file and return them as a flat list.
+        """
+        gradient = []
+        grad_found = False
+        for line in fchfile.split("\n"):
+            if "Cartesian Gradient" in line:
+                grad_found = True
+            elif grad_found:
+                if "Cartesian Force Constants" in line or "Dipole Moment" in line:
+                    grad_found = False
+                else:
+                    gradient.extend([float(grad) for grad in line.split()])
+        return gradient
+
+    @staticmethod
+    def parse_hessian(fchkfile: str) -> List[float]:
+        """
+        Parse the hessian from the fchk file and return as a flat list in atomic units.
+        Note gaussian gives us only half of the matrix so we have to fix this.
+        #TODO the reshaping or the array could be fragile is there a better way?
+        """
+        import numpy as np
+
+        hessian = []
+        hessian_found = False
+        for line in fchkfile.split("\n"):
+            if line.startswith("Cartesian Force Constants"):
+                hessian_found = True
+            elif hessian_found:
+                if line.startswith("Nonadiabatic coupling") or line.startswith(
+                    "Dipole Moment"
+                ):
+                    hessian_found = False
+                else:
+                    hessian.extend([float(hess) for hess in line.split()])
+        # now we can calculate the size of the hessian
+        # (2 * triangle length) + 0.25 = (x+0.5)^2
+        hess_size = int(np.sqrt(2 * len(hessian) + 0.25) - 0.5)
+        full_hessian = np.zeros((hess_size, hess_size))
+        m = 0
+        for i in range(hess_size):
+            for j in range(i + 1):
+                full_hessian[i, j] = hessian[m]
+                full_hessian[j, i] = hessian[m]
+                m += 1
+
+        return full_hessian.flatten().tolist()
+
+
+# register the gaussian harness this must be done so gaussian can be used!
+register_program(GaussianHarness())
diff --git a/QUBEKit/engines/geometry_optimiser.py b/QUBEKit/engines/geometry_optimiser.py
new file mode 100644
index 00000000..2c1f82f1
--- /dev/null
+++ b/QUBEKit/engines/geometry_optimiser.py
@@ -0,0 +1,243 @@
+"""
+A class which handles general geometry optimisation tasks.
+"""
+import copy
+from typing import Dict, Optional, Union
+
+import qcelemental as qcel
+import qcengine as qcng
+from pydantic import Field, PositiveInt, validator
+from typing_extensions import Literal
+
+from QUBEKit.engines.base_engine import BaseEngine
+from QUBEKit.molecules import Ligand
+from QUBEKit.utils import constants
+from QUBEKit.utils.exceptions import SpecificationError
+
+
+class GeometryOptimiser(BaseEngine):
+    """
+    A general geometry optimiser class which can dispatch the optimisation call to the correct program and method via qcengine.
+    #TODO do we want to expose more optimiser settings?
+    """
+
+    optimiser: str = Field(
+        "geometric",
+        description="The name of the optimisation engine which should be used note only gaussian supports native optimisation.",
+    )
+    program: str = "rdkit"
+    basis: Optional[str] = None
+    method: str = "mmff94"
+    maxiter: PositiveInt = Field(
+        350, description="The maximum number of optimisation steps."
+    )
+    convergence: Literal["GAU", "GAU_TIGHT", "GAU_VERYTIGHT"] = Field(
+        "GAU_TIGHT",
+        description="The convergence critera for the geometry optimisation.",
+    )
+    extras: Optional[Dict] = Field(
+        None,
+        description="Any extra arguments that should be passed to the geometry optimiser, like scf maxiter.",
+    )
+
+    @validator("optimiser")
+    def validate_optimiser(cls, optimiser: str) -> str:
+        """
+        Make sure the chosen optimiser is available.
+        """
+        procedures = qcng.list_available_procedures()
+        if optimiser.lower() not in procedures:
+            raise SpecificationError(
+                f"The optimiser {optimiser} is not available, available optimisers are {procedures}"
+            )
+        return optimiser.lower()
+
+    def _validate_specification(self) -> None:
+        """
+        Validate the specification this is called before running an optimisation to catch errors before run time.
+        """
+        from openff.toolkit.typing.engines.smirnoff import get_available_force_fields
+
+        openff_forcefields = [
+            ff.split(".offxml")[0].lower() for ff in get_available_force_fields()
+        ]
+        # set up some models
+        ani_methods = {"ani1x", "ani1ccx", "ani2x"}
+        xtb_methods = {
+            "gfn0-xtb",
+            "gfn0xtb",
+            "gfn1-xtb",
+            "gfn1xtb",
+            "gfn2-xtb",
+            "gfn2xtb",
+            "gfn-ff",
+            "gfnff",
+        }
+        rdkit_methods = {"uff", "mmff94", "mmff94s"}
+        gaff_forcefields = {
+            "gaff-1.4",
+            "gaff-1.8",
+            "gaff-1.81",
+            "gaff-2.1",
+            "gaff-2.11",
+        }
+        settings = {
+            "openmm": {"antechamber": gaff_forcefields, "smirnoff": openff_forcefields},
+            "torchani": {None: ani_methods},
+            "xtb": {None: xtb_methods},
+            "rdkit": {None: rdkit_methods},
+        }
+        # now check these settings
+        # TODO do we raise an error or just change at run time with a warning?
+        if self.program.lower() != "psi4" and self.optimiser == "optking":
+            raise SpecificationError(
+                f"The optimiser optking currently only supports psi4 as the engine."
+            )
+
+        # we do not validate QM as there are so many options
+        if self.program.lower() in settings:
+            program_settings = settings[self.program.lower()]
+
+            allowed_methods = program_settings.get(self.basis, None)
+            if allowed_methods is None:
+                raise SpecificationError(
+                    f"The Basis {self.basis} is not supported for the program {self.program} please chose from {program_settings.keys()}"
+                )
+            # now check the method
+            method = self.method.split(".offxml")[0].lower()
+            if method not in allowed_methods:
+                raise SpecificationError(
+                    f"The method {method} is not available for the program {self.program}  with basis {self.basis}, please chose from {allowed_methods}"
+                )
+
+    def __init__(self, **data):
+        """
+        Validate.
+        """
+        super().__init__(**data)
+        self._validate_specification()
+
+    def build_optimiser_keywords(self) -> Dict[str, Union[str, float]]:
+        """
+        Based on the selected optimiser return the optimisation keywords.
+        Are there any other options we want to expose?
+        """
+        opt_keywords = {"program": self.program}
+        if self.optimiser == "geometric":
+            opt_keywords["coordsys"] = "dlc"
+            opt_keywords["maxiter"] = self.maxiter
+            opt_keywords["convergence_set"] = self.convergence
+        elif self.optimiser == "optking":
+            opt_keywords["geom_maxiter"] = self.maxiter
+            opt_keywords["G_CONVERGENCE"] = self.convergence
+        return opt_keywords
+
+    @staticmethod
+    def check_available(program: str, optimiser: str) -> bool:
+        """
+        A wrapper around the general harness check in qcengine,
+        here we can check that a optimiser and program are installed before run time.
+        """
+        _ = qcng.get_program(name=program, check=True)
+        _ = qcng.get_procedure(name=optimiser)
+        return True
+
+    def optimise(self, molecule: Ligand, allow_fail: bool = False) -> Ligand:
+        """
+        For the given specification in the class run an optimisation on the ligand.
+
+        Run the specified optimisation on the ligand the final coordinates are extracted and stored in the ligand.
+        The optimisation schema is dumped to file along with the optimised geometry and the trajectory.
+
+        Args:
+            molecule:
+                The molecule which should be optimised
+            allow_fail:
+                If we should not raise an error if the molecule fails to be optimised, this will extract the last geometry
+                from the trajectory and return it.
+
+        Returns:
+            A new copy of the molecule at the optimised coordinates.
+        """
+        # first validate the settings
+        self._validate_specification()
+        # now validate that the programs are installed
+        self.check_available(program=self.program, optimiser=self.optimiser)
+
+        # now we need to distribute the job
+        model = self.qc_model
+        specification = qcel.models.procedures.QCInputSpecification(model=model)
+        initial_mol = molecule.to_qcschema()
+        opt_task = qcel.models.OptimizationInput(
+            initial_molecule=initial_mol,
+            input_specification=specification,
+            keywords=self.build_optimiser_keywords(),
+        )
+        opt_result = qcng.compute_procedure(
+            input_data=opt_task,
+            procedure=self.optimiser,
+            raise_error=False,
+            local_options=self.local_options,
+        )
+        # dump info to file
+        result_mol = self.handle_output(molecule=molecule, opt_output=opt_result)
+        # check if we can/have failed and raise the error
+        if not opt_result.success and not allow_fail:
+            raise RuntimeError(
+                f"{opt_result.error.error_type}: {opt_result.error.error_message}"
+            )
+
+        return result_mol
+
+    def handle_output(
+        self,
+        molecule: Ligand,
+        opt_output: Union[
+            qcel.models.OptimizationResult, qcel.models.common_models.FailedOperation
+        ],
+    ) -> Ligand:
+        """
+        Sort the output of the optimisation depending on success.
+
+        Take a molecule and an optimisation result or failed operation and unpack the trajectory and final geometry into
+        a copy of the molecule. Also save the result to file.
+
+        Args:
+            molecule:
+                The molecule that the optimisation was ran on.
+            opt_output:
+                The output of an optimisation which may be a valid result or a FailedOperation.
+
+        Returns:
+            A copy of the molecule with the final geometry in the trajectory saved.
+        """
+        result_mol = copy.deepcopy(molecule)
+        # Store the entire result to file
+        with open("result.json", "w") as out:
+            out.write(opt_output.json())
+
+        # Now sort the result
+        if opt_output.success:
+            # passed operation so act normal
+            trajectory = [result.molecule for result in opt_output.trajectory]
+        else:
+            try:
+                trajectory = [
+                    qcel.models.Molecule.from_data(mol["molecule"])
+                    for mol in opt_output.input_data["trajectory"]
+                ]
+            except KeyError:
+                # this is a fatal error as no compute was performed so exit early with no traj
+                return result_mol
+
+        # TODO how do we get this in the log
+        print(f"Optimisation finished in {len(trajectory)} steps.")
+        traj = [mol.geometry * constants.BOHR_TO_ANGS for mol in trajectory]
+        result_mol.coordinates = traj[-1]
+        # write to file
+        result_mol.to_file(file_name="opt.xyz")
+        result_mol.to_multiconformer_file(
+            file_name="opt_trajectory.xyz", positions=traj
+        )
+
+        return result_mol
diff --git a/QUBEKit/engines/psi4.py b/QUBEKit/engines/psi4.py
deleted file mode 100644
index fa3b7ccb..00000000
--- a/QUBEKit/engines/psi4.py
+++ /dev/null
@@ -1,308 +0,0 @@
-#!/usr/bin/env python3
-
-import subprocess as sp
-
-import numpy as np
-
-from QUBEKit.engines.base_engine import Engines
-from QUBEKit.utils.exceptions import PSI4Error
-from QUBEKit.utils.helpers import append_to_log
-
-
-class PSI4(Engines):
-    """
-    Writes and executes input files for psi4.
-    Also used to extract Hessian matrices; optimised structures; frequencies; etc.
-    """
-
-    def __init__(self, molecule):
-
-        super().__init__(molecule)
-
-        self.functional_dict = {"pbepbe": "PBE", "wb97xd": "wB97X-D"}
-        # Search for functional in dict, if it's not there, just leave the theory as it is.
-        self.molecule.theory = self.functional_dict.get(
-            self.molecule.theory.lower(), self.molecule.theory
-        )
-
-        # Test if PSI4 is callable
-        try:
-            sp.run("psi4 -h", shell=True, check=True, stdout=sp.PIPE)
-        except sp.CalledProcessError as exc:
-            raise ModuleNotFoundError(
-                "PSI4 not working. Please ensure PSI4 is installed and can be called with the command: psi4"
-            ) from exc
-
-        if self.molecule.geometric:
-            try:
-                sp.run("geometric-optimize -h", shell=True, check=True, stdout=sp.PIPE)
-            except sp.CalledProcessError as exc:
-                raise ModuleNotFoundError(
-                    "Geometric not working. Please ensure geometric is installed and can be called "
-                    "with the command: geometric-optimize"
-                ) from exc
-
-    # TODO add restart from log method
-    def generate_input(
-        self,
-        optimise=False,
-        hessian=False,
-        density=False,
-        energy=False,
-        fchk=False,
-        restart=False,
-        execute=True,
-    ):
-        """
-        Converts to psi4 input format to be run in psi4 without using geometric.
-        :param optimise: Optimise the molecule to the desired convergence criterion within the iteration limit
-        :param hessian: Calculate the hessian matrix
-        :param density: Calculate the electron density
-        :param energy: Calculate the single point energy of the molecule
-        :param fchk: Write out a gaussian style Fchk file
-        :param restart: Restart the calculation from a log point (required but unused to match g09's generate_input())
-        :param execute: Run the desired Psi4 job
-        :return: The completion status of the job True if successful False if not run or failed
-        """
-
-        setters = ""
-        tasks = ""
-
-        if energy:
-            append_to_log(self.molecule.home, "Writing psi4 energy calculation input")
-            tasks += f"\nenergy('{self.molecule.theory}')"
-
-        if optimise:
-            append_to_log(
-                self.molecule.home, "Writing PSI4 optimisation input", "minor"
-            )
-            setters += f" g_convergence {self.molecule.convergence}\n GEOM_MAXITER {self.molecule.iterations}\n"
-            tasks += f"\noptimize('{self.molecule.theory.lower()}')"
-
-        if hessian:
-            append_to_log(
-                self.molecule.home,
-                "Writing PSI4 Hessian matrix calculation input",
-            )
-            setters += " hessian_write on\n"
-
-            tasks += f"\nenergy, wfn = frequency('{self.molecule.theory.lower()}', return_wfn=True)"
-
-            tasks += "\nwfn.hessian().print_out()\n\n"
-
-        if density:
-            raise NotImplementedError(
-                "Due to PSI4 requiring a box size which cannot be automatically generated, "
-                "PSI4 cannot currently be used for density calculations. Please use Gaussian "
-                "instead."
-            )
-
-        if fchk:
-            append_to_log(
-                self.molecule.home, "Writing PSI4 input file to generate fchk file"
-            )
-            tasks += f"\ngrad, wfn = gradient('{self.molecule.theory.lower()}', return_wfn=True)"
-            tasks += "\nfchk_writer = psi4.core.FCHKWriter(wfn)"
-            tasks += f'\nfchk_writer.write("{self.molecule.name}_psi4.fchk")\n'
-
-        setters += "}\n"
-
-        if not execute:
-            setters += f"set_num_threads({self.molecule.threads})\n"
-
-        # input.dat is the PSI4 input file.
-        with open("input.dat", "w+") as input_file:
-            # opening tag is always writen
-            input_file.write(
-                f"memory {self.molecule.memory} GB\n\nmolecule {self.molecule.name} {{\n"
-                f"{self.molecule.charge} {self.molecule.multiplicity} \n"
-            )
-            # molecule is always printed
-            for i, atom in enumerate(self.molecule.coordinates):
-                input_file.write(
-                    f" {self.molecule.atoms[i].atomic_symbol}    "
-                    f"{float(atom[0]): .10f}  {float(atom[1]): .10f}  {float(atom[2]): .10f} \n"
-                )
-
-            input_file.write(
-                f" units angstrom\n no_reorient\n}}\n\nset {{\n basis {self.molecule.basis}\n"
-            )
-
-            input_file.write(setters)
-            input_file.write(tasks)
-
-        if execute:
-            with open("log.txt", "w+") as log:
-                try:
-                    sp.run(
-                        f"psi4 input.dat -n {self.molecule.threads}",
-                        shell=True,
-                        stdout=log,
-                        stderr=log,
-                        check=True,
-                    )
-                except sp.CalledProcessError as exc:
-                    raise PSI4Error(
-                        "PSI4 did not execute successfully check log file for details."
-                    ) from exc
-
-            # Now check the exit status of the job
-            return self.check_for_errors()
-
-        else:
-            return {"success": False, "error": "Not run"}
-
-    def check_for_errors(self):
-        """
-        Read the output file from the job and check for normal termination and any errors
-        :return: A dictionary of the success status and any errors.
-        """
-
-        with open("output.dat", "r") as log:
-            for line in log:
-                if "*** Psi4 exiting successfully." in line:
-                    return {"success": True}
-
-                elif "*** Psi4 encountered an error." in line:
-                    return {"success": False, "error": "Not known"}
-
-            return {"success": False, "error": "Segfault"}
-
-    def optimised_structure(self):
-        """
-        Parses the final optimised structure from the output.dat file (from psi4) to a numpy array.
-        Also returns the energy of the optimized structure.
-        """
-
-        # Run through the file and find all lines containing '==> Geometry', add these lines to a list.
-        # Reverse the list
-        # from the start of this list, jump down to the first atom and set this as the start point
-        # Split the row into 4 columns: centre, x, y, z.
-        # Add each row to a matrix.
-        # Return the matrix.
-
-        # output.dat is the psi4 output file.
-        with open("output.dat", "r") as file:
-            lines = file.readlines()
-            # Will contain index of all the lines containing '==> Geometry'.
-            geo_pos_list = []
-            for count, line in enumerate(lines):
-                if "==> Geometry" in line:
-                    geo_pos_list.append(count)
-
-                elif "**** Optimization is complete!" in line:
-                    opt_pos = count
-                    opt_steps = int(line.split()[5])
-
-            if not (opt_pos and opt_steps):
-                raise EOFError(
-                    "According to the output.dat file, optimisation has not completed."
-                )
-
-            # now get the final opt_energy
-            opt_energy = float(lines[opt_pos + opt_steps + 7].split()[1])
-
-            # Set the start as the last instance of '==> Geometry'.
-            start_of_vals = geo_pos_list[-1] + 9
-
-            opt_struct = []
-
-            for row in range(len(self.molecule.atoms)):
-
-                # Append the first 4 columns of each row, converting to float as necessary.
-                struct_row = []
-                for indx in range(3):
-                    struct_row.append(
-                        float(lines[start_of_vals + row].split()[indx + 1])
-                    )
-
-                opt_struct.append(struct_row)
-
-        return np.array(opt_struct), opt_energy
-
-    @staticmethod
-    def get_energy():
-        """Get the energy of a single point calculation."""
-
-        # open the psi4 log file
-        with open("output.dat", "r") as log:
-            for line in log:
-                if "Total Energy =" in line:
-                    return float(line.split()[3])
-
-        raise EOFError("Cannot find energy in output.dat file.")
-
-    def all_modes(self):
-        """Extract all modes from the psi4 output file."""
-
-        # Find "post-proj  all modes"
-        # Jump to first value, ignoring text.
-        # Move through data, adding it to a list
-        # continue onto next line.
-        # Repeat until the following line is known to be empty.
-
-        # output.dat is the psi4 output file.
-        with open("output.dat", "r") as file:
-            lines = file.readlines()
-            for count, line in enumerate(lines):
-                if "post-proj  all modes" in line:
-                    start_of_vals = count
-                    break
-            else:
-                raise EOFError("Cannot locate modes in output.dat file.")
-
-            # Barring the first (and sometimes last) line, dat file has 6 values per row.
-            end_of_vals = start_of_vals + (3 * len(self.molecule.atoms)) // 6
-
-            structures = lines[start_of_vals][24:].replace("'", "").split()
-            structures = structures[6:]
-
-            for row in range(1, end_of_vals - start_of_vals):
-                # Remove double strings and weird formatting.
-                structures += (
-                    lines[start_of_vals + row].replace("'", "").replace("]", "").split()
-                )
-
-            all_modes = [float(val) for val in structures]
-
-            return np.array(all_modes)
-
-    def geo_gradient(self, threads=False, execute=True):
-        """
-        Write the psi4 style input file to get the gradient for geometric
-        and run geometric optimisation.
-        """
-
-        molecule = self.molecule.coordinates
-
-        with open(f"{self.molecule.name}.psi4in", "w+") as file:
-
-            file.write(
-                f"memory {self.molecule.memory} GB\n\nmolecule {self.molecule.name} {{\n {self.molecule.charge} "
-                f"{self.molecule.multiplicity} \n"
-            )
-
-            for i, atom in enumerate(molecule):
-                file.write(
-                    f"  {self.molecule.atoms[i].atomic_symbol:2}    "
-                    f"{float(atom[0]): .10f}  {float(atom[1]): .10f}  {float(atom[2]): .10f}\n"
-                )
-
-            file.write(
-                f" units angstrom\n no_reorient\n}}\nset basis {self.molecule.basis}\n"
-            )
-
-            if threads:
-                file.write(f"set_num_threads({self.molecule.threads})")
-
-            file.write(f"\n\ngradient('{self.molecule.theory}')\n")
-
-        if execute:
-            with open("log.txt", "w+") as log:
-                sp.run(
-                    f"geometric-optimize --psi4 {self.molecule.name}.psi4in {self.molecule.constraints_file} "
-                    f"--nt {self.molecule.threads}",
-                    shell=True,
-                    stdout=log,
-                    stderr=log,
-                )
diff --git a/QUBEKit/engines/qcengine.py b/QUBEKit/engines/qcengine.py
index 68e7b38c..86cf4c3b 100644
--- a/QUBEKit/engines/qcengine.py
+++ b/QUBEKit/engines/qcengine.py
@@ -1,113 +1,50 @@
-#!/usr/bin/env python3
-
-import numpy as np
 import qcelemental as qcel
 import qcengine as qcng
+from pydantic import Field
+from typing_extensions import Literal
 
-from QUBEKit.engines.base_engine import Engines
-from QUBEKit.utils import constants
-from QUBEKit.utils.helpers import check_symmetry
-
-
-class QCEngine(Engines):
-    def __init__(self, molecule):
-
-        super().__init__(molecule)
+from QUBEKit.engines.base_engine import BaseEngine
+from QUBEKit.molecules import Ligand
 
-    def generate_qschema(self):
-        """
-        Using the molecule object, generate a QCEngine schema. This can then
-        be fed into the various QCEngine procedures.
-        :param input_type: The part of the molecule object that should be used when making the schema
-        :return: The qcelemental qschema
-        """
 
-        mol_data = f"{self.molecule.charge} {self.molecule.multiplicity}\n"
+class QCEngine(BaseEngine):
+    """
+    A wrapper around qcengine which calls single point calculations, this automatically converts the ligand and inputs
+    to the correct format.
+    TODO how do we want to support keywords.
+    """
 
-        for i, coord in enumerate(self.molecule.coordinates):
-            mol_data += f"{self.molecule.atoms[i].atomic_symbol} "
-            for item in coord:
-                mol_data += f"{item} "
-            mol_data += "\n"
+    driver: Literal["energy", "gradient", "hessian"] = Field(
+        "hessian", description="The single point job type that should be done."
+    )
 
-        return qcel.models.Molecule.from_data(mol_data)
-
-    def call_qcengine(self, engine, driver):
+    def call_qcengine(
+        self,
+        molecule: Ligand,
+    ) -> qcel.models.AtomicResult:
         """
-        Using the created schema, run a particular engine, specifying the driver (job type).
-        e.g. engine: geo, driver: energies.
-        :param engine: The engine to be used psi4 geometric
-        :param driver: The calculation type to be done e.g. energy, gradient, hessian, properties
-        :return: The required driver information
+        Calculate the requested property using qcengine for the given molecule.
+
+        Parameters
+        ----------
+        molecule:
+            The QUBEKit ligand that the calculation should be ran on.
+        input_type:
+            The coords the calculation should be ran on.
+
+        Returns
+        -------
+        qcel.models.AtomicResult
+            The full qcelemental atomic result so any required information can be extracted.
         """
+        qc_mol = molecule.to_qcschema()
+        task = qcel.models.AtomicInput(
+            molecule=qc_mol,
+            driver=self.driver,
+            model=self.qc_model,
+            keywords={"scf_type": "df"},
+        )
 
-        mol = self.generate_qschema()
-        options = {"memory": self.molecule.memory, "ncores": self.molecule.threads}
-
-        # Call psi4 for energy, gradient, hessian or property calculations
-        if engine == "psi4":
-            psi4_task = qcel.models.AtomicInput(
-                molecule=mol,
-                driver=driver,
-                model={"method": self.molecule.theory, "basis": self.molecule.basis},
-                keywords={"scf_type": "df"},
-            )
-
-            ret = qcng.compute(psi4_task, "psi4", local_options=options)
-
-            if driver == "hessian":
-                hess_size = 3 * len(self.molecule.atoms)
-                conversion = constants.HA_TO_KCAL_P_MOL / (constants.BOHR_TO_ANGS ** 2)
-                hessian = (
-                    np.reshape(ret.return_result, (hess_size, hess_size)) * conversion
-                )
-                check_symmetry(hessian)
-
-                return hessian
-
-            else:
-                return ret.return_result
-
-        # Call geometric with psi4 to optimise a molecule
-        elif engine == "geometric":
-            geo_task = {
-                "schema_name": "qcschema_optimization_input",
-                "schema_version": 1,
-                "keywords": {
-                    "coordsys": "dlc",
-                    "maxiter": self.molecule.iterations,
-                    "program": "psi4",
-                    "convergence_set": self.molecule.convergence,
-                },
-                "input_specification": {
-                    "schema_name": "qcschema_input",
-                    "schema_version": 1,
-                    "driver": "gradient",
-                    "model": {
-                        "method": self.molecule.theory,
-                        "basis": self.molecule.basis,
-                    },
-                    "keywords": {},
-                },
-                "initial_molecule": mol,
-            }
-            return qcng.compute_procedure(
-                geo_task, "geometric", return_dict=True, local_options=options
-            )
-
-        elif engine == "torchani":
-            ani_task = qcel.models.AtomicInput(
-                molecule=mol,
-                driver=driver,
-                model={"method": "ANI1x", "basis": None},
-                keywords={"scf_type": "df"},
-            )
-
-            return qcng.compute(
-                ani_task, "torchani", local_options=options
-            ).return_result
+        result = qcng.compute(task, self.program, local_options=self.local_options)
 
-        else:
-            raise KeyError(
-                'Invalid engine type provided. Please use "geometric", "psi4" or "torchani".'
-            )
+        return result
diff --git a/QUBEKit/lennard_jones.py b/QUBEKit/lennard_jones.py
index e435d974..6e4553a9 100644
--- a/QUBEKit/lennard_jones.py
+++ b/QUBEKit/lennard_jones.py
@@ -4,7 +4,7 @@
 from collections import namedtuple
 from typing import Dict, List, Optional
 
-from QUBEKit.ligand import Ligand
+from QUBEKit.molecules import Ligand
 from QUBEKit.utils import constants
 
 
diff --git a/QUBEKit/mod_seminario.py b/QUBEKit/mod_seminario.py
index 3cb094f5..bec304b6 100644
--- a/QUBEKit/mod_seminario.py
+++ b/QUBEKit/mod_seminario.py
@@ -11,7 +11,7 @@
 
 import numpy as np
 
-from QUBEKit.ligand import Ligand
+from QUBEKit.molecules import Ligand
 from QUBEKit.utils import constants
 
 
@@ -199,6 +199,9 @@ def __init__(self, molecule: Ligand):
         # reset the ligand data
         self.molecule.HarmonicBondForce = {}
         self.molecule.HarmonicAngleForce = {}
+        # convert the hessian from atomic units
+        conversion = constants.HA_TO_KCAL_P_MOL / (constants.BOHR_TO_ANGS ** 2)
+        self.hessian = molecule.hessian * conversion
 
     def __repr__(self):
         return f"{self.__class__.__name__}({self.__dict__!r})"
diff --git a/QUBEKit/molecules/__init__.py b/QUBEKit/molecules/__init__.py
new file mode 100644
index 00000000..a294694e
--- /dev/null
+++ b/QUBEKit/molecules/__init__.py
@@ -0,0 +1,3 @@
+from QUBEKit.molecules.components import Atom, Bond, Element, ExtraSite
+from QUBEKit.molecules.ligand import Ligand
+from QUBEKit.molecules.protein import Protein
diff --git a/QUBEKit/molecules/components.py b/QUBEKit/molecules/components.py
new file mode 100644
index 00000000..51a64c2e
--- /dev/null
+++ b/QUBEKit/molecules/components.py
@@ -0,0 +1,267 @@
+from enum import Enum
+from typing import List, Optional
+
+from pydantic import BaseModel, Field, PositiveInt
+from rdkit import Chem
+from rdkit.Chem.rdchem import GetPeriodicTable, PeriodicTable
+
+
+class Element:
+    """
+    Simple wrapper class for getting element info using RDKit.
+    """
+
+    @staticmethod
+    def p_table() -> PeriodicTable:
+        return GetPeriodicTable()
+
+    @staticmethod
+    def mass(identifier):
+        pt = Element.p_table()
+        return pt.GetAtomicWeight(identifier)
+
+    @staticmethod
+    def number(identifier):
+        pt = Element.p_table()
+        return pt.GetAtomicNumber(identifier)
+
+    @staticmethod
+    def name(identifier):
+        pt = Element.p_table()
+        return pt.GetElementSymbol(identifier)
+
+
+class AtomStereoChemistry(str, Enum):
+    """
+    Atom stereochemistry types.
+    """
+
+    R = "R"
+    S = "S"
+    U = "Unknown"
+
+
+class BondStereoChemistry(str, Enum):
+    """
+    Bond stereochemistry types.
+    """
+
+    E = "E"
+    Z = "Z"
+    U = "Unknown"
+
+
+class Atom(BaseModel):
+    """
+    Class to hold all of the "per atom" information.
+    All atoms in Molecule will have an instance of this Atom class to describe their properties.
+    """
+
+    class Config:
+        validate_assignment = True
+        json_encoders = {Enum: lambda v: v.value}
+
+    atomic_number: PositiveInt = Field(
+        ...,
+        description="The atomic number of the atom all other properties are based on this number.",
+    )
+    atom_index: int = Field(
+        ..., description="The index this atom has in the molecule object", gt=-1
+    )
+    atom_name: Optional[str] = Field(
+        None,
+        description="An optional unqiue atom name that should be assigned to the atom, the ligand object will make sure all atoms have unique names.",
+    )
+    formal_charge: int = Field(
+        ...,
+        description="The formal charge of the atom, used to calculate the molecule total charge",
+    )
+    aromatic: bool = Field(
+        ...,
+        description="If the atom should be considered aromatic `True` or not `False`.",
+    )
+    stereochemistry: Optional[AtomStereoChemistry] = Field(
+        None,
+        description="The stereochemistry of the atom where None means not stereogenic and U is unknown or ambiguous.",
+    )
+    bonds: Optional[List[int]] = Field(
+        None, description="The list of atom indices which are bonded to this atom."
+    )
+
+    @classmethod
+    def from_rdkit(cls, rd_atom: Chem.Atom) -> "Atom":
+        """
+        Build a QUBEKit atom from an rdkit atom instance.
+        """
+        atomic_number = rd_atom.GetAtomicNum()
+        index = rd_atom.GetIdx()
+        formal_charge = rd_atom.GetFormalCharge()
+        aromatic = rd_atom.GetIsAromatic()
+        bonds = [a.GetIdx() for a in rd_atom.GetNeighbors()]
+        # check for names in the normal places pdb, mol2 and mol
+        if rd_atom.HasProp("_Name"):
+            name = rd_atom.GetProp("_Name")
+        elif rd_atom.HasProp("_TriposAtomName"):
+            name = rd_atom.GetProp("_TriposAtomName")
+        else:
+            try:
+                name = rd_atom.GetMonomerInfo().GetName().strip()
+            except AttributeError:
+                name = None
+        # stereochem
+        if rd_atom.HasProp("_CIPCode"):
+            stereo_code = rd_atom.GetProp("_CIPCode")
+        else:
+            stereo_code = None
+        return cls(
+            atomic_number=atomic_number,
+            atom_index=index,
+            atom_name=name,
+            formal_charge=formal_charge,
+            aromatic=aromatic,
+            stereochemistry=stereo_code,
+            bonds=bonds,
+        )
+
+    @property
+    def atomic_mass(self) -> float:
+        """Convert the atomic number to mass."""
+        return Element.mass(self.atomic_number)
+
+    @property
+    def atomic_symbol(self) -> str:
+        """Convert the atomic number to the atomic symbol as per the periodic table."""
+        return Element.name(self.atomic_number).title()
+
+    def to_rdkit(self) -> Chem.Atom:
+        """
+        Convert the QUBEKit atom an RDKit atom.
+        """
+        # build the atom from atomic number
+        rd_atom = Chem.Atom(self.atomic_number)
+        rd_atom.SetFormalCharge(self.formal_charge)
+        rd_atom.SetIsAromatic(self.aromatic)
+        rd_atom.SetProp("_Name", self.atom_name)
+        # left is counter clockwise
+        if self.stereochemistry == "S":
+            rd_atom.SetChiralTag(Chem.CHI_TETRAHEDRAL_CCW)
+        # right is clockwise
+        elif self.stereochemistry == "R":
+            rd_atom.SetChiralTag(Chem.CHI_TETRAHEDRAL_CW)
+
+        return rd_atom
+
+    # def __repr__(self):
+    #     return f"{self.__class__.__name__}({self.__dict__!r})"
+    #
+    # def __str__(self):
+    #     """
+    #     Prints the Atom class objects' names and values one after another with new lines between each.
+    #     """
+    #
+    #     return_str = ""
+    #     for key, val in self.__dict__.items():
+    #         # Return all objects as {atom object name} = {atom object value(s)}.
+    #         return_str += f"\n{key} = {val}\n"
+    #
+    #     return return_str
+
+
+class Bond(BaseModel):
+    """
+    A basic bond class.
+    """
+
+    class Config:
+        validate_assignment = True
+        json_encoders = {Enum: lambda v: v.value}
+
+    atom1_index: int = Field(
+        ..., description="The index of the first atom in the bond."
+    )
+    atom2_index: int = Field(
+        ..., description="The index of the second atom in the bond."
+    )
+    bond_order: float = Field(..., description="The float value of the bond order.")
+    aromatic: bool = Field(
+        ..., description="If the bond should be considered aromatic."
+    )
+    stereochemistry: Optional[BondStereoChemistry] = Field(
+        None,
+        description="The stereochemistry of the bond, where None means not stereogenic.",
+    )
+
+    @classmethod
+    def from_rdkit(cls, rd_bond: Chem.Bond) -> "Bond":
+        """
+        Build a QUBEKit bond class from an rdkit reference.
+        """
+        atom1_index = rd_bond.GetBeginAtomIdx()
+        atom2_index = rd_bond.GetEndAtomIdx()
+        aromatic = rd_bond.GetIsAromatic()
+        order = rd_bond.GetBondTypeAsDouble()
+        stereo_tag = rd_bond.GetStereo()
+        if stereo_tag == Chem.BondStereo.STEREOZ:
+            stereo = "Z"
+        elif stereo_tag == Chem.BondStereo.STEREOE:
+            stereo = "E"
+        else:
+            stereo = None
+        return cls(
+            atom1_index=atom1_index,
+            atom2_index=atom2_index,
+            aromatic=aromatic,
+            bond_order=order,
+            stereochemistry=stereo,
+        )
+
+    @property
+    def rdkit_type(self) -> Chem.BondType:
+        """
+        Convert the bond order float to a bond type.
+        """
+        conversion = {
+            1: Chem.BondType.SINGLE,
+            1.5: Chem.BondType.AROMATIC,
+            2: Chem.BondType.DOUBLE,
+            3: Chem.BondType.TRIPLE,
+            4: Chem.BondType.QUADRUPLE,
+            5: Chem.BondType.QUINTUPLE,
+            6: Chem.BondType.HEXTUPLE,
+            7: Chem.BondType.ONEANDAHALF,
+        }
+        return conversion[self.bond_order]
+
+    @property
+    def rdkit_stereo(self) -> Optional[Chem.BondStereo]:
+        """
+        Return the rdkit style stereo enum.
+        """
+        if self.stereochemistry == "E":
+            return Chem.BondStereo.STEREOE
+        elif self.stereochemistry == "Z":
+            return Chem.BondStereo.STEREOZ
+        return None
+
+
+class ExtraSite:
+    """
+    Used to store extra sites for xml writer in ligand.
+    This class is used by both internal v-sites fitting and the ONETEP reader.
+    """
+
+    def __init__(self):
+        self.parent_index: Optional[int] = None
+        self.closest_a_index: Optional[int] = None
+        self.closest_b_index: Optional[int] = None
+        # Optional: Used for Nitrogen only.
+        self.closest_c_index: Optional[int] = None
+
+        self.o_weights: Optional[List[float]] = None
+        self.x_weights: Optional[List[float]] = None
+        self.y_weights: Optional[List[float]] = None
+
+        self.p1: Optional[float] = None
+        self.p2: Optional[float] = None
+        self.p3: Optional[float] = None
+        self.charge: Optional[float] = None
diff --git a/QUBEKit/ligand.py b/QUBEKit/molecules/ligand.py
similarity index 89%
rename from QUBEKit/ligand.py
rename to QUBEKit/molecules/ligand.py
index dfdb402d..ac1e0cad 100644
--- a/QUBEKit/ligand.py
+++ b/QUBEKit/molecules/ligand.py
@@ -27,22 +27,28 @@
 import pickle
 import xml.etree.ElementTree as ET
 from collections import OrderedDict
-from datetime import datetime
 from typing import Any, Dict, List, Optional, Set, Tuple
 from xml.dom.minidom import parseString
 
 import networkx as nx
 import numpy as np
+import qcelemental as qcel
 from rdkit import Chem
 from simtk import unit
 from simtk.openmm.app import Aromatic, Double, Single, Topology, Triple
 from simtk.openmm.app.element import Element
 
 import QUBEKit
-from QUBEKit.engines import RDKit
-from QUBEKit.utils.datastructures import Atom, Bond, ExtraSite
-from QUBEKit.utils.exceptions import FileTypeError, TopologyMismatch
-from QUBEKit.utils.file_handling import ReadInput, ReadInputProtein
+from QUBEKit.molecules.components import Atom, Bond, ExtraSite
+from QUBEKit.molecules.utils import RDKit, ReadInput
+from QUBEKit.utils import constants
+from QUBEKit.utils.exceptions import (
+    ConformerError,
+    FileTypeError,
+    StereoChemistryError,
+    TopologyMismatch,
+)
+from QUBEKit.utils.helpers import _assert_wrapper
 
 
 class DefaultsMixin:
@@ -98,7 +104,7 @@ def __init__(self, *args, **kwargs):
         self.div_index = 0
         self.parameter_engine = "antechamber"
         self.l_pen = 0.0
-        self.mm_opt_method = "openmm"
+        self.pre_opt_method: str = "rdkit_uff"
         self.relative_to_global = False
 
         self.excited_state = False
@@ -370,24 +376,6 @@ def get_bond_between(self, atom1_index: int, atom2_index: int) -> Bond:
             f"There is no bond between atoms {atom1_index} and {atom2_index} in this molecule."
         )
 
-    # def read_geometric_traj(self, trajectory):
-    #     """
-    #     Read in the molecule coordinates to the traj holder from a geometric optimisation using qcengine.
-    #     :param trajectory: The qcengine trajectory
-    #
-    #     TODO Move to QCEngine()
-    #     """
-    #
-    #     for frame in trajectory:
-    #         opt_traj = []
-    #         # Convert coordinates from bohr to angstroms
-    #         geometry = np.array(frame["molecule"]["geometry"]) * constants.BOHR_TO_ANGS
-    #         for i, atom in enumerate(frame["molecule"]["symbols"]):
-    #             opt_traj.append(
-    #                 [geometry[0 + i * 3], geometry[1 + i * 3], geometry[2 + i * 3]]
-    #             )
-    #         self.coords["traj"].append(np.array(opt_traj))
-
     @property
     def has_unique_atom_names(self) -> bool:
         """
@@ -976,11 +964,14 @@ def to_rdkit(self) -> Chem.Mol:
         """
         Generate an rdkit representation of the QUBEKit ligand object.
 
-        #TODO what properties should be put in the rdkit molecule? Multiplicity?
+        Here we build the molecule and assign the stereochemistry using the coordinates as we should always have a set of coordinates in the model.
+        This allows us to skip complicated local vs global stereo chemistry checks however this could break in future.
 
         Returns:
             An rdkit representation of the molecule.
         """
+        # TODO what properties should be put in the rdkit molecule? Multiplicity?
+
         # make an editable molecule
         rd_mol = Chem.RWMol()
         if self.name is not None:
@@ -1008,33 +999,39 @@ def to_rdkit(self) -> Chem.Mol:
         # must use openff MDL model for compatibility
         Chem.SetAromaticity(rd_mol, Chem.AromaticityModel.AROMATICITY_MDL)
 
-        Chem.AssignStereochemistry(
-            rd_mol, force=True, cleanIt=True, flagPossibleStereoCenters=True
+        # conformers
+        rd_mol = RDKit.add_conformer(
+            rdkit_mol=rd_mol, conformer_coordinates=self.coordinates
         )
+        Chem.AssignStereochemistryFrom3D(rd_mol)
 
         # now we should check that the stereo has not been broken
         for rd_atom in rd_mol.GetAtoms():
             index = rd_atom.GetIdx()
             qb_atom = self.atoms[index]
             if qb_atom.stereochemistry is not None:
-                assert qb_atom.stereochemistry == rd_atom.GetProp("_CIPCode")
+                with _assert_wrapper(StereoChemistryError):
+                    assert qb_atom.stereochemistry == rd_atom.GetProp(
+                        "_CIPCode"
+                    ), f"StereoChemistry incorrect expected {qb_atom.stereochemistry} got {rd_atom.GetProp('_CIPCode')} for atom {qb_atom}"
 
         for rd_bond in rd_mol.GetBonds():
             index = rd_bond.GetIdx()
             qb_bond = self.bonds[index]
-            assert rd_bond.GetBondTypeAsDouble() == qb_bond.bond_order
             if qb_bond.stereochemistry is not None:
                 rd_bond.SetStereo(qb_bond.rdkit_stereo)
             rd_stereo = rd_bond.GetStereo()
             if qb_bond.stereochemistry == "E":
-                assert rd_stereo == Chem.BondStereo.STEREOE
+                with _assert_wrapper(StereoChemistryError):
+                    assert (
+                        rd_stereo == Chem.BondStereo.STEREOE
+                    ), f"StereoChemistry incorrect expected E got {rd_stereo}"
             elif qb_bond.stereochemistry == "Z":
-                assert rd_stereo == Chem.BondStereo.STEREOZ
+                with _assert_wrapper(StereoChemistryError):
+                    assert (
+                        rd_stereo == Chem.BondStereo.STEREOZ
+                    ), f"StereoChemistry incorrect expected Z got {rd_stereo}"
 
-        # conformers
-        rd_mol = RDKit.add_conformer(
-            rdkit_mol=rd_mol, conformer_coordinates=self.coordinates
-        )
         return Chem.Mol(rd_mol)
 
     def get_smarts_matches(self, smirks: str) -> Optional[List[Tuple[int, ...]]]:
@@ -1411,6 +1408,33 @@ def _validate_atom_names(self) -> None:
         if not self.has_unique_atom_names:
             self.generate_atom_names()
 
+    def to_qcschema(self, extras: Optional[Dict] = None) -> qcel.models.Molecule:
+        """
+        build a qcschema molecule from the ligand object, this is useful to interface with QCEngine and QCArchive.
+        """
+        # make sure we have a conformer
+        coords = self.coordinates
+        if coords == [] or coords is None:
+            raise ConformerError(
+                "The molecule must have a conformation to make a qcschema molecule."
+            )
+        # input must be in bohr
+        coords *= constants.ANGS_TO_BOHR
+        # we do not store explicit bond order so guess at 1
+        bonds = [
+            (bond.atom1_index, bond.atom2_index, bond.bond_order) for bond in self.bonds
+        ]
+        symbols = [atom.atomic_symbol for atom in self.atoms]
+        schema_info = {
+            "symbols": symbols,
+            "geometry": coords,
+            "connectivity": bonds,
+            "molecular_charge": self.charge,
+            "molecular_multiplicity": self.multiplicity,
+            "extras": extras,
+        }
+        return qcel.models.Molecule.from_data(schema_info, validate=True)
+
     def add_conformer(self, file_name: str) -> None:
         """
         Read the given input file extract  the conformers and save them to the ligand.
@@ -1426,126 +1450,3 @@ def add_conformer(self, file_name: str) -> None:
             else:
                 coords = input_data.coords
         self.coordinates = coords
-
-
-class Protein(DefaultsMixin, Molecule):
-    """
-    This class handles the protein input to make the QUBEKit xml files and rewrite the pdb so we can use it.
-    """
-
-    def __init__(self, mol_input, name=None):
-        """
-        is_protein      Bool; True for Protein class
-        home            Current working directory (location for QUBEKit execution).
-        residues        List of all residues in the molecule in order e.g. ['ARG', 'HIS', ... ]
-        Residues        List of residue names for each atom e.g. ['ARG', 'ARG', 'ARG', ... 'HIS', 'HIS', ... ]
-        pdb_names       List
-        """
-
-        super().__init__(mol_input, name)
-
-        self.home = os.getcwd()
-        self.residues = None
-        self.Residues = None
-        self.pdb_names = None
-
-        self.combination = "opls"
-
-        if not isinstance(mol_input, ReadInputProtein):
-            self._check_file_type(file_name=mol_input)
-            input_data = ReadInputProtein.from_pdb(file_name=mol_input)
-        else:
-            input_data = mol_input
-        self._save_to_protein(input_data)
-
-    @classmethod
-    def from_file(cls, file_name: str, name: Optional[str] = None) -> "Protein":
-        """
-        Instance the protein class from a pdb file.
-        """
-        cls._check_file_type(file_name=file_name)
-        input_data = ReadInputProtein.from_pdb(file_name=file_name, name=name)
-        return cls(input_data)
-
-    @staticmethod
-    def _check_file_type(file_name: str) -> None:
-        """
-        Make sure the protien is being read from a pdb file.
-        """
-        if ".pdb" not in file_name:
-            raise FileTypeError("Proteins can only be read from pdb.")
-
-    def _save_to_protein(self, mol_input: ReadInputProtein):
-        """
-        Public access to private file_handlers.py file.
-        Users shouldn't ever need to interface with file_handlers.py directly.
-        All parameters will be set from a file (or other input) via this public method.
-            * Don't bother updating name, topology or atoms if they are already stored.
-            * Do bother updating coords, rdkit_mol, residues, Residues, pdb_names
-        """
-
-        if mol_input.name is not None:
-            self.name = mol_input.name
-        if mol_input.atoms is not None:
-            self.atoms = mol_input.atoms
-        if mol_input.coords is not None:
-            self.coordinates = mol_input.coords
-        if mol_input.residues is not None:
-            self.residues = mol_input.residues
-        if mol_input.pdb_names is not None:
-            self.pdb_names = mol_input.pdb_names
-        if mol_input.bonds is not None:
-            self.bonds = mol_input.bonds
-
-        if not self.bonds:
-            print(
-                "No connections found in pdb file; topology will be inferred by OpenMM."
-            )
-            return
-
-        self.symmetrise_from_topology()
-
-    def write_pdb(self, name=None):
-        """
-        This method replaces the ligand method as all of the atom names and residue names have to be replaced.
-        """
-
-        with open(f"{name if name is not None else self.name}.pdb", "w+") as pdb_file:
-
-            pdb_file.write(f"REMARK   1 CREATED WITH QUBEKit {datetime.now()}\n")
-            # Write out the atomic xyz coordinates
-            for i, (coord, atom) in enumerate(zip(self.coordinates, self.atoms)):
-                x, y, z = coord
-                # May cause issues if protein contains more than 10,000 atoms.
-                pdb_file.write(
-                    f"HETATM {i+1:>4}{atom.atom_name:>5} QUP     1{x:12.3f}{y:8.3f}{z:8.3f}"
-                    f"  1.00  0.00         {atom.atomic_symbol.upper():>3}\n"
-                )
-
-            # Add the connection terms based on the molecule topology.
-            for bond in self.bonds:
-                pdb_file.write(
-                    f"CONECT{bond.atom1_index + 1:5} {bond.atom2_index + 1:5}\n"
-                )
-
-            pdb_file.write("END\n")
-
-    def update(self):
-        """
-        After the protein has been passed to the parametrisation class we get back the bond info
-        use this to update all missing terms.
-        """
-        if not self.bonds:
-            self.bonds = []
-            # using the new harmonic bond force dict we can add the bond edges to the topology graph
-            for bond in self.HarmonicBondForce:
-                self.bonds.append(
-                    Bond(
-                        atom1_indx=bond[0],
-                        atom2_index=bond[1],
-                        bond_order=1,
-                        aromatic=False,
-                    )
-                )
-
-        self.symmetrise_from_topology()
diff --git a/QUBEKit/molecules/protein.py b/QUBEKit/molecules/protein.py
new file mode 100644
index 00000000..d366d38c
--- /dev/null
+++ b/QUBEKit/molecules/protein.py
@@ -0,0 +1,131 @@
+import os
+from datetime import datetime
+from typing import Optional
+
+from QUBEKit.molecules.components import Bond
+from QUBEKit.molecules.ligand import DefaultsMixin, Molecule
+from QUBEKit.molecules.utils import ReadInputProtein
+from QUBEKit.utils.exceptions import FileTypeError
+
+
+class Protein(DefaultsMixin, Molecule):
+    """
+    This class handles the protein input to make the QUBEKit xml files and rewrite the pdb so we can use it.
+    """
+
+    def __init__(self, mol_input, name=None):
+        """
+        is_protein      Bool; True for Protein class
+        home            Current working directory (location for QUBEKit execution).
+        residues        List of all residues in the molecule in order e.g. ['ARG', 'HIS', ... ]
+        Residues        List of residue names for each atom e.g. ['ARG', 'ARG', 'ARG', ... 'HIS', 'HIS', ... ]
+        pdb_names       List
+        """
+
+        super().__init__(mol_input, name)
+
+        self.home = os.getcwd()
+        self.residues = None
+        self.Residues = None
+        self.pdb_names = None
+
+        self.combination = "opls"
+
+        if not isinstance(mol_input, ReadInputProtein):
+            self._check_file_type(file_name=mol_input)
+            input_data = ReadInputProtein.from_pdb(file_name=mol_input)
+        else:
+            input_data = mol_input
+        self._save_to_protein(input_data)
+
+    @classmethod
+    def from_file(cls, file_name: str, name: Optional[str] = None) -> "Protein":
+        """
+        Instance the protein class from a pdb file.
+        """
+        cls._check_file_type(file_name=file_name)
+        input_data = ReadInputProtein.from_pdb(file_name=file_name, name=name)
+        return cls(input_data)
+
+    @staticmethod
+    def _check_file_type(file_name: str) -> None:
+        """
+        Make sure the protien is being read from a pdb file.
+        """
+        if ".pdb" not in file_name:
+            raise FileTypeError("Proteins can only be read from pdb.")
+
+    def _save_to_protein(self, mol_input: ReadInputProtein):
+        """
+        Public access to private file_handlers.py file.
+        Users shouldn't ever need to interface with file_handlers.py directly.
+        All parameters will be set from a file (or other input) via this public method.
+            * Don't bother updating name, topology or atoms if they are already stored.
+            * Do bother updating coords, rdkit_mol, residues, Residues, pdb_names
+        """
+
+        if mol_input.name is not None:
+            self.name = mol_input.name
+        if mol_input.atoms is not None:
+            self.atoms = mol_input.atoms
+        if mol_input.coords is not None:
+            self.coordinates = mol_input.coords
+        if mol_input.residues is not None:
+            self.residues = mol_input.residues
+        if mol_input.pdb_names is not None:
+            self.pdb_names = mol_input.pdb_names
+        if mol_input.bonds is not None:
+            self.bonds = mol_input.bonds
+
+        if not self.bonds:
+            print(
+                "No connections found in pdb file; topology will be inferred by OpenMM."
+            )
+            return
+
+        self.symmetrise_from_topology()
+
+    def write_pdb(self, name=None):
+        """
+        This method replaces the ligand method as all of the atom names and residue names have to be replaced.
+        """
+
+        with open(f"{name if name is not None else self.name}.pdb", "w+") as pdb_file:
+
+            pdb_file.write(f"REMARK   1 CREATED WITH QUBEKit {datetime.now()}\n")
+            # Write out the atomic xyz coordinates
+            for i, (coord, atom) in enumerate(zip(self.coordinates, self.atoms)):
+                x, y, z = coord
+                # May cause issues if protein contains more than 10,000 atoms.
+                pdb_file.write(
+                    f"HETATM {i+1:>4}{atom.atom_name:>5} QUP     1{x:12.3f}{y:8.3f}{z:8.3f}"
+                    f"  1.00  0.00         {atom.atomic_symbol.upper():>3}\n"
+                )
+
+            # Add the connection terms based on the molecule topology.
+            for bond in self.bonds:
+                pdb_file.write(
+                    f"CONECT{bond.atom1_index + 1:5} {bond.atom2_index + 1:5}\n"
+                )
+
+            pdb_file.write("END\n")
+
+    def update(self):
+        """
+        After the protein has been passed to the parametrisation class we get back the bond info
+        use this to update all missing terms.
+        """
+        if not self.bonds:
+            self.bonds = []
+            # using the new harmonic bond force dict we can add the bond edges to the topology graph
+            for bond in self.HarmonicBondForce:
+                self.bonds.append(
+                    Bond(
+                        atom1_indx=bond[0],
+                        atom2_index=bond[1],
+                        bond_order=1,
+                        aromatic=False,
+                    )
+                )
+
+        self.symmetrise_from_topology()
diff --git a/QUBEKit/engines/rdkit.py b/QUBEKit/molecules/utils.py
similarity index 63%
rename from QUBEKit/engines/rdkit.py
rename to QUBEKit/molecules/utils.py
index 74ac46e4..8e53b4c2 100644
--- a/QUBEKit/engines/rdkit.py
+++ b/QUBEKit/molecules/utils.py
@@ -1,15 +1,17 @@
-#!/usr/bin/env python3
-
 import copy
+import re
+from itertools import groupby
 from pathlib import Path
 from typing import Dict, List, Optional, Tuple
 
+import networkx as nx
 import numpy as np
 from rdkit import Chem
 from rdkit.Chem import AllChem, Descriptors
-from rdkit.Chem.rdForceFieldHelpers import MMFFOptimizeMolecule, UFFOptimizeMolecule
 from rdkit.Geometry.rdGeometry import Point3D
 
+from QUBEKit.molecules.components import Atom, Bond, Element
+from QUBEKit.utils.constants import BOHR_TO_ANGS
 from QUBEKit.utils.exceptions import FileTypeError, SmartsError
 
 
@@ -125,22 +127,6 @@ def smiles_to_rdkit_mol(smiles_string: str, name: Optional[str] = None):
 
         return mol_hydrogens
 
-    @staticmethod
-    def mm_optimise(rdkit_mol: Chem.Mol, ff="MMF") -> Chem.Mol:
-        """
-        Perform rough preliminary optimisation to speed up later optimisations.
-        Args:
-            rdkit_mol:
-                A complete Chem.Mol instance of a molecule.
-            ff:
-                The Force field to be used either MMF or UFF
-        return:
-            The name of the optimised pdb file that is made
-        """
-
-        {"MMF": MMFFOptimizeMolecule, "UFF": UFFOptimizeMolecule}[ff](rdkit_mol)
-        return rdkit_mol
-
     @staticmethod
     def rdkit_descriptors(rdkit_mol: Chem.Mol) -> Dict[str, float]:
         """
@@ -350,3 +336,219 @@ def add_conformer(
         rdkit_mol.AddConformer(conformer, assignId=True)
 
         return rdkit_mol
+
+
+class ReadInput:
+    """
+    Called inside Ligand; used to handle reading any kind of input valid in QUBEKit
+        QC JSON object
+        SMILES string
+        PDB, MOL2, XYZ file
+    """
+
+    def __init__(
+        self,
+        coords: Optional[np.ndarray] = None,
+        rdkit_mol: Optional = None,
+        name: Optional[str] = None,
+    ):
+
+        self.coords = coords
+        self.rdkit_mol = rdkit_mol
+        self.name = name
+
+    @classmethod
+    def from_smiles(cls, smiles: str, name: Optional[str] = None) -> "ReadInput":
+        """
+        Make a ReadInput object which can be taken by the Ligand class to make the model.
+
+        Note
+        ----
+        This method will generate a conformer for the molecule.
+
+        Parameters
+        ----------
+        smiles:
+            The smiles string which should be parsed by rdkit.
+        name:
+            The name that should be given to the molecule.
+        """
+        # Smiles string input
+        rdkit_mol = RDKit.smiles_to_rdkit_mol(smiles_string=smiles, name=name)
+        return cls(name=name, coords=None, rdkit_mol=rdkit_mol)
+
+    @classmethod
+    def from_file(cls, file_name: str) -> "ReadInput":
+        """
+        Read the input file using RDKit and return the molecule data.
+        """
+        input_file = Path(file_name)
+        # if the file is not there raise an error
+        if not input_file.exists():
+            raise FileNotFoundError(
+                f"{input_file.as_posix()} could not be found is this path correct?"
+            )
+        # xyz is a special case of file only internal readers catch
+        if input_file.suffix == ".xyz":
+            return cls.from_xyz(file_name=input_file.as_posix())
+        # read the input with rdkit
+        rdkit_mol = RDKit.file_to_rdkit_mol(file_path=input_file)
+        return cls(rdkit_mol=rdkit_mol, coords=None, name=rdkit_mol.GetProp("_Name"))
+
+    @classmethod
+    def from_qc_json(cls, qc_json) -> "ReadInput":
+        """
+        Given a QC JSON object, extracts the topology, atoms and coords of the molecule.
+        #TODO we need to be absle to read mapped smiles for this to work with stereochem and aromaticity
+        """
+
+        topology = nx.Graph()
+        atoms = []
+
+        for i, atom in enumerate(qc_json.symbols):
+            atoms.append(
+                Atom(
+                    atomic_number=Element().number(atom),
+                    atom_index=i,
+                    atom_name=f"{atom}{i}",
+                )
+            )
+            topology.add_node(i)
+
+        for bond in qc_json.connectivity:
+            topology.add_edge(*bond[:2])
+
+        coords = np.array(qc_json.geometry).reshape((len(atoms), 3)) * BOHR_TO_ANGS
+        return cls(name=None, rdkit_mol=None, coords=coords)
+
+    @classmethod
+    def from_xyz(cls, file_name: str) -> "ReadInput":
+        """
+        Internal xyz reader.
+        Extracts the coords of the molecule.
+        """
+
+        traj_molecules = []
+        coords = []
+
+        with open(file_name) as xyz_file:
+            lines = xyz_file.readlines()
+
+            n_atoms = float(lines[0])
+
+            for line in lines:
+                line = line.split()
+                # skip frame heading lines
+                if len(line) <= 1 or "Iteration" in line:
+                    continue
+
+                coords.append([float(line[1]), float(line[2]), float(line[3])])
+
+                if len(coords) == n_atoms:
+                    # we have collected the molecule now store the frame
+                    traj_molecules.append(np.array(coords))
+                    coords = []
+
+        coords = traj_molecules[0] if len(traj_molecules) == 1 else traj_molecules
+        return cls(coords=coords, name=None, rdkit_mol=None)
+
+
+class ReadInputProtein:
+    """
+    A class that specialises in reading Protein input files.
+    #TODO are we better doing this with openmm or another tool?
+    """
+
+    def __init__(
+        self,
+        atoms: List[Atom],
+        bonds: Optional[List[Bond]] = None,
+        coords: Optional[np.ndarray] = None,
+        pdb_names: Optional[List[str]] = None,
+        residues: Optional[List[str]] = None,
+        name: Optional[str] = None,
+    ):
+        self.atoms = atoms
+        self.bonds = bonds
+        self.coords = coords
+        self.name = name
+        self.residues = residues
+        self.pdb_names = pdb_names
+
+    @classmethod
+    def from_pdb(cls, file_name: str, name: Optional[str] = None):
+        """
+        Read the protein input pdb file.
+        :return:
+        """
+        with open(file_name, "r") as pdb:
+            lines = pdb.readlines()
+
+        coords = []
+        atoms = []
+        bonds = []
+        Residues = []
+        pdb_names = []
+
+        # atom counter used for graph node generation
+        atom_count = 0
+        for line in lines:
+            if "ATOM" in line or "HETATM" in line:
+                atomic_symbol = str(line[76:78])
+                atomic_symbol = re.sub("[0-9]+", "", atomic_symbol).strip()
+
+                # If the element column is missing from the pdb, extract the atomic_symbol from the atom name.
+                if not atomic_symbol:
+                    atomic_symbol = str(line.split()[2])
+                    atomic_symbol = re.sub("[0-9]+", "", atomic_symbol)
+
+                # now make sure we have a valid element
+                if atomic_symbol.lower() != "cl" and atomic_symbol.lower() != "br":
+                    atomic_symbol = atomic_symbol[0]
+
+                atom_name = f"{atomic_symbol}{atom_count}"
+                # TODO should we use a protein pdb package for this?
+                qube_atom = Atom(
+                    atomic_number=Element().number(atomic_symbol),
+                    atom_index=atom_count,
+                    atom_name=atom_name,
+                    formal_charge=0,
+                    aromatic=False,
+                )
+
+                atoms.append(qube_atom)
+
+                pdb_names.append(str(line.split()[2]))
+
+                # also get the residue order from the pdb file so we can rewrite the file
+                Residues.append(str(line.split()[3]))
+
+                atom_count += 1
+                coords.append(
+                    [float(line[30:38]), float(line[38:46]), float(line[46:54])]
+                )
+
+            elif "CONECT" in line:
+                conect_terms = line.split()
+                for atom in conect_terms[2:]:
+                    if int(atom):
+                        bond = Bond(
+                            atom1_index=int(conect_terms[1]) - 1,
+                            atom2_index=int(atom) - 1,
+                            bond_order=1,
+                            aromatic=False,
+                        )
+                        bonds.append(bond)
+
+        coords = np.array(coords)
+        residues = [res for res, group in groupby(Residues)]
+        if name is None:
+            name = Path(file_name).stem
+        return cls(
+            atoms=atoms,
+            bonds=bonds,
+            coords=coords,
+            pdb_names=pdb_names,
+            residues=residues,
+            name=name,
+        )
diff --git a/QUBEKit/parametrisation/__init__.py b/QUBEKit/parametrisation/__init__.py
index dd615d66..dc1e56b8 100644
--- a/QUBEKit/parametrisation/__init__.py
+++ b/QUBEKit/parametrisation/__init__.py
@@ -1,10 +1,5 @@
-#!/usr/bin/env python3
-
-from QUBEKit.utils.helpers import try_load
-
-from .base_parametrisation import Parametrisation
-from .openff import OpenFF
-from .xml import XML
-from .xml_protein import XMLProtein
-
-AnteChamber = try_load("AnteChamber", "QUBEKit.parametrisation.antechamber")
+from QUBEKit.parametrisation.antechamber import AnteChamber
+from QUBEKit.parametrisation.base_parametrisation import Parametrisation
+from QUBEKit.parametrisation.openff import OpenFF
+from QUBEKit.parametrisation.xml import XML
+from QUBEKit.parametrisation.xml_protein import XMLProtein
diff --git a/QUBEKit/parametrisation/base_parametrisation.py b/QUBEKit/parametrisation/base_parametrisation.py
index f0b9bfdb..69492b3f 100644
--- a/QUBEKit/parametrisation/base_parametrisation.py
+++ b/QUBEKit/parametrisation/base_parametrisation.py
@@ -1,12 +1,9 @@
-#!/usr/bin/env python3
-
 import xml.etree.ElementTree as ET
 from collections import OrderedDict
 from copy import deepcopy
 
-from QUBEKit.ligand import Ligand
+from QUBEKit.molecules.ligand import ExtraSite, Ligand
 from QUBEKit.utils import constants
-from QUBEKit.utils.datastructures import ExtraSite
 
 
 class Parametrisation:
diff --git a/QUBEKit/proteins/protein_run.py b/QUBEKit/proteins/protein_run.py
index 16b13be0..6ed48c96 100644
--- a/QUBEKit/proteins/protein_run.py
+++ b/QUBEKit/proteins/protein_run.py
@@ -6,7 +6,7 @@
 from functools import partial
 
 from QUBEKit.lennard_jones import LennardJones
-from QUBEKit.ligand import Protein
+from QUBEKit.molecules import Protein
 from QUBEKit.parametrisation import XMLProtein
 from QUBEKit.proteins.protein_tools import get_water, pdb_reformat, qube_general
 from QUBEKit.utils.file_handling import ExtractChargeData
@@ -32,7 +32,7 @@ def __call__(self, pars, namespace, values, option_string=None):
             protein = Protein(values)
             printf(f"{len(protein.residues)} residues found!")
             # TODO find the magic numbers for the box for onetep
-            protein.write_xyz(name="protein")
+            protein.to_file(file_name="protein.xyz")
             printf(f"protein.xyz file made for ONETEP\n Run this file")
             sys.exit()
 
@@ -49,7 +49,6 @@ def __call__(self, pars, namespace, values, option_string=None):
 
             # finally we need the non-bonded parameters from onetep
             # TODO should we also have the ability to get DDEC6 charges from the cube file?
-            pro.charge = 0
             pro.charges_engine = "onetep"
             pro.density_engine = "onetep"
 
diff --git a/QUBEKit/proteins/protein_tools.py b/QUBEKit/proteins/protein_tools.py
index 6685cdb2..8d91f51a 100644
--- a/QUBEKit/proteins/protein_tools.py
+++ b/QUBEKit/proteins/protein_tools.py
@@ -325,7 +325,7 @@ def get_water(water=None):
 def pdb_reformat(reference, target):
     """Rewrite a QUBEKit system pdb traj file to have the correct atom and residue names."""
 
-    from QUBEKit.ligand import Protein
+    from QUBEKit.molecules.protein import Protein
 
     # load the original file as a reference
     pro = Protein(reference)
diff --git a/QUBEKit/run.py b/QUBEKit/run.py
index 2d494590..9daba210 100644
--- a/QUBEKit/run.py
+++ b/QUBEKit/run.py
@@ -12,7 +12,6 @@
 
 import argparse
 import os
-import subprocess as sp
 import sys
 from collections import OrderedDict
 from datetime import datetime
@@ -23,12 +22,11 @@
 
 import QUBEKit
 from QUBEKit.dihedrals import TorsionOptimiser, TorsionScan
-from QUBEKit.engines import PSI4, Chargemol, Gaussian, QCEngine, RDKit
+from QUBEKit.engines import Chargemol, Gaussian, GeometryOptimiser, QCEngine
 from QUBEKit.lennard_jones import LennardJones
-from QUBEKit.ligand import Ligand
 from QUBEKit.mod_seminario import ModSeminario
+from QUBEKit.molecules import Ligand
 from QUBEKit.parametrisation import XML, AnteChamber, OpenFF, Parametrisation
-from QUBEKit.utils import constants
 from QUBEKit.utils.configs import Configure
 from QUBEKit.utils.constants import COLOURS
 from QUBEKit.utils.decorators import exception_logger
@@ -37,7 +35,7 @@
     pretty_print,
     pretty_progress,
 )
-from QUBEKit.utils.exceptions import HessianCalculationFailed, OptimisationFailed
+from QUBEKit.utils.exceptions import HessianCalculationFailed, SpecificationError
 from QUBEKit.utils.file_handling import (
     ExtractChargeData,
     extract_extra_sites_onetep,
@@ -304,17 +302,17 @@ def __call__(self, pars, namespace, values, option_string=None):
             type=int,
             help="Enter the ddec version for charge partitioning, does not effect ONETEP partitioning.",
         )
-        parser.add_argument(
-            "-geo",
-            "--geometric",
-            choices=[True, False],
-            type=string_to_bool,
-            help="Turn on geometric to use this during the qm optimisations, recommended.",
-        )
+        # parser.add_argument(
+        #     "-geo",
+        #     "--geometric",
+        #     choices=[True, False],
+        #     type=string_to_bool,
+        #     help="Turn on geometric to use this during the qm optimisations, recommended.",
+        # )
         parser.add_argument(
             "-bonds",
             "--bonds_engine",
-            choices=["psi4", "g09", "g16"],
+            choices=["psi4", "gaussian"],
             help="Choose the QM code to calculate the bonded terms.",
         )
         parser.add_argument(
@@ -352,10 +350,22 @@ def __call__(self, pars, namespace, values, option_string=None):
             "if xml make sure the xml has the same name as the pdb file.",
         )
         parser.add_argument(
-            "-mm",
-            "--mm_opt_method",
-            choices=["openmm", "rdkit_mff", "rdkit_uff", "none"],
-            help="Enter the mm optimisation method for pre qm optimisation.",
+            "-pre-opt",
+            "--pre_opt_method",
+            choices=[
+                "mmff94",
+                "mmff94s",
+                "uff",
+                "gfn1xtb",
+                "gfn2xtb",
+                "gfn0xtb",
+                "gaff-2.11",
+                "ani1x",
+                "ani1ccx",
+                "ani2x",
+                "openff-1.3.0",
+            ],
+            help="Enter the optimisation method for pre qm optimisation.",
         )
         parser.add_argument(
             "-config",
@@ -377,7 +387,7 @@ def __call__(self, pars, namespace, values, option_string=None):
             "--restart",
             choices=[
                 "parametrise",
-                "mm_optimise",
+                "pre_optimise",
                 "qm_optimise",
                 "hessian",
                 "mod_sem",
@@ -394,7 +404,7 @@ def __call__(self, pars, namespace, values, option_string=None):
             "-end",
             choices=[
                 "parametrise",
-                "mm_optimise",
+                "pre_optimise",
                 "qm_optimise",
                 "hessian",
                 "mod_sem",
@@ -548,7 +558,9 @@ def __call__(self, pars, namespace, values, option_string=None):
         groups.add_argument(
             "-i", "--input", help="Enter the molecule input pdb file (only pdb so far!)"
         )
-        groups.add_argument("-version", "--version", action="version", version="2.6.3")
+        groups.add_argument(
+            "-version", "--version", action="version", version=QUBEKit.__version__
+        )
 
         # Ensures help is shown (rather than an error) if no arguments are provided.
         return parser.parse_args(args=None if sys.argv[1:] else ["--help"])
@@ -641,7 +653,7 @@ def __init__(self, molecule):
         self.order = OrderedDict(
             [
                 ("parametrise", self.parametrise),
-                ("mm_optimise", self.mm_optimise),
+                ("pre_optimise", self.pre_optimise),
                 ("qm_optimise", self.qm_optimise),
                 ("hessian", self.hessian),
                 ("mod_sem", self.mod_sem),
@@ -657,7 +669,7 @@ def __init__(self, molecule):
         # Keep this for reference (used for numbering folders correctly)
         self.immutable_order = tuple(self.order)
 
-        self.engine_dict = {"psi4": PSI4, "g09": Gaussian, "g16": Gaussian}
+        self.engine_dict = {"g09": Gaussian, "g16": Gaussian}
 
         printf(self.start_up_msg)
 
@@ -822,8 +834,8 @@ def run(self, torsion_options=None):
                 f"Parametrising molecule with {self.molecule.parameter_engine}",
                 "Molecule parametrised",
             ],
-            "mm_optimise": [
-                "Partially optimising with MM",
+            "pre_optimise": [
+                f"Partially optimising with {self.molecule.pre_opt_method}",
                 "Partial optimisation complete",
             ],
             "qm_optimise": [
@@ -979,208 +991,111 @@ def parametrise(molecule: Ligand, verbose: bool = True) -> Ligand:
         return molecule
 
     @staticmethod
-    def mm_optimise(molecule: Ligand) -> Ligand:
+    def pre_optimise(molecule: Ligand) -> Ligand:
         """
-        Use an mm force field to get the initial optimisation of a molecule
+        Do a pre optimisation of the molecule using the specified program.
 
         options
         ---------
-        RDKit MFF or UFF force fields can have strange effects on the geometry of molecules
+        "mmff94", "uff", "mmff94s", "gfn1xtb", "gfn2xtb", "fgn0xtb", "gaff-2.11", "ani1x", "ani1ccx", "ani2x", "openff-1.3.0
 
-        Geometric / OpenMM depends on the force field the molecule was parameterised with gaff/2, OPLS smirnoff.
-        #TODO replace with a general optimiser using QCEngine.
         """
-
-        append_to_log(molecule.home, "Starting mm_optimisation", major=True)
-        # Check which method we want then do the optimisation
-        if molecule.mm_opt_method == "openmm":
-            if molecule.parameter_engine == "none":
-                raise OptimisationFailed(
-                    "You cannot optimise a molecule with OpenMM and no initial parameters; "
-                    "consider parametrising or using UFF/MFF in RDKit"
-                )
-            else:
-                # Make the inputs
-                molecule.to_file(file_name=f"{molecule.name}.pdb")
-                molecule.write_parameters()
-                # Run geometric
-                # TODO Should this be moved to a function? Seems like a likely point of failure
-                with open("log.txt", "w+") as log:
-                    sp.run(
-                        f"geometric-optimize --epsilon 0.0 --maxiter {molecule.iterations} --pdb "
-                        f"{molecule.name}.pdb --openmm {molecule.name}.xml "
-                        f'{molecule.constraints_file if molecule.constraints_file is not None else ""}',
-                        shell=True,
-                        stdout=log,
-                        stderr=log,
-                    )
-
-                molecule.add_conformer(f"{molecule.name}_optim.xyz")
-
+        # TODO drop all of this once we change configs
+
+        # now we want to build the optimiser from the inputs
+        method = molecule.pre_opt_method.lower()
+        if method in ["mmff94", "mmff94s", "uff"]:
+            program = "rdkit"
+            basis = None
+        elif method in ["gfn1xtb", "gfn2xtb", "gfn0xtb"]:
+            program = "xtb"
+            basis = None
+        elif method in ["ani1x", "ani1ccx", "ani2x"]:
+            program = "torchani"
+            basis = None
+        elif method == "gaff-2.11":
+            program = "openmm"
+            basis = "antechamber"
+        elif method == "openff-1.3.0":
+            program = "openmm"
+            basis = "smirnoff"
         else:
-            # TODO change to qcengine as this can already be done
-            # Run an rdkit optimisation with the right FF
-            rdkit_ff = {"rdkit_mff": "MFF", "rdkit_uff": "UFF"}[molecule.mm_opt_method]
-            rdkit_mol = RDKit.mm_optimise(molecule.to_rdkit(), ff=rdkit_ff)
-            molecule.coordinates = rdkit_mol.GetConformer().GetPositions()
+            raise SpecificationError(
+                f"The pre optimisation method {method} is not supported please chose from "
+                f"mmff94, mmff94s, uff, gfn1xtb, gfn2xtb, gfn0xtb, gaff-2.11, ani1x, ani1ccx, ani2x, openff-1.3.0"
+            )
 
         append_to_log(
             molecule.home,
-            f"Finishing mm_optimisation of the molecule with {molecule.mm_opt_method}",
+            f"Starting pre_optimisation with program: {program} basis: {basis} method: {method}",
             major=True,
         )
-        return molecule
-
-    def qm_optimise(self, molecule: Ligand) -> Ligand:
-        """Optimise the molecule coords. Can be through PSI4 (with(out) geometric) or through Gaussian."""
-
-        append_to_log(molecule.home, "Starting qm_optimisation", major=True)
-        MAX_RESTARTS = 3
-
-        if molecule.geometric and (molecule.bonds_engine == "psi4"):
-            qceng = QCEngine(molecule)
-            result = qceng.call_qcengine(
-                engine="geometric",
-                driver="gradient",
-            )
-
-            restart_count = 0
-
-            while (not result["success"]) and (restart_count < MAX_RESTARTS):
-                append_to_log(
-                    molecule.home,
-                    f'{molecule.bonds_engine} optimisation failed with error {result["error"]}; restarting',
-                )
-
-                try:
-                    molecule.coordinates = np.array(
-                        result["input_data"]["final_molecule"]["geometry"]
-                    ).reshape((len(molecule.atoms), 3))
-                    molecule.coordinates *= constants.BOHR_TO_ANGS
-
-                    result = qceng.call_qcengine(engine="geometric", driver="gradient")
-
-                except (KeyError, TypeError):
-                    result = qceng.call_qcengine(
-                        engine="geometric",
-                        driver="gradient",
-                    )
-
-                restart_count += 1
-
-            if not result["success"]:
-                raise OptimisationFailed("The optimisation did not converge")
-            trajectory = [
-                np.array(mol["molecule"]["geometry"]).reshape((molecule.n_atoms, 3))
-                * constants.BOHR_TO_ANGS
-                for mol in result["trajectory"]
-            ]
-
-            # store the final molecule as the qm optimised structure
-            molecule.coordinates = np.array(
-                result["final_molecule"]["geometry"]
-            ).reshape((molecule.n_atoms, 3))
-            molecule.coordinates *= constants.BOHR_TO_ANGS
-
-            molecule.qm_energy = result["energies"][-1]
-
-            # Write out the trajectory file
-            molecule.to_multiconformer_file(
-                file_name=f"{molecule.name}_opt.xyz", positions=trajectory
-            )
-            molecule.to_file(file_name="opt.xyz")
 
-        # Using Gaussian or geometric off
-        else:
-            qm_engine = self.engine_dict[molecule.bonds_engine](molecule)
-            result = qm_engine.generate_input(
-                optimise=True,
-                execute=molecule.bonds_engine,
-            )
-
-            restart_count = 0
-            while (not result["success"]) and (restart_count < MAX_RESTARTS):
-                append_to_log(
-                    molecule.home,
-                    f'{molecule.bonds_engine} optimisation failed with error {result["error"]}; restarting',
-                )
-
-                if result["error"] == "FileIO":
-                    result = qm_engine.generate_input(
-                        optimise=True, restart=True, execute=molecule.bonds_engine
-                    )
-                elif result["error"] == "Max iterations":
-                    result = qm_engine.generate_input(
-                        optimise=True,
-                        restart=True,
-                        execute=molecule.bonds_engine,
-                    )
-
-                else:
-                    molecule.coordinates = RDKit.generate_conformers(
-                        molecule.to_rdkit()
-                    )[-1]
-                    result = qm_engine.generate_input(
-                        optimise=True, execute=molecule.bonds_engine
-                    )
-
-                restart_count += 1
-
-            if not result["success"]:
-                raise OptimisationFailed(
-                    f"{molecule.bonds_engine} optimisation did not converge after 3 restarts; "
-                    f"last error {result['error']}"
-                )
+        g_opt = GeometryOptimiser(
+            program=program,
+            method=method,
+            basis=basis,
+            convergence="GAU",
+            maxiter=molecule.iterations,
+        )
+        # errors are auto raised from the class so catch the result, and write to file
+        result_mol = g_opt.optimise(molecule=molecule, allow_fail=True)
 
-            molecule.coordinates, molecule.qm_energy = qm_engine.optimised_structure()
-            molecule.to_file(file_name="opt.xyz")
+        append_to_log(
+            molecule.home,
+            f"Finishing mm_optimisation of the molecule with {molecule.pre_opt_method}",
+            major=True,
+        )
+        return result_mol
 
+    @staticmethod
+    def qm_optimise(molecule: Ligand) -> Ligand:
+        """
+        Optimise the molecule using qm via qcengine.
+        """
         append_to_log(
             molecule.home,
-            f'Finishing qm_optimisation of molecule{" using geometric" if molecule.geometric else ""}',
+            f"Starting qm_optimisation with program: {molecule.bonds_engine} basis: {molecule.basis} method: {molecule.theory}",
             major=True,
         )
+        # TODO do we want the geometry optimiser to handle restarts?
+        g_opt = GeometryOptimiser(
+            program=molecule.bonds_engine,
+            method=molecule.theory,
+            basis=molecule.basis,
+            convergence=molecule.convergence,
+            maxiter=molecule.iterations,
+        )
+        # errors are auto raised from the class output is always dumped to file
+        qm_result = g_opt.optimise(molecule=molecule, allow_fail=False)
+        append_to_log(molecule.home, f"QM optimisation finished", major=True)
 
-        return molecule
+        return qm_result
 
     @staticmethod
     def hessian(molecule: Ligand) -> Ligand:
-        """Using the assigned bonds engine, calculate and extract the Hessian matrix."""
+        """Using the assigned bonds engine, calculate the Hessian matrix and store in atomic units."""
+        from QUBEKit.utils.helpers import check_symmetry
 
         append_to_log(molecule.home, "Starting hessian calculation", major=True)
+        # build the QM engine
+        qm_engine = QCEngine(
+            program=molecule.bonds_engine,
+            method=molecule.theory,
+            basis=molecule.basis,
+            driver="hessian",
+        )
+        result = qm_engine.call_qcengine(molecule=molecule)
 
-        if molecule.bonds_engine in ["g09", "g16"]:
-            qm_engine = Gaussian(molecule)
-
-            # Use the checkpoint file as this has higher xyz precision
-            try:
-                copy(
-                    os.path.join(molecule.home, "03_qm_optimise", "lig.chk"), "lig.chk"
-                )
-                result = qm_engine.generate_input(
-                    "qm", hessian=True, restart=True, execute=molecule.bonds_engine
-                )
-            except FileNotFoundError:
-                append_to_log(
-                    molecule.home,
-                    "qm_optimise checkpoint not found, optimising first to refine atomic coordinates",
-                )
-                result = qm_engine.generate_input(
-                    optimise=True, hessian=True, execute=molecule.bonds_engine
-                )
-
-            if not result["success"]:
-                raise HessianCalculationFailed(
-                    "The hessian was not calculated check the log file."
-                )
-
-            hessian = qm_engine.hessian()
+        if not result.success:
+            raise HessianCalculationFailed(
+                "The hessian was not calculated check the log file."
+            )
 
-        else:
-            hessian = QCEngine(molecule).call_qcengine(engine="psi4", driver="hessian")
-            np.savetxt("hessian.txt", hessian)
+        np.savetxt("hessian.txt", result.return_result)
 
-        molecule.hessian = hessian
+        molecule.hessian = result.return_result
+        check_symmetry(molecule.hessian)
 
         append_to_log(
             molecule.home,
@@ -1214,7 +1129,9 @@ def density(self, molecule):
         if molecule.density_engine == "onetep":
             molecule.write_xyz(input_type="qm")
             # If using ONETEP, stop after this step
-            append_to_log(molecule.home, "Density analysis file made for ONETEP")
+            append_to_log(
+                molecule.home, "Density analysis file made for ONETEP", major=True
+            )
 
             # Edit the order to end here
             self.order = OrderedDict(
diff --git a/QUBEKit/utils/datastructures.py b/QUBEKit/utils/datastructures.py
index b970b1bc..4c71dd2b 100644
--- a/QUBEKit/utils/datastructures.py
+++ b/QUBEKit/utils/datastructures.py
@@ -1,12 +1,4 @@
-#!/usr/bin/env python3
-
-from enum import Enum
 from types import SimpleNamespace
-from typing import List, Optional
-
-from pydantic import BaseModel, Field, PositiveInt
-from rdkit import Chem
-from rdkit.Chem.rdchem import GetPeriodicTable, PeriodicTable
 
 
 class CustomNamespace(SimpleNamespace):
@@ -29,264 +21,3 @@ def items(self):
 
     def __iter__(self):
         return self.items()
-
-
-class Element:
-    """
-    Simple wrapper class for getting element info using RDKit.
-    """
-
-    @staticmethod
-    def p_table() -> PeriodicTable:
-        return GetPeriodicTable()
-
-    @staticmethod
-    def mass(identifier):
-        pt = Element.p_table()
-        return pt.GetAtomicWeight(identifier)
-
-    @staticmethod
-    def number(identifier):
-        pt = Element.p_table()
-        return pt.GetAtomicNumber(identifier)
-
-    @staticmethod
-    def name(identifier):
-        pt = Element.p_table()
-        return pt.GetElementSymbol(identifier)
-
-
-class AtomStereoChemistry(str, Enum):
-    """
-    Atom stereochemistry types.
-    """
-
-    R = "R"
-    S = "S"
-    U = "Unknown"
-
-
-class BondStereoChemistry(str, Enum):
-    """
-    Bond stereochemistry types.
-    """
-
-    E = "E"
-    Z = "Z"
-    U = "Unknown"
-
-
-class Atom(BaseModel):
-    """
-    Class to hold all of the "per atom" information.
-    All atoms in Molecule will have an instance of this Atom class to describe their properties.
-    """
-
-    class Config:
-        validate_assignment = True
-        json_encoders = {Enum: lambda v: v.value}
-
-    atomic_number: PositiveInt = Field(
-        ...,
-        description="The atomic number of the atom all other properties are based on this number.",
-    )
-    atom_index: int = Field(
-        ..., description="The index this atom has in the molecule object", gt=-1
-    )
-    atom_name: Optional[str] = Field(
-        None,
-        description="An optional unqiue atom name that should be assigned to the atom, the ligand object will make sure all atoms have unique names.",
-    )
-    formal_charge: int = Field(
-        ...,
-        description="The formal charge of the atom, used to calculate the molecule total charge",
-    )
-    aromatic: bool = Field(
-        ...,
-        description="If the atom should be considered aromatic `True` or not `False`.",
-    )
-    stereochemistry: Optional[AtomStereoChemistry] = Field(
-        None,
-        description="The stereochemistry of the atom where None means not stereogenic and U is unknown or ambiguous.",
-    )
-    bonds: Optional[List[int]] = Field(
-        None, description="The list of atom indices which are bonded to this atom."
-    )
-
-    @classmethod
-    def from_rdkit(cls, rd_atom: Chem.Atom) -> "Atom":
-        """
-        Build a QUBEKit atom from an rdkit atom instance.
-        """
-        atomic_number = rd_atom.GetAtomicNum()
-        index = rd_atom.GetIdx()
-        formal_charge = rd_atom.GetFormalCharge()
-        aromatic = rd_atom.GetIsAromatic()
-        bonds = [a.GetIdx() for a in rd_atom.GetNeighbors()]
-        # check for names in the normal places pdb, mol2 and mol
-        if rd_atom.HasProp("_Name"):
-            name = rd_atom.GetProp("_Name")
-        elif rd_atom.HasProp("_TriposAtomName"):
-            name = rd_atom.GetProp("_TriposAtomName")
-        else:
-            try:
-                name = rd_atom.GetMonomerInfo().GetName().strip()
-            except AttributeError:
-                name = None
-        # stereochem
-        if rd_atom.HasProp("_CIPCode"):
-            stereo_code = rd_atom.GetProp("_CIPCode")
-        else:
-            stereo_code = None
-        return cls(
-            atomic_number=atomic_number,
-            atom_index=index,
-            atom_name=name,
-            formal_charge=formal_charge,
-            aromatic=aromatic,
-            stereochemistry=stereo_code,
-            bonds=bonds,
-        )
-
-    @property
-    def atomic_mass(self) -> float:
-        """Convert the atomic number to mass."""
-        return Element.mass(self.atomic_number)
-
-    @property
-    def atomic_symbol(self) -> str:
-        """Convert the atomic number to the atomic symbol as per the periodic table."""
-        return Element.name(self.atomic_number).title()
-
-    def to_rdkit(self) -> Chem.Atom:
-        """
-        Convert the QUBEKit atom an RDKit atom.
-        """
-        # build the atom from atomic number
-        rd_atom = Chem.Atom(self.atomic_number)
-        rd_atom.SetFormalCharge(self.formal_charge)
-        rd_atom.SetIsAromatic(self.aromatic)
-        rd_atom.SetProp("_Name", self.atom_name)
-        # left is counter clockwise
-        if self.stereochemistry == "S":
-            rd_atom.SetChiralTag(Chem.CHI_TETRAHEDRAL_CCW)
-        # right is clockwise
-        elif self.stereochemistry == "R":
-            rd_atom.SetChiralTag(Chem.CHI_TETRAHEDRAL_CW)
-
-        return rd_atom
-
-    # def __repr__(self):
-    #     return f"{self.__class__.__name__}({self.__dict__!r})"
-    #
-    # def __str__(self):
-    #     """
-    #     Prints the Atom class objects' names and values one after another with new lines between each.
-    #     """
-    #
-    #     return_str = ""
-    #     for key, val in self.__dict__.items():
-    #         # Return all objects as {atom object name} = {atom object value(s)}.
-    #         return_str += f"\n{key} = {val}\n"
-    #
-    #     return return_str
-
-
-class Bond(BaseModel):
-    """
-    A basic bond class.
-    """
-
-    class Config:
-        validate_assignment = True
-        json_encoders = {Enum: lambda v: v.value}
-
-    atom1_index: int = Field(
-        ..., description="The index of the first atom in the bond."
-    )
-    atom2_index: int = Field(
-        ..., description="The index of the second atom in the bond."
-    )
-    bond_order: float = Field(..., description="The float value of the bond order.")
-    aromatic: bool = Field(
-        ..., description="If the bond should be considered aromatic."
-    )
-    stereochemistry: Optional[BondStereoChemistry] = Field(
-        None,
-        description="The stereochemistry of the bond, where None means not stereogenic.",
-    )
-
-    @classmethod
-    def from_rdkit(cls, rd_bond: Chem.Bond) -> "Bond":
-        """
-        Build a QUBEKit bond class from an rdkit reference.
-        """
-        atom1_index = rd_bond.GetBeginAtomIdx()
-        atom2_index = rd_bond.GetEndAtomIdx()
-        aromatic = rd_bond.GetIsAromatic()
-        order = rd_bond.GetBondTypeAsDouble()
-        stereo_tag = rd_bond.GetStereo()
-        if stereo_tag == Chem.BondStereo.STEREOZ:
-            stereo = "Z"
-        elif stereo_tag == Chem.BondStereo.STEREOE:
-            stereo = "E"
-        else:
-            stereo = None
-        return cls(
-            atom1_index=atom1_index,
-            atom2_index=atom2_index,
-            aromatic=aromatic,
-            bond_order=order,
-            stereochemistry=stereo,
-        )
-
-    @property
-    def rdkit_type(self) -> Chem.BondType:
-        """
-        Convert the bond order float to a bond type.
-        """
-        conversion = {
-            1: Chem.BondType.SINGLE,
-            1.5: Chem.BondType.AROMATIC,
-            2: Chem.BondType.DOUBLE,
-            3: Chem.BondType.TRIPLE,
-            4: Chem.BondType.QUADRUPLE,
-            5: Chem.BondType.QUINTUPLE,
-            6: Chem.BondType.HEXTUPLE,
-            7: Chem.BondType.ONEANDAHALF,
-        }
-        return conversion[self.bond_order]
-
-    @property
-    def rdkit_stereo(self) -> Optional[Chem.BondStereo]:
-        """
-        Return the rdkit style stereo enum.
-        """
-        if self.stereochemistry == "E":
-            return Chem.BondStereo.STEREOE
-        elif self.stereochemistry == "Z":
-            return Chem.BondStereo.STEREOZ
-        return None
-
-
-class ExtraSite:
-    """
-    Used to store extra sites for xml writer in ligand.
-    This class is used by both internal v-sites fitting and the ONETEP reader.
-    """
-
-    def __init__(self):
-        self.parent_index: Optional[int] = None
-        self.closest_a_index: Optional[int] = None
-        self.closest_b_index: Optional[int] = None
-        # Optional: Used for Nitrogen only.
-        self.closest_c_index: Optional[int] = None
-
-        self.o_weights: Optional[List[float]] = None
-        self.x_weights: Optional[List[float]] = None
-        self.y_weights: Optional[List[float]] = None
-
-        self.p1: Optional[float] = None
-        self.p2: Optional[float] = None
-        self.p3: Optional[float] = None
-        self.charge: Optional[float] = None
diff --git a/QUBEKit/utils/exceptions.py b/QUBEKit/utils/exceptions.py
index 48937acc..0f6c4a7a 100644
--- a/QUBEKit/utils/exceptions.py
+++ b/QUBEKit/utils/exceptions.py
@@ -56,3 +56,21 @@ class SmartsError(Exception):
     """
     The requested smarts could not be interpreted by rdkit
     """
+
+
+class SpecificationError(Exception):
+    """
+    When a requested specification is not possible.
+    """
+
+
+class ConformerError(Exception):
+    """
+    When a conformation is incorrect or missing.
+    """
+
+
+class StereoChemistryError(Exception):
+    """
+    When the stereochemistry is incorrect in an rdkit instance created from a qubekit molecule.
+    """
diff --git a/QUBEKit/utils/file_handling.py b/QUBEKit/utils/file_handling.py
index 236f0881..5bca9407 100644
--- a/QUBEKit/utils/file_handling.py
+++ b/QUBEKit/utils/file_handling.py
@@ -1,240 +1,15 @@
-#!/usr/bin/env python3
-
 """
 Purpose of this file is to read various inputs and produce the info required for
     Ligand() or Protein()
 Should be very little calculation here, simply file reading and some small validations / checks
 """
-
 import os
-import re
-from itertools import groupby
-from pathlib import Path
-from typing import List, Optional
 
-import networkx as nx
 import numpy as np
 
-from QUBEKit.engines import RDKit
-from QUBEKit.utils.constants import ANGS_TO_NM, BOHR_TO_ANGS
-from QUBEKit.utils.datastructures import Atom, Bond, CustomNamespace, Element, ExtraSite
-
-
-class ReadInput:
-    """
-    Called inside Ligand; used to handle reading any kind of input valid in QUBEKit
-        QC JSON object
-        SMILES string
-        PDB, MOL2, XYZ file
-    """
-
-    def __init__(
-        self,
-        coords: Optional[np.ndarray] = None,
-        rdkit_mol: Optional = None,
-        name: Optional[str] = None,
-    ):
-
-        self.coords = coords
-        self.rdkit_mol = rdkit_mol
-        self.name = name
-
-    @classmethod
-    def from_smiles(cls, smiles: str, name: Optional[str] = None) -> "ReadInput":
-        """
-        Make a ReadInput object which can be taken by the Ligand class to make the model.
-
-        Note
-        ----
-        This method will generate a conformer for the molecule.
-
-        Parameters
-        ----------
-        smiles:
-            The smiles string which should be parsed by rdkit.
-        name:
-            The name that should be given to the molecule.
-        """
-        # Smiles string input
-        rdkit_mol = RDKit.smiles_to_rdkit_mol(smiles_string=smiles, name=name)
-        return cls(name=name, coords=None, rdkit_mol=rdkit_mol)
-
-    @classmethod
-    def from_file(cls, file_name: str) -> "ReadInput":
-        """
-        Read the input file using RDKit and return the molecule data.
-        """
-        input_file = Path(file_name)
-        # if the file is not there raise an error
-        if not input_file.exists():
-            raise FileNotFoundError(
-                f"{input_file.as_posix()} could not be found is this path correct?"
-            )
-        # xyz is a special case of file only internal readers catch
-        if input_file.suffix == ".xyz":
-            return cls.from_xyz(file_name=input_file.as_posix())
-        # read the input with rdkit
-        rdkit_mol = RDKit.file_to_rdkit_mol(file_path=input_file)
-        return cls(rdkit_mol=rdkit_mol, coords=None, name=rdkit_mol.GetProp("_Name"))
-
-    @classmethod
-    def from_qc_json(cls, qc_json) -> "ReadInput":
-        """
-        Given a QC JSON object, extracts the topology, atoms and coords of the molecule.
-        #TODO we need to be absle to read mapped smiles for this to work with stereochem and aromaticity
-        """
-
-        topology = nx.Graph()
-        atoms = []
-
-        for i, atom in enumerate(qc_json.symbols):
-            atoms.append(
-                Atom(
-                    atomic_number=Element().number(atom),
-                    atom_index=i,
-                    atom_name=f"{atom}{i}",
-                )
-            )
-            topology.add_node(i)
-
-        for bond in qc_json.connectivity:
-            topology.add_edge(*bond[:2])
-
-        coords = np.array(qc_json.geometry).reshape((len(atoms), 3)) * BOHR_TO_ANGS
-        atoms = atoms or None
-        return cls(name=None, rdkit_mol=None, coords=coords)
-
-    @classmethod
-    def from_xyz(cls, file_name: str) -> "ReadInput":
-        """
-        Internal xyz reader.
-        Extracts the coords of the molecule.
-        """
-
-        traj_molecules = []
-        coords = []
-
-        with open(file_name) as xyz_file:
-            lines = xyz_file.readlines()
-
-            n_atoms = float(lines[0])
-
-            for line in lines:
-                line = line.split()
-                # skip frame heading lines
-                if len(line) <= 1 or "Iteration" in line:
-                    continue
-
-                coords.append([float(line[1]), float(line[2]), float(line[3])])
-
-                if len(coords) == n_atoms:
-                    # we have collected the molecule now store the frame
-                    traj_molecules.append(np.array(coords))
-                    coords = []
-
-        coords = traj_molecules[0] if len(traj_molecules) == 1 else traj_molecules
-        return cls(coords=coords, name=None, rdkit_mol=None)
-
-
-class ReadInputProtein:
-    """
-    A class that specialises in reading Protein input files.
-    #TODO are we better doing this with openmm or another tool?
-    """
-
-    def __init__(
-        self,
-        atoms: List[Atom],
-        bonds: Optional[List[Bond]] = None,
-        coords: Optional[np.ndarray] = None,
-        pdb_names: Optional[List[str]] = None,
-        residues: Optional[List[str]] = None,
-        name: Optional[str] = None,
-    ):
-        self.atoms = atoms
-        self.bonds = bonds
-        self.coords = coords
-        self.name = name
-        self.residues = residues
-        self.pdb_names = pdb_names
-
-    @classmethod
-    def from_pdb(cls, file_name: str, name: Optional[str] = None):
-        """
-        Read the protein input pdb file.
-        :return:
-        """
-        with open(file_name, "r") as pdb:
-            lines = pdb.readlines()
-
-        coords = []
-        atoms = []
-        bonds = []
-        Residues = []
-        pdb_names = []
-
-        # atom counter used for graph node generation
-        atom_count = 0
-        for line in lines:
-            if "ATOM" in line or "HETATM" in line:
-                atomic_symbol = str(line[76:78])
-                atomic_symbol = re.sub("[0-9]+", "", atomic_symbol).strip()
-
-                # If the element column is missing from the pdb, extract the atomic_symbol from the atom name.
-                if not atomic_symbol:
-                    atomic_symbol = str(line.split()[2])
-                    atomic_symbol = re.sub("[0-9]+", "", atomic_symbol)
-
-                # now make sure we have a valid element
-                if atomic_symbol.lower() != "cl" and atomic_symbol.lower() != "br":
-                    atomic_symbol = atomic_symbol[0]
-
-                atom_name = f"{atomic_symbol}{atom_count}"
-                # TODO should we use a protein pdb package for this?
-                qube_atom = Atom(
-                    atomic_number=Element().number(atomic_symbol),
-                    atom_index=atom_count,
-                    atom_name=atom_name,
-                    formal_charge=0,
-                    aromatic=False,
-                )
-
-                atoms.append(qube_atom)
-
-                pdb_names.append(str(line.split()[2]))
-
-                # also get the residue order from the pdb file so we can rewrite the file
-                Residues.append(str(line.split()[3]))
-
-                atom_count += 1
-                coords.append(
-                    [float(line[30:38]), float(line[38:46]), float(line[46:54])]
-                )
-
-            elif "CONECT" in line:
-                conect_terms = line.split()
-                for atom in conect_terms[2:]:
-                    if int(atom):
-                        bond = Bond(
-                            atom1_index=int(conect_terms[1]) - 1,
-                            atom2_index=int(atom) - 1,
-                            bond_order=1,
-                            aromatic=False,
-                        )
-                        bonds.append(bond)
-
-        coords = np.array(coords)
-        residues = [res for res, group in groupby(Residues)]
-        if name is None:
-            name = Path(file_name).stem
-        return cls(
-            atoms=atoms,
-            bonds=bonds,
-            coords=coords,
-            pdb_names=pdb_names,
-            residues=residues,
-            name=name,
-        )
+from QUBEKit.molecules import ExtraSite
+from QUBEKit.utils.constants import ANGS_TO_NM
+from QUBEKit.utils.datastructures import CustomNamespace
 
 
 class ExtractChargeData:
diff --git a/QUBEKit/virtual_sites.py b/QUBEKit/virtual_sites.py
index 56274fda..3a827ef3 100644
--- a/QUBEKit/virtual_sites.py
+++ b/QUBEKit/virtual_sites.py
@@ -18,7 +18,7 @@
 from mpl_toolkits.mplot3d import Axes3D
 from scipy.optimize import minimize
 
-from QUBEKit.ligand import Ligand
+from QUBEKit.molecules import ExtraSite, Ligand
 from QUBEKit.utils.constants import (
     ANGS_TO_NM,
     BOHR_TO_ANGS,
@@ -28,7 +28,6 @@
     PI,
     VACUUM_PERMITTIVITY,
 )
-from QUBEKit.utils.datastructures import ExtraSite
 from QUBEKit.utils.helpers import append_to_log
 
 
diff --git a/devtools/conda-envs/basic.yaml b/devtools/conda-envs/basic.yaml
index 5baf26d7..acc4ff50 100644
--- a/devtools/conda-envs/basic.yaml
+++ b/devtools/conda-envs/basic.yaml
@@ -1,7 +1,7 @@
 name: test
 channels:
   - conda-forge
-  - omnia
+  - psi4/label/dev
 dependencies:
     # Core
   - python
@@ -11,11 +11,14 @@ dependencies:
   - openff-toolkit
   - rdkit
   - torsiondrive
+  - optking
   - numpy >=1.18.1
   - scipy
   - networkx >=2.4
   - matplotlib-base
   - qcengine >=0.18.0
+  - jinja2
   - typing-extensions
   - pytest
-  - pytest-cov
\ No newline at end of file
+  - pytest-cov
+  - torchani
\ No newline at end of file
diff --git a/devtools/conda-envs/everything.yaml b/devtools/conda-envs/everything.yaml
index 0a76b5a4..16185df9 100644
--- a/devtools/conda-envs/everything.yaml
+++ b/devtools/conda-envs/everything.yaml
@@ -1,21 +1,30 @@
 name: test
 channels:
   - conda-forge
-  - omnia
-  - psi4
+  - psi4/label/dev
 dependencies:
     # Core
   - python
+  - pip
+
+    # Testing
+  - openff-toolkit
   - rdkit
   - torsiondrive
-    # extras
-  - qcengine >=0.14.0
-  - openforcefield
-  - geometric
-  - psi4
+  - optking
   - numpy >=1.18.1
   - scipy
   - networkx >=2.4
   - matplotlib-base
-
-
+  - qcengine >=0.18.0
+  - jinja2
+  - typing-extensions
+  - pytest
+  - pytest-cov
+  # compute engines
+  - xtb-python
+  - torchani
+  - psi4 =1.4a2.dev1058+670a850
+  - gau2grid =2.0.3
+  - dftd3
+  - gcp
diff --git a/devtools/design_doc.md b/devtools/design_doc.md
index f13a7548..d49dffeb 100644
--- a/devtools/design_doc.md
+++ b/devtools/design_doc.md
@@ -10,14 +10,16 @@ Whenever a class or method takes an argument: `molecule`, this is what is meant.
 
 The molecule object is either a ligand or a protein, the code for which is stored in the `ligand.py` file.
 A molecule object is initialised by simply passing a filename, smiles string or json dict like so:
-    
-    from QUBEKit.ligand import Ligand
-    
-    methane = Ligand('methane.pdb')
-    ethane = Ligand('CC')
-    ethane_with_name = Ligand(['CC', 'ethane'])
-    propane = Ligand('propane.mol2')
-    ...
+```python
+from QUBEKit.molecules import Ligand
+
+methane = Ligand('methane.pdb')
+ethane = Ligand('CC')
+ethane_with_name = Ligand(['CC', 'ethane'])
+propane = Ligand('propane.mol2')
+...
+```
+
     
 Upon initialisation, QUBEKit automatically builds the topology for the object, infers charge, multiplicity
 and calculates all the bonds, angles, dihedrals etc within the molecule.
@@ -115,12 +117,11 @@ Say, for example, a user would like to apply the modified Seminario method to a
 All that would be needed is some input file containing the coordinates, and the hessian matrix itself.
 
 ```python
-from QUBEKit.ligand import Ligand
+from QUBEKit.molecules.ligand import Ligand
 from QUBEKit.mod_seminario import ModSeminario
 
 import numpy as np
 
-
 # Initialise molecule
 molecule = Ligand('molecule.pdb')
 # Load hessian from file


Test patch:
diff --git a/QUBEKit/tests/engines/__init__.py b/QUBEKit/tests/engines/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/QUBEKit/tests/engines/gaussian_harness_test.py b/QUBEKit/tests/engines/gaussian_harness_test.py
new file mode 100644
index 00000000..0214ea50
--- /dev/null
+++ b/QUBEKit/tests/engines/gaussian_harness_test.py
@@ -0,0 +1,232 @@
+"""
+Tests for the qcengine gaussian harness
+"""
+
+import numpy as np
+import pytest
+import qcelemental as qcel
+import qcengine as qcng
+
+from QUBEKit.engines.gaussian_harness import GaussianHarness
+from QUBEKit.molecules import Ligand
+from QUBEKit.utils.file_handling import get_data
+
+
+def test_gaussian_version():
+    """
+    Try and check the gaussian version this will raise an error if not installed.
+    """
+    g = GaussianHarness()
+    if g.found():
+        assert g.get_version() in ["g09", "g16"]
+    else:
+        with pytest.raises(ModuleNotFoundError):
+            _ = g.get_version()
+
+
+@pytest.mark.parametrize(
+    "function, result",
+    [
+        pytest.param("pbe", "PBEPBE", id="pbe"),
+        pytest.param("wb97x-d", "wB97XD", id="wb97xd"),
+        pytest.param("b3lyp", "b3lyp", id="b3lyp"),
+    ],
+)
+def test_functional_converter(function, result):
+    """
+    Make sure that problem functions are converted.
+    """
+    theory = GaussianHarness.functional_converter(method=function)
+    assert theory == result
+
+
+@pytest.mark.parametrize(
+    "functional, result",
+    [
+        pytest.param("b3lyp-d3bj", "EmpiricalDispersion=GD3BJ b3lyp", id="b3lyp-d3bj"),
+        pytest.param("pbe-d3", "EmpiricalDispersion=GD3 PBEPBE", id="pbe-d3"),
+    ],
+)
+def test_dispersion_converter(functional, result):
+    """
+    Make sure that dispersion corrections are converted.
+    """
+    theory = GaussianHarness.functional_converter(method=functional)
+    assert result == theory
+
+
+@pytest.mark.parametrize(
+    "driver, result",
+    [
+        pytest.param("energy", "SP", id="Energy"),
+        pytest.param("gradient", "Force=NoStep", id="gradient"),
+        pytest.param("hessian", "FREQ", id="hessian"),
+    ],
+)
+def test_driver_conversion(driver, result):
+    """
+    Make sure each of the drivers are converted correctly.
+    """
+    g_driver = GaussianHarness.driver_conversion(driver=driver)
+    assert g_driver == result
+
+
+@pytest.mark.parametrize(
+    "driver, result",
+    [
+        pytest.param("energy", "SP", id="energy"),
+        pytest.param("gradient", "Force=NoStep", id="gradient"),
+        pytest.param("hessian", "FREQ", id="hessian"),
+    ],
+)
+def test_build_input(tmpdir, driver, result):
+    """
+    Build a gaussian input file for the given input model.
+    """
+    with tmpdir.as_cwd():
+        mol = Ligand.from_file(file_name=get_data("acetone.pdb"))
+        # build the atomic model
+        qc_spec = qcel.models.common_models.Model(method="pbe", basis="6-31G")
+        # build a job for s specific driver
+        qc_task = qcel.models.AtomicInput(
+            molecule=mol.to_qcschema(), driver=driver, model=qc_spec
+        )
+        g = GaussianHarness()
+        local_options = qcng.config.TaskConfig(
+            ncores=10, memory=10, nnodes=1, retries=1
+        )
+        input_data = g.build_input(input_model=qc_task, config=local_options)
+        file_input = input_data["infiles"]["gaussian.com"]
+        # now make sure the driver conversion is in the file input
+        assert result in file_input
+        # make sure we can write the file with no issues
+        with open("gaussian.com", "w") as com:
+            com.write(file_input)
+
+
+def test_get_version_from_log():
+    """
+    Try and extract the version info from a mock log file.
+    """
+    # mock output file with version stamp
+    with open(get_data("gaussian.log")) as log:
+        logfile = log.read()
+
+    version = GaussianHarness.parse_version(logfile=logfile)
+    assert version == "Gaussian 09:  ES64L-G09RevD.01 24-Apr-2013"
+
+
+def test_parse_gradient():
+    """
+    Make sure we can parse the gradient from a fchk file.
+    """
+    with open(get_data("gaussian.fchk")) as fchk:
+        fchkfile = fchk.read()
+
+    gradient = GaussianHarness.parse_gradient(fchfile=fchkfile)
+    assert len(gradient) == 30
+    # check the first and last values
+    assert gradient[0] == -1.01049212e-02
+    assert gradient[-1] == -7.31912443e-03
+
+
+def test_parse_hessian():
+    """
+    Make sure we can parse the hessian from a g09 mock fchk file.
+    """
+    with open(get_data("gaussian.fchk")) as fchk:
+        fchkfile = fchk.read()
+
+    hessian = GaussianHarness.parse_hessian(fchkfile=fchkfile)
+    assert len(hessian) == 900
+
+
+@pytest.mark.parametrize(
+    "driver",
+    [
+        pytest.param("energy", id="energy"),
+        pytest.param("gradient", id="gradient"),
+        pytest.param("hessian", id="hessian"),
+    ],
+)
+def test_parse_output(driver):
+    """
+    Test reading gaussian outfiles and extracting the correct information based on the
+    driver type.
+    """
+    outfiles = {}
+    with open(get_data("gaussian.log")) as log:
+        outfiles["gaussian.log"] = log.read()
+    with open(get_data("gaussian.fchk")) as fchk:
+        outfiles["lig.fchk"] = fchk.read()
+
+    # build the input
+    mol = Ligand.from_file(file_name=get_data("acetone.pdb"))
+    # build the atomic model
+    qc_spec = qcel.models.common_models.Model(method="pbe", basis="6-31G")
+    # build a job for a specific driver
+    qc_task = qcel.models.AtomicInput(
+        molecule=mol.to_qcschema(), driver=driver, model=qc_spec
+    )
+    g = GaussianHarness()
+    result = g.parse_output(outfiles=outfiles, input_model=qc_task)
+    if driver == "energy":
+        assert result.return_result == -1.931393770857046e02
+    elif driver == "gradient":
+        assert result.return_result.shape == (10, 3)
+    elif driver == "hessian":
+        assert result.return_result.shape == (30, 30)
+
+
+def test_fail_termination():
+    """
+    If we think the gaussian job did not finish then we need to make sure an error is raised.
+    """
+    random_string = "test string\n test string"
+    with pytest.raises(qcng.exceptions.UnknownError):
+        GaussianHarness.check_convergence(logfile=random_string)
+
+
+@pytest.mark.parametrize(
+    "driver",
+    [
+        pytest.param("energy", id="energy"),
+        pytest.param("gradient", id="gradient"),
+        pytest.param("hessian", id="hessian"),
+    ],
+)
+def test_full_run(driver, tmpdir):
+    """
+    For the given driver try a full execution if the user has gaussian installed.
+    """
+    if not GaussianHarness.found():
+        pytest.skip("Gaussian 09/16 not available test skipped.")
+
+    with tmpdir.as_cwd():
+        # build the input
+        mol = Ligand.from_file(file_name=get_data("acetone.pdb"))
+        # build the atomic model
+        qc_spec = qcel.models.common_models.Model(
+            method="wB97XD", basis="6-311++G(d,p)"
+        )
+        # build a job for a specific driver
+        qc_task = qcel.models.AtomicInput(
+            molecule=mol.to_qcschema(), driver=driver, model=qc_spec
+        )
+        g = GaussianHarness()
+        # run locally with 2 cores and 2 GB memory
+        result = g.compute(
+            input_data=qc_task,
+            config=qcng.config.TaskConfig(
+                **{"memory": 2, "ncores": 2, "nnodes": 1, "retries": 1}
+            ),
+        )
+
+        outfiles = {}
+        with open(get_data("gaussian.log")) as log:
+            outfiles["gaussian.log"] = log.read()
+        with open(get_data("gaussian.fchk")) as fchk:
+            outfiles["lig.fchk"] = fchk.read()
+        ref_result = g.parse_output(outfiles=outfiles, input_model=qc_task)
+
+        assert np.allclose(ref_result.return_result, result.return_result)
diff --git a/QUBEKit/tests/engines/geometry_optimiser_test.py b/QUBEKit/tests/engines/geometry_optimiser_test.py
new file mode 100644
index 00000000..d0c9b005
--- /dev/null
+++ b/QUBEKit/tests/engines/geometry_optimiser_test.py
@@ -0,0 +1,179 @@
+import os
+
+import pytest
+import qcengine
+
+from QUBEKit.engines import GeometryOptimiser
+from QUBEKit.molecules import Ligand
+from QUBEKit.utils.exceptions import SpecificationError
+from QUBEKit.utils.file_handling import get_data
+
+
+def test_rdkit_available():
+    """
+    Make sure the geometry optimiser allows rdkit as this comes with QUBEKit.
+    """
+    g_opt = GeometryOptimiser()
+    g_opt.program = "rdkit"
+
+
+def test_program_not_installed():
+    """
+    Make sure an error is raised when we try and use a program that is not available.
+    """
+    g_opt = GeometryOptimiser()
+    with pytest.raises(SpecificationError):
+        g_opt.program = "test"
+
+
+def test_local_options():
+    """
+    Make sure the task config is correctly made and keywords are converted.
+    """
+    g_opt = GeometryOptimiser(cores=10, memory=2)
+    local_options = g_opt.local_options
+    assert g_opt.cores == local_options["ncores"]
+    assert g_opt.memory == local_options["memory"]
+
+
+@pytest.mark.parametrize(
+    "optimiser",
+    [
+        pytest.param("geometric", id="geometric"),
+        pytest.param("optking", id="optking"),
+    ],
+)
+def test_optimiser_keywords(optimiser):
+    """
+    For the given optimiser make sure the keywords are updated correctly.
+    """
+    if "psi4" not in qcengine.list_available_programs():
+        pytest.skip("Psi4 missing skipping.")
+    g = GeometryOptimiser(
+        optimiser=optimiser, maxiter=1, convergence="GAU", program="psi4"
+    )
+    keywords = g.build_optimiser_keywords()
+    assert 1 in keywords.values()
+    assert "GAU" in keywords.values()
+
+
+def test_missing_optimiser():
+    """
+    Make sure an error is raised when we try and set a missing optimiser.
+    """
+    with pytest.raises(SpecificationError):
+        _ = GeometryOptimiser(optimiser="bad_optimiser")
+
+
+@pytest.mark.parametrize(
+    "program, basis, method",
+    [
+        pytest.param("rdkit", None, "UFF", id="rdkit uff"),
+        pytest.param("torchani", None, "ani2x", id="ani2x"),
+        pytest.param("psi4", "6-311G", "b3lyp", id="psi4 b3lyp"),
+        pytest.param("openmm", "smirnoff", "openff_unconstrained-1.3.0", id="openff"),
+        pytest.param("openmm", "antechamber", "gaff-2.11", id="gaff"),
+    ],
+)
+def test_spec_validation_pass(program, basis, method):
+    """
+    Make sure we can correctly validate a specification.
+    """
+    if program not in qcengine.list_available_programs():
+        pytest.skip(f"{program} missing skipping test")
+    _ = GeometryOptimiser(program=program, basis=basis, method=method)
+
+
+@pytest.mark.parametrize(
+    "program, basis, method",
+    [
+        pytest.param("rdkit", "uff", "uff", id="rdkit wrong basis"),
+        pytest.param("openmm", "gaff", "gaff-2.11", id="openmm wrong basis"),
+        pytest.param("xtb", "gfn", "gfn2xtb", id="xtb wrong basis"),
+        pytest.param("torchani", None, "ani3y", id="torchani wrong method"),
+        pytest.param("openmm", "smirnoff", "openff-sage", id="openff wrong method."),
+    ],
+)
+def test_spec_validation_fail(program, basis, method):
+    """
+    Make sure than an invalid specification raises an error.
+    """
+    if program not in qcengine.list_available_programs():
+        pytest.skip(f"{program} missing skipping test")
+    with pytest.raises(SpecificationError):
+        _ = GeometryOptimiser(program=program, basis=basis, method=method)
+
+
+@pytest.mark.parametrize(
+    "program, basis, method",
+    [
+        pytest.param("rdkit", None, "uff", id="rdkit uff"),
+        pytest.param("openmm", "antechamber", "gaff-2.11", id="gaff-2.11"),
+        pytest.param(
+            "openmm", "smirnoff", "openff_unconstrained-1.3.0", id="openff 1.3.0"
+        ),
+        pytest.param("xtb", None, "gfn2xtb", id="xtb gfn2xtb"),
+        pytest.param("torchani", None, "ani1ccx", id="torchani ccx"),
+        pytest.param("psi4", "3-21g", "hf", id="psi4 hf"),
+        pytest.param("gaussian", "3-21g", "hf", id="gaussian hf"),
+    ],
+)
+def test_optimise(program, basis, method, tmpdir):
+    """
+    Test running different optimisers with different programs.
+    """
+    if program not in qcengine.list_available_programs():
+        pytest.skip(f"{program} missing skipping test.")
+
+    with tmpdir.as_cwd():
+        mol = Ligand.from_file(get_data("water.pdb"))
+        g = GeometryOptimiser(
+            program=program,
+            basis=basis,
+            method=method,
+            optimiser="geometric",
+            convergence="GAU",
+        )
+        result_mol = g.optimise(molecule=mol)
+        assert result_mol.coordinates.tolist() != mol.coordinates.tolist()
+
+
+def test_optimise_fail_output(tmpdir):
+    """
+    Make sure the optimised geometries and result is still wrote out if we fail the molecule and an error is rasied.
+    """
+    with tmpdir.as_cwd():
+        mol = Ligand.from_file(file_name=get_data("water.pdb"))
+        g = GeometryOptimiser(
+            program="torchani", method="ani1ccx", basis=None, maxiter=5
+        )
+        with pytest.raises(RuntimeError):
+            g.optimise(molecule=mol, allow_fail=False)
+        files = os.listdir()
+        assert "opt.xyz" in files
+        assert "opt_trajectory.xyz" in files
+        assert "result.json" in files
+
+
+def test_optmiser_fail_no_output(tmpdir):
+    """
+    Make sure we raise an error correctly when there is no output from a failed optimisation.
+    """
+    if "psi4" not in qcengine.list_available_programs():
+        pytest.skip("Psi4 missing skipping test.")
+    with tmpdir.as_cwd():
+        mol = Ligand.from_file(file_name=get_data("water.pdb"))
+        g = GeometryOptimiser(
+            program="psi4", method="wb97x-dbj", basis="dzvp", maxiter=10
+        )
+        with pytest.raises(RuntimeError):
+            g.optimise(molecule=mol, allow_fail=False)
+
+
+def test_optking_fail():
+    """
+    Optking currently only works with psi4 make sure we raise an error if we use
+    a different program.
+    """
+    with pytest.raises(SpecificationError):
+        GeometryOptimiser(optimiser="optking", program="rdkit")
diff --git a/QUBEKit/tests/engines/qcengine_test.py b/QUBEKit/tests/engines/qcengine_test.py
new file mode 100644
index 00000000..7c0fed02
--- /dev/null
+++ b/QUBEKit/tests/engines/qcengine_test.py
@@ -0,0 +1,42 @@
+import pytest
+import qcengine
+
+from QUBEKit.engines import QCEngine
+from QUBEKit.molecules import Ligand
+from QUBEKit.utils.file_handling import get_data
+
+
+@pytest.mark.parametrize(
+    "program, basis, method",
+    [
+        pytest.param("rdkit", None, "mmff94", id="rdkit mmff"),
+        pytest.param("xtb", None, "gfn2xtb", id="gfn2xtb"),
+        pytest.param("openmm", "smirnoff", "openff-1.0.0.offxml", id="parsley"),
+        pytest.param("openmm", "antechamber", "gaff-2.11", id="gaff-2.11"),
+        pytest.param("torchani", None, "ani2x", id="ani2x"),
+        pytest.param("psi4", "3-21g", "hf", id="psi4 hf"),
+        pytest.param("gaussian", "3-21g", "hf", id="gaussian hf"),
+    ],
+)
+def test_single_point_energy(program, basis, method, tmpdir):
+    """
+    Make sure our qcengine wrapper works correctly.
+    """
+    if program not in qcengine.list_available_programs():
+        pytest.skip(f"{program} missing skipping test.")
+
+    with tmpdir.as_cwd():
+        mol = Ligand.from_file(file_name=get_data("water.pdb"))
+        engine = QCEngine(
+            program=program,
+            basis=basis,
+            method=method,
+            memory=1,
+            cores=1,
+            driver="energy",
+        )
+        result = engine.call_qcengine(molecule=mol)
+        assert result.driver == "energy"
+        assert result.model.basis == basis
+        assert result.model.method == method
+        assert result.provenance.creator.lower() == program
diff --git a/QUBEKit/tests/ligand/ligand_test.py b/QUBEKit/tests/ligand/ligand_test.py
index 09b88137..74d0fec9 100644
--- a/QUBEKit/tests/ligand/ligand_test.py
+++ b/QUBEKit/tests/ligand/ligand_test.py
@@ -5,8 +5,13 @@
 from rdkit.Chem import rdMolTransforms
 from simtk import unit
 
-from QUBEKit.ligand import Ligand
-from QUBEKit.utils.exceptions import FileTypeError, SmartsError, TopologyMismatch
+from QUBEKit.molecules import Ligand
+from QUBEKit.utils.exceptions import (
+    ConformerError,
+    FileTypeError,
+    SmartsError,
+    TopologyMismatch,
+)
 from QUBEKit.utils.file_handling import get_data
 from QUBEKit.utils.helpers import unpickle
 
@@ -550,6 +555,17 @@ def test_to_rdkit(molecule):
                 assert bond.GetStereo() == Chem.BondStereo.STEREOZ
 
 
+def test_to_rdkit_complicated_stereo():
+    """
+    Make sure we can convert a complicated molecule with multiple stereo centres to rdkit.
+    """
+    mol = Ligand.from_smiles(
+        "[H][C@]1([C@@]([C@](O[C@@]1([H])C([H])([H])OP(=O)(O[H])O[H])([H])N2C(=C([N+](C2([H])[H])([H])[H])C(=O)N([H])[H])O[H])([H])O[H])O[H]",
+        name="complicated",
+    )
+    mol.to_rdkit()
+
+
 @pytest.mark.parametrize(
     "molecule, charge",
     [
@@ -658,3 +674,24 @@ def test_get_bond_error(acetone):
     """
     with pytest.raises(TopologyMismatch):
         acetone.get_bond_between(atom1_index=4, atom2_index=9)
+
+
+def test_to_qcschema_no_conformer(acetone):
+    """
+    Make sure we raise an error when trying to make a qcelemental molecule with no coordinates.
+    """
+    with pytest.raises(ConformerError):
+        acetone.coordinates = None
+        acetone.to_qcschema()
+
+
+def test_to_qcschema(acetone):
+    """
+    Make sure we can convert to a valid qcelemental molecule.
+    """
+    qcel_mol = acetone.to_qcschema()
+    assert qcel_mol.atomic_numbers.tolist() == [
+        atom.atomic_number for atom in acetone.atoms
+    ]
+    assert len(qcel_mol.symbols) == acetone.n_atoms
+    assert qcel_mol.geometry.shape == (acetone.n_atoms, 3)
diff --git a/QUBEKit/tests/ligand/parametrisation_test.py b/QUBEKit/tests/ligand/parametrisation_test.py
index ea0092db..ed2b557e 100644
--- a/QUBEKit/tests/ligand/parametrisation_test.py
+++ b/QUBEKit/tests/ligand/parametrisation_test.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from QUBEKit.ligand import Ligand
+from QUBEKit.molecules import Ligand
 from QUBEKit.parametrisation import XML, AnteChamber, OpenFF, Parametrisation
 from QUBEKit.utils.file_handling import get_data
 
diff --git a/QUBEKit/tests/protein/parametrisation_test.py b/QUBEKit/tests/protein/parametrisation_test.py
index b9e11811..bfd4a99f 100644
--- a/QUBEKit/tests/protein/parametrisation_test.py
+++ b/QUBEKit/tests/protein/parametrisation_test.py
@@ -2,7 +2,7 @@
 
 import pytest
 
-from QUBEKit.ligand import Protein
+from QUBEKit.molecules import Protein
 from QUBEKit.parametrisation import XMLProtein
 from QUBEKit.utils.file_handling import get_data
 
diff --git a/QUBEKit/tests/protein/protein_test.py b/QUBEKit/tests/protein/protein_test.py
index eabea1df..4271ddd3 100644
--- a/QUBEKit/tests/protein/protein_test.py
+++ b/QUBEKit/tests/protein/protein_test.py
@@ -1,6 +1,6 @@
 import pytest
 
-from QUBEKit.ligand import Protein
+from QUBEKit.molecules import Protein
 from QUBEKit.utils.exceptions import FileTypeError
 from QUBEKit.utils.file_handling import get_data
 
diff --git a/QUBEKit/tests/run_test.py b/QUBEKit/tests/run_test.py
index 22e8ea47..a47014ea 100644
--- a/QUBEKit/tests/run_test.py
+++ b/QUBEKit/tests/run_test.py
@@ -6,7 +6,7 @@
 
 import pytest
 
-from QUBEKit.ligand import Ligand
+from QUBEKit.molecules import Ligand
 from QUBEKit.run import Execute
 from QUBEKit.utils.file_handling import get_data
 
diff --git a/QUBEKit/tests/seminario_test.py b/QUBEKit/tests/seminario_test.py
index 34c58cc1..ae2798eb 100644
--- a/QUBEKit/tests/seminario_test.py
+++ b/QUBEKit/tests/seminario_test.py
@@ -5,8 +5,8 @@
 import numpy as np
 import pytest
 
-from QUBEKit.ligand import Ligand
 from QUBEKit.mod_seminario import ModSeminario
+from QUBEKit.molecules import Ligand
 from QUBEKit.utils.file_handling import get_data
 
 
diff --git a/QUBEKit/tests/virtual_sites_test.py b/QUBEKit/tests/virtual_sites_test.py
index b008e31d..0cb0263c 100644
--- a/QUBEKit/tests/virtual_sites_test.py
+++ b/QUBEKit/tests/virtual_sites_test.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from QUBEKit.ligand import Ligand
+from QUBEKit.molecules import Ligand
 from QUBEKit.utils.constants import BOHR_TO_ANGS
 from QUBEKit.utils.datastructures import CustomNamespace
 from QUBEKit.utils.file_handling import get_data