Repository: libAtoms/matscipy (language: python)

Issue:
scipy.optimize.NoConvergence no longer exists
Corrected the use of the `scipy.optimize.NoConvergence` exception, which is no longer part of `scipy.optimize`.


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/build-wheels.yml b/.github/workflows/build-wheels.yml
index 160bc421..59db9914 100644
--- a/.github/workflows/build-wheels.yml
+++ b/.github/workflows/build-wheels.yml
@@ -25,7 +25,7 @@ jobs:
       matrix:
         buildplat:
         - [ubuntu-20.04, manylinux, x86_64]
-        - [macos-10.15, macosx, x86_64]
+        - [macos-12, macosx, x86_64]
         - [windows-2019, win, AMD64]
 
         python: ["cp38", "cp39", "cp310", "cp311"]
diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
index 616d104f..2d758785 100644
--- a/.github/workflows/publish.yml
+++ b/.github/workflows/publish.yml
@@ -20,7 +20,9 @@ jobs:
       - name: install_python
         run: |
           sudo apt-get update -qy
-          sudo apt-get install -y python3-dev python3-pip libxml2-dev libxslt-dev zlib1g-dev meson
+          sudo apt-get install -y python3-dev python3-pip libxml2-dev libxslt-dev zlib1g-dev
+          # Upgrade to latest meson and ninja
+          sudo pip install --upgrade meson ninja
 
       - name: install_package
         run: |
diff --git a/ChangeLog.md b/ChangeLog.md
index cd414f78..7a0d161f 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -9,6 +9,9 @@ v0.8.0 (not yet released)
 - Writing and reading inhomogenous (and non numeric) data with savetbl/readtbl
 - Hessian-based preconditioner
 - Changed build system to Meson
+- Removed `contact_mechanics` submodule. Stable alternatives:
+  [ContactMechanics](https://github.com/ContactEngineering/ContactMechanics) or
+  [Tamaas](https://gitlab.com/tamaas/tamaas)
 
 v0.7.0 (29Jul21)
 ----------------
diff --git a/docs/api/calculators.rst b/docs/api/calculators.rst
new file mode 100644
index 00000000..c01112c2
--- /dev/null
+++ b/docs/api/calculators.rst
@@ -0,0 +1,4 @@
+Calculators
+-----------
+
+.. automodule:: matscipy.calculators
diff --git a/docs/api/dislocation.rst b/docs/api/dislocation.rst
new file mode 100644
index 00000000..effee202
--- /dev/null
+++ b/docs/api/dislocation.rst
@@ -0,0 +1,4 @@
+Dislocation
+----------
+
+.. automodule:: matscipy.dislocation
diff --git a/docs/api/elasticity.rst b/docs/api/elasticity.rst
new file mode 100644
index 00000000..bbb1f506
--- /dev/null
+++ b/docs/api/elasticity.rst
@@ -0,0 +1,4 @@
+Elasticity
+----------
+
+.. automodule:: matscipy.elasticity
diff --git a/docs/api/fracture_mechanics.rst b/docs/api/fracture_mechanics.rst
new file mode 100644
index 00000000..ce3ec50e
--- /dev/null
+++ b/docs/api/fracture_mechanics.rst
@@ -0,0 +1,4 @@
+Fracture Mechanics
+------------------
+
+.. automodule:: matscipy.fracture_mechanics
diff --git a/docs/api/neighbours.rst b/docs/api/neighbours.rst
new file mode 100644
index 00000000..f53dd84f
--- /dev/null
+++ b/docs/api/neighbours.rst
@@ -0,0 +1,4 @@
+Neighbour lists
+---------------
+
+.. automodule:: matscipy.neighbours
diff --git a/docs/api/numerical.rst b/docs/api/numerical.rst
new file mode 100644
index 00000000..66bdfcf3
--- /dev/null
+++ b/docs/api/numerical.rst
@@ -0,0 +1,4 @@
+Numerical
+---------
+
+.. automodule:: matscipy.numerical
diff --git a/docs/api/root.rst b/docs/api/root.rst
new file mode 100644
index 00000000..df4c2860
--- /dev/null
+++ b/docs/api/root.rst
@@ -0,0 +1,4 @@
+Matscipy root
+-------------
+
+.. automodule:: matscipy
diff --git a/docs/api_reference.rst b/docs/api_reference.rst
new file mode 100644
index 00000000..8cc1358c
--- /dev/null
+++ b/docs/api_reference.rst
@@ -0,0 +1,11 @@
+API Reference
+=============
+
+.. toctree::
+   ./api/root
+   ./api/calculators
+   ./api/elasticity
+   ./api/neighbours
+   ./api/numerical
+   ./api/fracture_mechanics
+   ./api/dislocation
diff --git a/docs/calculators.rst b/docs/calculators.rst
new file mode 100644
index 00000000..14ca09bd
--- /dev/null
+++ b/docs/calculators.rst
@@ -0,0 +1,5 @@
+Calculators
+===========
+
+Here's a description of the calculators available in `matscipy` and how they can
+be composed.
diff --git a/docs/conf.py b/docs/conf.py
index 6385f5f9..4fa60355 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -26,8 +26,9 @@
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
-sys.path.insert(0, os.path.abspath('.'))
-sys.path.insert(0, os.path.abspath('../'))
+
+# sys.path.insert(0, os.path.abspath('.'))
+# sys.path.insert(0, os.path.abspath('../'))
 
 # -- General configuration ------------------------------------------------
 
@@ -269,3 +270,14 @@
 
 # If true, do not generate a @detailmenu in the "Top" node's menu.
 #texinfo_no_detailmenu = False
+
+# -- Extension configuration -------------------------------------------------
+
+autodoc_default_options = {
+    'members': True,
+    'member-order': 'bysource',
+    'undoc-members': True,
+    'show-inheritance': True,
+    'inherited-members': True,
+    'special-members': '__init__',
+}
diff --git a/docs/elastic_constants.rst b/docs/elastic_constants.rst
new file mode 100644
index 00000000..6ff88f6c
--- /dev/null
+++ b/docs/elastic_constants.rst
@@ -0,0 +1,4 @@
+Elastic Constants
+=================
+
+Here's how to use `matscipy` to evaluate elastic constants of an RVE.
diff --git a/docs/fracture_mechanics.rst b/docs/fracture_mechanics.rst
new file mode 100644
index 00000000..8eeeca0f
--- /dev/null
+++ b/docs/fracture_mechanics.rst
@@ -0,0 +1,4 @@
+Fracture Mechanics
+==================
+
+Here's what you can do with Matscipy in the field of fracture mechanics.
diff --git a/docs/index.rst b/docs/index.rst
index 773d94f7..0eb8607e 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -6,37 +6,49 @@
 .. include:: ../README.rst
 
 
-Contents
-========
-
-.. autosummary::
-
-   matscipy.angle_distribution
-   matscipy.atomic_strain
-   matscipy.calculators
-   matscipy.calculators.eam
-   matscipy.calculators.pair_potential
-   matscipy.calculators.mcfm
-   matscipy.calculators.supercell_calculator
-   matscipy.calculators.manybody
-   matscipy.contact_mechanics
-   matscipy.contact_mechanics.DMT
-   matscipy.contact_mechanics.greens_function
-   matscipy.contact_mechanics.Hertz
-   matscipy.contact_mechanics.JKR
-   matscipy.elasticity
-   matscipy.fracture_mechanics
-   matscipy.fracture_mechanics.clusters
-   matscipy.fracture_mechanics.crack
-   matscipy.fracture_mechanics.crackpathsel
-   matscipy.fracture_mechanics.energy_release
-   matscipy.fracture_mechanics.idealbrittlesolid
-   matscipy.hydrogenate
-   matscipy.io
-   matscipy.logger
-   matscipy.neighbours
-   matscipy.rings
-   matscipy.surface
+..
+   Matscipy submodules
+   ===================
+
+   .. autosummary::
+
+      matscipy.angle_distribution
+      matscipy.atomic_strain
+      matscipy.calculators
+      matscipy.calculators.eam
+      matscipy.calculators.pair_potential
+      matscipy.calculators.mcfm
+      matscipy.calculators.supercell_calculator
+      matscipy.calculators.manybody
+      matscipy.contact_mechanics
+      matscipy.contact_mechanics.DMT
+      matscipy.contact_mechanics.greens_function
+      matscipy.contact_mechanics.Hertz
+      matscipy.contact_mechanics.JKR
+      matscipy.elasticity
+      matscipy.fracture_mechanics
+      matscipy.fracture_mechanics.clusters
+      matscipy.fracture_mechanics.crack
+      matscipy.fracture_mechanics.crackpathsel
+      matscipy.fracture_mechanics.energy_release
+      matscipy.fracture_mechanics.idealbrittlesolid
+      matscipy.hydrogenate
+      matscipy.io
+      matscipy.logger
+      matscipy.neighbours
+      matscipy.numerical
+      matscipy.rings
+      matscipy.surface
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Table of Contents:
+
+   ./calculators
+   ./fracture_mechanics
+   ./plasticity
+   ./elastic_constants
+   ./api_reference
 
 
 Indices and tables
diff --git a/docs/plasticity.rst b/docs/plasticity.rst
new file mode 100644
index 00000000..b710c2c0
--- /dev/null
+++ b/docs/plasticity.rst
@@ -0,0 +1,4 @@
+Plasticity
+==========
+
+Here's what you can do with Matscipy in the field of plasticity.
\ No newline at end of file
diff --git a/matscipy/calculators/committee/__init__.py b/matscipy/calculators/committee/__init__.py
new file mode 100644
index 00000000..cb39b015
--- /dev/null
+++ b/matscipy/calculators/committee/__init__.py
@@ -0,0 +1,5 @@
+"""Implements the Committee (of Models) approach."""
+
+from .committee import CommitteeUncertainty, Committee, CommitteeMember
+from .log import set_logging
+from .utils import subsample
diff --git a/matscipy/calculators/committee/committee.py b/matscipy/calculators/committee/committee.py
new file mode 100644
index 00000000..a9622704
--- /dev/null
+++ b/matscipy/calculators/committee/committee.py
@@ -0,0 +1,460 @@
+import warnings
+import logging
+from pathlib import Path
+from collections import Counter
+
+import numpy as np
+import ase.io
+from ase.calculators.calculator import Calculator, all_changes
+
+
+logger = logging.getLogger('matscipy.calculators.committee')
+
+
+class CommitteeUncertainty(Calculator):
+    """
+    Calculator for a committee of machine learned interatomic potentials (MLIP).
+
+    The class assumes individual members of the committee already exist (i.e. their
+    training is performed externally). Instances of this class are initialized with
+    these committee members and results (energy, forces) are calculated as average
+    over these members. In addition to these values, also the uncertainty (standard
+    deviation) is calculated.
+
+    The idea for this Calculator class is based on the following publication:
+    Musil et al., J. Chem. Theory Comput. 15, 906−915 (2019)
+    https://pubs.acs.org/doi/full/10.1021/acs.jctc.8b00959
+
+    Parameter:
+    ----------
+    committee: Committee-instance
+        Representation for a collection of Calculators.
+    atoms : ase-Atoms, optional default=None
+        Optional object to which the calculator will be attached.
+    """
+
+    def __init__(self, committee, atoms=None):
+
+        self.implemented_properties = ['energy', 'forces', 'stress']
+
+        self.committee = committee
+
+        super().__init__(atoms=atoms)
+
+        logger.info('Initialized committee uncertainty calculator')
+        for line_i in self.committee.__repr__().splitlines():
+            logger.debug(line_i)
+
+    def calculate(self, atoms=None, properties=['energy', 'forces', 'stress'], system_changes=all_changes):
+        """Calculates committee (mean) values and variances."""
+
+        logger.info(f'Calculating properties {properties} with committee')
+        super().calculate(atoms, properties, system_changes)
+
+        property_committee = {k_i: [] for k_i in properties}
+
+        for cm_i in self.committee.members:
+            cm_i.calculator.calculate(atoms=atoms, properties=properties, system_changes=system_changes)
+
+            for p_i in properties:
+                property_committee[p_i].append(cm_i.calculator.results[p_i])
+
+        for p_i in properties:
+            self.results[p_i] = np.mean(property_committee[p_i], axis=0)
+            self.results[f'{p_i}_uncertainty'] = np.sqrt(np.var(property_committee[p_i], ddof=1, axis=0))
+
+            if self.committee.is_calibrated_for(p_i):
+                self.results[f'{p_i}_uncertainty'] = self.committee.scale_uncertainty(self.results[f'{p_i}_uncertainty'], p_i)
+            else:
+                warnings.warn(f'Uncertainty estimation has not been calibrated for {p_i}.')
+
+
+class Committee:
+    """
+    Instances of this class represent a committee of models.
+
+    It's use is to store the ```CommitteeMember```s representing the committee model
+    and to calibrate the obtained uncertainties (required when sub-sampling is used
+    to create the training data of the committee members).
+
+    Parameter:
+    ----------
+    members: list(M)
+        List of ```CommitteeMember``` instances representing the committee of (here `M`) models.
+    """
+    def __init__(self, members=None):
+        self.members = [] if members is None else members
+
+        logger.info('Initialized committee')
+        for line_i in self.__repr__().splitlines():
+            logger.debug(line_i)
+
+    @property
+    def members(self):
+        """List with committee members."""
+        return self._members
+
+    @members.setter
+    def members(self, members):
+        """Set list with committee members."""
+        for member_i in members:
+            self._check_member_type(member_i)
+        self._members = members
+        logger.info(f'Set {len(self.members)} members to represent the committee')
+        self._update()
+
+    @property
+    def number(self):
+        """Number of committee members."""
+        return self._number
+
+    @property
+    def atoms(self):
+        """Combined Atoms/samples in the committee."""
+        return self._atoms
+
+    @property
+    def ids(self):
+        """Identifiers of atoms/samples in the committee."""
+        return self._ids
+
+    @property
+    def id_to_atoms(self):
+        """Dictionary to translate identifiers to Atoms-objects."""
+        return self._id_to_atoms
+
+    @property
+    def id_counter(self):
+        """Counter-object for identifier appearances in the committee."""
+        return self._id_counter
+
+    @property
+    def alphas(self):
+        """(Linear) scaling factors for committee uncertainties."""
+        return self._alphas
+
+    @property
+    def calibrated_for(self):
+        """Set of properties the committee has been calibrated for."""
+        return self.alphas.keys()
+
+    @property
+    def validation_set(self):
+        """List of Atoms-objects."""
+        if not self._validation_set:
+            msg = '`Committee.set_internal_validation_set()` has not been called or ' + \
+                  '`Committee`-instance has been altered since last call.'
+            logger.warning(msg)
+            warnings.warn(msg)
+        return self._validation_set
+
+    def _update(self):
+        """Update status when ```Committee```-instance has been altered."""
+        self._number = len(self.members)
+        self._atoms = [atoms_ij for cm_i in self.members for atoms_ij in cm_i.atoms]
+        self._ids = [id_ij for cm_i in self.members for id_ij in cm_i.ids]
+        self._id_to_atoms = {id_i: atoms_i for id_i, atoms_i in zip(self.ids, self.atoms)}
+        self._id_counter = Counter(self.ids)
+        self._validation_set = []
+        self._alphas = {}
+        logger.info('Updated committee status')
+
+    def add_member(self, member):
+        """Extend committee by new ```member``` (i.e. ```CommitteeMember```-instance)."""
+        self._check_member_type(member)
+        self.members.append(member)
+        logger.info('Added +1 member to the committee')
+        self._update()
+
+    def __add__(self, member):
+        """Extend committee by new ```member``` (i.e. ```CommitteeMember```-instance)."""
+        self._check_member_type(member)
+        self.add_member(member)
+        return self
+
+    @staticmethod
+    def _check_member_type(member):
+        """Make sure ```member``` is of type ```CommitteeMember```."""
+        assert isinstance(member, CommitteeMember), \
+            f'Members of `Committee` need to be of type `CommitteeMember`. Found {type(member)}'
+
+    def set_internal_validation_set(self, appearance_threshold):
+        """
+        Define a validation set based on the Atoms-objects of sub-sampled committee training sets.
+
+        Parameter:
+        ----------
+        appearance_threshold: int
+            Number of times a sample for the validation set
+            is maximally allowed to appear across the training sets
+            of committee members.
+        """
+        if self._alphas:
+            msg = '`alphas` will be reset to avoid inconsistencies with new validation set.'
+            logger.warning(msg)
+            warnings.warn(msg)
+            self._reset_calibration_parameters()
+
+        assert 0 < appearance_threshold <= self.number - 2
+
+        self._validation_set = []
+        for id_i, appearance_i in self.id_counter.most_common()[::-1]:
+            if appearance_i > appearance_threshold:
+                break
+            self._validation_set.append(self.id_to_atoms[id_i])
+
+        logger.info(f'Set internal validation set with {len(self.validation_set)} entries')
+
+    def _reset_calibration_parameters(self):
+        """Reset parameters obtained from calling ```self.calibrate()```."""
+        self._alphas = {}
+        logger.info('Reset calibration parameters')
+
+    def calibrate(self, prop, key, location, system_changes=all_changes):
+        """
+        Obtain parameters to properly scale committee uncertainties and make
+        them available as an attribute (```alphas```) with another associated
+        attribute (```calibrated_for```) providing information about the property
+        for which the uncertainty will be scaled by it.
+
+        Parameter:
+        ----------
+        properties: list(str)
+            Properties for which the calibration will determine scaling factors.
+        key: str
+            Key under which the reference values in the validation set are stored
+            (i.e. under Atoms.info[```key```] / Atoms.arrays[```key```]).
+        location: str
+            Either 'info' or 'arrays'.
+        """
+        assert location in ['info', 'arrays'], f'`location` must be \'info\' or \'arrays\', not \'{location}\'.'
+
+        validation_ref = [np.asarray(getattr(sample_i, location)[key]).flatten() for sample_i in self.validation_set]
+        validation_pred, validation_pred_var = [], []
+
+        for idx_i, sample_i in enumerate(self.validation_set):
+
+            sample_committee_pred = []
+
+            for cm_j in self.members:
+
+                if cm_j.is_sample_in_atoms(sample_i):
+                    continue
+
+                cm_j.calculator.calculate(atoms=sample_i, properties=[prop], system_changes=system_changes)
+                sample_committee_pred.append(cm_j.calculator.results[prop])
+
+            validation_pred.append(np.mean(sample_committee_pred, axis=0).flatten())
+            validation_pred_var.append(np.var(sample_committee_pred, ddof=1, axis=0).flatten())
+
+        # For symmetry-reasons it can happen that e.g. all values for a force component of an atom are equal.
+        # This would lead to a division-by-zero error in self._calculate_alpha() due to zero-variances.
+        validation_ref = np.concatenate(validation_ref)
+        validation_pred = np.concatenate(validation_pred)
+        validation_pred_var = np.concatenate(validation_pred_var)
+        ignore_indices = np.where(validation_pred_var == 0)[0]
+        validation_ref = np.delete(validation_ref, ignore_indices)
+        validation_pred = np.delete(validation_pred, ignore_indices)
+        validation_pred_var = np.delete(validation_pred_var, ignore_indices)
+
+        self._alphas.update(
+                {prop: self._calculate_alpha(
+                    vals_ref=validation_ref,
+                    vals_pred=validation_pred,
+                    vars_pred=validation_pred_var,
+                    )
+                 })
+
+        logger.info(f'Calibrated committee for property \'{prop}\'')
+        logger.debug(f'\talpha = {self.alphas[prop]}')
+
+    def is_calibrated_for(self, prop):
+        """Check whether committee has been calibrated for ```prop```."""
+        return prop in self.calibrated_for
+
+    def _calculate_alpha(self, vals_ref, vals_pred, vars_pred):
+        """
+        Get (linear) uncertainty scaling factor alpha.
+
+        This implementation is based on:
+        Imbalzano et al., J. Chem. Phys. 154, 074102 (2021)
+        https://doi.org/10.1063/5.0036522
+
+        Parameter:
+        ----------
+        vals_ref: ndarray(N)
+            Reference values for validation set samples.
+        vals_pred: ndarray(N)
+            Values predicted by the committee for validation set samples.
+        vars_pred: ndarray(N)
+            Variance predicted by the committee for validation set samples.
+
+        Returns:
+        --------
+        (Linear) uncertainty scaling factor alpha.
+        """
+        N_val = len(vals_ref)
+        M = self.number
+        alpha_squared = -1/M + (M - 3)/(M - 1) * 1/N_val * np.sum(np.power(vals_ref-vals_pred, 2) / vars_pred)
+        logger.info(f'Calculated alpha')
+        logger.debug(f'\tN_val          = {N_val}')
+        logger.debug(f'\tM              = {M}')
+        logger.debug(f'\talpha_squared  = {alpha_squared}')
+        assert alpha_squared > 0, f'Obtained negative value for `alpha_squared`: {alpha_squared}'
+        return np.sqrt(alpha_squared)
+
+    def scale_uncertainty(self, value, prop):
+        """
+        Scale uncertainty ```value``` obtained with the committee according to the calibration
+        for the corresponding property (```prop```).
+
+        Parameter:
+        ----------
+        value: float / ndarray
+            Represents the uncertainty values (e.g. energy, forces) to be scaled.
+        prop: str
+            The property associated with ```value``` (for which the committee needs to be calibrated).
+
+        Returns:
+        --------
+        Scaled input ```value```.
+        """
+        return self.alphas[prop] * value
+
+    def __repr__(self):
+        s = ''
+
+        s_i = f'Committee Status\n'
+        s += s_i
+        s += '='*len(s_i) + '\n\n'
+
+        s += f'# members:                    {self.number:>10d}\n'
+        s += f'# atoms:                      {len(self.atoms):>10d}\n'
+        s += f'# ids:                        {len(self.ids):>10d}\n'
+        s += f'# atoms validation set:       {len(self._validation_set):>10d}\n'
+        if not self.calibrated_for:
+            s += f'calibrated for:               {"-":>10}\n'
+        else:
+            s += f'calibrated for:\n'
+            for p_i in sorted(self.calibrated_for):
+                s += f'{"":>4s}{p_i:<18}{self.alphas[p_i]:>18}\n'
+
+        for idx_i, cm_i in enumerate(self.members):
+            s += '\n\n'
+            s_i = f'Committee Member {idx_i}:\n'
+            s += s_i
+            s += '-'*len(s_i) + '\n'
+            s += cm_i.__repr__()
+
+        return s
+
+
+class CommitteeMember:
+    """
+    Lightweight class defining a member (i.e. a sub-model) of a committee model.
+
+    Parameter:
+    ----------
+    calculator: Calculator
+        Instance of a Calculator-class (or heirs e.g. quippy.potential.Potential)
+        representing a machine-learned model.
+    training_data: str / Path / list(Atoms), optional default=None
+        Path to or Atoms of (sub-sampled) training set used to create the machine-learned model
+        defined by the ```calculator```.
+    """
+    def __init__(self, calculator, training_data=None):
+        self._calculator = calculator
+
+        self._filename = 'no filename'
+        self._atoms = []
+        self._ids = []
+
+        if training_data is not None:
+            self.set_training_data(training_data)
+
+        logger.info('Created committee member')
+        for line_i in self.__repr__().splitlines():
+            logger.debug(line_i)
+
+    @property
+    def calculator(self):
+        """Model of the committee member."""
+        return self._calculator
+
+    @property
+    def filename(self):
+        """Path to the atoms/samples in the committee member."""
+        return self._filename
+
+    @filename.setter
+    def filename(self, filename):
+        """Set path to the atoms/samples in the committee member."""
+        msg = 'Use `set_training_data()` to modify the committee member'
+        logger.error(msg)
+        raise RuntimeError(msg)
+
+    @property
+    def atoms(self):
+        """Atoms/samples in the committee member."""
+        return self._atoms
+
+    @atoms.setter
+    def atoms(self, atoms):
+        """Set Atoms/samples in the committee member."""
+        msg = 'Use `set_training_data()` to modify the committee member'
+        logger.error(msg)
+        raise RuntimeError(msg)
+
+    @property
+    def ids(self):
+        """Identifiers of atoms/samples in the committee member."""
+        return self._ids
+
+    @ids.setter
+    def ids(self, ids):
+        """Set identifiers of atoms/samples in the committee member."""
+        msg = 'Use `set_training_data()` to modify the committee member'
+        logger.error(msg)
+        raise RuntimeError(msg)
+
+    def set_training_data(self, training_data):
+        """
+        Read in and store the training data of this committee members from the passed ```filename```.
+
+        Parameter:
+        ----------
+        training_data: str / Path / list(Atoms), optional default=None
+            Path to or Atoms of (sub-sampled) training set used to create the machine-learned model
+            defined by the ```calculator```. Individual Atoms need an Atoms.info['_Index_FullTrainingset']
+            for unique identification.
+        """
+        if len(self.atoms) > 0:
+            msg = 'Overwriting current training data.'
+            logger.warning(msg)
+            warnings.warn(msg)
+
+        if isinstance(training_data, (str, Path)):
+            self._filename = Path(training_data)
+            self._atoms = ase.io.read(self.filename, ':')
+        elif isinstance(training_data, list):
+            self._filename = 'No Filename'
+            self._atoms = training_data
+        self._ids = [atoms_i.info['_Index_FullTrainingSet'] for atoms_i in self.atoms]
+
+    def is_sample_in_atoms(self, sample):
+        """Check if passed Atoms-object is part of this committee member (by comparing identifiers)."""
+        if '_Index_FullTrainingSet' not in sample.info:
+            msg = 'Can\'t test if `sample` is in `atoms`. `sample` has no Atoms.info[\'_Index_FullTrainingSet\']'
+            logger.error(msg)
+            raise RuntimeError(msg)
+        else:
+            return sample.info['_Index_FullTrainingSet'] in self.ids
+
+    def __repr__(self):
+        s = ''
+        s += f'calculator: {str(self.calculator.__class__):>60s}\n'
+        s += f'filename:   {str(self.filename):>60s}\n'
+        s += f'# Atoms:    {len(self.atoms):>60d}\n'
+        s += f'# IDs:      {len(self.ids):>60d}'
+        return s
+
diff --git a/matscipy/calculators/committee/log.py b/matscipy/calculators/committee/log.py
new file mode 100644
index 00000000..a7caefb6
--- /dev/null
+++ b/matscipy/calculators/committee/log.py
@@ -0,0 +1,65 @@
+import os
+import logging
+
+
+def create_logger(name, level='warning', log_file=None):
+    """
+    Create a logger with defined level, and some nice formatter and handler.
+
+    Parameter:
+    ----------
+    name: str
+        Name of the logger.
+
+    level: str, optional default='warning'
+        The sensitivity level of both the logger and the handler.
+        Allowed options are 'debug', 'info', 'warning', 'error' and 'critical'.
+
+    log_file: str, optional default=None
+        Path to the file, in case you want to log to a file rather than standard out.
+
+    Returns:
+    --------
+    logger: logging.Logger
+        The respective logger.
+    """
+
+    levels = {'debug':      logging.DEBUG,
+              'info':       logging.INFO,
+              'warning':    logging.WARNING,
+              'error':      logging.ERROR,
+              'critical':   logging.CRITICAL}
+
+    logger = logging.getLogger(name)
+    logger.setLevel(levels[level])
+
+    # create console handler and set level to debug
+    if not log_file:
+        ch = logging.StreamHandler()
+    else:
+        ch = logging.FileHandler(os.path.abspath(log_file))
+
+    ch.setLevel(levels[level])
+
+    # create formatter
+    formatter = logging.Formatter('%(asctime)s: [%(levelname)s] %(name)s : %(message)s')
+    ch.setFormatter(formatter)
+    logger.addHandler(ch)
+
+    return logger
+
+
+def set_logging(level='warning', log_file=None):
+    """
+    Set the global logging level (and corresponding file location).
+
+    Parameter:
+    ----------
+    level: str, optional default='warning'
+        The sensitivity level of both the logger and the handler.
+        Allowed options are 'debug', 'info', 'warning', 'error' and 'critical'.
+
+    log_file: str, optional default=None
+        Path to the file, in case you want to log to a file rather than standard out.
+    """
+    logger = create_logger('matscipy.calculators.committee', level=level, log_file=log_file)
diff --git a/matscipy/contact_mechanics/meson.build b/matscipy/calculators/committee/meson.build
similarity index 57%
rename from matscipy/contact_mechanics/meson.build
rename to matscipy/calculators/committee/meson.build
index c45a055b..151ff639 100644
--- a/matscipy/contact_mechanics/meson.build
+++ b/matscipy/calculators/committee/meson.build
@@ -3,15 +3,13 @@
 # Pure Python sources
 python_sources = [
     '__init__.py',
-    'analysis.py',
-    'DMT.py',
-    'greens_function.py',
-    'Hertz.py',
-    'JKR.py'
+    'committee.py',
+    'log.py',
+    'utils.py'
 ]
 
 # Install pure Python
 python.install_sources(
     python_sources,
-    subdir: 'matscipy/contact_mechanics'
+    subdir: 'matscipy/calculators/committee'
 )
diff --git a/matscipy/calculators/committee/utils.py b/matscipy/calculators/committee/utils.py
new file mode 100644
index 00000000..06fc4e75
--- /dev/null
+++ b/matscipy/calculators/committee/utils.py
@@ -0,0 +1,86 @@
+import numpy as np
+import ase.io
+
+
+def subsample(samples, output_files, num_subsets, num_subset_samples, keep_isolated_atoms=True):
+    """
+    Draw sub-samples (without replacement) from given configurations and write to files.
+
+    Parameter:
+    ----------
+    samples: list(Atoms)
+        List of configurations to draw sub-samples from (e.g. full training set).
+    output_files: list(str / Path)
+        Target locations for sub-sampled sets of configurations.
+    num_subsets: int
+        Number of sub-sets to be drawn.
+    num_subset_samples: int
+        Number of configurations per sub-sets.
+    keep_isolated_atoms: bool, default True
+        Make isolated atoms (if present) be part of each sub-set.
+
+    Returns:
+    --------
+    subsamples: list(list(Atoms))
+        Contains the lists of sub-sampled configurations.
+    """
+    # keep track of position in original set of configurations
+    sample_atoms = []
+    isolated_atoms = []  # keep isolated atoms for each sub-sample
+    for idx_i, atoms_i in enumerate(samples):
+        atoms_i.info['_Index_FullTrainingSet'] = idx_i
+        if keep_isolated_atoms and len(atoms_i) == 1:
+            isolated_atoms.append(atoms_i)
+        else:
+            sample_atoms.append(atoms_i)
+
+    num_subset_samples -= len(isolated_atoms)
+    assert 1 < num_subset_samples <= len(sample_atoms), 'Negative `num_subset_samples` (after reduction by number of isolated atoms)'
+    assert len(output_files) == num_subsets, f'`outputs` requires `num_subsets` files to be specified.'
+
+    subsample_indices = _get_subsample_indices(len(sample_atoms), num_subsets, num_subset_samples)
+    subsamples = [isolated_atoms + [sample_atoms[idx_i] for idx_i in idxs] for idxs in subsample_indices]
+
+    for output_file_i, subsample_i in zip(output_files, subsamples):
+        ase.io.write(output_file_i, subsample_i)
+
+    return subsamples
+
+
+def _get_subsample_indices(num_samples, num_subsets, num_subset_samples):
+    """
+    Draw indices for sub-samples (without replacement).
+
+    Parameter:
+    ----------
+    num_samples: int
+        Number of configurations to draw sub-samples from (e.g. size of full training set).
+    num_subsets: int
+        Number of sub-sets to be drawn.
+    num_subset_samples: int
+        Number of configurations per sub-sets.
+
+    Returns:
+    --------
+    subsample_indices: list(list(int))
+        Contains the lists of indices representing sub-sampled configurations.
+    """
+    indice_pool = np.arange(num_samples)
+    subsample_indices = []
+    for _ in range(num_subsets):
+        if num_subset_samples <= len(indice_pool):
+            selected_indices = np.random.choice(indice_pool, num_subset_samples, False)
+            indice_pool = indice_pool[~np.isin(indice_pool, selected_indices)]
+            subsample_indices.append(selected_indices)
+        else:
+            selected_indices_part_1 = indice_pool
+            # re-fill pool with indices, taking account of already selected ones,
+            # in order to avoid duplicate selections
+            indice_pool = np.arange(num_samples)
+            indice_pool = indice_pool[~np.isin(indice_pool, selected_indices_part_1)]
+            selected_indices_part_2 = np.random.choice(indice_pool, num_subset_samples - len(selected_indices_part_1), False)
+            indice_pool = indice_pool[~np.isin(indice_pool, selected_indices_part_2)]
+            selected_indices = np.concatenate((selected_indices_part_1, selected_indices_part_2))
+            subsample_indices.append(selected_indices)
+
+    return subsample_indices
diff --git a/matscipy/calculators/meson.build b/matscipy/calculators/meson.build
index c5eeed89..3479a15f 100644
--- a/matscipy/calculators/meson.build
+++ b/matscipy/calculators/meson.build
@@ -17,6 +17,7 @@ python.install_sources(
 
 subdir('eam')
 subdir('ewald')
+subdir('committee')
 subdir('manybody')
 subdir('mcfm')
 subdir('pair_potential')
diff --git a/matscipy/cauchy_born.py b/matscipy/cauchy_born.py
new file mode 100644
index 00000000..760fc4da
--- /dev/null
+++ b/matscipy/cauchy_born.py
@@ -0,0 +1,742 @@
+import numpy as np
+from ase.lattice.cubic import Diamond
+
+from ase.optimize import LBFGS
+
+from scipy.linalg import polar, sqrtm
+
+from itertools import permutations
+
+from sklearn.linear_model import Lasso, LinearRegression
+
+from scipy.stats.qmc import LatinHypercube, scale
+
+class CubicCauchyBorn:
+    '''Corrector for the cauchy-born prediction of atomistic positions in multi-lattices subject to continuum strain
+        This model exploits the symmetry triad down the 111 axis that all cubic crystals posess by making models that 
+        only fit a single component of the shift-correction vector. The other 2 components can be obtained by calling the same model
+        with strain states that are rotated about the 111 axis.'''
+    def __init__(self,el,a0,calc,lattice=Diamond):
+        """Parameters
+        ----------
+        el : string
+            ASE chemical element that constitutes lattice
+        a0 : float
+            Lattice constant for cubic lattice
+        calc : ASE calculator object
+            ASE calculator.
+        lattice : ASE lattice builder function
+            ASE function from ase.lattice.cubic
+        """
+        self.el = el
+        self.a0 = a0
+        self.calc = calc
+        self.lattice = lattice
+        self.Lasso_Fitted = False #status of model fit
+        self.lattice1mask = None #mask for the lattice 1 atoms
+        self.lattice2mask = None #mask for the lattice 2 atoms
+
+
+    def set_sublattices(self,atoms,A):
+        '''Apply a small strain to all atoms in the supplied atoms structure and determine which atoms belong to which sublattice using forces.
+        NOTE as this method is based on forces, it does not work in cases where atom forces are high already, 
+        This means it should be applied to an unstrained version of the atoms structure and any sublattices
+        assigned adjacent to free surfaces cannot be trusted. This function updates self.lattice1mask and self.lattice2mask with the
+        measured sublattices.
+        Parameters
+        ----------
+        atoms : ASE atoms object
+            atoms structure with atoms on seperate sublattices that need assigning
+        A : 3x3 numpy array, floats
+            rotation matrix of the form [x^T,y^T,z^T], where x^T,y^T,z^T are normalised column vectors of the lab frame directions expressed
+            in terms of the crystal lattice directions.
+        '''
+
+        #currently, the idea for this is to apply a small strain to atoms
+        #measure the force difference, and set all the forces which are positive to sublattice
+        #1, and all those which are negative to sublattice 2. 
+
+        #generate U (right hand stretch tensor) to apply
+        U_voigt = np.array([1.001, 1.003, 1.002, 0.006, 0.002, 0.004])
+        U = np.zeros([3,3])
+        U[0,0] = U_voigt[0]
+        U[1,1] = U_voigt[1]
+        U[2,2] = U_voigt[2]
+        U[1,2],U[2,1] = U_voigt[3],U_voigt[3]
+        U[0,2],U[2,0] = U_voigt[4],U_voigt[4]
+        U[0,1],U[1,0] = U_voigt[5],U_voigt[5]
+        atoms_copy = atoms.copy()
+
+        #generate a copy of the atoms and stretch the cell, scaling atoms.
+        #get forces before and after and compare these
+        atoms_copy.set_calculator(self.calc)
+        f_before = atoms_copy.get_forces()
+        cell = atoms_copy.get_cell()
+        scaled_cell = U@cell
+        atoms_copy.set_cell(scaled_cell,scale_atoms=True)
+        f_after = atoms_copy.get_forces()
+
+        force_diff = f_after-f_before
+
+        #transform force_diff to lattice frame
+        force_diff_lattice = np.zeros_like(force_diff)
+        for i in range(len(atoms)):
+            force_diff_lattice[i,:] = A@force_diff[i,:]
+
+        #generate the lattice masks, checking to ensure that if by chance no shift occurs in one direction due to the applied strains 
+        #that the other directions are used instead and things don't break
+        lattice1mask = force_diff_lattice[:,0]>0
+        lattice2mask = force_diff_lattice[:,0]<0
+        if all(element == lattice1mask[0] for element in lattice1mask): #if they're all false or true, try the y component
+            lattice1mask = force_diff_lattice[:,1]>0
+            lattice2mask = force_diff_lattice[:,1]<0
+            if all(element == lattice1mask[0] for element in lattice1mask): #if they're all false, try the z component
+                lattice1mask = force_diff_lattice[:,2]>0
+                lattice2mask = force_diff_lattice[:,2]<0
+        
+
+        #set the masks for which atoms are in each of the sublattices
+        self.lattice1mask = lattice1mask
+        self.lattice2mask = lattice2mask     
+
+        
+
+    
+
+    def fit_taylor(self,de=1e-4):
+        '''fit a simple taylor-expansion type model to predict the cauchy-born shift corrector at a given applied strain
+            using finite differences. Sets self.grad_f to the vector of first derivatives of the first shift component 
+            with each of the 6 strain components and self.hess_f to the hessian of the first shift component with the strain components.
+            NOTE this model is generally worse than the model provided later in initial_lasso_fit.
+        Parameters
+        ----------
+        de : float
+            tolerance for finite difference gradient calculations
+        '''
+        calc = self.calc
+        def f_gl(E_vec,calc,unitcell):
+            '''Function that calculates the first component of the cauchy-born shift corrector for a given green-lagrange strain tensor.
+            Parameters
+            ----------
+            E_vec : array_like
+                green-lagrange strain state to apply in voigt notation
+            calc : ASE calculator object
+                ASE calculator to use to find shift
+            unitcell : ASE atoms object
+                Single unit cell of crystal to apply strain vector to and to find shift
+
+            Returns
+            -------
+            shift_diff[0] : float
+                first component of cauchy-born shift correction for applied strain state
+            '''
+            #green lagrange version of function
+            #turn E into matrix
+            E = np.zeros([3,3])
+            E[0,0] = E_vec[0]
+            E[1,1] = E_vec[1]
+            E[2,2] = E_vec[2]
+            E[2,1],E[1,2] = E_vec[3],E_vec[3]
+            E[2,0],E[0,2] = E_vec[4],E_vec[4]
+            E[0,1],E[1,0] = E_vec[5],E_vec[5]
+            # get U^2
+            Usqr = 2*E + np.eye(3)
+            #square root matrix to get U
+            U = sqrtm(Usqr,disp=True)
+
+            #this is just the symmetric stretch tensor, exactly what we need.
+            x = U
+
+
+            initial_shift = np.zeros([3])
+            relaxed_shift = np.zeros([3])
+
+            #build cell
+            unitcell_copy = unitcell.copy()
+            unitcell_copy.calc = self.calc
+            cell = unitcell_copy.get_cell()
+            cell_rescale = x[:,:]@cell
+            unitcell_copy.set_cell(cell_rescale,scale_atoms=True)
+            initial_shift[:] = unitcell_copy.get_positions()[1] - unitcell_copy.get_positions()[0]
+
+            #relax cell
+            opt = LBFGS(unitcell_copy,logfile=None)
+            opt.run(fmax = 1e-10)
+            relaxed_shift[:] = unitcell_copy.get_positions()[1] - unitcell_copy.get_positions()[0]
+            
+            #get shift
+            shift_diff = relaxed_shift-initial_shift
+            #print(shift_diff)
+            # THIS CAN BE SPED UP
+            # TECHNICALLY, EACH CALL OF THIS FUNCTION CAN BE USED
+            # TO COMPUTE 3 DERIVATIVES RATHER THAN JUST 1
+            # IF YOU EXPLOIT THE SYMMETRY OF CUBIC SYSTEMS
+            # WOULD MAKE THINGS 3x FASTER, BUT WOULD REQUIRE 
+            # *SOME* CODING
+            return shift_diff[0] #return 1 component
+
+        
+        def find_grad(eps,i,*args):
+            #compute gradient of function by varying a single component of eps
+            eps_up = eps.copy()
+            eps_down = eps.copy()
+            eps_up[i]+=self.de
+            eps_down[i]-=self.de
+            return (f_gl(eps_up,*args) - f_gl(eps_down,*args))/(2*self.de)
+
+        def find_second_deriv(eps,i,j,*args):
+            #find the second derivative of f
+            #by varying 2 components of eps
+            eps_up = eps.copy()
+            eps_down = eps.copy()
+            eps_up[j]+=self.de
+            eps_down[j]-=self.de
+            return (find_grad(eps_up,i,calc,unitcell)-find_grad(eps_down,i,calc,unitcell))/(2*self.de)
+
+
+        self.de = de
+        unitcell = self.lattice(self.el, latticeconstant=self.a0, size=[1, 1, 1])
+        unitcell.set_pbc([True,True,True])
+
+        #compute all the terms we need:
+        grad_f = np.zeros([6])
+        hess_f = np.zeros([6,6])
+        eps_vec = np.zeros([6])
+        for i in range(6):
+            grad_f[i] = find_grad(eps_vec,i,calc,unitcell)
+            #finds the gradient of the gradient of f with e[i] with respect to each strain component e[j]
+            for j in range(6-i):
+                hess_f[i,j] = find_second_deriv(eps_vec,i,j,calc,unitcell)
+                #include relevant symmetries
+                hess_f[j,i] = hess_f[i,j]
+
+        self.grad_f = grad_f
+        self.hess_f = hess_f
+
+
+    def save_taylor(self):
+        '''Save the results of the taylor expansion in to numpy readable files rather than re-calculating each time'''
+        try:
+            np.save('grad_f.npy',self.grad_f)
+            np.save('hess_f.npy',self.hess_f)
+        except AttributeError:
+            print('Need to fit a model with fit taylor before it can be saved!')
+
+    def load_taylor(self):
+        '''Read the results of the taylor expansion rather than re-calculating each time'''
+        try:
+            self.grad_f = np.load('grad_f.npy')
+            self.hess_f = np.load('hess_f.npy')
+        except FileNotFoundError:
+            print('No files found containing model parameters - make sure to call save_taylor() to save the model.')
+        
+        
+    def F_3D_from_atoms(self,atoms,F_func,coordinates='cart3D',cell=None,*args,**kwargs):
+        '''Generate the 3D deformation gradient tensor field (F) for an atoms object in a given coordinate system.
+        Parameters
+        ----------
+        atoms : ASE atoms object
+            Atoms object where atoms sit at the coordinates needed to generate the deformation gradient tensor field
+            from F_func.
+        F_func : Function
+            Function to calculate the deformation gradient tensor field from a set of coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates' and the function also accepts anything extra specified in *args and
+            **kwargs. Function must return the deformation gradient field for the atoms in form of a numpy array shape [natoms,ndims,ndims],
+            where F is specified in cartesian coordinates.
+        coordinates: string
+            The coordinates system that F_func accepts. Must be 'cylind2D', 'cylind3D', 'spherical', 'cart2D', 'cart3D'
+        cell : array_like, optional
+            An optional argument allowing the user to specify the cell parameters used to define the coordinate system.
+
+        Returns
+        -------
+        F3D : array_like
+            The 3D deformation gradient tensor field for all the atoms in the system, expressed in the lab frame.
+        '''
+        natoms = len(atoms)
+        if cell is None:
+            sx, sy, sz = atoms.cell.diagonal()
+        else:
+            [sx,sy,sz] = cell
+
+        if coordinates=='cylind2D':
+            x, y = atoms.positions[:, 0], atoms.positions[:, 1]
+            cx, cy = sx/2, sy/2
+            r = np.sqrt((x - cx)**2 + (y - cy)**2)
+            theta = np.arctan2((y-cy),(x-cx))
+            F_2D = F_func(r,theta,*args,**kwargs)
+            #pad F into a 3x3 matrix
+            F_3D = (np.zeros([natoms,3,3]))
+            F_3D[:,0:2,0:2] = F_2D[:,:,:]
+            F_3D[:,2,2] = 1
+
+        elif coordinates=='cylind3D':
+            x, y, z = atoms.positions[:,0],atoms.positions[:,1],atoms.positions[:,2]
+            cx, cy, cz = sx/2, sy/2, sz/2
+            r = np.sqrt((x - cx)**2 + (y - cy)**2)
+            theta = np.arctan2((y-cy),(x-cx))
+            F_3D = F_func(r,theta,z,*args,**kwargs)
+
+        elif coordinates=='spherical':
+            raise NotImplementedError
+
+        elif coordinates=='cart2D':
+            x, y = atoms.positions[:, 0], atoms.positions[:, 1]
+            F_2D = F_func(x,y,*args,**kwargs)            
+            F_3D = (np.zeros([natoms,3,3]))
+            F_3D[:,0:2,0:2] = F_2D[:,:,:]
+            F_3D[:,2,2] = 1
+        
+        elif coordinates=='cart3D':
+            x, y, z = atoms.positions[:, 0], atoms.positions[:, 1], atoms.positions[:,2]
+            F_3D = F_func(x,y,z,*args,**kwargs)
+            
+        else:
+            print('Coordinate system should be cart, cylind2D, cylind3D, spherical')
+            raise NotImplementedError
+
+        return F_3D
+
+    
+    def evaluate_F_or_E(self,A,atoms,F_func=None,E_func=None,cell=None,coordinates='cart3D',*args,**kwargs):
+        '''Get the deformation gradient tensor field or Green-Lagrange strain tensor field for a system of atoms in the lab frame,
+           depending on which function is provided. From these, find the Green-Lagrange strain tensor field in the lattice frame.
+        ----------
+        A : 3x3 numpy array, floats
+            rotation matrix of the form [x^T,y^T,z^T], where x^T,y^T,z^T are normalised column vectors of the lab frame directions expressed
+            in terms of the crystal lattice directions.
+        atoms : ASE atoms object
+            Atoms object where atoms sit at the coordinates needed to generate the deformation gradient tensor field
+            from F_func.
+        F_func : Function
+            Function to calculate the deformation gradient tensor field from a set of atomic coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates' and the function also accepts anything extra specified in *args and
+            **kwargs. Function must return the deformation gradient field for the atoms in form of a numpy array shape [natoms,ndims,ndims],
+            where F is specified in cartesian coordinates. The returned tensor field must be in the lab frame
+            Optional, but one of F_func or E_func must be provided.
+        E_func : Function
+            Function to calculate the Green-Lagrange tensor field from a set of atomic coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates' and the function also accepts anything extra specified in *args and
+            **kwargs. Function must return the Green-Lagrange tensor field for the atoms in form of a numpy array shape [natoms,ndims,ndims],
+            where E is specified in cartesian coordinates, and the returned tensor field is in the lab frame.
+            Optional, but one of F_func or E_func must be provided.
+        cell : array_like, optional
+            An optional argument allowing the user to specify the cell parameters used to define the coordinate system.
+
+        Returns
+        -------
+        E : array_like
+            The 3D Green-Lagrange tensor field for all the atoms in the system, expressed in the lattice frame.
+        R : array_like
+            The rotation tensor component of the deformation gradient tensor field applied to the system, according to the 
+            polar decomposition F = RU. See https://en.wikipedia.org/wiki/Finite_strain_theory#Polar_decomposition_of_the_deformation_gradient_tensor.
+        '''
+        natoms = len(atoms)
+
+        #in the case that only the F function is provided
+        if F_func is not None:
+            if E_func is not None:
+                print('Need to only provide one of E or F, not both')
+                raise ValueError
+            F_3D = self.F_3D_from_atoms(atoms,F_func,cell=cell,coordinates=coordinates,*args,**kwargs)
+            
+
+            #transform F from lab frame to lattice frame,
+            #find right handed stretch tensor U and rotation matrix
+            #R via polar decomposition, then use these to find 
+            #green-lagrange strain tensors E.
+
+            Fprime = np.zeros_like(F_3D)
+            R = np.zeros_like(F_3D)
+            U = np.zeros_like(F_3D)
+            E = np.zeros_like(F_3D)
+            for i in range(natoms):
+                Fprime[i,:,:] = A@F_3D[i,:,:]@np.transpose(A)
+                R[i,:,:],U[i,:,:] = polar(Fprime[i,:,:])
+                E[i,:,:] = 0.5*(U[i,:,:]@(np.transpose(U[i,:,:])) - np.eye(3))
+            
+        #in the case that only the E function is provided, find and rotate the Green-Lagrange strain tensor field directly
+        elif E_func is not None:
+            x, y, z = atoms.positions[:, 0], atoms.positions[:, 1], atoms.positions[:,2]
+            E_3D_lab = E_func(x,y,z,*args,**kwargs)
+            E = np.zeros_like(E_3D_lab)
+            R = np.zeros_like(E_3D_lab)
+            #for a specified strain, R is just I for all atoms
+            for i in range(natoms):
+                E[i,:,:] = A@E_3D_lab[i,:,:]@np.transpose(A)
+                R[i,:,:] = np.eye(3)
+            
+        else:
+            print('Error, please provide one of E or F')
+            raise ValueError
+        
+        return E, R
+
+
+    def predict_shifts(self,A,atoms,F_func=None,E_func=None,coordinates='cart3D',method='taylor',*args,**kwargs):
+        '''
+        F_func - function to calculate the deformation field tensors at
+        a vector of x,y,z, r,theta,(z), or r, theta, phi based on what
+        is passed to coordinates. Note that F must return a vector of shape
+        [number of atoms, 2, 2] for 2D, and [number of atoms, 3, 3] for 3D,
+        and the returned matrix must be expressed in cartesian coordinates. Optional, but must be 
+        given if E is not specified.
+
+        E_func - function to directly calculate the get the strain tensors applied to each atom
+        in the case where this is easier than F. Strain tensors must be returned of the form
+        [number of atoms, 3, 3]. Optional, but must be given if F is not specified. 
+
+        A - a matrix which holds the axes of the lab frame described in
+        the lattice frame, of the form (x^T,y^T,z^T), where x,y and z are the
+        normalised lab frame axes defined with respect to the cubic lattice basis.
+        This acts as a rotation matrix. For example
+        if our lab frame axis were x, y, z, which were defined
+        by x = 1/sqrt(2)[110], y = 1/sqrt(2)[-110], z = [001] with respect to our crystal lattice
+        A would be:
+        [[1/sqrt(2)   -1/sqrt(2)     0]
+         [1/sqrt(2)    1/sqrt(2)     0]
+         [0            0             1]]
+
+        coordinates - the form of atomic coordinates that need to be passed
+        to the function F. Can be cart2D, cart3D, cylind2D, cylind3D, spherical.
+        
+        *args, **kwargs - anything to additionally be passed to F_func or E_func
+        '''
+
+        E,R = self.evaluate_F_or_E(A,atoms,F_func=F_func,E_func=E_func,coordinates=coordinates,*args,**kwargs)
+        natoms = len(atoms)
+
+
+        #get the cauchy born shifts unrotated
+        shifts_no_rr = self.evaluate_shift_model(E,method=method)
+
+        shifts = np.zeros_like(shifts_no_rr)
+
+        #rotate the cauchy shifts both by the rotation induced by F
+        #and to get them back into the lab frame
+        for i in range(natoms):
+            #might have to transpose R here? seems to greatly reduce errors - no idea why though.
+            shifts[i,:] = np.transpose(A)@R[i,:,:]@shifts_no_rr[i,:]
+            #shifts[i,:] = np.transpose(A)@shifts_no_rr[i,:]
+
+
+        return shifts
+
+    def evaluate_shift_model(self,E,method='taylor'):
+        #takes E in the form of a list of matrices [natoms, :, :]
+        #define function for evaluating model
+
+        def permutation(strain_vec,perm_shift):
+            strain_vec_perm = np.zeros_like(strain_vec)
+            for i in range(3):
+                perm_num = (i-perm_shift)%3
+                strain_vec_perm[:,perm_num] = strain_vec[:,i]
+                strain_vec_perm[:,perm_num+3] = strain_vec[:,i+3]
+            return strain_vec_perm
+        
+        def evaluate_shift_lasso(eps):
+            epsx = eps
+            epsy = permutation(eps,1)
+            epsz = permutation(eps,2)
+            predictions = np.zeros([np.shape(eps)[0],3])
+            predictions[:,0] = self.LM.predict(self.basis_function_evaluation(epsx))
+            predictions[:,1] = self.LM.predict(self.basis_function_evaluation(epsy))
+            predictions[:,2] = self.LM.predict(self.basis_function_evaluation(epsz))
+            return predictions
+
+        def eval_tay_model(self,eps):
+            grad_f = self.grad_f
+            hess_f = self.hess_f
+            #takes E as a 2D array of voigt vectors, with number of rows=number atoms
+            # Note that the permutation is computed from observing how the matrix components permute
+            # when transforming the strain matrix by the 120 degree rotation matrices about [111],
+            # such that the axis are permuted in the order
+            # 1 ---> 2 ---> 3
+
+            
+            def eval_eps(eps_vec,grad_f,hess_f):
+                term_1 = grad_f@eps_vec
+                term_2 = (1/2)*(np.transpose(eps_vec)@hess_f@eps_vec)
+
+                return term_1 + term_2
+
+
+            epsx = eps
+            epsy = permutation(eps,1)
+            epsz = permutation(eps,2)
+            predictions = np.zeros([np.shape(eps)[0],3])
+            for i in range(np.shape(eps)[0]):
+                predictions[i,0] = eval_eps(epsx[i,:],grad_f,hess_f)
+                predictions[i,1] = eval_eps(epsy[i,:],grad_f,hess_f)
+                predictions[i,2] = eval_eps(epsz[i,:],grad_f,hess_f)
+            
+            return predictions
+
+
+        #turn E into voigt vector
+        E_voigt = np.zeros([np.shape(E)[0],6])
+        E_voigt[:,0] = E[:,0,0]
+        E_voigt[:,1] = E[:,1,1]
+        E_voigt[:,2] = E[:,2,2]
+        E_voigt[:,3] = E[:,1,2]
+        E_voigt[:,4] = E[:,0,2]
+        E_voigt[:,5] = E[:,0,1]
+
+
+        #return the shift vectors
+        if method == 'taylor':
+            return eval_tay_model(self,E_voigt)
+        elif method == 'lasso':
+            return evaluate_shift_lasso(E_voigt)
+        else:
+            print('Error! Can only predict cauchy born shift with implemented errors - "taylor" or "lasso"')
+            raise NotImplementedError
+
+    def apply_shifts(self,atoms,shifts,mask=None):
+
+        #mask covers atoms that we want to be shifted
+        if mask is not None:
+            lattice1mask = np.logical_and(self.lattice1mask,mask)
+            lattice2mask = np.logical_and(self.lattice2mask,mask)
+        else:
+            lattice1mask = self.lattice1mask
+            lattice2mask = self.lattice2mask
+        atoms.positions[lattice1mask] += -0.5*shifts[lattice1mask]
+        atoms.positions[lattice2mask] += 0.5*shifts[lattice2mask]
+
+
+    def get_cb_error(self,atoms,mask=None,forces=None):
+        atoms.set_calculator = self.calc
+        if forces is None:
+            forces = atoms.get_forces()
+
+        if mask is not None:
+            lattice1mask = np.logical_and(self.lattice1mask,mask)
+            lattice2mask = np.logical_and(self.lattice2mask,mask)
+        else:
+            lattice1mask = self.lattice1mask
+            lattice2mask = self.lattice2mask
+        
+        lattice1_force_vec = forces[lattice1mask]
+        lattice2_force_vec = forces[lattice2mask]
+        lattice1_force = np.sum(lattice1_force_vec)
+        lattice2_force = np.sum(lattice2_force_vec)
+        
+        cb_err = np.abs(lattice1_force-lattice2_force)
+
+        return cb_err #we should be trying to minimise this
+
+
+    
+    def initial_lasso_fit(self,alpha=0.1,initial_samples=10):
+        #if no lasso model has yet been fitted
+        #initialise the Lasso model
+        #self.LM = Lasso(alpha=alpha,warm_start=True,fit_intercept=False)
+        self.LM = LinearRegression(fit_intercept=False)
+        #build the initial simple dataset to begin fitting with
+        sampler = LatinHypercube(d=6)
+        samples = sampler.random(n=initial_samples)
+        lbounds = [-0.01,-0.01,-0.01,-0.01,-0.01,-0.01]
+        ubounds = [0.01,0.01,0.01,0.01,0.01,0.01]
+        samples = scale(samples,lbounds,ubounds)
+        self.lasso_data = samples
+        self.lasso_phi = None
+        E_vecs,shifts = self.get_data_points(samples)
+        #get design matrix
+        phi = self.basis_function_evaluation(E_vecs)
+        self.lasso_phi = phi
+        self.lasso_shifts = shifts
+        #print('shifts',np.shape(shifts))
+        print('Performing initial lasso fit for CB corrector.....')
+        self.LM.fit(self.lasso_phi,self.lasso_shifts.flatten())# fit the Lasso model to the given data
+        print('Fit completed')
+
+    def check_for_refit(self,A,atoms,forces,tol=1e-3,mask=None,E_func=None,F_func=None,\
+        coordinates='cart3D',refit_points=10,err_vec=None,*args,**kwargs):
+        #tolerance is error per atom
+        #multiply by number of atoms to get full tolerance
+
+        tol = tol*len(atoms[mask])
+        print('tol',tol)
+        
+        cb_err = self.get_cb_error(atoms,forces=forces,mask=mask)
+        print('cb_err',cb_err)
+        if abs(cb_err)<tol: #do not refit in this case
+            return 0, cb_err/len(atoms)
+
+
+        #otherwise, find the lattice 1 and lattice 2 atoms with the highest contribution to the cauchy-born errors,
+        #add these to the training dataset and refit the model.
+        if mask is not None:
+            lattice1mask = np.logical_and(self.lattice1mask,mask)
+            lattice2mask = np.logical_and(self.lattice2mask,mask)
+        else:
+            lattice1mask = self.lattice1mask
+            lattice2mask = self.lattice2mask
+        
+        if err_vec is None:
+            #get force magnitude vectors, only do this when no vector of errors is provided
+            force_mag = np.linalg.norm(forces,axis=1)
+        else:
+            force_mag = err_vec
+        
+        atom_numbers = np.array([i for i in range(len(atoms))])
+        
+        l1points = int(refit_points/2)
+        l2points = refit_points-l1points
+        #get random lattice 1 force 
+        #FIXME at some point, would be great if these could select the 
+        # points that the model is describing worst rather than random points.
+
+        force_mag_lattice1 = force_mag[lattice1mask]
+        atom_nos_lattice1 = atom_numbers[lattice1mask]
+        top_idx = np.argpartition(force_mag_lattice1,-l1points)[-l1points:]
+        high_force_atoms1 = atom_nos_lattice1[top_idx]
+        #high_force_atoms1 = np.random.choice(atom_nos_lattice1,l1points)
+        #get 5 largest lattice 2 force indices
+        force_mag_lattice2 = force_mag[lattice2mask]
+        atom_nos_lattice2 = atom_numbers[lattice2mask]
+        top_idx = np.argpartition(force_mag_lattice2,-l2points)[-l2points:]
+        high_force_atoms2 = atom_nos_lattice2[top_idx]
+        #high_force_atoms2 = np.random.choice(atom_nos_lattice2,l2points)
+
+
+        high_force_atoms = np.concatenate((high_force_atoms1,high_force_atoms2))
+        #print('atoms selected:', high_force_atoms)
+        #print('forces on atoms:', err_vec[high_force_atoms])
+        reduced_atoms = atoms[high_force_atoms]
+
+        sx, sy, sz = atoms.cell.diagonal()
+        cell = [sx,sy,sz]
+
+        #evaluate E and R for reduced list of atoms
+        E,R = self.evaluate_F_or_E(A,reduced_atoms,F_func=F_func,E_func=E_func,cell=cell,coordinates=coordinates,*args,**kwargs)
+        #print('Es',E)
+        #we get back the strain vectors in the lattice frame 
+        #turn E into voigt vectors
+        E_voigt = np.zeros([np.shape(E)[0],6])
+        E_voigt[:,0] = E[:,0,0]
+        E_voigt[:,1] = E[:,1,1]
+        E_voigt[:,2] = E[:,2,2]
+        E_voigt[:,3] = E[:,1,2]
+        E_voigt[:,4] = E[:,0,2]
+        E_voigt[:,5] = E[:,0,1]
+        
+        #pass set of Es to evaluate to refit_lasso_regression
+        self.refit_lasso_regression(atoms,E_voigt)
+        return 1, cb_err/len(atoms)
+        
+
+    def refit_lasso_regression(self,atoms,E_voigt):
+        calc = self.calc
+        #take our given E voigt and evaluate the model
+        E_vecs,shifts = self.get_data_points(E_voigt)
+        #print('evecs,shifts',E_vecs,shifts)
+        #get design matrix
+        phi = self.basis_function_evaluation(E_vecs)
+        temp_predict = self.LM.predict(phi)
+        #print('predicted shifts', temp_predict)
+        #print('predicted shifts - true shifts', shifts.flatten()-temp_predict)
+        if self.lasso_phi is None:
+            self.lasso_phi = phi
+            self.lasso_shifts = shifts
+        else:
+            self.lasso_phi = np.concatenate((self.lasso_phi,phi),axis=0)
+            self.lasso_shifts = np.concatenate((self.lasso_shifts,shifts),axis=0)
+
+        print('High error detected, refitting CB corrector.....')
+        self.LM.fit(self.lasso_phi,self.lasso_shifts.flatten())# fit the Lasso model to the given data
+        print('Re-fit completed')
+    
+    def get_data_points(self,E_vec):
+        unitcell = self.lattice(self.el, latticeconstant=self.a0, size=[1, 1, 1])
+        unitcell.set_pbc([True,True,True])
+        def permutation(strain_vec,perm_shift):
+            strain_vec_perm = np.zeros_like(strain_vec)
+            for i in range(3):
+                perm_num = (i-perm_shift)%3
+                strain_vec_perm[:,perm_num] = strain_vec[:,i]
+                strain_vec_perm[:,perm_num+3] = strain_vec[:,i+3]
+            return strain_vec_perm
+        
+        ndp = np.shape(E_vec)[0] #number of data points
+        shift_vals = np.zeros([3*ndp,1])
+        for i,E in enumerate(E_vec):
+            shift_diff = self.eval_shift(E,unitcell)
+            shift_vals[i] = shift_diff[0]
+            shift_vals[i+ndp] = shift_diff[1]
+            shift_vals[i+2*ndp] = shift_diff[2]
+        
+        eps = E_vec
+        epsx = eps
+        epsy = permutation(eps,1)
+        epsz = permutation(eps,2)
+
+        strains = np.zeros([3*ndp,6])
+        strains[0:ndp,:] = epsx
+        strains[ndp:2*ndp,:] = epsy
+        strains[2*ndp:3*ndp,:] = epsz
+
+
+        return strains,shift_vals
+        
+    def eval_shift(self,E_vec,unitcell):
+        #green lagrange version of function
+        #turn E into matrix
+        E = np.zeros([3,3])
+        E[0,0] = E_vec[0]
+        E[1,1] = E_vec[1]
+        E[2,2] = E_vec[2]
+        E[2,1],E[1,2] = E_vec[3],E_vec[3]
+        E[2,0],E[0,2] = E_vec[4],E_vec[4]
+        E[0,1],E[1,0] = E_vec[5],E_vec[5]
+        # get U^2
+        Usqr = 2*E + np.eye(3)
+        #square root matrix
+        U = sqrtm(Usqr,disp=True)
+
+        #this is just the symmetric stretch tensor, exactly what we need.
+        x = U
+
+
+        initial_shift = np.zeros([3])
+        relaxed_shift = np.zeros([3])
+
+        #build cell
+        unitcell_copy = unitcell.copy()
+        unitcell_copy.calc = self.calc
+        cell = unitcell_copy.get_cell()
+        cell_rescale = x[:,:]@cell
+        unitcell_copy.set_cell(cell_rescale,scale_atoms=True)
+        initial_shift[:] = unitcell_copy.get_positions()[1] - unitcell_copy.get_positions()[0]
+
+        #relax cell
+        opt = LBFGS(unitcell_copy,logfile=None)
+        opt.run(fmax = 1e-10)
+        relaxed_shift[:] = unitcell_copy.get_positions()[1] - unitcell_copy.get_positions()[0]
+
+        shift_diff = relaxed_shift-initial_shift
+        #TRANSFER SHIFT_DIFF BACK INTO THE LATTICE FRAME, OUT OF THE DEFORMED FRAME
+        #back_transform = np.transpose(cell_rescale)
+        #shift_diff_transform = back_transform@shift_diff
+        #print('shift_diff',shift_diff,'shift_diff_transform',shift_diff_transform)
+        return shift_diff #return all 3 components
+
+
+
+    def basis_function_evaluation(self,E_vecs):
+        #first basis function is just the e vector components evaluated, 2nd is all of the cross terms
+        triu_indices = np.triu_indices(6) #get upper triangle indices
+        phi_2nd_term = np.zeros([np.shape(E_vecs)[0],np.shape(triu_indices)[1]])
+        for i,E in enumerate(E_vecs):
+            e_outer_p = np.outer(E,E)
+            phi_2nd_term[i,:] = e_outer_p[triu_indices]
+        
+        #return the design matrix
+
+        #np.savetxt('design_matrix.txt',np.concatenate((E_vecs,phi_2nd_term),axis=1))
+        return np.concatenate((E_vecs,phi_2nd_term),axis=1)
+        
+        
+            
+
diff --git a/matscipy/contact_mechanics/DMT.py b/matscipy/contact_mechanics/DMT.py
deleted file mode 100644
index abcacc15..00000000
--- a/matscipy/contact_mechanics/DMT.py
+++ /dev/null
@@ -1,78 +0,0 @@
-#
-# Copyright 2014, 2021 Lars Pastewka (U. Freiburg)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-import math
-
-import numpy as np
-
-import matscipy.contact_mechanics.Hertz as Hertz
-
-###
-
-def radius_and_pressure(N, R, Es, w):
-    """
-    Given normal load, sphere radius and contact modulus compute contact radius
-    and peak pressure.
-
-    Parameters
-    ----------
-    N : float
-        Normal force.
-    R : float
-        Sphere radius.
-    Es : float
-        Contact modulus: Es = E/(1-nu**2) with Young's modulus E and Poisson
-        number nu.
-    w : float
-        Work of adhesion.
-    """
-
-    return Hertz.radius_and_pressure(N+2*math.pi*w*R, R, Es)
-
-
-def surface_stress(r, a, nu):
-    """
-    Given distance from the center of the sphere, contact radius and Poisson
-    number contact, compute the stress at the surface.
-
-    Parameters
-    ----------
-    r : array_like
-        Array of distance (from the center of the sphere).
-    a : float
-        Contact radius.
-    nu : float
-        Poisson number.
-
-    Returns
-    -------
-    pz : array
-        Contact pressure.
-    sr : array
-        Radial stress.
-    stheta : array
-        Azimuthal stress.
-    """
-
-    return Hertz.surface_stress(r, a, nu)
-    
-    
-def surface_displacements(r, a):
-    return Hertz.surface_displacements(r, a)
\ No newline at end of file
diff --git a/matscipy/contact_mechanics/Hertz.py b/matscipy/contact_mechanics/Hertz.py
deleted file mode 100644
index 92cf0cc6..00000000
--- a/matscipy/contact_mechanics/Hertz.py
+++ /dev/null
@@ -1,405 +0,0 @@
-#
-# Copyright 2014-2015, 2021 Lars Pastewka (U. Freiburg)
-#           2015 Till Junge (EPFL)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-import math
-
-import numpy as np
-
-###
-
-def radius_and_pressure(N, R, Es):
-    """
-    Given normal load, sphere radius and contact modulus compute contact radius
-    and peak pressure.
-
-    Parameters
-    ----------
-    N : float
-        Normal force.
-    R : float
-        Sphere radius.
-    Es : float
-        Contact modulus: Es = E/(1-nu**2) with Young's modulus E and Poisson
-        number nu.
-
-    Returns
-    -------
-    a : float
-        Contact radius.
-    p0 : float
-        Maximum pressure inside the contacting area (right under the apex).
-    """
-
-    a = R*(3./4*( N/(Es*R**2) ))**(1./3)
-    p0 = 3*N/(2*math.pi*a*a)
-    
-    return a, p0
-
-
-def surface_stress(r, poisson=0.5):
-    """
-    Given distance from the center of the sphere, contact radius and Poisson
-    number contact, compute the stress at the surface.
-
-    Parameters
-    ----------
-    r : array_like
-        Array of distance (from the center of the sphere in units of contact
-        radius a).
-    poisson : float
-        Poisson number.
-
-    Returns
-    -------
-    pz : array
-        Contact pressure (in units of maximum pressure p0).
-    sr : array
-        Radial stress (in units of maximum pressure p0).
-    stheta : array
-        Azimuthal stress (in units of maximum pressure p0).
-    """
-
-    mask0 = np.abs(r) < 1e-6
-    maski = np.logical_and(r < 1.0, np.logical_not(mask0))
-    masko = np.logical_and(np.logical_not(maski), np.logical_not(mask0))
-    r_0 = r[mask0]
-    r_i = r[maski]
-    r_o = r[masko]
-    
-    # Initialize
-    pz = np.zeros_like(r)
-    sr = np.zeros_like(r)
-    stheta = np.zeros_like(r)
-    
-    # Solution at r=0
-    if mask0.sum() > 0:
-        pz[mask0] = np.ones_like(r_0)
-        sr[mask0] = -(1.+2*poisson)/2.*np.ones_like(r_0)
-        stheta[mask0] = -(1.+2*poisson)/2.*np.ones_like(r_0)
-
-    # Solution inside the contact radius
-    if maski.sum() > 0:
-        r_a_sq = r_i**2
-        pz[maski] = np.sqrt(1-r_a_sq)
-        sr[maski] = (1.-2.*poisson)/(3.*r_a_sq)*(1.-(1.-r_a_sq)**(3./2))- \
-            np.sqrt(1.-r_a_sq)
-        stheta[maski] = -(1.-2.*poisson)/(3.*r_a_sq)*(1.-(1.-r_a_sq)**(3./2))- \
-            2*poisson*np.sqrt(1.-r_a_sq)
-    
-    # Solution outside of the contact radius
-    if masko.sum() > 0:
-        r_a_sq = r_o**2
-        po = (1.-2.*poisson)/(3.*r_a_sq)
-        sr[masko] = po
-        stheta[masko] = -po
-    
-    return pz, sr, stheta
-
-
-def centerline_stress(z, poisson=0.5):
-    """
-    Given distance from the center of the sphere, contact radius and Poisson
-    number contact, compute the stress at the surface.
-
-    Parameters
-    ----------
-    z : array_like
-        Array of depths (from the center of the sphere in units of contact
-        radius a).
-    poisson : float
-        Poisson number.
-
-    Returns
-    -------
-    srr : array
-        Radial stress (in units of maximum pressure p0).
-    szz : array
-        Contact pressure (in units of maximum pressure p0).
-    """
-
-    srr = -(1.+poisson)*(1.-z*np.arctan(1./z)) + 1./(2.*(1.+z**2))
-    szz = -1./(1.+z**2)
-
-    return srr, szz
-
-    
-def surface_displacements(r):
-    """
-    Return the displacements at the surface due to an indenting sphere.
-    See: K.L. Johnson, Contact Mechanics, p. 61
-
-    Parameters
-    ----------
-    r : array_like
-        Radial position normalized by contact radius a.
-
-    Returns
-    -------
-    uz : array
-        Normal displacements at the surface of the contact (in units of
-        p0/Es * a where p0 is maximum pressure, Es contact modulus and a 
-        contact radius).
-    """
-
-    maski = r < 1.0
-    masko = np.logical_not(maski)
-    r_i = r[maski]
-    r_o = r[masko]
-    
-    # Initialize
-    uz = np.zeros_like(r)
-    
-    # Solution inside the contact circle
-    if maski.sum() > 0:
-        uz[maski] = -math.pi*(2.-r_i**2)/4.
-
-    # Solution outside the contact circle
-    if masko.sum() > 0:
-        uz[masko] = (-(2.-r_o**2)*np.arcsin(1./r_o) - 
-            r_o*np.sqrt(1.-(1./r_o)**2))/2.
-
-    return uz
-
-
-def stress(r, z, poisson=0.5):
-    """
-    Return components of the stress tensor in the interior of the Hertz solid.
-    This is the solution given by: M.T. Huber, Ann. Phys. 319, 153 (1904)
-
-    Note that the stress tensor at any point in the solid has the form below in
-    a cylindrical coordinate system centered at the tip apex. Some off-diagonal
-    components are zero by rotational symmetry. stt is the circumferential
-    component, srr the radial component and szz the normal component of the
-    stress tensor.
-
-            / stt  0   0  \ 
-        s = |  0  srr srz |
-            \  0  srz szz /
-
-    Parameters
-    ----------
-    r : array_like
-        Radial position (in units of the contact radius a).
-    z : array_like
-        Depth (in units of the contact radius a).
-    poisson: float
-        Poisson number.
-
-    Returns
-    -------
-    stt : array
-        Circumferential component of the stress tensor (in units of maximum
-        pressure p0).
-    srr : array
-        Radial component of the stress tensor (in units of maximum pressure p0).
-    szz : array
-        Normal component of the stress tensor (in units of maximum pressure p0).
-    srz : array
-        Shear component of the stress tensor (in units of maximum pressure p0).
-    """
-
-    p = r**2+z**2-1
-    u = p/2 + np.sqrt(p**2/4+z**2)
-    sqrtu = np.sqrt(u)
-
-    # Note: The equation for u is
-    #     r**2/(1+u) + z**2/u = 1
-    # With r' = r/sqrt(1+u) and z' = z/sqrt(u) this expression becomes
-    #     r'**2 + z'**2 = 1, hence z' = sqrt(1-r'**2)
-    # Variables r and z below are substituted by r->r' and z->z'.
-    r = r/np.sqrt(1+u)
-    z = np.sqrt(1-r**2) # equiv. to z=u/sqrt(u), but defined for z=0
-
-    # Precompute arctan
-    sqrtu_arctan_inv_sqrtu = sqrtu*np.arctan(1./sqrtu)
-
-    # The next two expressions give numerical problems at the tip center and the
-    # contact edge, regularize with the asymptotic value.
-
-    # Note: (1-z**3)/r**2->3/2 for r->0
-    one_minus_z3_div_r2 = 3./2.*np.ones_like(r)
-    mask = r > 0
-    one_minus_z3_div_r2[mask] = (1.-z[mask]**3)/r[mask]**2
-
-    # Note: z**2/(u+z**2)->1 for u+z**2->0
-    u_plus_z2 = u+z**2
-    z2_div_u_plus_z2 = np.where(u_plus_z2 > 0., z**2/u_plus_z2, np.ones_like(z))
-
-    # Compute stresses
-    # Note: the factor 1/(1+u) stems from the substitution r->r' and z->z' above
-    stt = (1.-2.*poisson)/3. * 1./(1.+u) * one_minus_z3_div_r2 + \
-        z*(2.*poisson+ (1.-poisson)*u/(1.+u) - (1.+poisson)*sqrtu_arctan_inv_sqrtu)
-    szz = z*z2_div_u_plus_z2
-    srr = -( (1.-2.*poisson)/3. * 1./(1+u) * one_minus_z3_div_r2 + \
-         z*z2_div_u_plus_z2 + \
-         z*((1.-poisson)*u/(1.+u) + (1.+poisson)*sqrtu_arctan_inv_sqrtu - 2.) )
-    srz = r*z2_div_u_plus_z2 * sqrtu/np.sqrt(1.+u)
-
-    return -stt, -srr, -szz, -srz
-
-
-def stress_Cartesian(x, y, z, poisson=0.5):
-    """
-    Return components of the stress tensor in the interior of solid due to
-    normal Hertz loading.
-    This is the solution given by:
-    G.M. Hamilton, Proc. Instn. Mech. Engrs. 197C, 53-59 (1983)
-
-    Parameters
-    ----------
-    x, y : array_like
-        In-plane positions (in units of the contact radius a).
-    z : array_like
-        Depth (in units of the contact radius a).
-    poisson : float
-        Poisson number.
-
-    Returns
-    -------
-    sxx, syy, szz, syz, sxz, sxy : array
-        Individual components of the Cartesian stress tensor.
-    """
-
-    def stress_offcenter(x, y, z, r_sq, poisson=0.5):
-        A = r_sq + z**2 - 1
-        S = np.sqrt(A**2 + 4*z**2)
-
-        M = np.sqrt((S+A)/2)
-        N = np.sqrt((S-A)/2)
-        phi = np.arctan2(1, M)
-
-        G = M**2 - N**2 + z*M - N
-        H = 2*M*N + M + z*N
-
-        sxx = (1+poisson)*z*phi+1/r_sq*((y**2-x**2)/r_sq*((1-poisson)*N*z**2-(1-2*poisson)/3*(N*S+2*A*N+1)-poisson*M*z)-N*(x**2+2*poisson*y**2)-M*x**2*z/S)
-        syy = (1+poisson)*z*phi+1/r_sq*((x**2-y**2)/r_sq*((1-poisson)*N*z**2-(1-2*poisson)/3*(N*S+2*A*N+1)-poisson*M*z)-N*(y**2+2*poisson*x**2)-M*y**2*z/S)
-        szz = -N+z*M/S
-        sxy = x*y*(1-2*poisson)/r_sq**2*(-N*r_sq+2/3*N*(S+2*A)-z*(z*N+M)+2/3)+x*y*z/r_sq**2*(-M*r_sq/S-z*N+M)
-        syz = -z*(y*N/S-y*z*H/(G**2+H**2))
-        sxz = -z*(x*N/S-x*z*H/(G**2+H**2))
-
-        return sxx, syy, szz, syz, sxz, sxy
-
-    x = np.asarray(x)
-    y = np.asarray(y)
-    z = np.asarray(z)
-
-    sxx = np.zeros_like(x)
-    syy = np.zeros_like(x)
-    szz = np.zeros_like(x)
-    syz = np.zeros_like(x)
-    sxz = np.zeros_like(x)
-    sxy = np.zeros_like(x)
-
-    r_sq = x**2 + y**2
-
-    mask = r_sq > 0
-    xx, yy, zz, yz, xz, xy = stress_offcenter(x[mask], y[mask], z[mask],
-                                              r_sq[mask], poisson=poisson)
-    sxx[mask] = xx
-    syy[mask] = yy
-    szz[mask] = zz
-    syz[mask] = yz
-    sxz[mask] = xz
-    sxy[mask] = xy
-
-    mask = np.logical_not(mask)
-    if mask.sum() > 0:
-        z = z[mask]
-        sxx[mask] = ((1+poisson)*(z*np.arctan2(1, z)-1)+1/(2*(1+z**2)))
-        syy[mask] = sxx[mask]
-        szz[mask] = -1/(1+z**2)
-
-    return sxx, syy, szz, syz, sxz, sxy
-
-
-def stress_for_tangential_loading(x, y, z, poisson=0.5):
-    """
-    Return components of the stress tensor in the interior of solid due to
-    tangential (Hertz) loading.
-    This is the solution given by:
-    G.M. Hamilton, Proc. Instn. Mech. Engrs. 197C, 53-59 (1983)
-
-    Parameters
-    ----------
-    x, y : array_like
-        In-plane positions (in units of the contact radius a).
-    z : array_like
-        Depth (in units of the contact radius a).
-    poisson : float
-        Poisson number.
-
-    Returns
-    -------
-    sxx, syy, szz, syz, sxz, sxy : array
-        Individual components of the Cartesian stress tensor.
-    """
-
-    def stress_offcenter(x, y, z, r_sq, poisson=0.5):
-        A = r_sq + z**2 - 1
-        S = np.sqrt(A**2 + 4*z**2)
-
-        M = np.sqrt((S+A)/2)
-        N = np.sqrt((S-A)/2)
-        phi = np.arctan2(1, M)
-
-        G = M**2 - N**2 + z*M - N
-        H = 2*M*N + M + z*N
-
-        sxx = -x*(poisson/4+1)*phi+x*M/r_sq**2*((3/2-2*x**2/r_sq)*(S*poisson-2*A*poisson+z**2)+x**2*z**2/S+7*poisson*r_sq/4-2*poisson*x**2+r_sq)+x*z*N/r_sq**2*((3/2-2*x**2/r_sq)*(-S/6*(1-2*poisson)-A/3*(1-2*poisson)-1/2*(z**2+3))+x**2/S-poisson*r_sq/4-7*r_sq/4)+4*x*z/(3*r_sq**2)*(3/2-2*x**2/r_sq)*(1-2*poisson)
-        syy = -3*poisson*x*phi/4+x*M/r_sq**2*((1/2-2*y**2/r_sq)*(poisson*(S-2*A+r_sq)+z**2)+y**2*z**2/S+3/4*poisson*r_sq)+z*x*N/r_sq**2*((1/2-2*y**2/r_sq)*(-S/6*(1-2*poisson)-A/3*(1-2*poisson)-z**2/2-3/2)+y**2/S-3/4*poisson*r_sq-r_sq/4)+4/3*z*x/r_sq**2*(1/2-2*y**2/r_sq)*(1-2*poisson)
-        szz = z*x*N/(2*r_sq)*(1-(r_sq+z**2+1)/S)
-        sxy = y/2*(poisson/2-1)*phi+y*M/r_sq**2*(x**2*z**2/S+poisson*((S-2*A)*(1/2-2*x**2/r_sq)-2*x**2+r_sq/4)+r_sq/2+z**2*(1/2-2*x**2/r_sq))+y*z*N/r_sq**2*((1/2-2*x**2/r_sq)*((2*poisson-1)*(S/6+A/3)-z**2/2-3/2-r_sq/2)+r_sq*poisson/4+x**2/S-y**2/2-3*x**2/2)+4*y*z/(3*r_sq**2)*(1/2-2*x**2/r_sq)*(1-2*poisson)
-        syz = x*y*z/(2*r_sq**2)*(M*(1/2+1/S*(z**2/2-3/2-r_sq/2))+z*N/2*(-3+1/S*(5+z**2+r_sq)))
-        sxz = 3*z*phi/2+z*M/r_sq*(1+x**2/r_sq-x**2/S)+N/r_sq*(-3/4*(S+2*A)+z**2-3/4-1/4*r_sq+z**2/2*(1/2-2*x**2/r_sq))
-
-        return sxx, syy, szz, syz, sxz, sxy
-
-    x = np.asarray(x)
-    y = np.asarray(y)
-    z = np.asarray(z)
-
-    sxx = np.zeros_like(x)
-    syy = np.zeros_like(x)
-    szz = np.zeros_like(x)
-    syz = np.zeros_like(x)
-    sxz = np.zeros_like(x)
-    sxy = np.zeros_like(x)
-
-    r_sq = x**2 + y**2
-
-    mask = r_sq > 0
-    xx, yy, zz, yz, xz, xy = stress_offcenter(x[mask], y[mask], z[mask],
-                                              r_sq[mask], poisson=poisson)
-    sxx[mask] = xx
-    syy[mask] = yy
-    szz[mask] = zz
-    syz[mask] = yz
-    sxz[mask] = xz
-    sxy[mask] = xy
-
-    mask = np.logical_not(mask)
-    if mask.sum() > 0:
-        z = z[mask]
-        sxz[mask] += -1+3/2*z*np.arctan2(1, z)-z**2/(2*(1+z**2))
-
-    return sxx, syy, szz, syz, sxz, sxy
-
diff --git a/matscipy/contact_mechanics/JKR.py b/matscipy/contact_mechanics/JKR.py
deleted file mode 100644
index 4b3326f8..00000000
--- a/matscipy/contact_mechanics/JKR.py
+++ /dev/null
@@ -1,48 +0,0 @@
-#
-# Copyright 2014-2015, 2021 Lars Pastewka (U. Freiburg)
-#           2014 James Kermode (Warwick U.)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-import math
-
-import numpy as np
-
-###
-
-def radius(N, R, Es, w):
-    """
-    Given normal load, sphere radius and contact modulus compute contact radius
-    and peak pressure.
-
-    Parameters
-    ----------
-    N : float
-        Normal force.
-    R : float
-        Sphere radius.
-    Es : float
-        Contact modulus: Es = E/(1-nu**2) with Young's modulus E and Poisson
-        number nu.
-    w : float
-        Work of adhesion.
-    """
-
-    return ( 3.*R/(4.*Es)*(N + 6*w*math.pi*R + np.sqrt(12*w*math.pi*R*N + 
-                           (6*w*math.pi*R)**2)) )**(1./3)
-
diff --git a/matscipy/contact_mechanics/__init__.py b/matscipy/contact_mechanics/__init__.py
deleted file mode 100644
index b65615b4..00000000
--- a/matscipy/contact_mechanics/__init__.py
+++ /dev/null
@@ -1,1 +0,0 @@
-# Nothing here yet.
diff --git a/matscipy/contact_mechanics/analysis.py b/matscipy/contact_mechanics/analysis.py
deleted file mode 100644
index 68655fa3..00000000
--- a/matscipy/contact_mechanics/analysis.py
+++ /dev/null
@@ -1,95 +0,0 @@
-#
-# Copyright 2015 Lars Pastewka (U. Freiburg)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-import numpy as np
-
-###
-
-from ..ffi import count_islands, count_segments, distance_map
-
-###
-
-# Stencils for determining nearest-neighbor relationships on a square grid
-nn_stencil = [(1,0),(0,1),(-1,0),(0,-1)]
-nnn_stencil = [(1,1),(1,0),(1,-1),(0,-1),(-1,-1),(-1,0),(-1,1),(0,1)]
-
-###
-
-def outer_perimeter(c, stencil=nn_stencil):
-    """
-    Return a map where the outer perimeter is marked with the patch number.
-
-    Parameters
-    ----------
-    c : array_like
-        2D map of islands.
-
-    Returns
-    -------
-    c : array
-        2D map of island perimeters.
-    """
-    
-    c_nearby = np.zeros_like(c, dtype=bool)
-    for dx, dy in stencil:
-        tmp = c.copy()
-        if dx != 0:
-            tmp = np.roll(tmp, dx, 0)
-        if dy != 0:
-            tmp = np.roll(tmp, dy, 1)
-        c_nearby = np.logical_or(c_nearby, tmp)
-    return np.logical_and(np.logical_not(c), c_nearby)
-
-
-def inner_perimeter(patch_ids, stencil=nn_stencil):
-    """
-    Return a map where the inner perimeter is marked with the patch number
-
-    Parameters
-    ----------
-    c : array_like
-        2D map of islands.
-
-    Returns
-    -------
-    c : array
-        2D map of island perimeters.
-    """
-    
-    c = outer_perimeter(patch_ids == 0, stencil)
-    return c*patch_ids
-
-
-def island_areas(island_ids):
-    """
-    Return a list containing island areas for each island number
-
-    Parameters
-    ----------
-    c : array_like
-        2D map of island ids.
-  
-    Returns
-    -------
-    areas : array
-        Array containing the area for each island id.
-    """
-    
-    return np.bincount(island_ids.reshape((-1,)))[1:]
diff --git a/matscipy/contact_mechanics/greens_function.py b/matscipy/contact_mechanics/greens_function.py
deleted file mode 100644
index 49a77b5a..00000000
--- a/matscipy/contact_mechanics/greens_function.py
+++ /dev/null
@@ -1,398 +0,0 @@
-#
-# Copyright 2014-2015, 2021 Lars Pastewka (U. Freiburg)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-"""
-Contact mechanics module.
-
-This module contains functions creating real and reciprocal space Green's
-function. Continuous real space Green's function can be converted to their
-discrete reciprocal space representation by `real_to_reciprocal_space`.
-Typically, all *nonperiodic* kernel return real space Green's functions,
-all *periodic* kernel reciprocal space ones.
-"""
-
-from math import isnan, pi, sqrt
-
-import numpy as np
-
-###
-
-def square_pressure__nonperiodic(x, y, a=0.5, b=0.5):
-    """
-    Real-space representation of Green's function for the normal displacements
-    of a non-periodic linear elastic half-space with contact modulus 2 and
-    Poisson number 1/2 in response to a uniform pressure applied to a
-    rectangular region.
-    See: K.L. Johnson, Contact Mechanics, p. 54
-
-    Parameters
-    ----------
-    x : array_like
-        x-coordinates.
-    y : array_like
-        y-coordinates.
-    a, b : float
-        Size of rectangle. Rectangle has corners at -a,-b and a,b, i.e.
-        rectangle has extend of 2*a x 2*b. Solution converges to
-        Boussinesq-Cerrutti form for a->0 and b->0.
-
-    Returns
-    -------
-    gf : array
-        Green's function for displacements at coordinates given by x,y as a
-        function of pressure applied to the rectangular region.
-    """
-
-    return ( (x+a)*np.log( ( (y+b)+np.sqrt((y+b)*(y+b)+(x+a)*(x+a)) )/
-                           ( (y-b)+np.sqrt((y-b)*(y-b)+(x+a)*(x+a)) ) )+
-             (y+b)*np.log( ( (x+a)+np.sqrt((y+b)*(y+b)+(x+a)*(x+a)) ) /
-                           ( (x-a)+np.sqrt((y+b)*(y+b)+(x-a)*(x-a)) ) )+
-             (x-a)*np.log( ( (y-b)+np.sqrt((y-b)*(y-b)+(x-a)*(x-a)) ) /
-                           ( (y+b)+np.sqrt((y+b)*(y+b)+(x-a)*(x-a)) ) )+
-             (y-b)*np.log( ( (x-a)+np.sqrt((y-b)*(y-b)+(x-a)*(x-a)) ) /
-                           ( (x+a)+np.sqrt((y-b)*(y-b)+(x+a)*(x+a)) ) ) )/(2*pi);
-
-
-def point_traction__nonperiodic(quantities, x, y, z, G=1.0, poisson=0.5):
-    """
-    Real-space representation of Green's function for the displacement and
-    stress in the bulk of a non-periodic linear elastic half-space in response
-    to a concentrated surface force. This is the Boussinesq-Cerrutti solution.
-    See: K.L. Johnson, Contact Mechanics, p. 51 and p. 69
-    Sign convention is as in Johnson!
-
-    Parameters
-    ----------
-    quantities : str
-        Each character in this string defines a return quantity. They are
-        returned in a tuple of the same order. Possible quantities are
-            'x' : Displacement for a concentrated surface traction in
-                  x-direction.
-            'X' : Stress for a concentrated surface traction in x-direction.
-            'Y' : Stress for a concentrated surface traction in y-direction.
-            'Z' : Stress for a concentrated surface pressure
-                  (i.e. "traction" in z-direction).
-    x : array_like
-        x-coordinates.
-    y : array_like
-        y-coordinates.
-    z : array_like
-        z-coordinates. Into the solid is positive.
-    G : float
-        Shear modulus.
-    poisson : float
-        Poisson number.
-
-    Returns
-    -------
-    sxx : array
-        Green's function xx-component of the stress tensor at coordinates given
-        by x,y,z as a function of pressure applied at point 0,0.
-    syy : array
-        yy-component of the stress tensor.
-    szz : array
-        zz-component of the stress tensor.
-    syz : array
-        yz-component of the stress tensor.
-    sxz : array
-        xz-component of the stress tensor.
-    sxy : array
-        xy-component of the stress tensor.
-    """
-
-    r_sq = x**2 + y**2
-    r_sq = np.array(r_sq, dtype=float)
-    rho = np.sqrt(r_sq + z**2)
-
-    r_sq[r_sq <= 0.0] = 1e-9
-    rho[rho <= 0.0] = 1e-9
-
-    retvals = []
-    for q in quantities:
-        if q == 'x':
-            ux = (1/rho+x**2/rho**3+(1-2*poisson)*(1/(rho+z)-x**2/(rho*(rho+z)**2)))/(4*pi*G)
-            uy = (x*y/rho**3-(1-2*poisson)*x*y/(rho*(rho+z)**2))/(4*pi*G)
-            uz = (x*z/rho**3+(1-2*poisson)*x/(rho*(rho+z)))/(4*pi*G)
-            retvals += [np.array([ux, uy, uz])]
-        elif q == 'y':
-            raise NotImplementedError()
-        elif q == 'z':
-            ux = (x*z/rho**3-(1-2*poisson)*x/(rho*(rho+z)))/(4*pi*G)
-            uy = (y*z/rho**3-(1-2*poisson)*y/(rho*(rho+z)))/(4*pi*G)
-            uz = (z**2/rho**3+2*(1-poisson)/rho)/(4*pi*G)
-            retvals += [np.array([ux, uy, uz])]
-        elif q == 'X':
-            sxx = ( -3*x**3/rho**5 + (1-2*poisson)*(x/rho**3 - 3*x/(rho*(rho+z)**2) + x**3/(rho**3*(rho+z)**2) + 2*x**3/(rho**2*(rho+z)**3)) )/(2*pi)
-            syy = ( -3*x*y**2/rho**5 + (1-2*poisson)*(x/rho**3 - x/(rho*(rho+z)**2) + x*y**2/(rho**3*(rho+z)**2) + 2*x*y**2/(rho**2*(rho+z)**3)) )/(2*pi)
-            szz = ( -3*x*z**2/rho**5 )/(2*pi)
-            sxy = ( -3*x**2*y/rho**5 + (1-2*poisson)*(-y/(rho*(rho+z)**2) + x**2*y/(rho**3*(rho+z)**2) + 2*x**2*y/(rho**2*(rho+z)**3)) )/(2*pi)
-            sxz = ( -3*x**2*z/rho**5 )/(2*pi)
-            syz = ( -3*x*y*z/rho**5 )/(2*pi)
-            retvals += [np.array([sxx, syy, szz, syz, sxz, sxy])]
-        elif q == 'Y':
-            raise NotImplementedError()
-        elif q == 'Z':
-            sxx = ( (1-2*poisson)/r_sq * ((1 - z/rho) * (x**2 - y**2)/r_sq + z*y**2/rho**3) - 3*z*x**2/rho**5 )/(2*pi)
-            syy = ( (1-2*poisson)/r_sq * ((1 - z/rho) * (y**2 - x**2)/r_sq + z*x**2/rho**3) - 3*z*y**2/rho**5 )/(2*pi)
-            szz = -3*z**3/(2*pi*rho**5)
-            # Note: Johnson is lacking a factor of 2 here!!!
-            sxy = ( (1-2*poisson)/r_sq * (2*(1 - z/rho) * x*y/r_sq - x*y*z/rho**3) - 3*x*y*z/rho**5 )/(2*pi)
-            sxz = -3*x*z**2/(2*pi*rho**5)
-            syz = -3*y*z**2/(2*pi*rho**5)
-            retvals += [np.array([sxx, syy, szz, syz, sxz, sxy])]
-        else:
-            raise ValueError("Unknown quantity '{0}' requested.".format(q))
-
-    if len(quantities) == 1:
-        return retvals[0]
-    else:
-        return retvals
-
-
-def real_to_reciprocal_space(nx, ny=None, gf=square_pressure__nonperiodic,
-                             coordinates=False):
-    """
-    Return the reciprocal space representation of a real-space Green's function
-    on an FFT grid.
-    Note: If the Green's function is for the non-periodic (free-bounday)
-    problem, then only a section of 1/2 nx by 1/2 ny of the grid can have
-    non-zero pressure. The other region is a padding region. See R.W. Hockney,
-    Methods in Computational Physics, Vol. 9, pp. 135-211 (1970) for an
-    application of this method to the electrostatic problem.
-
-    Parameters
-    ----------
-    nx : int
-        Number of grid points in x-direction.
-    ny : int
-        Number of grid points in y-direction.
-    gf : function
-        Function returning the real-space Green's function.
-    coordinates : bool
-        If True, return grid coordinates.
-
-    Returns
-    -------
-    G : array
-        Reciprocal space representation of the stiffness on an nx,ny grid.
-    x : array
-        x-coordinates of grid.
-    y : array
-        y-coordinates of grid.
-    """
-
-    if ny is None:
-        nx, ny = nx
-
-    x = np.arange(nx)
-    x = np.where(x <= nx//2, x, x-nx)
-    x.shape = (-1,1)
-    y = np.arange(ny)
-    y = np.where(y <= ny//2, y, y-ny)
-    y.shape = (1,-1)
-
-    G = gf(x, y)
-    if isinstance(G, tuple) or isinstance(G, list):
-        r = [ np.fft.fft2(np.real(_G)) for _G in G ]
-    else:
-        r = np.fft.fft2(np.real(G))
-
-    if coordinates:
-        return r, x, y
-    else:
-        return r
-
-
-def point_displacement__periodic(nx, ny=None, phi0=None, size=None):
-    """
-    Return reciprocal space stiffness coefficients (i.e. inverse of the Green's
-    function) for a periodic system with contact modulus 2 and Poisson number
-    1/2. This gives force as a function of displacement.
-
-    Parameters
-    ----------
-    nx : int
-        Number of grid points in x-direction.
-    ny : int
-        Number of grid points in y-direction.
-    phi0 : float
-        Stiffness at q=0 (gamma point).
-    size : tuple
-        System size. Assumed to be nx,ny if parameter is omitted.
-
-    Returns
-    -------
-    phi : array
-        Reciprocal space representation of the stiffness on an nx,ny grid.
-    """
-
-    if ny is None:
-        nx, ny = nx
-
-    qx = np.arange(nx, dtype=np.float64)
-    qx = np.where(qx <= nx//2, 2*pi*qx/sx, 2*pi*(nx-qx)/sx)
-    qy = np.arange(ny, dtype=np.float64)
-    qy = np.where(qy <= ny//2, 2*pi*qy/sy, 2*pi*(ny-qy)/sy)
-    phi  = np.sqrt( (qx*qx).reshape(-1, 1) + (qy*qy).reshape(1, -1) )
-    if phi0 is None:
-        phi[0, 0] = (phi[1, 0].real + phi[0, 1].real)/2
-    else:
-        phi[0, 0] = phi0
-
-    return phi
-
-
-def min_ccg(h_r, gf_q, u_r=None, pentol=1e-6, maxiter=100000,
-            logger=None):
-    """
-    Use a constrained conjugate gradient optimization to find the equilibrium
-    configuration deflection of an elastic manifold. The conjugate gradient
-    iteration is reset using the steepest descent direction whenever the contact
-    area changes.
-    Method is described in I.A. Polonsky, L.M. Keer, Wear 231, 206 (1999)
-
-    Parameters
-    ----------
-    h_r : array_like
-        Height profile of the rigid counterbody.
-    gf_q : array_like
-        Green's function (in reciprocal space).
-    u_r : array
-        Array used for initial displacements. A new array is created if omitted.
-    pentol : float
-        Maximum penetration of contacting regions required for convergence.
-    maxiter : float
-        Maximum number of iterations.
-
-    Returns
-    -------
-    u : array
-        2d-array of displacements.
-    p : array
-        2d-array of pressure.
-    """
-
-    # Note: Suffix _r deontes real-space _q reciprocal space 2d-arrays
-
-    if logger is not None:
-        logger.pr('maxiter = {0}'.format(maxiter))
-        logger.pr('pentol = {0}'.format(pentol))
-
-    if u_r is None:
-        u_r = np.zeros_like(h_r)
-
-    u_r[:, :] = np.where(u_r > h_r, h_r, u_r)
-
-    # Compute forces
-    p_r = -np.fft.ifft2(np.fft.fft2(u_r)/gf_q).real
-
-    # iteration
-    delta = 0
-    delta_str = 'reset'
-    G_old = 1.0
-    for it in range(1, maxiter+1):
-        # Reset contact area
-        c_r = p_r > 0.0
-
-        # Compute total contact area (area with repulsive force)
-        A = np.sum(c_r)
-
-        # Compute G = sum(g*g) (over contact area only)
-        g_r = h_r-u_r
-        G = np.sum(c_r*g_r*g_r)
-
-        # t = (g + delta*(G/G_old)*t) inside contact area and 0 outside
-        if delta > 0:
-            t_r = c_r*(g_r + delta*(G/G_old)*t_r)
-        else:
-            t_r = c_r*g_r
-
-        # Compute elastic displacement that belong to t_r
-        # (Note: r_r is negative of Polonsky, Kerr's r)
-        r_r = -np.fft.ifft2(gf_q*np.fft.fft2(t_r)).real
-
-        # Note: Sign reversed from Polonsky, Keer because this r_r is negative
-        # of theirs.
-        tau = 0.0
-        if A > 0:
-            # tau = -sum(g*t)/sum(r*t) where sum is only over contact region
-            x = -np.sum(c_r*r_r*t_r)
-            if x > 0.0:
-                tau = np.sum(c_r*g_r*t_r)/x
-            else:
-                G = 0.0
-
-        # Save forces for later
-        fold_r = p_r.copy()
-
-        p_r -= tau*c_r*t_r
-
-        # Find area with negative forces and negative gap
-        # (i.e. penetration of the two surfaces)
-        nc_r = np.logical_and(p_r <= 0.0, g_r < 0.0)
-
-        # Set all negative forces to zero
-        p_r *= p_r > 0.0
-
-        if np.sum(nc_r) > 0:
-            # nc_r contains area that just jumped into contact. Update their
-            # forces.
-            p_r -= tau*nc_r*g_r
-        
-            delta = 0
-            delta_str = 'sd'
-        else:
-            delta = 1
-            delta_str = 'cg'
-
-        # Compute new displacements from updated forces
-        u_r = -np.fft.ifft2(gf_q*np.fft.fft2(p_r)).real
-       
-        # Store G for next step
-        G_old = G
-        
-        # Compute root-mean square penetration, max penetration and max force
-        # difference between the steps
-        if A > 0:
-            rms_pen = sqrt(G/A)
-        else:
-            rms_pen = sqrt(G)
-        max_pen = max(0.0, np.max(c_r*(u_r-h_r)))
-
-        # Elastic energy would be
-        # e_el = -0.5*np.sum(p_r*u_r)
-
-        if rms_pen < pentol and max_pen < pentol:
-            if logger is not None:
-                logger.st(['status', 'it', 'A', 'tau', 'rms_pen', 'max_pen'],
-                          ['CONVERGED', it, A, tau, rms_pen, max_pen],
-                          force_print=True)
-            return u_r, p_r
-
-        if logger is not None:
-            logger.st(['status', 'it', 'A', 'tau', 'rms_pen', 'max_pen'],
-                      [delta_str, it, A, tau, rms_pen, max_pen])
-
-        if isnan(G) or isnan(rms_pen):
-            raise RuntimeError('nan encountered.')
-
-    raise RuntimeError('Maximum number of iterations ({0}) exceeded.' \
-                           .format(maxiter))
-
diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index d3915620..e3802749 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -1,6 +1,6 @@
 #
 # Copyright 2019, 2021 Lars Pastewka (U. Freiburg)
-#           2018-2021 Petr Grigorev (Warwick U.)
+#           2018-2023 Petr Grigorev (Warwick U.)
 #           2020 James Kermode (Warwick U.)
 #           2019 Arnaud Allera (U. Lyon 1)
 #           2019 Wolfram G. Nöhring (U. Freiburg)
@@ -25,16 +25,15 @@
 
 from scipy.optimize import minimize
 
-from ase.lattice.cubic import BodyCenteredCubic
+from ase.lattice.cubic import (BodyCenteredCubic, FaceCenteredCubic,
+                               Diamond, SimpleCubicFactory)
 from ase.constraints import FixAtoms, StrainFilter
 #from ase.optimize import FIRE
 from ase.optimize.precon import PreconLBFGS
 from ase.build import bulk
 from ase.calculators.lammpslib import LAMMPSlib
 from ase.units import GPa  # unit conversion
-from ase.lattice.cubic import SimpleCubicFactory, Diamond
 from ase.io import read
-from ase.geometry import get_distances
 
 from matscipy.neighbours import neighbour_list, mic
 from matscipy.elasticity import fit_elastic_constants
@@ -63,7 +62,7 @@ def make_screw_cyl(alat, C11, C12, C44,
         radius of cylinder of unconstrained atoms around the
         dislocation  in angstrom
     cutoff : float
-        Potential cutoff for marinica potentials for FS cutoff = 4.4
+        Potential cutoff for Marinica potentials for FS cutoff = 4.4
     hard_core : bool
         Description of parameter `hard_core`.
     center : type
@@ -209,17 +208,17 @@ def make_edge_cyl(alat, C11, C12, C44,
     cylinder_r - radius of cylinder of unconstrained atoms around the
                  dislocation  in angstrom
 
-    cutoff - potential cutoff for marinica potentials for FS cutoff = 4.4
+    cutoff - potential cutoff for Marinica potentials for FS cutoff = 4.4
 
     symbol : string
-        Symbol of the element to pass to ase.lattuce.cubic.SimpleCubicFactory
+        Symbol of the element to pass to ase.lattice.cubic.SimpleCubicFactory
         default is "W" for tungsten
     '''
     from atomman import ElasticConstants
     from atomman.defect import Stroh
-    # Create a Stroh ojbect with junk data
+    # Create a Stroh object with junk data
     stroh = Stroh(ElasticConstants(C11=141, C12=110, C44=98),
-                                   np.array([0, 0, 1]))
+                  np.array([0, 0, 1]))
 
     axes = np.array([[1, 1, 1],
                      [1, -1, 0],
@@ -239,7 +238,7 @@ def make_edge_cyl(alat, C11, C12, C44,
 
     bulk = unit_cell.copy()
 
-    # shift to make the zeros of the cell betweem the atomic planes
+    # shift to make the zeros of the cell between the atomic planes
     # and under the midplane on Y axes
     X_midplane_shift = (1.0/3.0)*alat*np.sqrt(3.0)/2.0
     Y_midplane_shift = 0.25*alat*np.sqrt(2.0)
@@ -255,7 +254,7 @@ def make_edge_cyl(alat, C11, C12, C44,
     Lx = int(round(tot_r/(alat*np.sqrt(3.0)/2.0)))
     Ly = int(round(tot_r/(alat*np.sqrt(2.))))
 
-    # factor 2 to make shure odd number of images is translated
+    # factor 2 to make sure odd number of images is translated
     # it is important for the correct centering of the dislocation core
     bulk = bulk * (2*Lx, 2*Ly, 1)
 
@@ -291,9 +290,8 @@ def make_edge_cyl(alat, C11, C12, C44,
 
     x, y, z = bulk.positions.T
     # move lower left segment
-    bulk.positions[(y < 0.0) & (x < X_midplane_shift)] -= \
-        [alat * np.sqrt(3.0) / 2.0, 0.0, 0.0]
-    # make the doslocation extraplane center
+    bulk.positions[(y < 0.0) & (x < X_midplane_shift)] -= [alat * np.sqrt(3.0) / 2.0, 0.0, 0.0]
+    # make the dislocation extra half plane center
     bulk.positions += [(1.0/3.0)*alat*np.sqrt(3.0)/2.0, 0.0, 0.0]
 
     return ED, bulk
@@ -310,7 +308,7 @@ def plot_vitek(dislo, bulk,
     bulk : ase.Atoms
         Corresponding bulk configuration for calculation of displacements.
     alat : float
-        Lattice parameter for calculation of neghbour list cutoff.
+        Lattice parameter for calculation of neighbour list cutoff.
     plot_axes : matplotlib.Axes.axes object
         Existing axes to plot on, allows to pass existing matplotlib axes
         have full control of the graph outside the function.
@@ -337,8 +335,8 @@ def plot_vitek(dislo, bulk,
 
     neighborListCutoff = 0.95 * alat
 
-    # plot window is +-10 angstroms from center in x,y directions, and one Burgers vector
-    # thickness along z direction
+    # plot window is +-10 angstroms from center in x,y directions,
+    # and one Burgers vector thickness along z direction
     x, y, _ = bulk.positions.T
 
     plot_range = np.array([[x.mean() - xyscale, x.mean() + xyscale],
@@ -349,14 +347,15 @@ def plot_vitek(dislo, bulk,
     # distance between atoms on the plot
     plot_scale = 1.885618083
 
-    fig = differential_displacement(base_system, disl_system,
-                                    burgers,
-                                    cutoff=neighborListCutoff,
-                                    xlim=plot_range[0],
-                                    ylim=plot_range[1],
-                                    zlim=plot_range[2],
-                                    matplotlib_axes=plot_axes,
-                                    plot_scale=plot_scale)
+    _ = differential_displacement(base_system, disl_system,
+                                  burgers,
+                                  cutoff=neighborListCutoff,
+                                  xlim=plot_range[0],
+                                  ylim=plot_range[1],
+                                  zlim=plot_range[2],
+                                  matplotlib_axes=plot_axes,
+                                  plot_scale=plot_scale)
+
 
 def show_NEB_configurations(images, bulk, xyscale=7,
                             show=True, core_positions=None):
@@ -453,11 +452,11 @@ def get_elastic_constants(pot_path=None,
     pot_path - path to the potential
 
     symbol : string
-        Symbol of the element to pass to ase.lattuce.cubic.SimpleCubicFactory
+        Symbol of the element to pass to ase.lattice.cubic.SimpleCubicFactory
         default is "W" for tungsten
     """
 
-    unit_cell = bulk(symbol)
+    unit_cell = bulk(symbol, cubic=True)
 
     if (pot_path is not None) and (calculator is None):
         # create lammps calculator with the potential
@@ -468,16 +467,17 @@ def get_elastic_constants(pot_path=None,
 
     unit_cell.calc = calculator
 
-#   simple calculation to get the lattice constant and cohesive energy
-#    alat0 = W.cell[0][1] - W.cell[0][0]
-    sf = StrainFilter(unit_cell)  # or UnitCellFilter(W) -> to minimise wrt pos, cell
-    #opt = FIRE(sf)
-    opt = PreconLBFGS(sf)
+    #   simple calculation to get the lattice constant and cohesive energy
+    #    alat0 = W.cell[0][1] - W.cell[0][0]
+    sf = StrainFilter(unit_cell)
+    # or UnitCellFilter(W)
+    # -> to minimise wrt pos, cell
+    opt = FIRE(sf)
     opt.run(fmax=1e-4)  # max force in eV/A
-    alat = unit_cell.cell[0][1] - unit_cell.cell[0][0]
-#    print("a0 relaxation %.4f --> %.4f" % (a0, a))
-#    e_coh = W.get_potential_energy()
-#    print("Cohesive energy %.4f" % e_coh)
+    alat = unit_cell.cell.lengths()[0]
+    #    print("a0 relaxation %.4f --> %.4f" % (a0, a))
+    #    e_coh = W.get_potential_energy()
+    #    print("Cohesive energy %.4f" % e_coh)
 
     Cij, Cij_err = fit_elastic_constants(unit_cell,
                                          symmetry="cubic",
@@ -593,11 +593,15 @@ def make_barrier_configurations(elastic_param=None,
 
     return disloc_ini, disloc_fin, bulk_ini
 
-def make_screw_cyl_kink(alat, C11, C12, C44,
-                        cylinder_r=40, kink_length=26, kind="double", **kwargs):
-    """Function to create kink configuration based on make_screw_cyl() function.
-        Double kink configuration is in agreement with quadrupoles in terms of formation energy.
-        Single kink configurations provide correct and stable structure, but formation energy is not accessible?
+
+def make_screw_cyl_kink(alat, C11, C12, C44, cylinder_r=40,
+                        kink_length=26, kind="double", **kwargs):
+    """
+    Function to create kink configuration based on make_screw_cyl() function.
+    Double kink configuration is in agreement with
+    quadrupoles in terms of formation energy.
+    Single kink configurations provide correct and stable structure,
+    but formation energy is not accessible?
 
     Parameters
     ----------
@@ -631,8 +635,11 @@ def make_screw_cyl_kink(alat, C11, C12, C44,
     b = np.sqrt(3.0) * alat / 2.0
     cent_x = np.sqrt(6.0) * alat / 3.0
 
-    disloc_ini, disloc_fin, bulk_ini = make_barrier_configurations((alat, C11, C12, C44),
-                                                                    cylinder_r=cylinder_r, **kwargs)
+    (disloc_ini,
+     disloc_fin,
+     bulk_ini) = make_barrier_configurations((alat, C11, C12, C44),
+                                             cylinder_r=cylinder_r,
+                                             **kwargs)
 
     if kind == "double":
 
@@ -690,7 +697,8 @@ def make_screw_cyl_kink(alat, C11, C12, C44,
         cell[2][2] -= 2.0 * b / 3.0
         kink.set_cell(cell, scale_atoms=False)
 
-        # make sure all the atoms are removed and cell is modified for the bulk as well.
+        # make sure all the atoms are removed and cell is modified
+        # for the bulk as well.
         large_bulk.cell[2][0] += cent_x
         large_bulk.cell[2][2] -= 2.0 * b / 3.0
         large_bulk = large_bulk[right_kink_mask]
@@ -727,7 +735,8 @@ def make_screw_cyl_kink(alat, C11, C12, C44,
         cell[2][2] -= 1.0 * b / 3.0
         kink.set_cell(cell, scale_atoms=False)
 
-        # make sure all the atoms are removed and cell is modified for the bulk as well.
+        # make sure all the atoms are removed and cell is modified
+        # for the bulk as well.
         large_bulk.cell[2][0] -= cent_x
         large_bulk.cell[2][2] -= 1.0 * b / 3.0
         large_bulk = large_bulk[left_kink_mask]
@@ -739,8 +748,10 @@ def make_screw_cyl_kink(alat, C11, C12, C44,
 
     return kink, reference_straight_disloc, large_bulk
 
+
 def slice_long_dislo(kink, kink_bulk, b):
-    """Function to slice a long dislocation configuration to perform dislocation structure and core position analysis
+    """Function to slice a long dislocation configuration to perform
+       dislocation structure and core position analysis
 
     Parameters
     ----------
@@ -1284,7 +1295,7 @@ def screw_cyl_octahedral(alat, C11, C12, C44,
 
     # Create a Stroh object with junk data
     stroh = Stroh(ElasticConstants(C11=141, C12=110, C44=98),
-                                   np.array([0, 0, 1]))
+                  np.array([0, 0, 1]))
 
     c = ElasticConstants(C11=C11, C12=C12, C44=C44)
     burgers = alat * np.array([1., 1., 1.])/2.
@@ -1338,8 +1349,6 @@ def dipole_displacement_angle(W_bulk, dislo_coord_left, dislo_coord_right,
         Generates a simple displacement field for two dislocations in a dipole
         configuration uding simple Voltera solution as u = b/2 * angle
     """
-
-
     burgers = W_bulk.cell[2][2]
 
     shifted_positions = W_bulk.positions + shift - dislo_coord_left
@@ -1381,7 +1390,7 @@ def get_u_img(W_bulk,
                 u_img += dipole_displacement_angle(W_bulk,
                                                    dislo_coord_left + shift,
                                                    dislo_coord_right + shift,
-                                                   shift=n2_shift*C2_quadrupole + n1_shift*C1_quadrupole)
+                                                   shift=n2_shift * C2_quadrupole + n1_shift * C1_quadrupole)
 
     return u_img
 
@@ -1521,8 +1530,8 @@ def make_screw_quadrupole(alat,
     """
 
     unit_cell = BodyCenteredCubic(directions=[[1, -2, 1],
-                                      [2, -1, -1],
-                                      [1, 1, 1]],
+                                              [2, -1, -1],
+                                              [1, 1, 1]],
                                   symbol=symbol,
                                   pbc=(True, True, True),
                                   latticeconstant=alat, debug=0)
@@ -1693,10 +1702,12 @@ def make_screw_quadrupole(alat,
     return disloc_quadrupole, bulk, dislo_coord_left, dislo_coord_right
 
 
-def make_screw_quadrupole_kink(alat, kind="double", n1u=5, kink_length=20, symbol="W"):
+def make_screw_quadrupole_kink(alat, kind="double",
+                               n1u=5, kink_length=20, symbol="W"):
     """Generates kink configuration using make_screw_quadrupole() function
        works for BCC structure.
-       The method is based on paper https://doi.org/10.1016/j.jnucmat.2008.12.053
+       The method is based on paper
+       https://doi.org/10.1016/j.jnucmat.2008.12.053
 
     Parameters
     ----------
@@ -1705,7 +1716,8 @@ def make_screw_quadrupole_kink(alat, kind="double", n1u=5, kink_length=20, symbo
     kind : string
         kind of the kink: right, left or double
     n1u : int
-        Number of lattice vectors for the quadrupole cell (make_screw_quadrupole() function)
+        Number of lattice vectors for the quadrupole cell
+        (make_screw_quadrupole() function)
     kink_length : int
         Length of the cell per kink along b in unit of b, must be even.
     symbol : string
@@ -1727,15 +1739,17 @@ def make_screw_quadrupole_kink(alat, kind="double", n1u=5, kink_length=20, symbo
     b = np.sqrt(3.0) * alat / 2.0
     cent_x = np.sqrt(6.0) * alat / 3.0
 
-    ini_disloc_quadrupole, W_bulk, _, _ = make_screw_quadrupole(alat, n1u=n1u,
-                                                                   left_shift=0.0,
-                                                                   right_shift=0.0,
-                                                                   symbol=symbol)
+    (ini_disloc_quadrupole,
+     W_bulk, _, _) = make_screw_quadrupole(alat, n1u=n1u,
+                                           left_shift=0.0,
+                                           right_shift=0.0,
+                                           symbol=symbol)
 
-    fin_disloc_quadrupole, W_bulk, _, _ = make_screw_quadrupole(alat, n1u=n1u,
-                                                                   left_shift=1.0,
-                                                                   right_shift=1.0,
-                                                                   symbol=symbol)
+    (fin_disloc_quadrupole,
+     W_bulk, _, _) = make_screw_quadrupole(alat, n1u=n1u,
+                                           left_shift=1.0,
+                                           right_shift=1.0,
+                                           symbol=symbol)
 
     reference_straight_disloc = ini_disloc_quadrupole * [1, 1, kink_length]
     large_bulk = W_bulk * [1, 1, kink_length]
@@ -1756,7 +1770,8 @@ def make_screw_quadrupole_kink(alat, kind="double", n1u=5, kink_length=20, symbo
         kink.cell[2][2] += upper_kink.cell[2][2]
         kink.extend(upper_kink)
 
-        # left kink is created the kink vector is in negative x direction assuming (x, y, z) is right group of vectors
+        # left kink is created the kink vector is in negative x direction
+        # assuming (x, y, z) is right group of vectors
         kink = kink[left_kink_mask]
         kink.cell[2][0] -= cent_x
         kink.cell[2][2] -= 1.0 * b / 3.0
@@ -1798,7 +1813,8 @@ def make_screw_quadrupole_kink(alat, kind="double", n1u=5, kink_length=20, symbo
 
         kink.cell[2][2] += upper_kink.cell[2][2]
 
-        large_bulk = W_bulk * [1, 1, 2 * kink_length]  # double kink is double length
+        # double kink is double length
+        large_bulk = W_bulk * [1, 1, 2 * kink_length]
 
     else:
         raise ValueError('Kind must be "right", "left" or "double"')
@@ -1904,8 +1920,8 @@ def make_edge_cyl_001_100(a0, C11, C12, C44,
         print('disloc SCF', i, '|d1-d2|_inf =', res)
         i += 1
         if i > 10:
-            raise RuntimeError('Self-consistency did \
-                                not converge in 10 cycles')
+            raise RuntimeError('Self-consistency did ' +
+                               'not converge in 10 cycles')
     disp = disp2
 
     x, y, z = disloc.positions.T
@@ -2117,7 +2133,8 @@ def get_centering_mask(atoms, radius,
 
 def check_duplicates(atoms, distance=0.1):
     """
-    Returns a mask of atoms that have at least one other atom closer than distance
+    Returns a mask of atoms that have at least
+    one other atom closer than distance
     """
     mask = atoms.get_all_distances() < distance
     duplicates = np.full_like(atoms, False)
@@ -2195,13 +2212,11 @@ def init_stroh(self):
         c = ElasticConstants(C11=self.C11, C12=self.C12, C44=self.C44)
         self.stroh = Stroh(c, burgers=self.burgers, axes=self.axes)
 
-
     def set_burgers(self, burgers):
         self.burgers = burgers
         if self.stroh is None:
             self.init_stroh()
 
-
     def plot_unit_cell(self, ms=250, ax=None):
         import matplotlib.pyplot as plt
 
@@ -2222,7 +2237,7 @@ def plot_unit_cell(self, ms=250, ax=None):
                 [0.0, y0, y0, 0.0, 0.0], color="black", zorder=0)
 
         bulk_atoms = ax.scatter([], [], color="w", edgecolor="k",
-                             label="lattice atoms")
+                                label="lattice atoms")
         core1 = ax.scatter([], [], marker="x", label="initial core position",
                            c="r")
         core2 = ax.scatter([], [], marker="x", label="glide core position",
@@ -2232,8 +2247,7 @@ def plot_unit_cell(self, ms=250, ax=None):
         ax.set_xlabel(r"$\AA$")
         ax.set_ylabel(r"$\AA$")
 
-
-    def displacements(self, bulk_positions, center, self_consistent=True, 
+    def displacements(self, bulk_positions, center, self_consistent=True,
                       tol=1e-6, max_iter=100, verbose=True):
         if self.stroh is None:
             self.init_stroh()
@@ -2253,7 +2267,8 @@ def displacements(self, bulk_positions, center, self_consistent=True,
                 print('disloc SCF', i, '|d1-d2|_inf =', res)
             i += 1
             if i > max_iter:
-                raise RuntimeError(f'Self-consistency did not converge in {max_iter} cycles')
+                raise RuntimeError('Self-consistency' +
+                                   f'did not converge in {max_iter} cycles')
         return disp2
 
     def build_cylinder(self, radius,
@@ -2273,7 +2288,7 @@ def build_cylinder(self, radius,
         repeat_extension = np.floor(2.0 * extension /
                                     np.diag(self.unit_cell.cell)).astype(int)
         repeat_core_position = np.floor(2.0 * core_position /
-                                    np.diag(self.unit_cell.cell)).astype(int)
+                                        np.diag(self.unit_cell.cell)).astype(int)
 
         extra_repeat = np.stack((repeat_core_position,
                                  repeat_extension)).max(axis=0)
@@ -2296,7 +2311,8 @@ def build_cylinder(self, radius,
         center = np.diag(bulk.cell) / 2
         shifted_center = center + core_position
 
-        cylinder_mask = get_centering_mask(bulk, radius, core_position, extension)
+        cylinder_mask = get_centering_mask(bulk, radius,
+                                           core_position, extension)
 
         # add square borders for the case of large extension or core position
         x, y, _ = bulk.positions.T
@@ -2408,7 +2424,7 @@ def build_impurity_cylinder(self, disloc, impurity, radius,
         repeat_extension = np.floor(extension /
                                     np.diag(self.unit_cell.cell)).astype(int)
         repeat_core_position = np.floor(core_position /
-                                    np.diag(self.unit_cell.cell)).astype(int)
+                                        np.diag(self.unit_cell.cell)).astype(int)
 
         extra_repeat = np.stack((repeat_core_position,
                                  repeat_extension)).max(axis=0)
@@ -2461,8 +2477,8 @@ def build_impurity_cylinder(self, disloc, impurity, radius,
         non_core_mask = np.logical_not(core_mask)
 
         displacemets = self.displacements(impurities_bulk.positions[non_core_mask],
-                                                          shifted_center,
-                                                          self_consistent=self_consistent)
+                                          shifted_center,
+                                          self_consistent=self_consistent)
 
         impurities_disloc.positions[non_core_mask] += displacemets
 
@@ -2489,6 +2505,7 @@ def build_impurity_cylinder(self, disloc, impurity, radius,
 
         return impurities_disloc
 
+
 class BCCScrew111Dislocation(CubicCrystalDislocation):
     def __init__(self, alat, C11, C12, C44, symbol='W'):
         axes = np.array([[1, 1, -2],
@@ -2507,7 +2524,7 @@ def __init__(self, alat, C11, C12, C44, symbol='W'):
                          axes, burgers, unit_cell_core_position, parity,
                          glide_distance)
 
-        
+
 class BCCEdge111Dislocation(CubicCrystalDislocation):
     def __init__(self, alat, C11, C12, C44, symbol='W'):
         axes = np.array([[1, 1, 1],
@@ -2515,7 +2532,7 @@ def __init__(self, alat, C11, C12, C44, symbol='W'):
                          [1, 1, -2]])
         burgers = alat * np.array([1, 1, 1]) / 2.0
         unit_cell_core_position = alat * np.array([(1.0/3.0) * np.sqrt(3.0)/2.0,
-                                  0.25 * np.sqrt(2.0), 0])
+                                                  0.25 * np.sqrt(2.0), 0])
         parity = [0, 0]
         unit_cell = BodyCenteredCubic(directions=axes.tolist(),
                                       size=(1, 1, 1), symbol=symbol,
@@ -2668,20 +2685,82 @@ def __init__(self, alat, C11, C12, C44, symbol='C'):
 
 
 class CubicCrystalDissociatedDislocation(CubicCrystalDislocation):
-    """
-        This class represents a dissociated dislocation in a cubic crystal
-        with b = b_left + b_right.
-        'left_dislocation' and 'right_dislocations' are expected
-        to be instances of classes derived from CubicCrystalDislocation class.
-    """
-    def __init__(self, left_dislocation, right_dislocation,
-                 *args, **kwargs):
+    def __init__(self, left_dislocation, right_dislocation, burgers):
+        """This class represents a dissociated dislocation in a cubic crystal
+        with burgers vercor b = b_left + b_right.
+
+        Args:
+            left_dislocation (CubicCrystalDislocation): dislocation with b_left
+            right_dislocation (CubicCrystalDislocation): dislocation with b_right
+            burgers (ndarray of 3 floats): resulting burgers vector
+
+        Raises:
+            ValueError: If resulting burgers vector
+                        burgers is not a sum of burgers vectors of
+                        left and right dislocations.
+            ValueError: If one of the properties of
+                        left and righ dislocations are not the same.
+        """
+        try:
+            np.testing.assert_almost_equal(left_dislocation.burgers +
+                                           right_dislocation.burgers,
+                                           burgers)
+        except AssertionError as error:
+            print(error)
+            raise ValueError("Burgers vectors of left and right disloctions" +
+                             "do not add up to the desired vector")
+
+        # checking that parameters of
+        # left and right dislocations are the same
+        try:
+            assert left_dislocation.alat == right_dislocation.alat
+            assert left_dislocation.C11 == right_dislocation.C11
+            assert left_dislocation.C12 == right_dislocation.C12
+            assert left_dislocation.C44 == right_dislocation.C44
+
+            np.testing.assert_equal(left_dislocation.unit_cell.get_chemical_symbols(),
+                                    right_dislocation.unit_cell.get_chemical_symbols())
+            np.testing.assert_equal(left_dislocation.unit_cell.cell.cellpar(),
+                                    right_dislocation.unit_cell.cell.cellpar())
+            np.testing.assert_equal(left_dislocation.unit_cell.positions,
+                                    right_dislocation.unit_cell.positions)
+
+            np.testing.assert_equal(left_dislocation.axes,
+                                    right_dislocation.axes)
+
+            np.testing.assert_equal(left_dislocation.unit_cell_core_position,
+                                    right_dislocation.unit_cell_core_position)
+
+            np.testing.assert_equal(left_dislocation.parity,
+                                    right_dislocation.parity)
+
+            np.testing.assert_equal(left_dislocation.glide_distance,
+                                    right_dislocation.glide_distance)
+
+            assert left_dislocation.n_planes == right_dislocation.n_planes
+            assert left_dislocation.self_consistent == right_dislocation.self_consistent
+
+        except AssertionError as error:
+            print("Parameters of left and right partials are not the same!")
+            print(error)
+            raise ValueError("Parameters of left and right" +
+                             "partials must be the same")
 
         self.left_dislocation = left_dislocation
         self.right_dislocation = right_dislocation
 
-        super().__init__(*args, **kwargs)
-
+        super().__init__(left_dislocation.unit_cell,
+                         left_dislocation.alat,
+                         left_dislocation.C11,
+                         left_dislocation.C12,
+                         left_dislocation.C44,
+                         left_dislocation.axes,
+                         burgers,
+                         unit_cell_core_position=left_dislocation.unit_cell_core_position,
+                         parity=left_dislocation.parity,
+                         glide_distance=right_dislocation.glide_distance,
+                         n_planes=left_dislocation.n_planes,
+                         self_consistent=left_dislocation.self_consistent)
 
     def build_cylinder(self, radius, partial_distance=0,
                        core_position=np.array([0., 0., 0.]),
@@ -2709,104 +2788,186 @@ def build_cylinder(self, radius, partial_distance=0,
             [self.glide_distance * partial_distance, 0.0, 0.0])
 
         bulk, disloc = self.left_dislocation.build_cylinder(radius,
-                                                  extension=extension + partial_distance_Angstrom,
-                                                  core_position=core_position,
-                                                  fix_width=fix_width,
-                                                  self_consistent=self_consistent)
+                                                            extension=extension + partial_distance_Angstrom,
+                                                            core_position=core_position,
+                                                            fix_width=fix_width,
+                                                            self_consistent=self_consistent)
 
         _, disloc_right = self.right_dislocation.build_cylinder(radius,
-                                                      core_position=core_position + partial_distance_Angstrom,
-                                                      extension=extension,
-                                                      fix_width=fix_width,
-                                                      self_consistent=self_consistent)
+                                                                core_position=core_position + partial_distance_Angstrom,
+                                                                extension=extension,
+                                                                fix_width=fix_width,
+                                                                self_consistent=self_consistent)
 
         u_right = disloc_right.positions - bulk.positions
         disloc.positions += u_right
 
         return bulk, disloc
 
+    def displacements(self, bulk_positions, center,
+                      partial_distance=0, **kwargs):
+        """Overloaded function to provide correct displacements
+        for the dissociated dislocation.
+        Partial distance is provided as an integer to define number
+        of glide distances between two partials.
+        """
+
+        partial_distance_Angstrom = np.array(
+            [self.glide_distance * partial_distance, 0.0, 0.0])
+
+        left_u = self.left_dislocation.displacements(bulk_positions, center,
+                                                     **kwargs)
+        right_u = self.right_dislocation.displacements(bulk_positions,
+                                                       center + partial_distance_Angstrom,
+                                                       **kwargs)
+
+        return left_u + right_u
+
 
 class DiamondGlideScrew(CubicCrystalDissociatedDislocation):
     def __init__(self, alat, C11, C12, C44, symbol='C'):
 
-        axes = np.array([[1, 1, -2],
-                        [1, 1, 1],
-                        [1, -1, 0]])
-
         # aiming for the resulting burgers vector
         burgers = alat * np.array([1, -1, 0]) / 2.
 
-        disloCenterX = 0.5 * (alat * np.linalg.norm(axes[0])) / 6.0
-        # 1/4 + 1/2 (1/3 - 1/4) - to be in the middle of the glide set
-        disloCenterY = 7.0 * (alat * np.linalg.norm(axes[1])) / 24.0
+        # 30 degree
+        burgers_left = alat * np.array([2., -1., -1.]) / 6.
 
-        unit_cell_core_position = np.array([disloCenterX,
-                                            disloCenterY, 0])
+        left30 = DiamondGlide30degreePartial(alat, C11, C12, C44,
+                                             symbol=symbol)
+        left30.set_burgers(burgers_left)
+        # another 30 degree
+        # burgers_right = alat * np.array([1, -2, 1.]) / 6. - default value
+        right30 = DiamondGlide30degreePartial(alat, C11, C12, C44,
+                                              symbol=symbol)
 
-        parity = [0, 0]
+        super().__init__(left30, right30, burgers)
 
-        unit_cell = Diamond(symbol, directions=axes.tolist(),
-                            pbc=(False, False, True),
-                            latticeconstant=alat)
 
-        glide_distance = alat * np.linalg.norm(axes[0]) / 4.0
+class DiamondGlide60Degree(CubicCrystalDissociatedDislocation):
+    def __init__(self, alat, C11, C12, C44, symbol='C'):
 
-        n_planes = 2
+        # aiming for the resulting burgers vector
+        burgers = alat * np.array([1, 0, -1]) / 2.
 
         # 30 degree
         burgers_left = alat * np.array([2., -1., -1.]) / 6.
-        left30 = DiamondGlide30degreePartial(alat, C11, C12, C44)
+        left30 = DiamondGlide30degreePartial(alat, C11, C12, C44,
+                                             symbol=symbol)
         left30.set_burgers(burgers_left)
-        # another 30 degree
-        burgers_right = alat * np.array([1, -2, 1.]) / 6.
-        right30 = DiamondGlide30degreePartial(alat, C11, C12, C44)
-        self_consistent = False
-        super().__init__(left30, right30, unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance, n_planes=n_planes,
-                         self_consistent=self_consistent)
+        # 90 degree
+        # burgers_right = alat * np.array([1, 1, -2.]) / 6. - default value
+        right90 = DiamondGlide90degreePartial(alat, C11, C12, C44,
+                                              symbol=symbol)
 
+        super().__init__(left30, right90, burgers)
 
-class DiamondGlide60Degree(CubicCrystalDissociatedDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='C'):
+
+class FCCScrewShockleyPartial(CubicCrystalDislocation):
+    def __init__(self, alat, C11, C12, C44, symbol='Fe'):
         axes = np.array([[1, 1, -2],
                          [1, 1, 1],
                          [1, -1, 0]])
 
-        # aiming for the resulting burgers vector
-        burgers = alat * np.array([1, 0, -1]) / 2.
+        burgers = alat * np.array([1, -2, 1.]) / 6.
 
-        disloCenterX = 0.5 * (alat * np.linalg.norm(axes[0])) / 6.0
-        # 1/4 + 1/2 (1/3 - 1/4) - to be in the middle of the glide set
-        disloCenterY = 7.0 * (alat * np.linalg.norm(axes[1])) / 24.0
+        parity = [0, 0]
 
-        unit_cell_core_position = np.array([disloCenterX,
-                                            disloCenterY, 0])
+        unit_cell = FaceCenteredCubic(symbol, directions=axes.tolist(),
+                                      pbc=(False, False, True),
+                                      latticeconstant=alat)
 
-        parity = [0, 0]
+        # put the dislocation in the centroid of the first triangle in xy plane
+        # get sorting of the atoms in the distance of atoms in xy plane
+        sorted_indices = np.argsort(np.linalg.norm(unit_cell.positions[:, :2],
+                                                   axis=1))
+        # centroid coordinates are simply
+        # mean of x y coordinates of the trianlge
+        (disloCenterX,
+         disloCenterY,
+         _) = unit_cell.positions[sorted_indices[:3]].mean(axis=0)
 
-        unit_cell = Diamond(symbol, directions=axes.tolist(),
-                            pbc=(False, False, True),
-                            latticeconstant=alat)
+        unit_cell_core_position = np.array([disloCenterX,
+                                            disloCenterY, 0])
 
         glide_distance = alat * np.linalg.norm(axes[0]) / 4.0
 
         n_planes = 2
+        self_consistent = True
+        super().__init__(unit_cell, alat, C11, C12, C44,
+                         axes, burgers, unit_cell_core_position, parity,
+                         glide_distance, n_planes=n_planes,
+                         self_consistent=self_consistent)
 
-        # 30 degree
+
+class FCCScrew110Dislocation(CubicCrystalDissociatedDislocation):
+    def __init__(self, alat, C11, C12, C44, symbol='Fe'):
+
+        # aiming for the resulting burgers vector
+        burgers = alat * np.array([1, -1, 0]) / 2.
+
+        # Shockley partial
         burgers_left = alat * np.array([2., -1., -1.]) / 6.
-        left30 = DiamondGlide30degreePartial(alat, C11, C12, C44)
-        left30.set_burgers(burgers_left)
-        # 90 degree
-        burgers_right = alat * np.array([1, 1, -2.]) / 6.
-        right90 = DiamondGlide90degreePartial(alat, C11, C12, C44)
-        self_consistent = False
-        super().__init__(left30, right90, unit_cell, alat, C11, C12, C44,
+        left_shockley = FCCScrewShockleyPartial(alat, C11, C12, C44,
+                                                symbol=symbol)
+        left_shockley.set_burgers(burgers_left)
+        # another Shockley partial
+        # burgers_right = alat * np.array([1, -2, 1.]) / 6. - default value
+        right_shockley = FCCScrewShockleyPartial(alat, C11, C12, C44,
+                                                 symbol=symbol)
+        super().__init__(left_shockley, right_shockley, burgers)
+
+
+class FCCEdgeShockleyPartial(CubicCrystalDislocation):
+    def __init__(self, alat, C11, C12, C44, symbol='Fe'):
+        axes = np.array([[1, -1, 0],
+                         [1, 1, 1],
+                         [-1, -1, 2]])
+
+        burgers = alat * np.array([1, -2, 1.]) / 6.
+
+        parity = [0, 0]
+
+        unit_cell = FaceCenteredCubic(symbol, directions=axes.tolist(),
+                                      pbc=(False, False, True),
+                                      latticeconstant=alat)
+
+        disloCenterX = 0.0
+        # middle between two (111) planes
+        disloCenterY = unit_cell.cell[1][1] / 6.0
+
+        unit_cell_core_position = np.array([disloCenterX,
+                                            disloCenterY, 0])
+
+        glide_distance = alat * np.linalg.norm(axes[0]) / 4.0
+
+        n_planes = 6
+        self_consistent = True
+        super().__init__(unit_cell, alat, C11, C12, C44,
                          axes, burgers, unit_cell_core_position, parity,
                          glide_distance, n_planes=n_planes,
                          self_consistent=self_consistent)
 
 
+class FCCEdge110Dislocation(CubicCrystalDissociatedDislocation):
+    def __init__(self, alat, C11, C12, C44, symbol='Fe'):
+
+        # aiming for the resulting burgers vector
+        burgers = alat * np.array([1, -1, 0]) / 2.
+
+        # Shockley partial
+        burgers_left = alat * np.array([2., -1., -1.]) / 6.
+        left_shockley = FCCEdgeShockleyPartial(alat, C11, C12, C44,
+                                               symbol=symbol)
+        left_shockley.set_burgers(burgers_left)
+        # another Shockley partial
+        # burgers_right = alat * np.array([1, -2, 1.]) / 6. - default value
+        right_shockley = FCCEdgeShockleyPartial(alat, C11, C12, C44,
+                                                symbol=symbol)
+
+        super().__init__(left_shockley, right_shockley, burgers)
+
+
 class FixedLineAtoms:
     """Constrain atoms to move along a given direction only."""
     def __init__(self, a, direction):
@@ -2912,7 +3073,8 @@ def gamma_line(unit_cell, calc=None, shift_dir=0, surface=2,
 
         if image.get_forces().max() < fmax:
             raise RuntimeError(
-                "Initial max force is smaller than fmax! Check surface direction")
+                "Initial max force is smaller than fmax!" +
+                "Check surface direction")
 
         if relax:
             opt = LBFGSLineSearch(image)
diff --git a/matscipy/elasticity.py b/matscipy/elasticity.py
index ada2a2da..426045d2 100644
--- a/matscipy/elasticity.py
+++ b/matscipy/elasticity.py
@@ -1096,11 +1096,11 @@ def elastic_moduli(C, l=np.array([1, 0, 0]), R=None, tol=1e-6):
                               np.eye(3, dtype=float)) > tol):
             raise RuntimeError('Matrix *R* does not describe a rotation.')
     else:
-        u_a = np.array([1, 0, 0])
-        u_b = l/norm(l)  # Normalise directions
+        u_a = l/norm(l)  # Normalise directions
+        u_b = np.array([1, 0, 0])
         R = np.eye(3)
 
-        if not np.allclose(l, u_a, rtol=tol, atol=tol):
+        if np.abs(np.dot(u_a, u_b)) < 1.0 - tol:
             u_v = np.cross(u_a, u_b)
             u_v_mat = np.array([[ 0,      -u_v[2], u_v[1]],
                                 [ u_v[2], 0,      -u_v[0]],
diff --git a/matscipy/fracture_mechanics/clusters.py b/matscipy/fracture_mechanics/clusters.py
index a603758d..5af1ad40 100644
--- a/matscipy/fracture_mechanics/clusters.py
+++ b/matscipy/fracture_mechanics/clusters.py
@@ -59,7 +59,7 @@ def set_groups(a, n, skin_x, skin_y, central_x=-1./2, central_y=-1./2,
 
     a.set_array('groups', g)
 
-def set_regions(cryst, r_I, cutoff, r_III):
+def set_regions(cryst, r_I, cutoff, r_III, extended_far_field=False,extended_region_I=False,exclude_surface=False):
     sx, sy, sz = cryst.cell.diagonal()
     x, y = cryst.positions[:, 0], cryst.positions[:, 1]
     cx, cy = sx/2, sy/2
@@ -67,10 +67,12 @@ def set_regions(cryst, r_I, cutoff, r_III):
 
     # Regions I and III defined by radial distance from center
     regionI = r < r_I
-    regionIII = (r >= r_I) & (r < r_III)
+    regionII = (r >= r_I) & (r < (r_I + cutoff))
+    regionIII = (r >= r_I+cutoff) & (r < r_III)
+    #regionIII = (r >= r_I) & (r < r_III)
     regionIV = (r >= r_III) & (r < (r_III + cutoff))
 
-    regionII = np.zeros(len(cryst), bool)
+    """    regionII = np.zeros(len(cryst), bool)
     i, j = neighbour_list('ij', cryst, cutoff)
     for idx in regionI.nonzero()[0]:
         neighbs = j[i == idx] 
@@ -79,8 +81,35 @@ def set_regions(cryst, r_I, cutoff, r_III):
         # print(f'adding {mask.sum()} neigbours of atom {idx} to regionII')
         mask[regionI] = False # exclude those in region I already
         regionII[mask] = True # add to region I
-        regionIII[mask] = False # remove from region III
-
+        regionIII[mask] = False # remove from region III"""
+
+    if exclude_surface or extended_region_I:
+        #build a mask of the material surface based on the following criteria:
+        # - the atoms x coordinate should lie between rI and 
+        #  rI + 2*cutoff (if non extended) and rIII + cutoff if extended.
+        # - the atoms y coordinate should lie between +- cutoff
+        if extended_far_field:
+            x_criteria = ((x-cx)<-(r_I-cutoff))&\
+                    (((x-cx)>-((r_III)+(cutoff))))
+        else:
+            x_criteria = ((x-cx)<-((r_I-cutoff)))&\
+                (((x-cx)>-((r_I)+((2*cutoff)))))
+        y_criteria = ((y-cy)>-(cutoff))&\
+                ((y-cy)<(cutoff))
+        surface_mask = np.logical_and(x_criteria,y_criteria)
+    if extended_region_I:
+        #re-do all the region logic so they are correct
+        #(and followed by an xor to exclude new region I atoms.)
+        regionI = regionI|surface_mask
+        regionII = regionII^np.logical_and(regionII,surface_mask)
+        regionIII = regionIII^np.logical_and(regionIII,surface_mask)
+        regionIV = regionIV^np.logical_and(regionIV,surface_mask)
+    
+    elif exclude_surface:
+        regionII = regionII^np.logical_and(regionII,surface_mask)
+        regionIII = regionIII^np.logical_and(regionIII,surface_mask)
+        regionIV = regionIV|(surface_mask^np.logical_and(surface_mask,regionI)) #add the surface mask to region IV
+    
     cryst.new_array('region', np.zeros(len(cryst), dtype=int))
     region = cryst.arrays['region']
     region[regionI]  = 1
diff --git a/matscipy/fracture_mechanics/crack.py b/matscipy/fracture_mechanics/crack.py
index 05d3b6e1..56106f35 100644
--- a/matscipy/fracture_mechanics/crack.py
+++ b/matscipy/fracture_mechanics/crack.py
@@ -28,7 +28,6 @@
 from numpy.linalg import inv
 
 try:
-    from scipy.optimize.nonlin import NoConvergence
     from scipy.optimize import brentq, leastsq, minimize, root
     from scipy.sparse import csc_matrix, spdiags
     from scipy.sparse.linalg import spsolve, spilu, LinearOperator
@@ -47,6 +46,12 @@
 
 from matscipy.neighbours import neighbour_list
 
+import matplotlib.pyplot as plt
+from ase.optimize.precon import Exp, PreconLBFGS
+from ase.optimize import LBFGS
+from ase.constraints import FixAtoms
+from ase import Atom
+import ase.io
 
 ###
 
@@ -252,6 +257,13 @@ def converged(F, X):
 
 
 ###
+#wrapper to count function calls for writing
+def counted(f):
+    def wrapped(*args):
+        wrapped.calls += 1
+        return f(*args,str(wrapped.calls))
+    wrapped.calls = 0
+    return wrapped
 
 class RectilinearAnisotropicCrack:
     """
@@ -557,7 +569,8 @@ class CubicCrystalCrack:
     """
 
     def __init__(self, crack_surface, crack_front, C11=None, C12=None,
-                 C44=None, stress_state=PLANE_STRAIN, C=None, Crot=None):
+                 C44=None, stress_state=PLANE_STRAIN, C=None, Crot=None,
+                 cb=None):
         """
         Initialize a crack in a cubic crystal with elastic constants C11, C12
         and C44 (or optionally a full 6x6 elastic constant matrix C).
@@ -603,6 +616,12 @@ def __init__(self, crack_surface, crack_front, C11=None, C12=None,
                                         S6[0, 5], S6[1, 5], S6[2, 5],
                                         S6[5, 5])
 
+        self.third_dir = third_dir
+        self.crack_surface = crack_surface
+        self.crack_front = crack_front
+        self.RotationMatrix = A
+        self.cb = cb
+
 
     def k1g(self, surface_energy):
         """
@@ -913,7 +932,7 @@ class SinclairCrack:
     def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
                  variable_alpha=True, variable_k=False,
                  alpha_scale=None, k_scale=None,
-                 extended_far_field=False):
+                 extended_far_field=False,rI=0.0,rIII=0.0,cutoff=0.0):
         """
 
         Parameters
@@ -953,7 +972,10 @@ def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
         self.u = np.zeros((self.N1, 3))
         self.alpha = alpha
         self.k = k
-
+        self.rI = rI
+        self.rIII = rIII
+        self.cutoff = cutoff
+        self.shiftmask = None
         # check atoms are sorted by distance from centre so we can use N1,N2,N3
         tip_x = cryst.cell.diagonal()[0] / 2.0
         tip_y = cryst.cell.diagonal()[1] / 2.0
@@ -975,6 +997,8 @@ def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
 
         self.precon = None
         self.precon_count = 0
+        self.f_alpha_correction = 0
+        self.alpha0 = alpha
 
     def pack(self, u, alpha, k):
         dofs = list(u.reshape(-1))
@@ -1008,6 +1032,10 @@ def get_dofs(self):
 
     def set_dofs(self, x):
         self.u[:], self.alpha, self.k = self.unpack(x, reshape=True)
+        print('alpha',self.alpha)
+        if self.variable_k:
+            print('k',self.k)
+        #self.k = self.k-(np.abs(self.alpha-self.alpha0)*self.k1g)
         self.update_atoms()
 
     def __len__(self):
@@ -1034,6 +1062,98 @@ def u_cle(self, alpha=None):
         u = np.c_[ux, uy, np.zeros_like(ux)] # convert to 3D field
         return u
 
+
+    #-----------------delete this part from main branch-------------------
+
+    def relax_crack_surface(self):
+        """Relaxes the free surface of the crack in regions II-IV. This is crucial
+        as the unrelaxed surface leads to huge f_{alpha} components which leads to a 
+        strong finite size effect. """
+        self.atoms.set_pbc([False, False, True])
+        self.atoms.calc = self.calc
+
+        self.atoms.info['k'] = self.k
+        self.atoms.info['alpha'] = self.alpha
+        
+        # x = x_cryst + K * u_cle + u
+        self.atoms.positions[:, :] = self.cryst.positions
+        #self.u_cle()[:,1].tofile('uCLEsolution.csv',sep=',')
+        u_cle_solution = self.k * self.u_cle()
+        self.atoms.positions[:, :] += u_cle_solution
+        self.atoms.positions[self.regionI, :] += self.u
+
+        # add vacuum
+        self.atoms.cell = self.cryst.cell
+        self.atoms.cell[0, 0] += self.vacuum
+        self.atoms.cell[1, 1] += self.vacuum
+
+        #clone atoms for relaxing
+        atoms = self.atoms.copy()
+        i = neighbour_list('i', atoms, 12.0) #large cutoff to get the first 2 layers
+        coord = np.bincount(i)
+        np.savetxt('coords.txt',coord)
+        #build a mask for relaxing based on the following criteria:
+        # - the coordination of the atoms must be lower than 0.8*the max coordination
+        # - the atoms x coordinate should lie between rI - cutoff and 
+        #  rI + 2*cutoff (if non extended) and rIII + cutoff/2 if extended.
+        # - the atoms y coordinate should lie between the positive and negative
+        # maximum y Ucle displacement + cutoff
+        coordination_criteria = (coord<(0.8*np.max(coord)))
+        #get the coords of the centre of the cell.
+        mid_cell_x = self.cryst.cell.diagonal()[0] / 2.0
+        mid_cell_y = self.cryst.cell.diagonal()[1] / 2.0
+        if self.extended_far_field:
+            x_criteria = ((atoms.positions[:,0]-mid_cell_x)<-((self.rI)-(2*self.cutoff)))&\
+                (((atoms.positions[:,0]-mid_cell_x)>-((self.rIII)+(self.cutoff/2))))
+        else:
+            x_criteria = ((atoms.positions[:,0]-mid_cell_x)<-((self.rI)-(2*self.cutoff)))&\
+                (((atoms.positions[:,0]-mid_cell_x)>-((self.rI)+((2*self.cutoff)))))
+        y_criteria = ((atoms.positions[:,1]-mid_cell_y)<(np.max(u_cle_solution[:,1])+(2*self.cutoff)))&\
+            ((atoms.positions[:,1]-mid_cell_y)>(-(np.max(u_cle_solution[:,1])+(2*self.cutoff))))
+        
+        
+        #relax_mask = np.logical_and(np.logical_and(x_criteria,y_criteria),coordination_criteria)
+        relax_mask = (np.logical_and(x_criteria,y_criteria))
+
+        #--------output for testing------------
+        atoms.new_array('relax_region', np.zeros(len(atoms), dtype=int))
+        atoms.new_array('coord_criteria', np.zeros(len(atoms), dtype=int))
+        atoms.new_array('x_criteria', np.zeros(len(atoms), dtype=int))
+        atoms.new_array('y_criteria', np.zeros(len(atoms), dtype=int))
+        relax_region = atoms.arrays['relax_region']
+        coord_criteria_arr = atoms.arrays['coord_criteria']
+        x_criteria_arr = atoms.arrays['x_criteria']
+        y_criteria_arr = atoms.arrays['y_criteria']
+        relax_region[relax_mask] = 1
+        coord_criteria_arr[coordination_criteria] = 1
+        x_criteria_arr[x_criteria] = 1
+        y_criteria_arr[y_criteria] = 1
+        print('number of atoms to relax =', len(relax_region[relax_mask]))
+        ase.io.write('relaxed_atom_mask.xyz',atoms)
+        #-----------------------------------------
+
+        print('Relaxing', len(relax_region[relax_mask]),'surface atoms.')
+        #perform relaxation
+        atoms.calc = self.calc
+        atoms.set_constraint(FixAtoms(mask=~relax_mask))
+        opt = PreconLBFGS(atoms)
+        opt.run(fmax=1e-3,steps=25)
+
+        #subtract u_cle and u from the solution
+        atoms.positions -= u_cle_solution
+        atoms.positions[self.regionI, :] -= self.u
+
+        #set the cryst.atom.positions to the solution
+        self.cryst.positions = atoms.positions
+
+        #update the positions of the atoms
+        self.update_atoms()
+
+        #write to file for testing
+        self.write_atoms_to_file()
+
+#------------------------------------------------------------------------
+
     def fit_cle(self, r_fit=20.0, variable_alpha=True, variable_k=True, x0=None,
                 grid=None):
         def residuals(x, mask):
@@ -1074,6 +1194,12 @@ def residuals(x, mask):
                 raise RuntimeError('CLE fit failed')
         return res
 
+    def get_deformation_gradient(self,r,theta,k):
+        return self.crk.crack.deformation_gradient(r, theta, k)
+
+    def set_shiftmask(self,radial_dist):
+        self.shiftmask = self.r>radial_dist
+    
     def update_atoms(self):
         """
         Update self.atoms from degrees of freedom (self.u, self.alpha, self.k)
@@ -1088,12 +1214,57 @@ def update_atoms(self):
         self.atoms.positions[:, :] = self.cryst.positions
         self.atoms.positions[:, :] += self.k * self.u_cle()
         self.atoms.positions[self.regionI, :] += self.u
+        if self.crk.cb is not None: #if the crack has a multilattice cauchy-born object
+            #get rotation matrix 
+            A = np.transpose(self.crk.RotationMatrix)
+            #find shifts
+            print('finding shifts.....')
+            # very important to pass cryst rather than atoms here as the displacement gradient field
+            # is found from the positions of the original atoms, not the deformed atoms
+            shifts = self.crk.cb.predict_shifts(A,self.cryst,\
+                F_func=self.get_deformation_gradient, coordinates='cylind2D',method='lasso',k=self.k)
+            print('done!')
+            #apply shifts
+            self.crk.cb.apply_shifts(self.atoms,shifts,mask=self.shiftmask)
+
 
         # add vacuum
         self.atoms.cell = self.cryst.cell
         self.atoms.cell[0, 0] += self.vacuum
         self.atoms.cell[1, 1] += self.vacuum
 
+    def get_f_alpha_correction(self):
+        #function which gets the f_alpha contribution from regions II and III without u being applied
+        #'update atoms' but don't add u
+        self.atoms.set_pbc([False, False, True])
+        self.atoms.calc = self.calc
+
+        self.atoms.info['k'] = self.k
+        self.atoms.info['alpha'] = self.alpha
+
+        # x = x_cryst + K * u_cle + u
+        self.atoms.positions[:, :] = self.cryst.positions
+        self.atoms.positions[:, :] += self.k * self.u_cle()
+        if self.crk.cb is not None: #if the crack has a multilattice cauchy-born object
+            #get rotation matrix 
+            A = np.transpose(self.crk.RotationMatrix)
+            #find shifts
+            print('finding shifts.....')
+            # very important to pass cryst rather than atoms here as the displacement gradient field
+            # is found from the positions of the original atoms, not the deformed atoms
+            shifts = self.crk.cb.predict_shifts(A,self.cryst,\
+                F_func=self.get_deformation_gradient, coordinates='cylind2D',method='lasso',k=self.k)
+            print('done!')
+            #apply shifts
+            self.crk.cb.apply_shifts(self.atoms,shifts,mask=self.shiftmask)
+        if self.extended_far_field:
+            mask = self.regionII|self.regionIII
+        else:
+            mask = self.regionII
+        f_alpha = self.get_crack_tip_force(mask=mask)
+        self.update_atoms()
+        return f_alpha
+
     def set_atoms(self, atoms):
         N1_in = (atoms.arrays['region'] == 1).sum()
         if 'alpha' in atoms.info:
@@ -1105,13 +1276,13 @@ def set_atoms(self, atoms):
         self.update_atoms()  # now we have same u_cle in atoms and self.atoms
         min_len = min(N1_in, self.N1)
         # FIXME this assumes stable sort order for atoms and self.atoms
-        u = atoms.positions[:min_len] - self.atoms.positions[:min_len]
+        u = atoms.positions - self.atoms.positions
         shift = np.diag(self.atoms.cell)/2 - np.diag(atoms.cell)/2
         u += shift
-        self.u[:min_len] = u
+        self.u = u[self.regionI]
         self.update_atoms()
 
-    def get_crack_tip_force(self, forces=None, mask=None):
+    def get_crack_tip_force(self, forces=None, mask=None, full_array_output=False):
         # V_alpha = -\nabla_1 U_CLE(alpha)
         tip_x = self.cryst.cell.diagonal()[0] / 2.0 + self.alpha
         tip_y = self.cryst.cell.diagonal()[1] / 2.0
@@ -1145,6 +1316,10 @@ def get_crack_tip_force(self, forces=None, mask=None):
             mask = self.regionII
             if self.extended_far_field:
                 mask = self.regionII | self.regionIII
+        if full_array_output is True:
+            reduced_forces = forces[mask,:]
+            reduced_V = V[mask,:]
+            return np.tensordot(forces[mask, :], V[mask, :]), np.array([np.dot(reduced_forces[i, :], reduced_V[i, :]) for i in range(np.shape(reduced_forces)[0])])
         return np.tensordot(forces[mask, :], V[mask, :])
 
     def get_xdot(self, x1, x2, ds=None):
@@ -1187,10 +1362,18 @@ def get_forces(self, x1=None, xdot1=None, ds=None, forces=None, mask=None):
         F = list(forces[self.regionI, :].reshape(-1))
         if self.variable_alpha:
             f_alpha = self.get_crack_tip_force(forces, mask=mask)
-            F.append(f_alpha)
+            self.f_alpha_correction = 0.0 #self.get_f_alpha_correction()
+            F.append((f_alpha-self.f_alpha_correction)) #ninsert some kind of botched preconditioning.
+            print('f_alpha',f_alpha)
+            print('corrected f_alpha',(f_alpha-self.f_alpha_correction))
         if self.variable_k:
             f_k = self.get_k_force(x1, xdot1, ds)
             F.append(f_k)
+            print('f_k', f_k)
+
+        #self.norm_F_vals.append(np.linalg.norm(F,np.inf))
+        #self.F_alpha_vals.append(f_alpha)
+        #self.alpha_vals.append(self.alpha)
         return np.array(F)
 
     def update_precon(self, x, F=None):
@@ -1285,6 +1468,11 @@ def log(x, f=None):
 
         def residuals(x, *args):
             self.set_dofs(x)
+            #self.write_atoms_to_file() #comment out at some point
+            #if self.extended_far_field:
+            #    mask = self.regionII|self.regionIII
+            #else:
+            #    mask = self.regionII
             return self.get_forces(*args)
 
         def cg_objective(x):
@@ -1432,7 +1620,7 @@ def odeprecon(f,x):
             return res.nit
         else:
             self.atoms.write('no_convergence.xyz')
-            raise NoConvergence
+            raise RuntimeError(f"no convergence of scipy optimizer {method}")
 
     def get_potential_energy(self):
         # E1: energy of region I and II atoms
@@ -1514,9 +1702,9 @@ def arc_length_continuation(self, x0, x1, N=10, ds=0.01, ftol=1e-2,
             self.set_dofs(x2)
             try:
                 if opt_method == 'krylov':
-                    num_steps = self.optimize(ftol, max_steps, args=(x1, xdot1, ds), precon=precon, method=opt_method)
+                    num_steps = self.optimize(ftol, max_steps, args=(x1, xdot1, ds), precon=precon, method=opt_method,verbose=2)
                 elif opt_method == 'ode12r':
-                    num_steps = self.optimize(ftol, max_steps, args=[x1, xdot1, ds], precon=precon, method=opt_method)
+                    num_steps = self.optimize(ftol, max_steps, args=[x1, xdot1, ds], precon=precon, method=opt_method,verbose=2)
                 print(f'Corrector converged in {num_steps}/{max_steps} steps')
             except NoConvergence:
                 if ds < ds_min:
@@ -1687,6 +1875,201 @@ def frame(idx):
 
         return FuncAnimation(fig, frame, frames)
 
+    @counted
+    def write_atoms_to_file(self, fname):
+        alpha_at = Atom('Au')
+        crack_atoms = self.atoms.copy()
+
+        r = crack_atoms.arrays['region']
+        r[self.regionI] = 1
+        r[self.regionII] = 2
+        r[self.regionIII] = 3
+        r[self.regionIV] = 4
+
+        crack_atoms.calc = self.calc
+        forces = crack_atoms.get_forces()
+        crack_atoms.append(alpha_at)
+        crack_atoms[-1].position = [self.cryst.cell[0, 0] / 2.0 + self.alpha,
+                             self.cryst.cell[1, 1] / 2.0,0]
+        crack_atoms.new_array('fx', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('fy', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('fz', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('ftot', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('logabsftot', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('falphacomponents',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('logabsfalphacomp',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('exx', np.zeros(len(crack_atoms),dtype=float))
+        crack_atoms.new_array('eyy',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('ezz',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('eyz',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('exz',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('exy',np.zeros(len(crack_atoms)),dtype =float)
+
+        crack_atoms.new_array('ux', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('uy', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('uz', np.zeros(len(crack_atoms), dtype=float))
+        
+
+        fx = crack_atoms.arrays['fx']
+        fy = crack_atoms.arrays['fy']
+        fz = crack_atoms.arrays['fz']
+        ux = crack_atoms.arrays['ux']
+        uy = crack_atoms.arrays['uy']
+        uz = crack_atoms.arrays['uz']
+        ftot = crack_atoms.arrays['ftot']
+        logabsftot = crack_atoms.arrays['logabsftot']
+        falphacomponents = crack_atoms.arrays['falphacomponents']
+        logabsfalphacomp = crack_atoms.arrays['logabsfalphacomp']
+        store_strain=False
+        if store_strain:
+            exx = crack_atoms.arrays['exx']
+            eyy = crack_atoms.arrays['eyy']
+            ezz = crack_atoms.arrays['ezz']
+            eyz = crack_atoms.arrays['eyz']
+            exz = crack_atoms.arrays['exz']
+            exy = crack_atoms.arrays['exy']
+
+            #get strains applied
+            A = np.transpose(self.crk.RotationMatrix)
+            E,R = self.crk.cb.evaluate_F_or_E(A,self.cryst,F_func=self.get_deformation_gradient, coordinates='cylind2D',k=self.k)
+            for i in range(len(crack_atoms)-1):
+                E[i,:,:] = np.transpose(A)@E[i,:,:]@A #transform back to lab frame
+            #store values of E
+            exx[0:len(crack_atoms)-1] = E[:,0,0]
+            eyy[0:len(crack_atoms)-1] = E[:,1,1]
+            ezz[0:len(crack_atoms)-1] = E[:,2,2]
+            eyz[0:len(crack_atoms)-1] = E[:,1,2]
+            exz[0:len(crack_atoms)-1] = E[:,0,2]
+            exy[0:len(crack_atoms)-1] = E[:,0,1]
+        #store other things
+        fx[0:len(crack_atoms)-1] = (forces[:,0])
+        fy[0:len(crack_atoms)-1] = (forces[:,1])
+        fz[0:len(crack_atoms)-1] = (forces[:,2])
+        ftot[0:len(crack_atoms)-1] = (np.linalg.norm(forces,ord=2,axis=1))
+        ux[0:len(crack_atoms)-1][self.regionI] = (self.u[:,0])
+        uy[0:len(crack_atoms)-1][self.regionI] = (self.u[:,1])
+        uz[0:len(crack_atoms)-1][self.regionI] = (self.u[:,2])
+
+        mask = np.append(self.regionII,np.array([False],dtype=bool)) #extra false added as we have added alpha to array
+        crack_tip_force_mask = self.regionII
+        if self.extended_far_field:
+            mask = np.append((self.regionII | self.regionIII),np.array([False],dtype=bool))
+            crack_tip_force_mask = self.regionII | self.regionIII
+        
+        falpha, falphas = self.get_crack_tip_force(mask = crack_tip_force_mask,full_array_output=True)
+        logfalphas = np.log10(np.abs(falphas))
+        falphacomponents[mask] = falphas
+        logabsfalphacomp[mask] = logfalphas
+        fx[len(crack_atoms)-1] = falpha
+        ftot[len(crack_atoms)-1] = falpha
+        logabsftot[:] = np.log10(np.abs(ftot))
+        ase.io.write(fname+'.xyz',crack_atoms)
+    
+    def strain_err(self,cutoff,seperate_surface=False):
+        '''Function that returns the atomistic corrector strain error Dv for each atom using the norm of the difference
+        between the corrector on each atom and all those around it within some cutoff. Also has an option to return states
+        adjacent to the free surface of the crack seperately for comparison.
+
+        Parameters
+        ----------
+        cutoff : float
+            cutoff around each atom to use to find the strain error norm.
+        seperate_surface : bool
+            whether or not to return the surface atoms in a seperate array.
+
+        Returns
+        -------
+        r : array of atom radial distance from the centre of the crack system
+
+        dv : norm strain error for each atom in r.
+        """
+        
+        '''
+        #want to get all neighbours in region I
+        I, J = neighbour_list('ij', self.atoms[self.regionI], cutoff)
+        #print(I)
+        #print(J)
+        v = self.u
+        dv = np.linalg.norm(v[I, :] - v[J, :], axis=1)
+        r = self.r[self.regionI][I]
+        mask = r<(self.rI-cutoff)
+        if seperate_surface:
+            mid_cell_x = self.cryst.cell.diagonal()[0] / 2.0
+            mid_cell_y = self.cryst.cell.diagonal()[1] / 2.0
+            I_positions = self.cryst.get_positions()[I]
+            x_criteria = ((I_positions[:,0]-mid_cell_x)<0.5)
+            y_criteria = ((I_positions[:,1]-mid_cell_y)<((1.5*self.cutoff)))&\
+            ((I_positions[:,1]-mid_cell_y)>(-1.5*self.cutoff))
+            surface_mask_full = np.logical_and(x_criteria,y_criteria)
+            bulk_mask_full = np.logical_not(surface_mask_full)
+            surface_mask = np.logical_and(surface_mask_full,mask)
+            bulk_mask = np.logical_and(bulk_mask_full,mask)
+            return r[surface_mask],dv[surface_mask],r[bulk_mask],dv[bulk_mask]
+        else:
+            return r[mask],dv[mask]
+
+        
+    def convergence_line_plot(self,num=0):
+        def gen_mask(r,rval,dx,dr,cx,cy):
+            #for a line going vertically up
+            #so r is essentiall y here
+            rcriteria = (r > rval)&(r<(rval+dr))
+            xcriteria = (x-cx > 0)&(x-cx<dx)
+            return np.logical_and(rcriteria,xcriteria)
+        crack_atoms = self.atoms.copy()
+        crack_atoms.calc = self.calc
+        forces = crack_atoms.get_forces()
+
+        sx, sy, sz = crack_atoms.cell.diagonal()
+        x, y = crack_atoms.positions[:, 0], crack_atoms.positions[:, 1]
+        cx, cy = sx/2, sy/2
+        r = np.sqrt((x - cx)**2 + (y - cy)**2)
+
+        dr = 0.2 #Angstroms
+        dx = 1 #angstroms
+        regionIvals = np.arange(0,self.rI,dr)
+        other_region_vals = np.arange(self.rI,self.rIII,dr)
+        full_list = np.concatenate((regionIvals,other_region_vals))
+        U = np.zeros([len(regionIvals),3])
+        F = np.zeros([len(full_list),3])
+        for i,rval in enumerate(regionIvals):
+            mask = gen_mask(r,rval,dx,dr,cx,cy)
+            if len(forces[mask])==0:
+                U[i,:] = np.zeros([1,3])
+                F[i,:] = np.zeros([1,3])
+            else:
+                print(f'non 0 at r={rval}')
+                print('vals',self.u[mask[self.regionI]])
+                print('mean',np.mean(self.u[mask[self.regionI]],axis=0))
+                U[i,:] = np.mean(self.u[mask[self.regionI]],axis=0)
+                F[i,:] = np.mean(forces[mask],axis=0)
+
+        for i,rval in enumerate(other_region_vals):
+            mask = gen_mask(r,rval,dx,dr,cx,cy)
+            k = i+len(regionIvals)
+            print(forces[mask])
+            if len(forces[mask])==0:
+                F[k,:] = np.zeros([1,3])
+            else:
+                print(f'non 0 at r={rval}')
+                F[k,:] = np.mean(forces[mask],axis=0)
+
+        plt.figure()
+        for i in range(3):
+            plt.plot(regionIvals,U[:,i])
+        plt.xlabel('r')
+        plt.ylabel('atomistic corrector U at r')
+        plt.legend(['Ux','Uy','Uz'])
+        plt.savefig(f'AtomisticCorrector{num}.png')
+        
+        plt.figure()
+        for i in range(3):
+            plt.plot(full_list,F[:,i])
+        plt.xlabel('r')
+        plt.ylabel('Force on atoms at r')
+        plt.legend(['Fx','Fy','Fz'])
+        plt.savefig(f'ForcesOnAtoms{num}.png')
+
 
 
 def isotropic_modeI_crack_tip_stress_field(K, r, t, xy_only=True,
@@ -2249,7 +2632,7 @@ def find_tip_stress_field(atoms, r_range=None, initial_params=None, fix_params=N
 
 def plot_stress_fields(atoms, r_range=None, initial_params=None, fix_params=None,
                        sigma=None, avg_sigma=None, avg_decay=0.005, calc=None):
-    """
+    r"""
     Fit and plot atomistic and continuum stress fields
 
     Firstly a fit to the Irwin `K`-field solution is carried out using
diff --git a/matscipy/io/tbl.py b/matscipy/io/tbl.py
index c0f99e5d..ce8e04dc 100644
--- a/matscipy/io/tbl.py
+++ b/matscipy/io/tbl.py
@@ -85,7 +85,7 @@ def loadtbl(fn, usecols=None, types=None, fromfile=False, **kwargs):
     column_labels = None
     while line.startswith('#'):
         line = line[1:].strip()
-        column_labels = [s.strip() for s in re.split('[\s,]+', line)]
+        column_labels = [s.strip() for s in re.split(r'[\s,]+', line)]
         pos = f.tell()
         line = f.readline()
     f.seek(pos)
diff --git a/matscipy/meson.build b/matscipy/meson.build
index ddeb4644..bd5a2ba3 100644
--- a/matscipy/meson.build
+++ b/matscipy/meson.build
@@ -44,7 +44,6 @@ python.install_sources(
 )
 
 subdir('calculators')
-subdir('contact_mechanics')
 subdir('electrochemistry')
 subdir('fracture_mechanics')
 subdir('io')
diff --git a/matscipy/rings.py b/matscipy/rings.py
index 02d843ab..7217bdef 100644
--- a/matscipy/rings.py
+++ b/matscipy/rings.py
@@ -50,4 +50,12 @@ def ring_statistics(a, cutoff, maxlength=-1):
     """
     i, j, r = neighbour_list('ijD', a, cutoff)
     d = distances_on_graph(i, j)
+
+    if maxlength > 0:
+        ringstat = np.zeros(maxlength)
+        rs = find_sp_rings(i, j, r, d, maxlength)
+        ringstat[:len(rs)] += rs
+    else:
+        ringstat = find_sp_rings(i, j, r, d, maxlength)
+
     return find_sp_rings(i, j, r, d, maxlength)
diff --git a/matscipy/surface_reconstruction.py b/matscipy/surface_reconstruction.py
new file mode 100644
index 00000000..0a1614de
--- /dev/null
+++ b/matscipy/surface_reconstruction.py
@@ -0,0 +1,211 @@
+import numpy as np
+from ase.lattice.cubic import Diamond, FaceCenteredCubic, SimpleCubic, BodyCenteredCubic
+
+from ase.optimize import LBFGS
+
+from ase.constraints import FixAtoms
+
+from matscipy.cauchy_born import CubicCauchyBorn
+from ase.optimize.precon import PreconLBFGS
+import ase.io
+
+
+
+class SurfaceReconstruction:
+    def __init__(self,el,a0,calc,directions,surf_dir,lattice='diamond'):
+
+        #directions is a list of 3 vectors specifiying the surface of interest
+        #in the orientation specified by the problem. surf_dir says which axis is free surface
+        #e.g directions=[[1,1,1], [0,-1,1], [-2,1,1]]
+        self.el = el
+        self.a0 = a0
+        self.calc = calc
+        self.directions = directions
+        self.surf_dir = surf_dir
+        A = np.zeros([3,3])
+        for i,direction in enumerate(directions):
+            direction = np.array(direction)
+            A[:,i] = direction/np.linalg.norm(direction)
+        self.A = A #rotation matrix
+        self.eval_cb = None
+        if lattice == 'diamond':
+            self.lattice = Diamond
+            self.cb = CubicCauchyBorn(self.el,self.a0,self.calc,lattice=Diamond)
+
+    def map_surface(self,fmax=0.0001,layers=6):
+        #set the number of atomic layers to model
+        height = 30
+        size = [1,1,1]
+        size[self.surf_dir]*=height
+        self.layers = layers
+
+        #build the bulk structure, and the slab structure which has vacuum in surface direction
+        bulk = self.lattice(directions=self.directions, size=size, symbol=self.el, latticeconstant=self.a0, pbc=(1,1,1))
+        cell = bulk.get_cell() ; cell[self.surf_dir,:] *=2 # vacuum along self.surf_dir axis (surface normal)
+        slab = bulk.copy() ; slab.set_cell(cell,scale_atoms=False)
+        
+        #get the mask for the bottom surface and fix it
+        pos_before = slab.get_positions()
+        mask = pos_before[:,self.surf_dir]<30 #fix bottom atoms to stop lower surface reconstruction
+        slab.set_constraint(FixAtoms(mask=mask))
+        ase.io.write('initial_slab.xyz',slab)
+        slab.calc = self.calc
+
+        #run an optimisation to relax the surface
+        opt_slab = PreconLBFGS(slab) ; opt_slab.run(fmax=fmax)
+
+        pos_after = slab.get_positions()
+
+
+        #measure the shift between the final position and the intial position
+        pos_shift = pos_after-pos_before
+
+        #for a multi-lattice:
+        if self.cb is not None:
+            #if this is a multi-lattice
+
+            #set which atoms are in each multilattice based on the bulk structure. This is important because
+            #different relaxations happen for each multilattice atom in each atomic layer near the surface
+            self.cb.set_sublattices(bulk,self.A)
+            #in a multilattice, there are 2 atomic layers per unit cell, so set the height to double the atomic layers
+            total_layers = 2*height
+            n_atoms_per_layer = int(len(slab)/total_layers)
+
+            #get the indices of the atoms based on how close they are to the surface in order, for each sub lattice
+            sorted_surf_indices_lattice1 = np.argsort(pos_before[self.cb.lattice1mask][:,self.surf_dir])
+            sorted_surf_indices_lattice2 = np.argsort(pos_before[self.cb.lattice2mask][:,self.surf_dir])
+            
+            #split these into layers and save them to an array.
+            #each layer holds the indices of the atoms it contains [note - these are the indices of the array
+            # pos_shift[self.cb.latticemask] - not pos_shift]
+            surf_layer_list_lattice1 = np.zeros([total_layers,int(n_atoms_per_layer/2)],dtype=int)
+            surf_layer_list_lattice2 = np.zeros([total_layers,int(n_atoms_per_layer/2)],dtype=int)
+            for i in range(total_layers):
+                surf_layer_list_lattice1[i,:] = sorted_surf_indices_lattice1[i*int(n_atoms_per_layer/2):(i+1)*int(n_atoms_per_layer/2)]
+                surf_layer_list_lattice2[i,:] = sorted_surf_indices_lattice2[i*int(n_atoms_per_layer/2):(i+1)*int(n_atoms_per_layer/2)]
+
+            #get the inter-layer distance in the direction we care about.
+            self.inter_surf_dist = np.abs(pos_before[self.cb.lattice1mask][surf_layer_list_lattice1[1,:]][0,self.surf_dir]\
+                 - pos_before[self.cb.lattice1mask][surf_layer_list_lattice1[0,:]][0,self.surf_dir])
+            
+            #now go through the top layers that we are applying a surface relaxation to and save the shifts of 
+            # the lattice 1 and lattice 2 atoms in these layers in the form [x,y,z] for each layer
+        
+            self.relaxation_array_lattice1 = np.zeros([layers,3])
+            self.relaxation_array_lattice2 = np.zeros([layers,3])
+            for layer in range(layers):
+                pos_shifts_lattice1 = pos_shift[self.cb.lattice1mask][surf_layer_list_lattice1[total_layers-layer-1,:]]
+                pos_shifts_lattice2 = pos_shift[self.cb.lattice2mask][surf_layer_list_lattice2[total_layers-layer-1,:]]
+
+                self.relaxation_array_lattice1[layer,:] = pos_shifts_lattice1[0,:]
+                self.relaxation_array_lattice2[layer,:] = pos_shifts_lattice2[0,:]
+        
+        else:
+            #for a non-multilattice #need to look at this ngl
+            total_layers = height
+            n_atoms_per_layer = int(len(slab)/height)
+            sorted_surf_indices = np.argsort(pos_before[:,self.surf_dir])
+            surf_layer_list = np.zeros([height,n_atoms_per_layer],dtype=int)
+            for i in range(height):
+                surf_layer_list[i,:] = sorted_surf_indices[i*n_atoms_per_layer:((i+1)*n_atoms_per_layer)]
+            
+            self.inter_surf_dist = np.abs(pos_before[surf_layer_list[1,:]][0,surf_dir] - pos_before[surf_layer_list[0,:]][0,surf_dir])
+            self.relaxation_array = np.zeros([layers,3])
+            for layer in range(layers):
+                pos_shifts_by_layer = pos_shift[surf_layer_list[layer,:]]
+                self.relaxation_array[layer,:] = pos_shifts_by_layer[0,:]
+            
+        
+
+    def identify_layers(self,atoms,surface_coords,xlim=None,ylim=None,zlim=None,search_dir=-1,atoms_for_cb=None):
+        #function for identifying the layers adjacent to a surface, where the top of the surface
+        #is found near coordinates surface_coords
+        pos = atoms.get_positions()
+
+        #get the masks for the limits imposed by xlim,ylim and zlim
+        if xlim is not None:
+            xmask = (pos[:,0]<xlim[1])&(pos[:,0]>xlim[0])
+        else:
+            xmask = np.ones([len(atoms)],dtype=bool)
+        
+        if ylim is not None:
+            ymask = (pos[:,1]<ylim[1])&(pos[:,1]>ylim[0])
+        else:
+            ymask = np.ones([len(atoms)],dtype=bool)
+
+        if zlim is not None:
+            zmask = (pos[:,2]<zlim[1])&(pos[:,2]>zlim[0])
+        else:
+            zmask = np.ones([len(atoms)],dtype=bool)
+
+        lim_mask = xmask & ymask & zmask
+        
+        #now, get the mask for the provided atoms structure
+        #for each individual surface layer
+        layer_mask_set = np.zeros([len(atoms),self.layers],dtype=bool)
+        for layer in range(self.layers):
+            #this is the crucial line, essentially finds a mask for the
+            #atoms that lie within a +- 1 layer distance from the starting point.
+            #we then shift down or up a layer each iteration, depending on search_dir
+            if search_dir == -1:
+                #this corresponds to searching downwards for surface layers
+                layer_mask = np.logical_and((pos[:,self.surf_dir]<=surface_coords[self.surf_dir]\
+                    -(layer*self.inter_surf_dist)),(pos[:,self.surf_dir]\
+                    >surface_coords[self.surf_dir]-((layer+1)*self.inter_surf_dist)))
+            elif search_dir == 1:
+                #this corresponds to searching upward for surface layers
+                layer_mask = np.logical_and((pos[:,self.surf_dir]>=surface_coords[self.surf_dir]\
+                    +(layer*self.inter_surf_dist)),(pos[:,self.surf_dir]\
+                    <surface_coords[self.surf_dir]+((layer+1)*self.inter_surf_dist)))
+            else:
+                print('search_dir should be -1 for a downward surface, 1 for an upward surface!')
+                raise ValueError
+            layer_mask = layer_mask&lim_mask
+            
+            layer_mask_set[:,layer] = layer_mask #add the layer mask found here to the full set
+        
+        
+        if self.eval_cb is None:
+            return layer_mask_set
+        else:
+            #for a multi-lattice, also have to account for shifts between multilattice atoms
+            #in layers, so need to refine the masks further to be the lattice1 and lattice2 atoms in each layer
+            if atoms_for_cb is None:
+                atoms_for_cb = atoms
+            #note that to find the sublattices, need to use a version of the atoms structure with no free surface
+            if self.eval_cb.lattice1mask is None:
+                self.eval_cb.set_sublattices(atoms_for_cb,self.A)
+            layer_mask_set_lattice1 = np.zeros_like(layer_mask_set,dtype=bool)
+            layer_mask_set_lattice2 = np.zeros_like(layer_mask_set,dtype=bool)
+            for layer in range(self.layers):
+                layer_mask_set_lattice1[:,layer] = layer_mask_set[:,layer]&self.eval_cb.lattice1mask
+                layer_mask_set_lattice2[:,layer] = layer_mask_set[:,layer]&self.eval_cb.lattice2mask
+            return layer_mask_set_lattice1, layer_mask_set_lattice2
+
+        
+
+    def apply_surface_shift(self,atoms,surface_coords,cb=None,xlim=None,ylim=None,zlim=None,search_dir=-1,atoms_for_cb=None):
+        #This finds the layers of atoms that need shifting and applies the surface shift accordingly.
+        self.eval_cb = cb
+        if self.eval_cb is not None:
+            layer_mask_set_lattice1,layer_mask_set_lattice2 = self.identify_layers(atoms,surface_coords,xlim=xlim,ylim=ylim,zlim=zlim,search_dir=search_dir,atoms_for_cb=atoms_for_cb)
+            for layer in range(self.layers):
+                #print('diff',self.relaxation_array_lattice1[layer,:])
+                #print('initial',atoms.positions[:, :][layer_mask_set_lattice1[:,layer]])
+                lattice_1_shift = self.relaxation_array_lattice1[layer,:]
+                lattice_2_shift = self.relaxation_array_lattice2[layer,:]
+
+                #if the surface is going upward, invert surf_dir component of applied shift
+                lattice_1_shift[self.surf_dir] = (-search_dir)*lattice_1_shift[self.surf_dir]
+                lattice_2_shift[self.surf_dir] = (-search_dir)*lattice_2_shift[self.surf_dir]
+                atoms.positions[:, :][layer_mask_set_lattice1[:,layer]] += lattice_1_shift
+                atoms.positions[:, :][layer_mask_set_lattice2[:,layer]] += lattice_2_shift
+                #print('after',atoms.positions[:, :][layer_mask_set_lattice1[:,layer]])
+
+        else:
+            layer_mask_set = self.identify_layers(atoms,surface_coords,xlim=xlim,ylim=ylim,zlim=zlim,search_dir=search_dir)
+            for layer in range(self.layers):
+                shift = self.relaxation_array[layer,:]
+                shift[self.surf_dir] = (-search_dir)*shift[self.surf_dir]
+                atoms.positions[:,:][layer_mask_set[:,layer]] += (-search_dir)*self.relaxation_array      
+        
\ No newline at end of file
diff --git a/meson.build b/meson.build
index 167b3c77..97371eef 100644
--- a/meson.build
+++ b/meson.build
@@ -9,6 +9,7 @@ project(
 pymod = import('python')
 python = pymod.find_installation('python3',
     required: true,
+    pure: false
 )
 
 subdir('c')
diff --git a/pyproject.toml b/pyproject.toml
index 04d2d35a..dcb05288 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,5 +1,5 @@
 [build-system]
-requires = ["meson>=0.46.0", "meson-python>=0.11.0", "ninja", "oldest-supported-numpy"]
+requires = ["meson>=1.0.0", "meson-python>=0.13.0", "ninja", "oldest-supported-numpy"]
 build-backend = "mesonpy"
 
 [project]


Test patch:
diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml
index ab8905f0..1aa952b5 100644
--- a/.github/workflows/run_tests.yml
+++ b/.github/workflows/run_tests.yml
@@ -23,7 +23,9 @@ jobs:
       - name: install_python
         run: |
           sudo apt-get update -qy
-          sudo apt-get install -y python3-dev python3-pip python3-venv libxml2-dev libxslt-dev zlib1g-dev meson
+          sudo apt-get install -y python3-dev python3-pip python3-venv libxml2-dev libxslt-dev zlib1g-dev
+          # Upgrade to latest meson and ninja
+          sudo pip install --upgrade meson ninja
 
       - name: build_c
         run: |
diff --git a/tests/committee_data/test_data.xyz b/tests/committee_data/test_data.xyz
new file mode 100644
index 00000000..58188aee
--- /dev/null
+++ b/tests/committee_data/test_data.xyz
@@ -0,0 +1,10 @@
+3
+Properties=species:S:1:pos:R:3:forces_uncertainty:R:3:forces:R:3:F_lj:R:3 energy_uncertainty=0.006221848541648684 energy=-0.9743916502614107 E_lj=-0.9792880907149858 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       0.01095547       0.00000000       0.00000000       1.74328636       0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000       0.01030488       0.00000000       0.00000000      -1.63976171       0.00000000       0.00000000      -1.64800172       0.00000000       0.00000000
+O        2.16954488       0.00000000       0.00000000       0.00065059       0.00000000       0.00000000      -0.10352465       0.00000000       0.00000000      -0.10404487       0.00000000       0.00000000
+3
+Properties=species:S:1:pos:R:3:forces_uncertainty:R:3:forces:R:3:F_lj:R:3 energy_uncertainty=0.00619385497395605 energy=-0.9700076318401718 E_lj=-0.9748820420504242 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       0.01095547       0.00000000       0.00000000       1.74328636       0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000       0.01039149       0.00000000       0.00000000      -1.65354324       0.00000000       0.00000000      -1.66185250       0.00000000       0.00000000
+O        2.21499943       0.00000000       0.00000000       0.00056398       0.00000000       0.00000000      -0.08974312       0.00000000       0.00000000      -0.09019409       0.00000000       0.00000000
diff --git a/tests/committee_data/training_data.xyz b/tests/committee_data/training_data.xyz
new file mode 100644
index 00000000..cddf0ab6
--- /dev/null
+++ b/tests/committee_data/training_data.xyz
@@ -0,0 +1,150 @@
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=0 E_lj=447.40100301767563 appears_in_committee="1 3 5 6 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       2.07189248       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000   -8436.80972524       0.00000000       0.00000000
+O        0.66954488       0.00000000      -0.00000000    8434.73783277       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=1 E_lj=37.6864504997281 appears_in_committee="1 2 5 6 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.94959641       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000    -686.51940142       0.00000000       0.00000000
+O        0.80590852       0.00000000      -0.00000000     684.56980502       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=2 E_lj=-0.6438122260115658 appears_in_committee="1 2 5 6 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.86023314       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000     -31.94401496       0.00000000       0.00000000
+O        0.98772670       0.00000000      -0.00000000      30.08378183       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=3 E_lj=-1.9630289630844797 appears_in_committee="1 2 5 6 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.82292822       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.66040137       0.00000000       0.00000000
+O        1.12409033       0.00000000      -0.00000000      -0.16252684       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=4 E_lj=-1.596993146139846 appears_in_committee="1 3 5 7 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.79381392       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000       0.45915589       0.00000000       0.00000000
+O        1.30590852       0.00000000       0.00000000      -2.25296981       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=5 E_lj=-1.3434096765111547 appears_in_committee="0 3 4 6 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.78081177       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -0.31408295       0.00000000       0.00000000
+O        1.44227215       0.00000000       0.00000000      -1.46672882       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=6 E_lj=-1.1868293849353155 appears_in_committee="1 3 4 7 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.77222437       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -0.89669474       0.00000000       0.00000000
+O        1.57863579       0.00000000       0.00000000      -0.87552963       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=7 E_lj=-1.0719901327713648 appears_in_committee="1 2 5 6 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.76495177       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.32484452       0.00000000       0.00000000
+O        1.76045397       0.00000000       0.00000000      -0.44010725       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=8 E_lj=-1.0253375039936834 appears_in_committee="0 3 4 6 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.49204815       0.00000000       0.00000000
+O        1.89681761       0.00000000       0.00000000      -0.25999844       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=9 E_lj=-0.9902733045162482 appears_in_committee="1 3 4 6 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.61245525       0.00000000       0.00000000
+O        2.07863579       0.00000000       0.00000000      -0.13959134       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=10 E_lj=-0.9748820420504242 appears_in_committee="0 2 4 6 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.66185250       0.00000000       0.00000000
+O        2.21499943       0.00000000       0.00000000      -0.09019409       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=11 E_lj=-0.9622857407222333 appears_in_committee="0 2 4 7 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.69978766       0.00000000       0.00000000
+O        2.39681761       0.00000000       0.00000000      -0.05225893       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=12 E_lj=-0.9563790996251943 appears_in_committee="1 3 4 7 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.71646316       0.00000000       0.00000000
+O        2.53318124       0.00000000       0.00000000      -0.03558343       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=13 E_lj=-0.952319955906694 appears_in_committee="1 2 4 7 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.72734391       0.00000000       0.00000000
+O        2.66954488       0.00000000       0.00000000      -0.02470268       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=14 E_lj=-0.9487276967963708 appears_in_committee="0 2 4 7 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.73644301       0.00000000       0.00000000
+O        2.85136306       0.00000000       0.00000000      -0.01560358       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=15 E_lj=-0.9469142984200198 appears_in_committee="0 2 5 7 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.74078493       0.00000000       0.00000000
+O        2.98772670       0.00000000       0.00000000      -0.01126166       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=16 E_lj=-0.946778041737247 appears_in_committee="0 3 5 7 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.16954488       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=17 E_lj=-0.946778041737247 appears_in_committee="1 2 5 6 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.30590852       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=18 E_lj=-0.946778041737247 appears_in_committee="0 2 4 6 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.44227215       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=19 E_lj=-0.946778041737247 appears_in_committee="0 2 4 7 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.62409033       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=20 E_lj=-0.946778041737247 appears_in_committee="0 3 5 7 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.76045397       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=21 E_lj=-0.946778041737247 appears_in_committee="1 2 5 7 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.94227215       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=22 E_lj=-0.946778041737247 appears_in_committee="0 3 4 7 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.07863579       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=23 E_lj=-0.946778041737247 appears_in_committee="1 3 5 7 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.26045397       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=24 E_lj=-0.946778041737247 appears_in_committee="0 3 4 7 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.39681761       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=25 E_lj=-0.946778041737247 appears_in_committee="1 3 5 6 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.53318124       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=26 E_lj=-0.946778041737247 appears_in_committee="0 2 4 6 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.71499943       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=27 E_lj=-0.946778041737247 appears_in_committee="1 2 4 7 9" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.85136306       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=28 E_lj=-0.946778041737247 appears_in_committee="0 3 5 6 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        5.03318124       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=29 E_lj=-0.946778041737247 appears_in_committee="0 3 5 6 8" pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        5.16954488       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
diff --git a/tests/committee_data/training_data_minimal.xyz b/tests/committee_data/training_data_minimal.xyz
new file mode 100644
index 00000000..811a5c23
--- /dev/null
+++ b/tests/committee_data/training_data_minimal.xyz
@@ -0,0 +1,10 @@
+3
+Properties=species:S:1:pos:R:3 pbc="F F F" _Index_FullTrainingSet=3
+O       -1.16954488      -0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000
+O        0.66954488       0.00000000      -0.00000000
+3
+Properties=species:S:1:pos:R:3 pbc="F F F" _Index_FullTrainingSet=4
+O       -1.16954488      -0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000
+O        0.71499943       0.00000000      -0.00000000
diff --git a/tests/committee_data/validation_set.xyz b/tests/committee_data/validation_set.xyz
new file mode 100644
index 00000000..a27b217e
--- /dev/null
+++ b/tests/committee_data/validation_set.xyz
@@ -0,0 +1,150 @@
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=8 E_lj=-1.0253375039936834 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.49204815       0.00000000       0.00000000
+O        1.89681761       0.00000000       0.00000000      -0.25999844       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=4 E_lj=-1.596993146139846 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.79381392       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000       0.45915589       0.00000000       0.00000000
+O        1.30590852       0.00000000       0.00000000      -2.25296981       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=0 E_lj=447.40100301767563 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       2.07189248       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000   -8436.80972524       0.00000000       0.00000000
+O        0.66954488       0.00000000      -0.00000000    8434.73783277       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=27 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.85136306       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=7 E_lj=-1.0719901327713648 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.76495177       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.32484452       0.00000000       0.00000000
+O        1.76045397       0.00000000       0.00000000      -0.44010725       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=16 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.16954488       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=13 E_lj=-0.952319955906694 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.72734391       0.00000000       0.00000000
+O        2.66954488       0.00000000       0.00000000      -0.02470268       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=2 E_lj=-0.6438122260115658 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.86023314       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000     -31.94401496       0.00000000       0.00000000
+O        0.98772670       0.00000000      -0.00000000      30.08378183       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=19 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.62409033       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=12 E_lj=-0.9563790996251943 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.71646316       0.00000000       0.00000000
+O        2.53318124       0.00000000       0.00000000      -0.03558343       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=10 E_lj=-0.9748820420504242 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.66185250       0.00000000       0.00000000
+O        2.21499943       0.00000000       0.00000000      -0.09019409       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=11 E_lj=-0.9622857407222333 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.69978766       0.00000000       0.00000000
+O        2.39681761       0.00000000       0.00000000      -0.05225893       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=9 E_lj=-0.9902733045162482 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.61245525       0.00000000       0.00000000
+O        2.07863579       0.00000000       0.00000000      -0.13959134       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=22 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.07863579       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=15 E_lj=-0.9469142984200198 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.74078493       0.00000000       0.00000000
+O        2.98772670       0.00000000       0.00000000      -0.01126166       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=6 E_lj=-1.1868293849353155 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.77222437       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -0.89669474       0.00000000       0.00000000
+O        1.57863579       0.00000000       0.00000000      -0.87552963       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=3 E_lj=-1.9630289630844797 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.82292822       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.66040137       0.00000000       0.00000000
+O        1.12409033       0.00000000      -0.00000000      -0.16252684       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=23 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.26045397       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=5 E_lj=-1.3434096765111547 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.78081177       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -0.31408295       0.00000000       0.00000000
+O        1.44227215       0.00000000       0.00000000      -1.46672882       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=1 E_lj=37.6864504997281 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.94959641       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000    -686.51940142       0.00000000       0.00000000
+O        0.80590852       0.00000000      -0.00000000     684.56980502       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=14 E_lj=-0.9487276967963708 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.73644301       0.00000000       0.00000000
+O        2.85136306       0.00000000       0.00000000      -0.01560358       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=25 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.53318124       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=26 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.71499943       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=28 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        5.03318124       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=21 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.94227215       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=20 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.76045397       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=17 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.30590852       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=18 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        3.44227215       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=29 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        5.16954488       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
+3
+Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:F_lj:R:3 _Index_FullTrainingSet=24 E_lj=-0.946778041737247 pbc="F F F"
+O       -1.16954488      -0.00000000       0.00000000       1.75204659       0.00000000       0.00000000
+C        0.00000000       0.00000000      -0.00000000      -1.75204659       0.00000000       0.00000000
+O        4.39681761       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000
diff --git a/tests/test_analysis.py b/tests/test_analysis.py
deleted file mode 100644
index 69b7268d..00000000
--- a/tests/test_analysis.py
+++ /dev/null
@@ -1,210 +0,0 @@
-#
-# Copyright 2015, 2020 Lars Pastewka (U. Freiburg)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-from math import sqrt
-import unittest
-
-import numpy as np
-
-from matscipy.contact_mechanics.analysis import (count_islands,
-                                                 count_segments,
-                                                 distance_map,
-                                                 inner_perimeter,
-                                                 outer_perimeter)
-
-###
-
-class TestAnalysis(unittest.TestCase):
-
-    def test_count_islands(self):
-        m_xy = np.zeros([3,3], dtype=bool)
-        m_xy[1,1] = True
-
-        nump, p_xy = count_islands(m_xy)
-       
-        self.assertEqual(nump, 1)
-        self.assertTrue(np.all(
-            p_xy == np.array([[0,0,0],
-                              [0,1,0],
-                              [0,0,0]])))
-                              
-        m_xy = np.zeros([5,3], dtype=bool)
-        m_xy[1,1] = True
-        m_xy[3,1] = True
-
-        nump, p_xy = count_islands(m_xy)
-
-        self.assertEqual(nump, 2)
-        self.assertTrue(np.all(
-            p_xy == np.array([[0,0,0],
-                              [0,1,0],
-                              [0,0,0],
-                              [0,2,0],
-                              [0,0,0]])))
-
-        m_xy = np.zeros([6,3], dtype=bool)
-        m_xy[1,1] = True
-        m_xy[3,1] = True
-        m_xy[3,2] = True
-
-        nump, p_xy = count_islands(m_xy)
-
-        self.assertEqual(nump, 2)
-        self.assertTrue(np.all(
-            p_xy == np.array([[0,0,0],
-                              [0,1,0],
-                              [0,0,0],
-                              [0,2,2],
-                              [0,0,0],
-                              [0,0,0]])))
-
-        m_xy = np.zeros([6,3], dtype=bool)
-        m_xy[1,1] = True
-        m_xy[2,1] = True
-        m_xy[3,1] = True
-        m_xy[3,2] = True
-
-        nump, p_xy = count_islands(m_xy)
-
-        self.assertEqual(nump, 1)
-        self.assertTrue(np.all(
-            p_xy == np.array([[0,0,0],
-                              [0,1,0],
-                              [0,1,0],
-                              [0,1,1],
-                              [0,0,0],
-                              [0,0,0]])))
-
-
-    def test_count_segments(self):
-        m_xy = np.zeros([3,3], dtype=bool)
-        m_xy[1,1] = True
-        
-        nump, p_xy = count_segments(m_xy)
-       
-        self.assertEqual(nump, 1)
-        self.assertTrue(np.all(
-            p_xy == np.array([[0,0,0],
-                              [0,1,0],
-                              [0,0,0]])))
-                              
-        m_xy = np.zeros([5,3], dtype=bool)
-        m_xy[1,1] = True
-        m_xy[3,1] = True
-
-        nump, p_xy = count_segments(m_xy)
-
-        self.assertEqual(nump, 2)
-        self.assertTrue(np.all(
-            p_xy == np.array([[0,0,0],
-                              [0,1,0],
-                              [0,0,0],
-                              [0,2,0],
-                              [0,0,0]])))
-
-        m_xy = np.zeros([6,3], dtype=bool)
-        m_xy[1,1] = True
-        m_xy[2,1] = True
-        m_xy[3,1] = True
-        m_xy[3,2] = True
-
-        nump, p_xy = count_segments(m_xy)
-
-        self.assertEqual(nump, 3)
-        self.assertTrue(np.all(
-            p_xy == np.array([[0,0,0],
-                              [0,1,0],
-                              [0,2,0],
-                              [0,3,3],
-                              [0,0,0],
-                              [0,0,0]])))
-
-
-    def test_distance_map(self):
-        m_xy = np.zeros([3,3], dtype=bool)
-        m_xy[1,1] = True
-        
-        d_xy = distance_map(m_xy)
-        
-        sqrt_2 = sqrt(2.0)
-        self.assertTrue(np.all(np.abs(
-            d_xy-np.array([[sqrt_2,1.0,sqrt_2],[1.0,0.0,1.0],[sqrt_2,1.0,sqrt_2]])
-            <1e-12)))
-            
-        m_xy = np.zeros([5,3], dtype=bool)
-        m_xy[1,1] = True
-        m_xy[3,1] = True
-        
-        d_xy = distance_map(m_xy)
-       
-        self.assertTrue(np.all(np.abs(
-            d_xy-np.array([[sqrt_2,1.0,sqrt_2],[1.0,0.0,1.0],
-                           [sqrt_2,1.0,sqrt_2],[1.0,0.0,1.0],
-                           [sqrt_2,1.0,sqrt_2]])
-            <1e-12)))
-
-        m_xy = np.zeros([6,3], dtype=bool)
-        m_xy[1,1] = True
-        m_xy[3,1] = True
-        
-        d_xy = distance_map(m_xy)
-        
-        sqrt_5 = sqrt(5.0)
-        self.assertTrue(np.all(np.abs(
-            d_xy-np.array([[sqrt_2,1.0,sqrt_2],[1.0,0.0,1.0],
-                           [sqrt_2,1.0,sqrt_2],[1.0,0.0,1.0],
-                           [sqrt_2,1.0,sqrt_2],[sqrt_5,2.0,sqrt_5]])
-            <1e-12)))
-            
-            
-    def test_perimeter(self):
-        m_xy = np.zeros([3,3], dtype=bool)
-        m_xy[1,1] = True
-        
-        i_xy = inner_perimeter(m_xy)
-        o_xy = outer_perimeter(m_xy)
-        
-        self.assertTrue(np.array_equal(i_xy, m_xy))
-        self.assertTrue(np.array_equal(o_xy, np.array([[False,True, False],
-                                                       [True, False,True ],
-                                                       [False,True, False]])))
-            
-        m_xy = np.zeros([5,3], dtype=bool)
-        m_xy[1,1] = True
-        m_xy[3,1] = True
-        
-        i_xy = inner_perimeter(m_xy)
-        o_xy = outer_perimeter(m_xy)
-        
-        self.assertTrue(np.array_equal(i_xy, np.array([[False,False,False],
-                                                       [False,True, False],
-                                                       [False,False,False],
-                                                       [False,True, False],
-                                                       [False,False,False]])))
-        self.assertTrue(np.array_equal(o_xy, np.array([[False,True, False],
-                                                       [True, False,True ],
-                                                       [False,True, False],
-                                                       [True, False,True ],
-                                                       [False,True, False]])))
-
-###
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/tests/test_committee.py b/tests/test_committee.py
new file mode 100644
index 00000000..5c681e7d
--- /dev/null
+++ b/tests/test_committee.py
@@ -0,0 +1,219 @@
+import os
+import re
+import copy
+import pathlib
+import pytest
+import numpy as np
+
+import ase.io
+import ase.calculators.emt
+import ase.calculators.lj
+
+import matscipy.calculators.committee
+
+
+@pytest.fixture
+def committeemember():
+    member = matscipy.calculators.committee.CommitteeMember(calculator=ase.calculators.emt.EMT())
+    training_data = os.path.join(f'{os.path.dirname(__file__)}/committee_data/training_data_minimal.xyz')
+    member.set_training_data(training_data)
+    return member
+
+
+@pytest.fixture
+def committee_minimal(committeemember):
+    committee = matscipy.calculators.committee.Committee(
+        members=[copy.deepcopy(committeemember), copy.deepcopy(committeemember)]
+    )
+    return committee
+
+
+@pytest.fixture
+def committee():
+    committee = matscipy.calculators.committee.Committee()
+    committee_training_data = ase.io.read(f'{os.path.dirname(__file__)}/committee_data/training_data.xyz', ':')
+    num_members = 10
+    epsilons = np.linspace(0.98, 1.01, num_members)
+    np.random.seed(123)
+    np.random.shuffle(epsilons)
+    for idx_i in range(num_members):
+        committee += matscipy.calculators.committee.CommitteeMember(
+            calculator=ase.calculators.lj.LennardJones(sigma=1, epsilon=epsilons[idx_i]),
+            training_data=[atoms_i for atoms_i in committee_training_data
+                           if idx_i in atoms_i.info['appears_in_committee']]
+        )
+    return committee
+
+
+@pytest.fixture
+def committee_calibrated(committee):
+    committee.set_internal_validation_set(appearance_threshold=5)
+    committee.calibrate(prop='energy', key='E_lj', location='info')
+    committee.calibrate(prop='forces', key='F_lj', location='arrays')
+    return committee
+
+
+def test_committeemember_initialize():
+    matscipy.calculators.committee.CommitteeMember(calculator=ase.calculators.emt.EMT())
+
+    training_data = os.path.join(f'{os.path.dirname(__file__)}/committee_data/training_data_minimal.xyz')
+    matscipy.calculators.committee.CommitteeMember(calculator=ase.calculators.emt.EMT(),
+                                                   training_data=training_data)
+    matscipy.calculators.committee.CommitteeMember(calculator=ase.calculators.emt.EMT(),
+                                                   training_data=ase.io.read(training_data, ':'))
+
+
+def test_committeemember_set_training_data(committeemember):
+    training_data = os.path.join(f'{os.path.dirname(__file__)}/committee_data/training_data_minimal.xyz')
+
+    with pytest.warns(Warning, match=re.escape('Overwriting current training data.')):
+        committeemember.set_training_data(training_data)
+    with pytest.warns(Warning, match=re.escape('Overwriting current training data.')):
+        committeemember.set_training_data(pathlib.Path(training_data))
+    with pytest.warns(Warning, match=re.escape('Overwriting current training data.')):
+        committeemember.set_training_data(ase.io.read(training_data, ':'))
+
+
+def test_committeemember_is_sample_in_atoms(committeemember):
+    training_data = ase.io.read(os.path.join(f'{os.path.dirname(__file__)}/committee_data/training_data_minimal.xyz'), ':')
+    test_data = ase.io.read(os.path.join(f'{os.path.dirname(__file__)}/committee_data/test_data.xyz'), ':')
+
+    assert committeemember.is_sample_in_atoms(sample=training_data[0])
+    with pytest.raises(RuntimeError,
+                       match=re.escape('Can\'t test if `sample` is in `atoms`. '
+                                       '`sample` has no Atoms.info[\'_Index_FullTrainingSet\']')):
+        assert not committeemember.is_sample_in_atoms(sample=test_data[0])
+    test_data[0].info['_Index_FullTrainingSet'] = -1
+    assert not committeemember.is_sample_in_atoms(sample=test_data[0])
+
+
+def test_committeemember_setter(committeemember):
+
+    with pytest.raises(RuntimeError, match=re.escape('Use `set_training_data()` to modify the committee member')):
+        committeemember.filename = ''
+
+    with pytest.raises(RuntimeError, match=re.escape('Use `set_training_data()` to modify the committee member')):
+        committeemember.atoms = []
+
+    with pytest.raises(RuntimeError, match=re.escape('Use `set_training_data()` to modify the committee member')):
+        committeemember.ids = []
+
+
+def test_committee_initialize(committeemember):
+    committee = matscipy.calculators.committee.Committee()
+    expected_status = [
+        ('members', []),
+        ('number', 0),
+        ('atoms', []),
+        ('ids', []),
+        ('alphas', {}),
+        ('calibrated_for', set()),
+    ]
+    for attribute_i, value_i in expected_status:
+        assert getattr(committee, attribute_i) == value_i
+    with pytest.warns(Warning, match=re.escape('`Committee.set_internal_validation_set()` has not been called or '
+                                               '`Committee`-instance has been altered since last call.')):
+        assert getattr(committee, 'validation_set') == []
+
+    member_0 = copy.deepcopy(committeemember)
+    member_1 = copy.deepcopy(committeemember)
+    committee = matscipy.calculators.committee.Committee(
+        members=[member_0, member_1]
+    )
+    expected_status = [
+        ('members', [member_0, member_1]),
+        ('number', 2),
+        ('atoms', member_0.atoms + member_1.atoms),
+        ('ids', member_0.ids + member_1.ids),
+        ('alphas', {}),
+        ('calibrated_for', set()),
+    ]
+    for attribute_i, value_i in expected_status:
+        assert getattr(committee, attribute_i) == value_i
+    with pytest.warns(Warning, match=re.escape('`Committee.set_internal_validation_set()` has not been called or '
+                                               '`Committee`-instance has been altered since last call.')):
+        assert getattr(committee, 'validation_set') == []
+
+
+def test_committee_member(committee_minimal):
+
+    with pytest.raises(AssertionError,
+                       match=re.escape('Members of `Committee` need to be of type `CommitteeMember`. Found ')):
+        matscipy.calculators.committee.Committee(members=[0, 1])
+
+    with pytest.raises(AssertionError,
+                       match=re.escape('Members of `Committee` need to be of type `CommitteeMember`. Found ')):
+        committee_minimal.members = [0, 1]
+
+    with pytest.raises(AssertionError,
+                       match=re.escape('Members of `Committee` need to be of type `CommitteeMember`. Found ')):
+        committee_minimal.add_member(0)
+
+    with pytest.raises(AssertionError,
+                       match=re.escape('Members of `Committee` need to be of type `CommitteeMember`. Found ')):
+        committee_minimal += 0
+
+
+def test_committee_set_internal_validation_set(committee):
+
+    with pytest.raises(AssertionError):
+        committee.set_internal_validation_set(appearance_threshold=0)
+
+    with pytest.raises(AssertionError):
+        committee.set_internal_validation_set(appearance_threshold=committee.number - 1)
+
+    committee.set_internal_validation_set(appearance_threshold=5)
+    obtained = set([atoms_i.info['_Index_FullTrainingSet'] for atoms_i
+                    in committee.validation_set])
+    expected = set([atoms_i.info['_Index_FullTrainingSet'] for atoms_i
+                    in ase.io.read(os.path.join(f'{os.path.dirname(__file__)}/committee_data/validation_set.xyz'), ':')])
+    assert obtained == expected
+
+
+def test_committee_calibrate(committee):
+    committee.set_internal_validation_set(appearance_threshold=5)
+
+    committee.calibrate(prop='energy', key='E_lj', location='info')
+    assert committee.calibrated_for == set(['energy'])
+    np.testing.assert_array_almost_equal(committee.alphas['energy'], 0.6295416920992463, decimal=6)
+
+    committee.calibrate(prop='forces', key='F_lj', location='arrays')
+    assert committee.calibrated_for == set(['energy', 'forces'])
+    np.testing.assert_array_almost_equal(committee.alphas['forces'], 0.6195847443699875, decimal=6)
+
+    with pytest.warns(Warning,
+                      match=re.escape('`alphas` will be reset to avoid inconsistencies with new validation set.')):
+        committee.set_internal_validation_set(appearance_threshold=4)
+        assert committee.alphas == {}
+
+
+def test_committee__calculate_alpha(committee):
+    vals_ref = np.array([1.01, 1.02, 1.03])
+    vals_pred = np.array([2.01, 1.02, 1.03])
+    var_pred = np.array([1.01, 0.02, 0.03])
+
+    obtained = committee._calculate_alpha(vals_ref, vals_pred, var_pred)
+    np.testing.assert_array_almost_equal(obtained, 0.39584382766472004, decimal=6)
+
+
+def test_committee_scale_uncertainty(committee):
+    committee._alphas = {'energy': 2.5}
+
+    assert committee.scale_uncertainty(2, 'energy') == 5.0
+
+
+def test_committeeuncertainty_initialize(committee_calibrated):
+    matscipy.calculators.committee.CommitteeUncertainty(committee=committee_calibrated)
+
+
+def test_committeeuncertainty_calculate(committee_calibrated):
+    calculator = matscipy.calculators.committee.CommitteeUncertainty(committee=committee_calibrated)
+    test_data = ase.io.read(os.path.join(f'{os.path.dirname(__file__)}/committee_data/test_data.xyz'), ':')
+    for atoms_i in test_data:
+        calculator.calculate(atoms=atoms_i, properties=['energy', 'forces'])
+        for prop_j in ['energy', 'energy_uncertainty']:
+            np.testing.assert_array_almost_equal(calculator.results[prop_j], atoms_i.info[prop_j], decimal=6,
+                                                 err_msg=f'Missmatch in property \'{prop_j}\'')
+        for prop_j in ['forces', 'forces_uncertainty']:
+            np.testing.assert_array_almost_equal(calculator.results[prop_j], atoms_i.arrays[prop_j], decimal=6,
+                                                 err_msg=f'Missmatch in property \'{prop_j}\'')
diff --git a/tests/test_dislocation.py b/tests/test_dislocation.py
index 7bce5203..38e80823 100644
--- a/tests/test_dislocation.py
+++ b/tests/test_dislocation.py
@@ -28,11 +28,13 @@
 import numpy as np
 
 from ase.calculators.lammpslib import LAMMPSlib
-from scipy.optimize import minimize
 from matscipy.calculators.eam import EAM
 
+test_dir = os.path.dirname(os.path.realpath(__file__))
+
 try:
     import matplotlib
+    matplotlib.use("Agg")  # Activate 'agg' backend for off-screen plotting for testing.
 except ImportError:
     print("matplotlib not found: skipping some tests")
 
@@ -55,7 +57,8 @@
 class TestDislocation(matscipytest.MatSciPyTestCase):
     """Class to store test for dislocation.py module."""
 
-    @unittest.skipIf("atomman" not in sys.modules, 'requires Stroh solution from atomman to run')
+    @unittest.skipIf("atomman" not in sys.modules,
+                     'requires Stroh solution from atomman to run')
     def test_core_position(self):
         """Make screw dislocation and fit a core position to it.
 
@@ -91,20 +94,23 @@ def test_core_position(self):
         self.assertArrayAlmostEqual(core_pos, initial_guess + center[:2], tol=1e-4)
 
     def test_elastic_constants_EAM(self):
-        """Test the get_elastic_constants() function using matscipy EAM calculator."""
-        target_values = np.array([3.14339177996466,  # alat
-                                  523.0266819809012,  # C11
-                                  202.1786296941397,  # C12
-                                  160.88179872237012])  # C44 for eam4
+        """Test the get_elastic_constants()
+           function using matscipy EAM calculator."""
+        target_values = np.array([3.1433,  # alat
+                                  523.03,  # C11
+                                  202.18,  # C12
+                                  160.88])  # C44 for eam4
 
         pot_name = "w_eam4.fs"
-        calc_EAM = EAM(pot_name)
+        pot_path = os.path.join(test_dir, pot_name)
+        calc_EAM = EAM(pot_path)
         obtained_values = sd.get_elastic_constants(calculator=calc_EAM,
                                                    delta=1.0e-3)
 
-        self.assertArrayAlmostEqual(obtained_values, target_values, tol=1e-4)
+        self.assertArrayAlmostEqual(obtained_values, target_values, tol=1e-2)
 
-    # This function tests the lammpslib and LAMMPS installation and thus skipped during automated testing
+    # This function tests the lammpslib and LAMMPS installation
+    # skipped during automated testing
     @unittest.skipIf("lammps" not in sys.modules,
                      "LAMMPS installation is required and thus is not good for automated testing")
     def test_elastic_constants_lammpslib(self):
@@ -115,15 +121,15 @@ def test_elastic_constants_lammpslib(self):
         """
         print("WARNING: In case lammps crashes no error message is printed: ",
               "check 'lammps.log' file in test folder")
-        target_values = np.array([3.14339177996466,  # alat
-                                  523.0266819809012,  # C11
-                                  202.1786296941397,  # C12
-                                  160.88179872237012])  # C44 for eam4
+        target_values = np.array([3.1434,  # alat
+                                  523.03,  # C11
+                                  202.18,  # C12
+                                  160.882])  # C44 for eam4
 
         pot_name = "w_eam4.fs"
-
+        pot_path = os.path.join(test_dir, pot_name)
         lammps = LAMMPSlib(lmpcmds=["pair_style eam/fs",
-                                    "pair_coeff * * %s W" % pot_name],
+                                    "pair_coeff * * %s W" % pot_path],
                            atom_types={'W': 1}, keep_alive=True,
                            log_file="lammps.log")
 
@@ -132,9 +138,10 @@ def test_elastic_constants_lammpslib(self):
 
         os.remove("lammps.log")
 
-        self.assertArrayAlmostEqual(obtained_values, target_values, tol=1e-4)
+        self.assertArrayAlmostEqual(obtained_values, target_values, tol=1e-2)
 
-    # This function tests the lammpslib and LAMMPS installation and thus skipped during automated testing
+    # This function tests the lammpslib and LAMMPS installation
+    # skipped during automated testing
     @unittest.skipIf("lammps" not in sys.modules or
                      "atomman" not in sys.modules,
                      "LAMMPS installation and Stroh solution are required")
@@ -158,10 +165,11 @@ def test_screw_cyl_lammpslib(self):
         center = (cent_x, 0.0, 0.0)
 
         pot_name = "w_eam4.fs"
+        pot_path = os.path.join(test_dir, pot_name)
         target_toten = -13086.484626  # Target value for w_eam4
 
         lammps = LAMMPSlib(lmpcmds=["pair_style eam/fs",
-                                    "pair_coeff * * %s W" % pot_name],
+                                    "pair_coeff * * %s W" % pot_path],
                            atom_types={'W': 1}, keep_alive=True,
                            log_file="lammps.log")
 
@@ -211,7 +219,8 @@ def test_differential_displacement(self):
         print("'dd_test.png' will be created: check the displacement map")
         fig.savefig("dd_test.png")
 
-    @unittest.skipIf("atomman" not in sys.modules, 'requires Stroh solution from atomman to run')
+    @unittest.skipIf("atomman" not in sys.modules,
+                     'requires Stroh solution from atomman to run')
     def test_read_dislo_QMMM(self):
         """Test read_dislo_QMMM() function"""
 
@@ -254,13 +263,17 @@ def test_read_dislo_QMMM(self):
             total_Nat_type += Nat_type
             self.assertEqual(Nat_type, target_values[atom_type])
 
-        self.assertEqual(Nat, total_Nat_type)  # total number of atoms in region is equal to Nat (no atoms unmapped)
+        # total number of atoms in region is equal to Nat (no atoms unmapped)
+        self.assertEqual(Nat, total_Nat_type)
         # TODO
-        #  self.assertAtomsAlmostEqual(disloc, test_disloc) - gives an error of _cell attribute new ase version?
+        # self.assertAtomsAlmostEqual(disloc, test_disloc) -
+        # gives an error of _cell attribute new ase version?
 
-    @unittest.skipIf("atomman" not in sys.modules, 'requires Stroh solution from atomman to run')
+    @unittest.skipIf("atomman" not in sys.modules,
+                     'requires Stroh solution from atomman to run')
     def test_stroh_solution(self):
-        """Builds isotropic Stroh solution and compares it to Volterra solution"""
+        """Builds isotropic Stroh solution
+        and compares it to Volterra solution"""
 
         alat = 3.14
         C11 = 523.0
@@ -269,7 +282,8 @@ def test_stroh_solution(self):
 
         # A = 2. * C44 / (C11 - C12)
         # print(A) # A = 0.999937 very isotropic material.
-        # At values closer to 1.0 Stroh solution is numerically unstable and does not pass checks
+        # At values closer to 1.0 Stroh solution is numerically unstable
+        # and does not pass checks
         cylinder_r = 40
         burgers = alat * np.sqrt(3.0) / 2.
 
@@ -282,26 +296,32 @@ def test_stroh_solution(self):
         u_volterra = np.arctan2(y, x) * burgers / (2.0 * np.pi)
 
         # compare x and y components with zeros - isotropic solution
-        self.assertArrayAlmostEqual(np.zeros_like(u_volterra), u_stroh[:, 0], tol=1e-4)
-        self.assertArrayAlmostEqual(np.zeros_like(u_volterra), u_stroh[:, 1], tol=1e-4)
+        self.assertArrayAlmostEqual(np.zeros_like(u_volterra), u_stroh[:, 0],
+                                    tol=1e-4)
+        self.assertArrayAlmostEqual(np.zeros_like(u_volterra), u_stroh[:, 1],
+                                    tol=1e-4)
         #  compare z component with simple Volterra solution
         self.assertArrayAlmostEqual(u_volterra, u_stroh[:, 2])
 
     def test_make_screw_quadrupole_kink(self):
-        """Test the total number of atoms in the quadrupole double kink configuration"""
+        """Test the total number of atoms in the quadrupole
+            double kink configuration"""
 
         alat = 3.14
         n1u = 5
         kink_length = 20
 
-        kink, _, _ = sd.make_screw_quadrupole_kink(alat=alat, n1u=n1u, kink_length=kink_length)
+        kink, _, _ = sd.make_screw_quadrupole_kink(alat=alat, n1u=n1u,
+                                                   kink_length=kink_length)
         quadrupole_base, _, _, _ = sd.make_screw_quadrupole(alat=alat, n1u=n1u)
 
         self.assertEqual(len(kink), len(quadrupole_base) * 2 * kink_length)
 
-    @unittest.skipIf("atomman" not in sys.modules, 'requires Stroh solution from atomman to run')
+    @unittest.skipIf("atomman" not in sys.modules,
+                     'requires Stroh solution from atomman to run')
     def test_make_screw_cyl_kink(self):
-        """Test the total number of atoms and number of fixed atoms in the cylinder double kink configuration"""
+        """Test the total number of atoms and number of fixed atoms
+            in the cylinder double kink configuration"""
 
         alat = 3.14339177996466
         C11 = 523.0266819809012
@@ -314,11 +334,19 @@ def test_make_screw_cyl_kink(self):
         cylinder_r = 40
         kink_length = 26
 
-        kink, large_disloc, straight_bulk = sd.make_screw_cyl_kink(alat, C11, C12, C44, kink_length=kink_length,
-                                                                   cylinder_r=cylinder_r, kind="double")
+        (kink,
+         large_disloc, _) = sd.make_screw_cyl_kink(alat,
+                                                   C11,
+                                                   C12,
+                                                   C44,
+                                                   kink_length=kink_length,
+                                                   cylinder_r=cylinder_r,
+                                                   kind="double")
 
         # check the total number of atoms as compared to make_screw_cyl()
-        disloc, _, _ = sd.make_screw_cyl(alat, C11, C12, C12, cylinder_r=cylinder_r, l_extend=center)
+        disloc, _, _ = sd.make_screw_cyl(alat, C11, C12, C12,
+                                         cylinder_r=cylinder_r,
+                                         l_extend=center)
 
         self.assertEqual(len(kink), len(disloc) * 2 * kink_length)
 
@@ -340,12 +368,15 @@ def test_slice_long_dislo(self):
         n1u = 5
         kink_length = 20
 
-        kink, straight_dislo, kink_bulk = sd.make_screw_quadrupole_kink(alat=alat, n1u=n1u, kink_length=kink_length)
+        kink, _, kink_bulk = sd.make_screw_quadrupole_kink(alat=alat,
+                                                           n1u=n1u,
+                                                           kink_length=kink_length)
         quadrupole_base, _, _, _ = sd.make_screw_quadrupole(alat=alat, n1u=n1u)
 
         sliced_kink, core_positions = sd.slice_long_dislo(kink, kink_bulk, b)
 
-        # check the number of sliced configurations is equal to length of 2 * kink_length * 3 (for double kink)
+        # check the number of sliced configurations is equal
+        # to length of 2 * kink_length * 3 (for double kink)
         self.assertEqual(len(sliced_kink), kink_length * 3 * 2)
 
         # check that the bulk and kink slices are the same size
@@ -356,16 +387,26 @@ def test_slice_long_dislo(self):
         # check that the size of slices are the same as single b configuration
         self.assertArrayAlmostEqual(len(quadrupole_base), len(sliced_kink[0][0]))
 
-        right_kink, straight_dislo, kink_bulk = sd.make_screw_quadrupole_kink(alat=alat, n1u=n1u, kind="right",
-                                                                              kink_length=kink_length)
+        (right_kink, _,
+         kink_bulk) = sd.make_screw_quadrupole_kink(alat=alat,
+                                                    n1u=n1u,
+                                                    kind="right",
+                                                    kink_length=kink_length)
+
         sliced_right_kink, _ = sd.slice_long_dislo(right_kink, kink_bulk, b)
-        # check the number of sliced configurations is equal to length of kink_length * 3 - 2 (for right kink)
+        # check the number of sliced configurations is equal
+        # to length of kink_length * 3 - 2 (for right kink)
         self.assertEqual(len(sliced_right_kink), kink_length * 3 - 2)
 
-        left_kink, straight_dislo, kink_bulk = sd.make_screw_quadrupole_kink(alat=alat, n1u=n1u, kind="left",
-                                                                              kink_length=kink_length)
+        (left_kink, _,
+         kink_bulk) = sd.make_screw_quadrupole_kink(alat=alat,
+                                                    n1u=n1u,
+                                                    kind="left",
+                                                    kink_length=kink_length)
+
         sliced_left_kink, _ = sd.slice_long_dislo(left_kink, kink_bulk, b)
-        # check the number of sliced configurations is equal to length of kink_length * 3 - 1 (for left kink)
+        # check the number of sliced configurations is equal
+        # to length of kink_length * 3 - 1 (for left kink)
         self.assertEqual(len(sliced_left_kink), kink_length * 3 - 1)
 
     def check_disloc(self, cls, ref_angle, structure="BCC", test_u=True,
@@ -375,8 +416,13 @@ def check_disloc(self, cls, ref_angle, structure="BCC", test_u=True,
         C12 = 202.1786296941397
         C44 = 160.88179872237012
 
-        d = cls(alat, C11, C12, C44)
+        d = cls(alat, C11, C12, C44, symbol="Al")
         bulk, disloc = d.build_cylinder(20.0)
+
+        # test that assigning non default symbol worked
+        assert np.unique(bulk.get_chemical_symbols()) == "Al"
+        assert np.unique(disloc.get_chemical_symbols()) == "Al"
+
         assert len(bulk) == len(disloc)
 
         if test_u:
@@ -388,7 +434,8 @@ def check_disloc(self, cls, ref_angle, structure="BCC", test_u=True,
 
             displacement = disloc.positions - bulk.positions
 
-            np.testing.assert_array_almost_equal(displacement, stroh_displacement)
+            np.testing.assert_array_almost_equal(displacement,
+                                                 stroh_displacement)
 
         results = sd.ovito_dxa_straight_dislo_info(disloc, structure=structure)
         assert len(results) == 1
@@ -399,16 +446,16 @@ def check_disloc(self, cls, ref_angle, structure="BCC", test_u=True,
         print(f'angle = {angle} ref_angle = {ref_angle} err = {err}')
         assert abs(err) < tol
 
-    @unittest.skipIf("atomman" not in sys.modules or 
+    @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_screw_dislocation(self):
         self.check_disloc(sd.BCCScrew111Dislocation, 0.0)
 
-    @unittest.skipIf("atomman" not in sys.modules or 
+    @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
-    def test_edge_dislocation(self):        
+    def test_edge_dislocation(self):
         self.check_disloc(sd.BCCEdge111Dislocation, 90.0)
 
     @unittest.skipIf("atomman" not in sys.modules or
@@ -425,7 +472,7 @@ def test_edge100110_dislocation(self,):
         self.check_disloc(sd.BCCEdge100110Dislocation, 90.0,
                           burgers=np.array([1.0, 0.0, 0.0]))
 
-    @unittest.skipIf("atomman" not in sys.modules or 
+    @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_mixed_dislocation(self):
@@ -447,24 +494,54 @@ def test_90degree_diamond_partial_dislocation(self,):
                           structure="Diamond",
                           burgers=(1.0 / 6.0) * np.array([2.0, 1.0, 1.0]))
 
-
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_screw_diamond_dislocation(self,):
         self.check_disloc(sd.DiamondGlideScrew, 0.0,
-                          structure="Diamond", test_u=False,
+                          structure="Diamond",
                           burgers=(1.0 / 2.0) * np.array([0.0, 1.0, 1.0]))
 
-
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_60degree_diamond_dislocation(self,):
         self.check_disloc(sd.DiamondGlide60Degree, 60.0,
-                          structure="Diamond", test_u=False,
+                          structure="Diamond",
                           burgers=(1.0 / 2.0) * np.array([1.0, 0.0, 1.0]))
 
+    @unittest.skipIf("atomman" not in sys.modules or
+                     "ovito" not in sys.modules,
+                     "requires atomman and ovito")
+    def test_fcc_screw_shockley_partial_dislocation(self,):
+        self.check_disloc(sd.FCCScrewShockleyPartial, 30.0,
+                          structure="FCC",
+                          burgers=(1.0 / 6.0) * np.array([1.0, 2.0, 1.0]))
+
+    @unittest.skipIf("atomman" not in sys.modules or
+                     "ovito" not in sys.modules,
+                     "requires atomman and ovito")
+    def test_fcc_screw110_dislocation(self,):
+        self.check_disloc(sd.FCCScrew110Dislocation, 0.0,
+                          structure="FCC",
+                          burgers=(1.0 / 2.0) * np.array([0.0, 1.0, 1.0]))
+
+    @unittest.skipIf("atomman" not in sys.modules or
+                     "ovito" not in sys.modules,
+                     "requires atomman and ovito")
+    def test_fcc_edge_shockley_partial_dislocation(self,):
+        self.check_disloc(sd.FCCEdgeShockleyPartial, 60.0,
+                          structure="FCC",
+                          burgers=(1.0 / 6.0) * np.array([1.0, 2.0, 1.0]))
+
+    @unittest.skipIf("atomman" not in sys.modules or
+                     "ovito" not in sys.modules,
+                     "requires atomman and ovito")
+    def test_fcc_edge110_partial_dislocation(self,):
+        self.check_disloc(sd.FCCEdge110Dislocation, 90.0,
+                          structure="FCC",
+                          burgers=(1.0 / 2.0) * np.array([1.0, 1.0, 0.0]))
+
     def check_glide_configs(self, cls, structure="BCC"):
         alat = 3.14339177996466
         C11 = 523.0266819809012
@@ -492,8 +569,8 @@ def check_glide_configs(self, cls, structure="BCC"):
         # test that difference between initial and final positions are
         # roughly equal to glide distance.
         # Since the configurations are unrelaxed dxa gives
-        # a very rough estimation (especially for edge dislcoations)
-        # thus tolerance is taken to be ~1 Angstom
+        # a very rough estimation (especially for edge dislocations)
+        # thus tolerance is taken to be ~1 Angstrom
         np.testing.assert_almost_equal(fin_x_position - ini_x_position,
                                        d.glide_distance, decimal=0)
 
@@ -521,51 +598,75 @@ def test_mixed_glide(self):
     def test_edge100_glide(self):
         self.check_glide_configs(sd.BCCEdge100Dislocation)
 
-
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_edge100110_glide(self):
-            self.check_glide_configs(sd.BCCEdge100110Dislocation)
-
+        self.check_glide_configs(sd.BCCEdge100110Dislocation)
 
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_30degree_diamond_partial_glide(self):
-            self.check_glide_configs(sd.DiamondGlide30degreePartial,
-                                     structure="Diamond")
-
+        self.check_glide_configs(sd.DiamondGlide30degreePartial,
+                                 structure="Diamond")
 
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_90degree_diamond_partial_glide(self):
-            self.check_glide_configs(sd.DiamondGlide90degreePartial,
-                                     structure="Diamond")
+        self.check_glide_configs(sd.DiamondGlide90degreePartial,
+                                 structure="Diamond")
+
+    @unittest.skipIf("atomman" not in sys.modules or
+                     "ovito" not in sys.modules,
+                     "requires atomman and ovito")
+    def test_screw_diamond_glide(self):
+        self.check_glide_configs(sd.DiamondGlideScrew,
+                                 structure="Diamond")
 
+    @unittest.skipIf("atomman" not in sys.modules or
+                     "ovito" not in sys.modules,
+                     "requires atomman and ovito")
+    def test_60degree_diamond_glide(self):
+        self.check_glide_configs(sd.DiamondGlide60Degree,
+                                 structure="Diamond")
 
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
-    def test_screw_diamond_partial_glide(self):
-            self.check_glide_configs(sd.DiamondGlideScrew,
-                                     structure="Diamond")
+    def test_fcc_screw_shockley_partial_glide(self):
+        self.check_glide_configs(sd.FCCScrewShockleyPartial,
+                                 structure="FCC")
 
+    @unittest.skipIf("atomman" not in sys.modules or
+                     "ovito" not in sys.modules,
+                     "requires atomman and ovito")
+    def test_fcc_screw_110_glide(self):
+        self.check_glide_configs(sd.FCCScrew110Dislocation,
+                                 structure="FCC")
 
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
-    def test_60degree_diamond_partial_glide(self):
-            self.check_glide_configs(sd.DiamondGlide60Degree,
-                                     structure="Diamond")
+    def test_fcc_edge_shockley_partial_glide(self):
+        self.check_glide_configs(sd.FCCEdgeShockleyPartial,
+                                 structure="FCC")
+
+    @unittest.skipIf("atomman" not in sys.modules or
+                     "ovito" not in sys.modules,
+                     "requires atomman and ovito")
+    def test_fcc_edge_dislocation_glide(self):
+        self.check_glide_configs(sd.FCCEdge110Dislocation,
+                                 structure="FCC")
 
     def test_fixed_line_atoms(self):
 
         from ase.build import bulk
 
         pot_name = "w_eam4.fs"
-        calc_EAM = EAM(pot_name)
+        pot_path = os.path.join(test_dir, pot_name)
+        calc_EAM = EAM(pot_path)
         # slightly pressurised cell to avoid exactly zero forces
         W = bulk("W", a=0.9 * 3.143392, cubic=True)
 
@@ -591,9 +692,9 @@ def test_fixed_line_atoms(self):
             # forces on unconstrained atoms are non zero
             assert (W.get_forces()[~fixed_mask] != 0.0).all()
 
-
     @unittest.skipIf("lammps" not in sys.modules,
-                     "LAMMPS installation is required and thus is not good for automated testing")
+                     "LAMMPS installation is required and " +
+                     "thus is not good for automated testing")
     def test_gamma_line(self):
 
         # eam_4 parameters
@@ -602,9 +703,10 @@ def test_gamma_line(self):
         unit_cell = dislocation.unit_cell
 
         pot_name = "w_eam4.fs"
+        pot_path = os.path.join(test_dir, pot_name)
         # calc_EAM = EAM(pot_name) eam calculator is way too slow
         lammps = LAMMPSlib(lmpcmds=["pair_style eam/fs",
-                                    "pair_coeff * * %s W" % pot_name],
+                                    "pair_coeff * * %s W" % pot_path],
                            atom_types={'W': 1}, keep_alive=True,
                            log_file="lammps.log")
 
@@ -629,5 +731,6 @@ def test_gamma_line(self):
         np.testing.assert_almost_equal(E, target_E, decimal=3)
         np.testing.assert_almost_equal(shift, target_shift, decimal=3)
 
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_elastic_moduli.py b/tests/test_elastic_moduli.py
index 47444d89..37fc0b46 100644
--- a/tests/test_elastic_moduli.py
+++ b/tests/test_elastic_moduli.py
@@ -83,6 +83,24 @@ def test_rotation():
             np.testing.assert_array_almost_equal(G, EM[i]['G'], decimal=9)
 
 
+def test_monoclinic():
+    C = np.array([
+            [5,2,2,0,1,0],
+            [2,5,2,0,1,0],
+            [2,2,5,0,1,0],
+            [0,0,0,2,0,1],
+            [1,1,1,0,2,0],
+            [0,0,0,1,0,2]
+        ])
+    l = [np.array([1, 0, 1]), np.array([1, 0, -1])]
+    EM = np.array([5.4545,3.1579])
+
+    for i, directions in enumerate(l):
+        directions = directions / np.linalg.norm(directions)
+        E, nu, Gm, B, K = elastic_moduli(C, l=directions)
+        np.testing.assert_almost_equal(E[0], EM[i], decimal=1)
+
+
 ###
 
 if __name__ == '__main__':
diff --git a/tests/test_greens_function.py b/tests/test_greens_function.py
deleted file mode 100644
index 87e46f53..00000000
--- a/tests/test_greens_function.py
+++ /dev/null
@@ -1,192 +0,0 @@
-#
-# Copyright 2014-2015, 2020-2021 Lars Pastewka (U. Freiburg)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-# ======================================================================
-# matscipy - Python materials science tools
-# https://github.com/libAtoms/matscipy
-#
-# Copyright (2014) James Kermode, King's College London
-#                  Lars Pastewka, Karlsruhe Institute of Technology
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-# ======================================================================
-
-from math import pi, sqrt
-
-import unittest
-
-import numpy as np
-
-import matscipytest
-import matscipy.contact_mechanics.greens_function as gf
-import matscipy.contact_mechanics.Hertz as Hertz
-
-###
-
-class TestGreensFunction(matscipytest.MatSciPyTestCase):
-
-    def test_Hertz_displacements_square(self):
-        nx = 256 # Grid size
-        a = 32. # Contact radius
-        G, x, y = gf.real_to_reciprocal_space(
-            nx, nx, gf=gf.square_pressure__nonperiodic,
-            coordinates=True)
-
-        r_sq = (x**2 + y**2)/a**2
-        P = np.where(r_sq > 1., np.zeros_like(r_sq), np.sqrt(1.-r_sq))
-        u = -np.fft.ifft2(G*np.fft.fft2(P)).real
-        # Note: contact modulus Es = 2 for the GF
-        uref = a/2. * Hertz.surface_displacements(np.sqrt(r_sq))
-
-        # Only the middle section of displacements is correct in for the
-        # free-space Green's function.
-        u = u[:nx//4, :nx//4]
-        uref = uref[:nx//4, :nx//4]
-        #np.savetxt('u.out', u)
-        #np.savetxt('u0.out', np.transpose([y[0,:nx//4], u[0,:]]))
-        #np.savetxt('uref.out', uref)
-        #np.savetxt('uref0.out', np.transpose([y[0,:nx//4], uref[0,:]]))
-        #np.savetxt('u_uref.out', u/uref)
-        self.assertTrue(np.max(np.abs((u-uref)/uref)) < 0.01)
-
-    def test_Hertz_stress(self):
-        nx = 256 # Grid size
-        a = 32. # Contact radius
-        tol = 1e-2
-
-        # z: Depth at which to compute stress
-        # nu: Poisson
-        for z, nu in [ ( 16., 0.5 ), ( 32., 0.5 ), ( 16., 0.3 ), ( a/8, 1/3) ]:
-            ( Gxx, Gyy, Gzz, Gyz, Gxz, Gxy ), x, y = \
-                gf.real_to_reciprocal_space(
-                    nx, nx, gf=lambda x, y:
-                    gf.point_traction__nonperiodic('Z', x, y, z, poisson=nu),
-                    coordinates=True)
-
-            r_sq = (x**2 + y**2)/a**2
-            P = np.where(r_sq > 1., np.zeros_like(r_sq), np.sqrt(1.-r_sq))
-            sxx = np.fft.ifft2(Gxx*np.fft.fft2(P)).real
-            syy = np.fft.ifft2(Gyy*np.fft.fft2(P)).real
-            szz = np.fft.ifft2(Gzz*np.fft.fft2(P)).real
-            syz = np.fft.ifft2(Gyz*np.fft.fft2(P)).real
-            sxz = np.fft.ifft2(Gxz*np.fft.fft2(P)).real
-            sxy = np.fft.ifft2(Gxy*np.fft.fft2(P)).real
-
-            x, y, z = np.meshgrid(x, y, z, indexing='ij')
-
-            sxx2, syy2, szz2, syz2, sxz2, sxy2 = \
-                Hertz.stress_Cartesian(x/a, y/a, z/a, poisson=nu)
-
-            m = np.abs(szz2[:nx//4,:nx//4,0]).max()
-            self.assertArrayAlmostEqual(sxx[:nx//4,:nx//4]/m, sxx2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(syy[:nx//4,:nx//4]/m, syy2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(szz[:nx//4,:nx//4]/m, szz2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(syz[:nx//4,:nx//4]/m, syz2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(sxz[:nx//4,:nx//4]/m, sxz2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(sxy[:nx//4,:nx//4]/m, sxy2[:nx//4,:nx//4,0]/m, tol)
-
-    def test_Hertz_stress_tangential(self):
-        nx = 256 # Grid size
-        a = 32. # Contact radius
-        tol = 1e-2
-
-        # z: Depth at which to compute stress
-        # nu: Poisson
-        for z, nu in [ ( 16., 0.5 ), ( 32., 0.5 ), ( 16., 0.3 ), ( a/6, 1/3) ]:
-            ( Gxx, Gyy, Gzz, Gyz, Gxz, Gxy ), x, y = \
-                gf.real_to_reciprocal_space(
-                    nx, nx, gf=lambda x, y:
-                    gf.point_traction__nonperiodic('X', x, y, z, poisson=nu),
-                    coordinates=True)
-
-            r_sq = (x**2 + y**2)/a**2
-            P = np.where(r_sq > 1., np.zeros_like(r_sq), np.sqrt(1.-r_sq))
-            sxx = np.fft.ifft2(Gxx*np.fft.fft2(P)).real
-            syy = np.fft.ifft2(Gyy*np.fft.fft2(P)).real
-            szz = np.fft.ifft2(Gzz*np.fft.fft2(P)).real
-            syz = np.fft.ifft2(Gyz*np.fft.fft2(P)).real
-            sxz = np.fft.ifft2(Gxz*np.fft.fft2(P)).real
-            sxy = np.fft.ifft2(Gxy*np.fft.fft2(P)).real
-
-            x, y, z = np.meshgrid(x, y, z, indexing='ij')
-
-            sxx2, syy2, szz2, syz2, sxz2, sxy2 = \
-                Hertz.stress_for_tangential_loading(x/a, y/a, z/a, poisson=nu)
-
-            # DEBUG CODE
-            #plt.figure(figsize=[15,4])
-            #plt.subplot(1,3,1)
-            #plt.pcolormesh(np.fft.fftshift(sxx)[nx//4:3*nx//4,nx//4:3*nx//4])
-            #plt.colorbar()
-            #plt.subplot(1,3,2)
-            #plt.pcolormesh(np.fft.fftshift(sxx2)[nx//4:3*nx//4,nx//4:3*nx//4,0])
-            #plt.colorbar()
-            #plt.subplot(1,3,3)
-            #plt.pcolormesh(np.fft.fftshift(sxx-sxx2[:,:,0])[nx//4:3*nx//4,nx//4:3*nx//4])
-            #plt.colorbar()
-            #plt.show()
-
-            m = np.abs(szz2[:nx//4,:nx//4,0]).max()
-            self.assertArrayAlmostEqual(sxx[:nx//4,:nx//4]/m, sxx2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(syy[:nx//4,:nx//4]/m, syy2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(szz[:nx//4,:nx//4]/m, szz2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(syz[:nx//4,:nx//4]/m, syz2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(sxz[:nx//4,:nx//4]/m, sxz2[:nx//4,:nx//4,0]/m, tol)
-            self.assertArrayAlmostEqual(sxy[:nx//4,:nx//4]/m, sxy2[:nx//4,:nx//4,0]/m, tol)
-
-    def test_min_ccg(self):
-        nx = 256 # Grid size
-        R = 256. # sphere radius
-        d = R/20. # penetration depth
-        G, x, y = gf.real_to_reciprocal_space(
-            nx, nx, gf=gf.square_pressure__nonperiodic,
-            coordinates=True)
-
-        r_sq = (x**2 + y**2)/R**2
-        h = np.where(r_sq > 1., (R-d)*np.ones_like(r_sq), (R-d)-R*np.sqrt(1.-r_sq))
-
-        # Contact radius
-        a = R*sqrt(d/R)
-        p0 = 2/pi*sqrt(d/R)
-        r_sq = (x**2 + y**2)/a**2
-        p_analytic, sr, stheta = Hertz.surface_stress(np.sqrt(r_sq), poisson=0.5)
-        p_analytic *= p0
-
-        u, p = gf.min_ccg(h, G)
-        p /= 2 # E* of GF is 2
-
-        assert (p-p_analytic).max() < 1e-3
-
-###
-
-if __name__ == '__main__':
-    unittest.main()
-
diff --git a/tests/test_hertz.py b/tests/test_hertz.py
deleted file mode 100644
index 7624018c..00000000
--- a/tests/test_hertz.py
+++ /dev/null
@@ -1,103 +0,0 @@
-#
-# Copyright 2014-2015, 2020 Lars Pastewka (U. Freiburg)
-#           2014 James Kermode (Warwick U.)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-# ======================================================================
-# matscipy - Python materials science tools
-# https://github.com/libAtoms/matscipy
-#
-# Copyright (2014) James Kermode, King's College London
-#                  Lars Pastewka, Karlsruhe Institute of Technology
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-# ======================================================================
-
-import unittest
-
-import numpy as np
-
-import matscipytest
-import matscipy.contact_mechanics.Hertz as Hertz
-
-###
-
-class TestHertz(matscipytest.MatSciPyTestCase):
-
-    def test_Hertz_centerline_stress(self):
-        z = np.linspace(0.0, 5.0, 101)
-        for nu in [ 0.3, 0.5 ]:
-            srr1, szz1 = Hertz.centerline_stress(z, poisson=nu)
-            stt2, srr2, szz2, srz2 = Hertz.stress(np.zeros_like(z), z, poisson=nu)
-
-            self.assertTrue(np.max(np.abs(srr1-srr2)) < 1e-6)
-            self.assertTrue(np.max(np.abs(srr1-stt2)) < 1e-6)
-            self.assertTrue(np.max(np.abs(szz1-szz2)) < 1e-6)
-
-
-    def test_Hertz_surface_stress(self):
-        r = np.linspace(0.0, 5.0, 101)
-        for nu in [ 0.3, 0.5 ]:
-            pzz1, srr1, stt1 = Hertz.surface_stress(r, poisson=nu)
-            stt2, srr2, szz2, srz2 = Hertz.stress(r, np.zeros_like(r), poisson=nu)
-
-            self.assertTrue(np.max(np.abs(pzz1+szz2)) < 1e-6)
-            self.assertTrue(np.max(np.abs(srr1-srr2)) < 1e-6)
-            self.assertTrue(np.max(np.abs(stt1-stt2)) < 1e-6)
-
-    def test_Hertz_subsurface_stress(self):
-        nx = 51 # Grid size
-        a = 32. # Contact radius
-
-        y = np.linspace(0, 3*a, nx)
-        z = np.linspace(0, 3*a, nx)
-
-        y, z = np.meshgrid(y, z)
-
-        x = np.zeros_like(y)
-
-        # nu: Poisson
-        for nu in [0.3, 0.5]:
-            sxx, syy, szz, syz, sxz, sxy = \
-                Hertz.stress_Cartesian(x/a, y/a, z/a, poisson=nu)
-
-            r_sq = (x**2 + y**2)/a**2
-            stt2, srr2, szz2, srz2 = \
-                Hertz.stress(np.sqrt(r_sq), z/a, poisson=nu)
-
-            self.assertTrue(np.max(np.abs(sxx-stt2)) < 1e-12)
-            self.assertTrue(np.max(np.abs(syy-srr2)) < 1e-12)
-            self.assertTrue(np.max(np.abs(szz-szz2)) < 1e-12)
-            self.assertTrue(np.max(np.abs(syz-srz2)) < 1e-12)
-
-###
-
-if __name__ == '__main__':
-    unittest.main()