diff --git a/H-atom-4c.ipynb b/H-atom-4c.ipynb
new file mode 100644
index 0000000..8f9822b
--- /dev/null
+++ b/H-atom-4c.ipynb
@@ -0,0 +1,299 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "90e862bb-783a-4f97-a227-8f6e0b4012b8",
+   "metadata": {},
+   "source": [
+    "# The Dirac equation for the Hydrogen atom\n",
+    "\n",
+    "In this notebook we will illustrate how one can solve the Dirac equation for a Hydrogen atom using the Multiwavelet framework provided by VAMPyR\n",
+    "\n",
+    "The Dirac equation can be coincisely written as follows:\n",
+    "\n",
+    "\\begin{equation}\n",
+    "(\\beta m c^2+ c \\boldsymbol{\\alpha} \\cdot \\mathbf{p} + V) \\phi = \\epsilon \\phi \n",
+    "\\end{equation}\n",
+    "\n",
+    "where $\\phi = (\\phi^{L\\alpha}, \\phi^{L\\beta}, \\phi^{S\\alpha}, \\phi^{S\\beta})^t$ represents a 4-component spinor, $\\boldsymbol{\\alpha} = \n",
+    "\\begin{pmatrix}\n",
+    "0_2 & \\boldsymbol{\\sigma} \\\\\n",
+    "\\boldsymbol{\\sigma} & 0_2 & \n",
+    "\\end{pmatrix}\n",
+    "$ and\n",
+    "$\\beta = \n",
+    "\\begin{pmatrix}\n",
+    "1_2 & 0_2 \\\\\n",
+    "0_2 & -1_2\n",
+    "\\end{pmatrix}\n",
+    "$ are the 4x4 Dirac matrices, $\\boldsymbol{\\sigma}$ is a cartesian vector collecting the three 2x2 Pauli matrices, $\\mathbf{p}$ is the momentum operator, $c$ is the speed of light, $m$ is the electron mass and $V$ is the nuclear potential.\n",
+    "\n",
+    "As for the non-relativistic case,  the equation is first rewritten in its integral formulation:\n",
+    "$$\\phi = \\frac{1}{2mc^2}(\\beta m c^2+ c \\boldsymbol{\\alpha} \\cdot \\mathbf{p} + \\epsilon) \\left[ G_\\mu \\star (V \\psi) \\right]$$\n",
+    "\n",
+    "where $G_\\mu(x) = \\frac{e^{-\\mu |x|}}{4 \\pi |x|}$ is the Helmholtz Green's kernel and $\\mu = \\sqrt{\\frac{m^2c^4-\\epsilon}{mc^2}}$. An initial guess can be obtained by taking a Slater orbital or a Gaussian function for the $\\psi^{L\\alpha}$ component and then applying the restricted kinetic balance:\n",
+    "\n",
+    "$$\n",
+    "\\begin{pmatrix}\n",
+    "\\phi^{S\\alpha} \\\\\n",
+    "\\phi^{S\\beta}\n",
+    "\\end{pmatrix}\n",
+    "= \\frac{1}{2c}\\boldsymbol{\\sigma} \\cdot \\mathbf{p} \n",
+    "\\begin{pmatrix}\n",
+    "\\phi^{L\\alpha} \\\\\n",
+    "0\n",
+    "\\end{pmatrix}\n",
+    "$$\n",
+    "The guess obtained is then plugged into the integral form of the Dirac equation, which can then be iterated until convergence"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "54c75006-0358-42c5-9075-37a5b8a909ba",
+   "metadata": {},
+   "source": [
+    "We start by loading the relevant packages: the 3d version of `vampyr`, `numpy`, the `complex_function` class which deals with complex funtions and the `orbital` class which deals with 4-component spinors. Each complex function is handled as a pair of `function_tree`s and each spinor is handled as a 4c vector of complex functions. The `nuclear_potential` package is self-explanatory"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "37e0e6b2-886e-4415-8612-d2a652f24c4f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from vampyr import vampyr3d as vp\n",
+    "from orbital4c import orbital as orb\n",
+    "from orbital4c import nuclear_potential as nucpot\n",
+    "from orbital4c import complex_fcn as cf\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7fbe4360-9936-47d4-bb6f-d51bb3b34b69",
+   "metadata": {},
+   "source": [
+    "As a reference, the exact Dirac energy for the ground state Hydrogen atom is computed in the function below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "0ef3d827-810b-4411-8ba3-bd63271033e4",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Exact Energy -0.5000066565989982\n"
+     ]
+    }
+   ],
+   "source": [
+    "def analytic_1s(light_speed, n, k, Z):\n",
+    "    alpha = 1/light_speed\n",
+    "    gamma = orb.compute_gamma(k,Z,alpha)\n",
+    "    tmp1 = n - np.abs(k) + gamma\n",
+    "    tmp2 = Z * alpha / tmp1\n",
+    "    tmp3 = 1 + tmp2**2\n",
+    "    return light_speed**2 / np.sqrt(tmp3)\n",
+    "\n",
+    "light_speed = 137.03599913900001\n",
+    "alpha = 1/light_speed\n",
+    "k = -1\n",
+    "l = 0\n",
+    "n = 1\n",
+    "m = 0.5\n",
+    "Z = 1\n",
+    "atom = \"H\"\n",
+    "\n",
+    "energy_1s = analytic_1s(light_speed, n, k, Z)\n",
+    "print('Exact Energy',energy_1s - light_speed**2, flush = True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6277d776-2236-496a-bdd5-41929d11ac3d",
+   "metadata": {},
+   "source": [
+    "The `MultiResolutionAnalysis` object defining the simulation box is constructed and passed to the classes for complex functions and spinors"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "id": "5ac41574-6f43-411d-841f-f8fdf0749b66",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mra = vp.MultiResolutionAnalysis(box=[-30,30], order=6)\n",
+    "prec = 1.0e-4\n",
+    "origin = [0.1, 0.2, 0.3]  # origin moved to avoid placing the nuclar charge on a node\n",
+    "\n",
+    "orb.orbital4c.light_speed = light_speed\n",
+    "orb.orbital4c.mra = mra\n",
+    "cf.complex_fcn.mra = mra"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d5abdbe4-ba5b-49d3-ae7e-febab5ced4fe",
+   "metadata": {},
+   "source": [
+    "We construct a starting guess by taking a simple Gaussian function and initialize the real part of the $\\phi^{L\\alpha}$ component of the spinor with it. Thereafter the restricted kinetic balance is employed. This is implemented in the `init_small_components` method of the `orbital` class"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "id": "95826623-b759-4318-afc5-e0c225cb5f1d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "a_coeff = 3.0\n",
+    "b_coeff = np.sqrt(a_coeff/np.pi)**3\n",
+    "gauss = vp.GaussFunc(b_coeff, a_coeff, origin)\n",
+    "gauss_tree = vp.FunctionTree(mra)\n",
+    "vp.advanced.build_grid(out=gauss_tree, inp=gauss)\n",
+    "vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss)\n",
+    "gauss_tree.normalize()\n",
+    "\n",
+    "spinor_H = orb.orbital4c()\n",
+    "La_comp = cf.complex_fcn()\n",
+    "La_comp.copy_fcns(real = gauss_tree)\n",
+    "spinor_H.copy_components(La = La_comp)\n",
+    "spinor_H.init_small_components(prec/10)\n",
+    "spinor_H.normalize()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f6a40ddd-7933-46ee-8128-96caf89f3682",
+   "metadata": {},
+   "source": [
+    "The nuclear potential is defined and projected onto the `V_tree`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "id": "4f0cda6f-217c-4941-8ce5-fa52576c2d98",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Peps = vp.ScalingProjector(mra, prec)\n",
+    "f = lambda x: nucpot.coulomb_HFYGB(x, origin, Z, prec)\n",
+    "V_tree = Peps(f)\n",
+    "\n",
+    "default_der = 'BS'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1678bc97-d85d-4a48-9aea-13901631534d",
+   "metadata": {},
+   "source": [
+    "The orbital is optimized by iterating the integral version of the Dirac equation as follows:\n",
+    "1. Application of the Dirac Hamiltonian $f^n = \\hat{h}_D \\phi^n = (\\beta m c^2+ c \\boldsymbol{\\alpha} \\cdot \\mathbf{p}) \\phi^n$\n",
+    "2. Application of the potnetial operator $g^n = \\hat{V} \\phi^n$\n",
+    "3. Sum $h^n = f^n + g^n$\n",
+    "4. Expectation value of the energy $\\left\\langle \\phi^n | h^n \\right\\rangle$\n",
+    "5. Calculation of the Helmholtz parameter $\\mu$\n",
+    "6. Convolution with the Helmholtz kernel $t^n = G_\\mu \\star g^n$\n",
+    "7. application of the shifted Dirac Hamiltonian $\\tilde{\\phi}^{n+1} = (\\hat{h}_D + \\epsilon) t^n$\n",
+    "8. normalization of the new iterate\n",
+    "9. calculation of the change in the orbital $\\delta \\phi^n = \\phi^{n+1}-\\phi^n$\n",
+    "\n",
+    "Once the orbital error is below the requested threshold the iteration is interrupted and the final energy is computed."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "id": "40a443b6-33f4-4771-8b06-ccd12c9da35d",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Energy -0.14180720307558659\n",
+      "Error 0.3648353655184581\n",
+      "Energy -0.49611934473068686\n",
+      "Error 0.0024839139737050653\n",
+      "Energy -0.4997319277244969\n",
+      "Error 0.00023304055836475225\n",
+      "Energy -0.49998214000152075\n",
+      "Error 2.3987883028320866e-05\n",
+      "Final Energy -0.5000042447172746\n",
+      "Exact Energy -0.5000066565989982\n",
+      "Difference -2.411881723674014e-06\n"
+     ]
+    }
+   ],
+   "source": [
+    "orbital_error = 1\n",
+    "while orbital_error > prec:\n",
+    "    # 1. Application of the Dirac Hamiltonian\n",
+    "    hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)\n",
+    "    # 2. Application of the potnetial operator\n",
+    "    v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)\n",
+    "    # 3. Sum\n",
+    "    add_psi = hd_psi + v_psi\n",
+    "    # 4. Expectation value of the energy\n",
+    "    energy = (spinor_H.dot(add_psi)).real\n",
+    "    print('Energy',energy-light_speed**2)\n",
+    "    # 5. Calculation of the Helmholtz parameter\n",
+    "    mu = orb.calc_dirac_mu(energy, light_speed)\n",
+    "    # 6. Convolution with the Helmholtz kernel\n",
+    "    tmp = orb.apply_helmholtz(v_psi, mu, prec)\n",
+    "    tmp.crop(prec/10)\n",
+    "    # 7. application of the shifted Dirac Hamiltonian\n",
+    "    new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = default_der) \n",
+    "    new_orbital.crop(prec/10)\n",
+    "    # 8. normalization of the new iterate\n",
+    "    new_orbital.normalize()\n",
+    "    delta_psi = new_orbital - spinor_H\n",
+    "    # 9. calculation of the change in the orbital\n",
+    "    orbital_error = (delta_psi.dot(delta_psi)).real\n",
+    "    print('Error',orbital_error, flush = True)\n",
+    "    spinor_H = new_orbital\n",
+    "\n",
+    "# Computing the final energy\n",
+    "hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)\n",
+    "v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)\n",
+    "add_psi = hd_psi + v_psi\n",
+    "energy = (spinor_H.dot(add_psi)).real\n",
+    "print('Final Energy',energy - light_speed**2)\n",
+    "\n",
+    "energy_1s = analytic_1s(light_speed, n, k, Z)\n",
+    "\n",
+    "print('Exact Energy',energy_1s - light_speed**2)\n",
+    "print('Difference',energy_1s - energy)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/orbital4c/__init__.py b/orbital4c/__init__.py
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/orbital4c/__init__.py
@@ -0,0 +1 @@
+
diff --git a/orbital4c/complex_fcn.py b/orbital4c/complex_fcn.py
new file mode 100644
index 0000000..6680b9d
--- /dev/null
+++ b/orbital4c/complex_fcn.py
@@ -0,0 +1,295 @@
+from vampyr import vampyr3d as vp
+import numpy as np
+
+class complex_fcn:
+    """Complex function trees as pairs of real and imaginary trees"""
+    mra = None
+    def __init__(self):
+        self.real = vp.FunctionTree(self.mra)
+        self.imag = vp.FunctionTree(self.mra)
+        self.setZero()
+        
+    def copy_fcns(self, real=None, imag=None):
+        if not real == None:
+            vp.advanced.copy_grid(self.real, real)
+            vp.advanced.copy_func(self.real, real)
+        if not imag == None:
+            vp.advanced.copy_grid(self.imag, imag)
+            vp.advanced.copy_func(self.imag, imag)
+            
+    def squaredNorm(self):
+        re = self.real.squaredNorm()
+        im = self.imag.squaredNorm()
+        return re + im
+    
+    def normalize(self):
+        norm = np.sqrt(self.squaredNorm())
+        factor = 1.0 / norm 
+        self.real *= factor
+        self.imag *= factor
+    
+    def setZero(self):
+        self.real.setZero()
+        self.imag.setZero()
+
+    def save(self, name):
+        self.real.saveTree(f"{name}_real")
+        self.imag.saveTree(f"{name}_imag")
+
+    def read(self, name):
+        self.real.loadTree(f"{name}_real")
+        self.imag.loadTree(f"{name}_imag")
+
+    def cropRealImag(self, prec):
+        norm = np.sqrt(self.squaredNorm())
+        self.crop(prec * norm, True)
+
+    def crop(self, prec, abs = False):
+        self.real.crop(prec, abs)
+        self.imag.crop(prec, abs)
+
+    def __call__(self, position):
+        rval = self.real(position)
+        ival = self.imag(position)
+        return rval + 1j * ival
+
+    def __add__(self, other):
+        output = complex_fcn()
+        output.real = self.real + other.real
+        output.imag = self.imag + other.imag
+        return output
+
+    def __sub__(self, other):
+        output = complex_fcn()
+        output.real = self.real - other.real
+        output.imag = self.imag - other.imag
+        return output
+
+    def real_mul(self, coeff):
+        output = complex_fcn()
+        output.real = self.real * coeff
+        output.imag = self.imag * coeff
+        return output
+
+    def imag_mul(self, coeff):
+        output = complex_fcn()
+        output.real = self.imag * (-1.0 * coeff)
+        output.imag = self.real * coeff
+        return output
+
+    def __rmul__(self, other):
+        output = complex_fcn()
+        output.real = self.real * np.real(other) - self.imag * np.imag(other)
+        output.imag = self.real * np.imag(other) + self.imag * np.real(other)
+        return output
+        
+    def __mul__(self, other):
+        output = complex_fcn()
+        output.real = self.real * np.real(other) - self.imag * np.imag(other)
+        output.imag = self.real * np.imag(other) + self.imag * np.real(other)
+        return output
+        
+    def __str__(self):
+        return ('Real part {}\n Imag part {}'.format(self.real, self.imag))
+    
+    def gradient(self, der = 'ABGV'):
+        if(der == 'ABGV'):
+            D = vp.ABGVDerivative(self.mra, 0.0, 0.0)
+        elif(der == 'PH'):
+            D = vp.PHDerivative(self.mra)
+        elif(der == 'BS'):
+            D = vp.BSDerivative(self.mra)
+        else:
+            exit("Derivative operator not found")
+        grad_re = vp.gradient(D, self.real)
+        grad_im = vp.gradient(D, self.imag)
+        grad = []
+        for i in range(3):
+            comp = complex_fcn()
+            comp.copy_fcns(real=grad_re[i], imag=grad_im[i])
+            grad.append(comp)
+        return grad
+
+    def derivative(self, dir = 0, der = 'ABGV'):
+        if(der == 'ABGV'):
+            D = vp.ABGVDerivative(self.mra, 0.0, 0.0)
+        elif(der == 'PH'):
+            D = vp.PHDerivative(self.mra)
+        elif(der == 'BS'):
+            D = vp.BSDerivative(self.mra)
+        else:
+            exit("Derivative operator not found")
+        re_der = D(self.real, dir)
+        im_der = D(self.imag, dir)
+        der_func = complex_fcn()
+        der_func.real = re_der
+        der_func.imag = im_der
+        return der_func
+
+    def complex_conj(self):
+        output = complex_fcn()
+        output.real = self.real 
+        output.imag = -1.0 * self.imag
+        return output
+       
+        
+    def density(self, prec):
+        density = vp.FunctionTree(self.mra)
+        add_vector = []
+        temp_r = vp.FunctionTree(self.mra)
+        temp_r.setZero()
+        temp_i = vp.FunctionTree(self.mra)
+        temp_i.setZero()
+        if(self.real.squaredNorm() > 0):
+            vp.advanced.multiply(prec, temp_r, 1.0, self.real, self.real)
+        if(self.imag.squaredNorm() > 0):
+            vp.advanced.multiply(prec, temp_i, 1.0, self.imag, self.imag)
+        vp.advanced.add(prec/10, density, [temp_r, temp_i])
+        return density
+
+    def dot(self, other, cc_first = True):
+        out_real = 0
+        out_imag = 0
+        func_a = self.real
+        func_b = self.imag
+        func_c = other.real
+        func_d = other.imag
+
+        fbd = 1.0
+        fbc = -1.0
+        if(not cc_first):
+            fbd = -1.0
+            fbc = 1.0
+        
+        if(func_a.squaredNorm() > 0 and func_c.squaredNorm() > 0):
+           out_real = out_real + vp.dot(func_a, func_c)
+        if(func_b.squaredNorm() > 0 and func_d.squaredNorm() > 0):
+           out_real = out_real + fbd * vp.dot(func_b, func_d)
+        if(func_a.squaredNorm() > 0 and func_d.squaredNorm() > 0):
+           out_imag = out_imag + vp.dot(func_a, func_d)
+        if(func_b.squaredNorm() > 0 and func_c.squaredNorm() > 0):
+           out_imag = out_imag + fbc * vp.dot(func_b, func_c)
+
+        return out_real + 1j * out_imag
+
+    def advanced_overlap_density(self, other, prec):
+        func_a = self.real
+        func_b = self.imag
+        func_c = other.real
+        func_d = other.imag
+
+        rr = vp.FunctionTree(self.mra)
+        ii = vp.FunctionTree(self.mra)
+        ri = vp.FunctionTree(self.mra)
+        ir = vp.FunctionTree(self.mra)
+
+        if(func_a.squaredNorm() > 0 and func_c.squaredNorm() > 0):
+            vp.advanced.multiply(prec, rr, 1.0, func_a, func_c)
+        if(func_b.squaredNorm() > 0 and func_d.squaredNorm() > 0):
+            vp.advanced.multiply(prec, ii, 1.0, func_b, func_d)
+        if(func_a.squaredNorm() > 0 and func_d.squaredNorm() > 0):
+            vp.advanced.multiply(prec, ri, 1.0, func_a, func_d)
+        if(func_b.squaredNorm() > 0 and func_c.squaredNorm() > 0):
+            vp.advanced.multiply(prec, ir, -1.0, func_b, func_c)
+
+        output = complex_fcn()
+        output.real = rr + ii
+        output.imag = ri + ir
+
+        return output
+
+def apply_potential(factor, potential, func, prec):
+    output = complex_fcn()
+    vp.advanced.multiply(prec, output.real, factor, potential, func.real)
+    vp.advanced.multiply(prec, output.imag, factor, potential, func.imag)
+    return output
+
+def apply_helmholtz(func, mu, light_speed, prec):
+    out_func = complex_fcn()
+    H = vp.HelmholtzOperator(func.mra, mu, prec)
+    if(func.real.squaredNorm() > 1e-12):
+        vp.advanced.apply(prec, out_func.real, H, func.real)
+    if(func.imag.squaredNorm() > 1e-12):
+        vp.advanced.apply(prec, out_func.imag, H, func.imag)
+    return out_func
+
+def apply_poisson(func, mra, P, prec, thresholdNorm = 0, factor = 1.0):
+    out_func = complex_fcn()
+    if(func.real.squaredNorm() > thresholdNorm):
+        vp.advanced.apply(prec, out_func.real, P, func.real)
+        out_func.real *= factor
+    if(func.imag.squaredNorm() > thresholdNorm):
+        vp.advanced.apply(prec, out_func.imag, P, func.imag)
+        out_func.imag *= factor
+    out_func.cropRealImag(prec)
+    return out_func
+
+def multiply(prec, lhs, rhs):
+    rr = vp.FunctionTree(lhs.mra)
+    ri = vp.FunctionTree(lhs.mra)
+    ir = vp.FunctionTree(lhs.mra)
+    ii = vp.FunctionTree(lhs.mra)
+    vp.advanced.multiply(prec, rr, 1.0, lhs.real, rhs.real, -1, True)
+    vp.advanced.multiply(prec, ri, 1.0, lhs.real, rhs.imag, -1, True)
+    vp.advanced.multiply(prec, ir, 1.0, lhs.imag, rhs.real, -1, True)
+    vp.advanced.multiply(prec, ii, 1.0, lhs.imag, rhs.imag, -1, True)
+    output = complex_fcn()
+    output.real = rr - ii
+    output.imag = ri + ri
+    output.crop(prec)
+    return output
+
+def divergence(vector, prec, der = "BS"):
+    out = complex_fcn()
+    der_vec = {}
+    for i in range(3):
+        der_vec[i] = vector[i].derivative(i, der)
+    out = der_vec[0] + der_vec[1] + der_vec[2]
+    out.cropRealImag(prec)
+    return out
+
+def vector_dot_r(vector, prec):
+    out = complex_fcn()
+    components = {}
+    projection_operator = vp.ScalingProjector(out.mra, prec)
+    for i in range(3):
+        r_i = projection_operator(lambda x : x[i])
+        components[i] = apply_potential(1.0, r_i, vector[i], prec)
+    out = components[0] + components[1] + components[2]
+    out.cropRealImag(prec)
+    return out
+
+def scalar_times_r(function, prec):
+    components = {}
+    projection_operator = vp.ScalingProjector(function.mra, prec)
+    for i in range(3):
+        r_i = projection_operator(lambda x : x[i])
+        components[i] = apply_potential(1.0, r_i, function, prec)
+        components[i].cropRealImag(prec)
+    return components
+
+def vector_tensor_r(vector, prec):
+    tensor = []
+    for i in range(len(vector)):
+        tensor.append(scalar_times_r(vector[i], prec))
+    return tensor
+
+def vector_gradient(vector, der = "BS"):
+    tensor = []
+    for i in range(len(vector)):
+        tensor.append(vector[i].gradient(der))
+    return tensor
+
+def add_vector(func_array, coeff_array, prec):
+    output = complex_fcn()
+    real_array = []
+    imag_array = []
+    for i in range(len(func_array)):
+        real_array.append(( coeff_array[i].real, func_array[i].real))
+        real_array.append((-coeff_array[i].imag, func_array[i].imag))
+        imag_array.append(( coeff_array[i].real, func_array[i].imag))
+        imag_array.append(( coeff_array[i].imag, func_array[i].real))
+    vp.advanced.add(prec, output.real, real_array)
+    vp.advanced.add(prec, output.imag, imag_array)
+    return output
+
diff --git a/orbital4c/nuclear_potential.py b/orbital4c/nuclear_potential.py
new file mode 100644
index 0000000..e29946d
--- /dev/null
+++ b/orbital4c/nuclear_potential.py
@@ -0,0 +1,27 @@
+from scipy.special import erf
+from vampyr import vampyr3d as vp
+from vampyr import vampyr1d as vp1
+import numpy as np
+
+def point_charge(position, center , charge):
+    d2 = ((position[0] - center[0])**2 +
+          (position[1] - center[1])**2 +
+          (position[2] - center[2])**2)
+    distance = np.sqrt(d2)
+    return charge / distance
+
+def coulomb_HFYGB(position, center, charge, precision):
+    d2 = ((position[0] - center[0])**2 +
+          (position[1] - center[1])**2 +
+          (position[2] - center[2])**2)
+    distance = np.sqrt(d2)
+    def smoothing_HFYGB(charge, prec):
+        factor = 0.00435 * prec / charge**5
+        return factor**(1./3.)
+    def uHFYGB(r):
+        u = erf(r)/r + (1/(3*np.sqrt(np.pi)))*(np.exp(-(r**2)) + 16*np.exp(-4*r**2))
+        return u
+    factor = smoothing_HFYGB(charge, precision)
+    value = uHFYGB(distance/factor)
+    return charge * value / factor
+
diff --git a/orbital4c/orbital.py b/orbital4c/orbital.py
new file mode 100644
index 0000000..a0fbf88
--- /dev/null
+++ b/orbital4c/orbital.py
@@ -0,0 +1,393 @@
+from vampyr import vampyr3d as vp
+import numpy as np
+import copy as cp
+from scipy.special import gamma
+from orbital4c import complex_fcn as cf
+
+class orbital4c:
+    """Four components orbital."""
+    mra = None
+    light_speed = -1.0
+    comp_dict = {'La': 0, 'Lb': 1, 'Sa': 2, 'Sb': 3}
+    def __init__(self):
+        self.comp_array = np.array([cf.complex_fcn(),
+                                    cf.complex_fcn(),
+                                    cf.complex_fcn(),
+                                    cf.complex_fcn()])
+        
+    def __getitem__(self, key):
+        return self.comp_array[self.comp_dict[key]]
+    
+    def __setitem__(self, key, val):
+        self.comp_array[self.comp_dict[key]] = val
+        
+    def __len__(self):
+        return 4
+
+    def __str__(self):
+        return ('Large components\n alpha\n{} beta\n{} Small components\n alpha\n{} beta\n{}'.format(self["La"],
+                                  self["Lb"],
+                                  self["Sa"],
+                                  self["Sb"]))
+    
+    def __add__(self, other):
+        output = orbital4c()
+        output.comp_array = self.comp_array + other.comp_array
+        return output
+
+    def __sub__(self, other):
+        output = orbital4c()
+        output.comp_array = self.comp_array - other.comp_array
+        return output
+
+    def __call__(self, position):
+        return [x(position) for x in self.comp_array]
+
+    def save(self, name):
+        self.comp_array[0].save(f"{name}_Large_alpha")
+        self.comp_array[1].save(f"{name}_Large_beta")
+        self.comp_array[2].save(f"{name}_Small_alpha")
+        self.comp_array[3].save(f"{name}_Small_beta")
+
+    def read(self, name):
+        self.comp_array[0].read(f"{name}_Large_alpha")
+        self.comp_array[1].read(f"{name}_Large_beta")
+        self.comp_array[2].read(f"{name}_Small_alpha")
+        self.comp_array[3].read(f"{name}_Small_beta")
+
+    def __rmul__(self, factor):
+        output = orbital4c()
+        output.comp_array =  factor * self.comp_array
+        return output
+
+    def __mul__(self, factor):
+        output = orbital4c()
+        output.comp_array =  factor * self.comp_array 
+        return output   
+
+    def norm(self):
+        out = 0
+        for comp in self.comp_dict.keys():
+            comp_norm = self[comp].squaredNorm()
+            out += comp_norm
+        out = np.sqrt(out)
+        return out
+
+    def squaredNorm(self):
+        out = 0
+        for comp in self.comp_dict.keys():
+            comp_norm = self[comp].squaredNorm()
+            out += comp_norm
+        return out
+
+    def squaredLargeNorm(self):
+        alpha_ns = self.squaredNormComp('La')
+        beta_ns = self.squaredNormComp('Lb')
+        return alpha_ns + beta_ns
+
+    def squaredSmallNorm(self):
+        alpha_ns = self.squaredNormComp('Sa')
+        beta_ns = self.squaredNormComp('Sb')
+        return alpha_ns + beta_ns
+    
+    def squaredNormComp(self, comp):
+        return self[comp].squaredNorm()
+
+    def crop(self, prec):
+        for func in self.comp_array:
+            func.crop(prec)
+
+    def cropLargeSmall(self, prec):
+        largeNorm = np.sqrt(self.squaredLargeNorm())
+        smallNorm = np.sqrt(self.squaredSmallNorm())
+        precLarge = prec * largeNorm
+        precSmall = prec * smallNorm
+#        print('precisions', precLarge, precSmall)
+        self['La'].crop(precLarge, True)
+        self['Lb'].crop(precLarge, True)
+        self['Sa'].crop(precSmall, True)
+        self['Sb'].crop(precSmall, True)        
+        
+    def setZero(self):
+        for func in self.comp_array:
+            func.setZero()
+
+    def rescale(self, factor):
+        for comp in self.comp_array:
+            comp.real *= factor
+            comp.imag *= factor
+            
+    def copy_component(self, func, component='La'):
+        self[component].copy_fcns(func.real, func.imag)
+        
+    def normalize(self):
+        norm_sq = 0
+        for comp in self.comp_array:
+            norm_sq += comp.squaredNorm()
+        norm = np.sqrt(norm_sq)
+        self.rescale(1.0/norm)
+    
+    def copy_components(self, La=None, Lb=None, Sa=None, Sb=None):
+        nr_of_functions = 0
+        if(La != None):
+            nr_of_functions += 1
+            self.copy_component(La, 'La')
+        if(Lb != None):
+            nr_of_functions += 1
+            self.copy_component(Lb, 'Lb')
+        if(Sa != None):
+            nr_of_functions += 1
+            self.copy_component(Sa, 'Sa')
+        if(Sb != None):
+            nr_of_functions += 1
+            self.copy_component(Sb, 'Sb')
+        if(nr_of_functions == 0):
+            print("WARNING: No component copied!")
+        
+    def init_small_components(self,prec):
+        grad_a = self['La'].gradient()
+        grad_b = self['Lb'].gradient()
+        plx = np.array([grad_a[0],grad_b[0]])
+        ply = np.array([grad_a[1],grad_b[1]])
+        plz = np.array([grad_a[2],grad_b[2]])
+        sigma_x = np.array([[0,1],  
+                            [1,0]])
+        sigma_y = np.array([[0,-1j],
+                            [1j,0]])
+        sigma_z = np.array([[1,0],  
+                            [0,-1]])
+
+        sLx = sigma_x@plx
+        sLy = sigma_y@ply
+        sLz = sigma_z@plz
+        
+        sigma_p_L = sLx + sLy + sLz
+
+        sigma_p_L *= -0.5j/orbital4c.light_speed
+        self['Sa'] = sigma_p_L[0]
+        self['Sb'] = sigma_p_L[1]
+        
+    def derivative(self, dir = 0, der = 'ABGV'):
+        orb_der = orbital4c()
+        for key in self.comp_dict:
+            orb_der[key] = self[key].derivative(dir, der) 
+        return orb_der
+    
+    def gradient(self, der = 'ABGV'):
+        orb_grad = {}
+        for key in self.comp_dict.keys():
+            orb_grad[key] = self[key].gradient(der)
+        grad = []
+        for i in range(3):
+            comp = orbital4c()
+            comp.copy_components(La = orb_grad['La'][i], 
+                          Lb = orb_grad['Lb'][i], 
+                          Sa = orb_grad['Sa'][i], 
+                          Sb = orb_grad['Sb'][i])
+            grad.append(comp)
+        return grad
+    
+    def complex_conj(self):
+        orb_out = orbital4c()
+        for key in self.comp_dict.keys():
+            orb_out[key] = self[key].complex_conj() 
+        return orb_out
+
+    def density(self, prec):
+        density = vp.FunctionTree(self.mra)
+        add_vector = []
+        for comp in self.comp_array:
+            temp = comp.density(prec).crop(prec)
+            if(temp.squaredNorm() > 0):
+                add_vector.append((1.0,temp))
+        vp.advanced.add(prec, density, add_vector)
+        return density    
+
+    def overlap_density(self, other, prec):
+        density = cf.complex_fcn()
+        add_vector_real = []
+        add_vector_imag = []
+        for comp in self.comp_dict.keys():
+            func_i = self[comp]
+            func_j = other[comp]
+            temp = cf.multiply(prec, func_i.complex_conj(), func_j)
+            if(temp.real.squaredNorm() > 0):
+                add_vector_real.append((1.0,temp.real))
+            if(temp.imag.squaredNorm() > 0):
+                add_vector_imag.append((1.0,temp.imag))
+        vp.advanced.add(prec, density.real, add_vector_real)
+        vp.advanced.add(prec, density.imag, add_vector_imag)
+        return density
+
+    def alpha(self, direction, prec):
+        out_orb = orbital4c()
+        alpha_order = np.array([[3, 2, 1, 0],
+                                [3, 2, 1, 0],
+                                [2, 3, 0, 1]])
+        
+        alpha_coeff = np.array([[ 1,  1,   1,  1],
+                                [-1j, 1j, -1j, 1j],
+                                [ 1, -1,   1, -1]])
+
+        for idx in range(4):
+            coeff = alpha_coeff[direction][idx]
+            comp = alpha_order[direction][idx]
+            out_orb.comp_array[idx] = coeff * self.comp_array[comp]
+            out_orb.comp_array[idx].crop(prec)
+        return out_orb
+
+    def alpha_p(self, prec, der = "ABGV"):
+        out_orb = orbital4c()
+        orb_grad = self.gradient(der)
+        apx = orb_grad[0].alpha(0, prec)
+        apy = orb_grad[1].alpha(1, prec)
+        apz = orb_grad[2].alpha(2, prec)
+        result = -1j * (apx + apy + apz)
+        result.cropLargeSmall(prec)
+        return result
+
+    def classicT(self, der = 'ABGV'):
+        orb_grad = self.gradient(der)
+        val = 0
+        for i in range(3):
+            val += 0.5 * orb_grad[i].squaredNorm()
+        return val
+    
+    def alpha_vector(self, prec):
+        return [self.alpha(0, prec), self.alpha(1, prec), self.alpha(2, prec)]
+    
+    def ktrs(self, prec):   #Kramers´ Time Reversal Symmetry
+        out_orb = orbital4c()
+        tmp = self.complex_conj()
+        ktrs_order = np.array([1, 0, 3, 2])
+        ktrs_coeff = np.array([-1,  1,  -1,  1])
+        for idx in range(4):
+            coeff = ktrs_coeff[idx]
+            comp = ktrs_order[idx]
+            out_orb.comp_array[idx] = tmp.comp_array[comp].real_mul(coeff)
+        return out_orb
+
+#Beta c**2
+    def beta(self, shift = 0):
+        out_orb = orbital4c()
+#        beta = np.array([[orbital4c.light_speed**2 + shift, 0, 0, 0  ],
+#                         [0, orbital4c.light_speed**2 + shift, 0, 0  ],
+#                         [0, 0, -orbital4c.light_speed**2 + shift, 0 ],
+#                         [0, 0,  0, -orbital4c.light_speed**2 + shift]])
+#        out_orb.comp_array = beta@self.comp_array
+        beta = np.array([orbital4c.light_speed**2 + shift,
+                         orbital4c.light_speed**2 + shift,
+                        -orbital4c.light_speed**2 + shift,
+                        -orbital4c.light_speed**2 + shift])
+        for idx in range(4):
+            out_orb.comp_array[idx] = beta[idx] * self.comp_array[idx]
+        return out_orb
+    
+    def beta2(self):
+        out_orb = orbital4c()
+        beta = np.array([1.0,
+                         1.0,
+                        -1.0,
+                        -1.0])
+        for idx in range(4):
+            out_orb.comp_array[idx] = beta[idx] * self.comp_array[idx]
+        return out_orb
+    
+    def dot(self, other):
+        result = 0
+        for comp in self.comp_dict.keys():
+            factor = 1
+#            if('S' in comp) factor = c**2
+            component = self[comp].dot(other[comp])
+            result += component
+        return result
+
+def apply_dirac_hamiltonian(orbital, prec, shift = 0.0, der = 'ABGV'):
+    beta_phi = orbital.beta(shift)
+    grad_phi = orbital.gradient(der)
+    alpx_phi = -1j * orbital4c.light_speed * grad_phi[0].alpha(0, prec)
+    alpy_phi = -1j * orbital4c.light_speed * grad_phi[1].alpha(1, prec)
+    alpz_phi = -1j * orbital4c.light_speed * grad_phi[2].alpha(2, prec)
+    return beta_phi + alpx_phi + alpy_phi + alpz_phi
+
+def apply_potential(factor, potential, orbital, prec):
+    out_orbital = orbital4c()
+    for comp in orbital.comp_dict:
+        if orbital[comp].squaredNorm() > 0:
+            out_orbital[comp] = cf.apply_potential(factor, potential, orbital[comp], prec)
+    return out_orbital
+
+def apply_complex_potential(factor, potential, orbital, prec):
+    out_orbital = orbital4c()
+    for comp in orbital.comp_dict:
+        if orbital[comp].squaredNorm() > 0:
+            out_orbital[comp] = potential * orbital[comp] 
+    return out_orbital
+
+def add_vector(orbital_array, coeff_array, prec):
+    output = orbital4c()
+    for comp in output.comp_dict:
+        func_array = []
+        for orbital in orbital_array:
+            func_array.append(orbital[comp])
+        output[comp] = cf.add_vector(func_array, coeff_array, prec)
+    return output
+        
+
+def apply_helmholtz(orbital, mu, prec):
+    out_orbital = orbital4c()
+    for comp in orbital.comp_dict.keys():
+        out_orbital[comp] = cf.apply_helmholtz(orbital[comp], mu, orbital4c.light_speed, prec)
+    out_orbital.rescale(-1.0/(2*np.pi))
+    return out_orbital
+
+def init_1s_orbital(orbital,k,Z,n,alpha,origin,prec):
+    gamma_factor = compute_gamma(k,Z,alpha)
+    norm_const = compute_norm_const(n, gamma_factor)
+    idx = 0
+    for comp in orbital.comp_array:
+        func_real = lambda x: one_s_alpha_comp([x[0]-origin[0], x[1]-origin[1], x[2]-origin[2]],
+                                                Z, alpha, gamma_factor, norm_const, idx)
+        func_imag = lambda x: one_s_alpha_comp([x[0]-origin[0], x[1]-origin[1], x[2]-origin[2]],
+                                                Z, alpha, gamma_factor, norm_const, idx+1 )
+        vp.advanced.project(prec, comp.real, func_real)
+        vp.advanced.project(prec, comp.imag, func_imag)
+        idx += 2
+    orbital.normalize()
+    return orbital
+
+def compute_gamma(k,Z,alpha):
+    return np.sqrt(k**2 - Z**2 * alpha**2)
+
+def one_s_alpha(x,Z,alpha,gamma_factor):
+    r = np.sqrt(x[0]**2 + x[1]**2 + x[2]**2)
+    tmp1 = 1.0 + gamma_factor
+    tmp4 = Z * alpha
+    u = x/r
+    lar =   tmp1
+    sai =   tmp4 * u[2]
+    sbr = - tmp4 * u[1]
+    sbi =   tmp4 * u[0]
+    return lar, 0, 0, 0, 0, sai, sbr, sbi
+
+def one_s_alpha_comp(x,Z,alpha,gamma_factor,norm_const,comp):
+    r = np.sqrt(x[0]**2 + x[1]**2 + x[2]**2)
+    tmp2 = r ** (gamma_factor - 1)
+    tmp3 = np.exp(-Z*r)
+    values = one_s_alpha(x,Z,alpha,gamma_factor)
+    return values[comp] * tmp2 * tmp3 * norm_const / np.sqrt(2*np.pi)
+
+def alpha_gradient(orbital, prec):
+    out = orbital4c()
+    grad_vec = orbital.gradient(der = "BS")
+    alpha_vec = {}
+    for i in range(3):
+        alpha_vec[i] = grad_vec[i].alpha(i, prec)
+    out = alpha_vec[0] + alpha_vec[1] + alpha_vec[2]
+    return out
+
+def calc_dirac_mu(energy, light_speed):
+    return np.sqrt((light_speed**4-energy**2)/light_speed**2)
+
+    
+