Repository: pace-neutrons/Euphonic (language: python)

Issue:
Python type-hinting
Raised by @DBJochym in a meeting: we are currently using Python's type-hinting in the Euphonic code-base, but not entirely consistently. It would be helpful to document the project's attitude towards this and outline a policy for contributed code if appropriate.

My current thoughts on this are that while type-hinting is an expressive way of writing and documenting Python functions, the static analysis tool Mypy does not work very usefully with Numpy arrays and Pint quantities. As these constitute a significant fraction of the data passed around Euphonic, I would be reluctant to impose an "all code should be annotated" or "mypy must pass" policy. But maybe if we think about this a little more there is another useful policy we can consider. At the very least we should mention our preferences in CONTRIBUTING.

@rebeccafair 


Discussion / hints:
Just updating some of the things we've touched on in other discussions here.

The generic use of `Quantity` and `np.ndarray` for typing is not ideal as it omits shape, type (float/int) and units information, but there are signs from the pint developers this may improve in the future in Pint issues 778, 1166 and currently there is a Pint PR in progress 1259. This would probably require bumping the versions of Pint/Numpy we depend on, but maybe Mantid's version restrictions will also have moved on by then so it might not be an issue. For now we should continue using `Quantity` and `np.ndarray` and putting shape, type and unit info in the docstring.

One other issue is that we currently have a bit of inconsistency in our use of `TypeVar`. `QpointFrequencies` and `QpointPhononmodes` use a generic `T` symbol to indicate 'whatever this class is', which works when subclassing because you get a generic `T` in the docs, rather than something contradictory. Whereas `Spectrum`, `Spectrum1D` etc. use `S`, `S1D` which means that in the `Spectrum1D` subclass you get a mix of `S` and `S1D` depending on which methods are inherited, which is a bit confusing. I would propose to move everything to use a generic `T` to avoid this.

In any case, Euphonic should be updated so all classes/methods/functions handle typing consistently, which would be an improvement, as currently some functions have no typing at all.

Whether we impose any hard restrictions on linting/type checking and what tools we'll use should maybe be another discussion. There are some brief Euphonic style guidelines here which could be added to https://github.com/pace-neutrons/pace-developers/blob/master/euphonic/development/03_style_guidelines.md.

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 1e25f9d0a..99a1e366c 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,10 +1,17 @@
 `Unreleased <https://github.com/pace-neutrons/Euphonic/compare/v0.6.1...HEAD>`_
 ----------
 
+- Improvements:
+
+  - Type hinting is now handled more consistently across different Euphonic
+    classes and functions
+
 - Bug Fixes:
 
   - Will no longer raise a KeyError reading from ``phonopy.yaml`` if
     ``physical_unit`` key is not present, instead will assume default units
+  - Can now read Phonopy BORN files where the (optional) NAC conversion
+    factor is not present
 
 `v0.6.1 <https://github.com/pace-neutrons/Euphonic/compare/v0.6.0...v0.6.1>`_
 ------
diff --git a/c/_euphonic.c b/c/_euphonic.c
index f8e645f31..3eb065517 100644
--- a/c/_euphonic.c
+++ b/c/_euphonic.c
@@ -29,7 +29,6 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
     PyArrayObject *py_modegs;
     PyArrayObject *py_all_ogs_cart;
     int dipole;
-    int reciprocal_asr;
     int splitting;
     int n_threads = 1;
 
@@ -84,6 +83,7 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
     int n_rqpts;
     int dmats_len;
     int modegs_len;
+    int asr_corr_len;
     int n_split_qpts;
     int q, i, qpos;
     int max_ims;
@@ -93,7 +93,7 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
     int n_gvecs;
 
     // Parse inputs
-    if (!PyArg_ParseTuple(args, "OO!O!O!O!O!O!O!O!O!iiiO!O!O!O!i",
+    if (!PyArg_ParseTuple(args, "OO!O!O!O!O!O!O!O!O!iiO!O!O!O!i",
                           &py_idata,
                           &PyArray_Type, &py_cell_vec,
                           &PyArray_Type, &py_recip_vec,
@@ -105,7 +105,6 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
                           &PyArray_Type, &py_asr_correction,
                           &PyArray_Type, &py_dmat_weighting,
                           &dipole,
-                          &reciprocal_asr,
                           &splitting,
                           &PyArray_Type, &py_evals,
                           &PyArray_Type, &py_dmats,
@@ -178,6 +177,7 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
     n_split_qpts = PyArray_DIMS(py_split_idx)[0];
     dmats_len = PyArray_DIMS(py_dmats)[0];
     modegs_len = PyArray_DIMS(py_modegs)[0];
+    asr_corr_len = PyArray_DIMS(py_asr_correction)[0];
     max_ims = PyArray_DIMS(py_sc_im_idx)[3];
     dmat_elems = 2*9*n_atoms*n_atoms;
     if (dipole) {
@@ -242,7 +242,7 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
                 add_arrays(dmat_elems, corr, dmat);
             }
 
-            if (reciprocal_asr) {
+            if (asr_corr_len > 0) {
                 add_arrays(dmat_elems, asr_correction, dmat);
             }
 
diff --git a/euphonic/cli/dispersion.py b/euphonic/cli/dispersion.py
index a4ac9c3b2..21406da47 100644
--- a/euphonic/cli/dispersion.py
+++ b/euphonic/cli/dispersion.py
@@ -1,13 +1,15 @@
-from typing import List
+from argparse import ArgumentParser
+from typing import List, Optional
 
 import euphonic
 from euphonic.plot import plot_1d
+from euphonic import Spectrum1D
 from .utils import (load_data_from_file, get_args, _bands_from_force_constants,
                     _get_q_distance, matplotlib_save_or_show, _get_cli_parser,
                     _calc_modes_kwargs)
 
 
-def main(params: List[str] = None):
+def main(params: Optional[List[str]] = None) -> None:
     args = get_args(get_parser(), params)
     data = load_data_from_file(args.filename)
 
@@ -46,7 +48,7 @@ def main(params: List[str] = None):
     if x_tick_labels:
         spectrum.x_tick_labels = x_tick_labels
 
-    spectra = spectrum.split(**split_args)
+    spectra = spectrum.split(**split_args)  # type: List[Spectrum1D]
 
     _ = plot_1d(spectra,
                 title=args.title,
@@ -57,7 +59,7 @@ def main(params: List[str] = None):
     matplotlib_save_or_show(save_filename=args.save_to)
 
 
-def get_parser():
+def get_parser() -> ArgumentParser:
     parser, _ = _get_cli_parser(features={'read-fc', 'read-modes', 'plotting',
                                           'q-e', 'btol'})
     parser.description = (
diff --git a/euphonic/cli/dos.py b/euphonic/cli/dos.py
index d9a47ad50..266b12a42 100644
--- a/euphonic/cli/dos.py
+++ b/euphonic/cli/dos.py
@@ -1,4 +1,5 @@
-from typing import List
+from argparse import ArgumentParser
+from typing import List, Optional
 
 from euphonic import (ureg, ForceConstants, QpointPhononModes,
                       Spectrum1DCollection)
@@ -11,7 +12,7 @@
                     _arrange_pdos_groups)
 
 
-def main(params: List[str] = None):
+def main(params: Optional[List[str]] = None) -> None:
     parser = get_parser()
     args = get_args(parser, params)
 
@@ -78,7 +79,7 @@ def main(params: List[str] = None):
     matplotlib_save_or_show(save_filename=args.save_to)
 
 
-def get_parser():
+def get_parser() -> ArgumentParser:
     parser, _ = _get_cli_parser(features={'read-fc', 'read-modes', 'mp-grid',
                                           'plotting', 'ebins',
                                           'adaptive-broadening',
diff --git a/euphonic/cli/intensity_map.py b/euphonic/cli/intensity_map.py
index fda682c5b..d8b1e9123 100755
--- a/euphonic/cli/intensity_map.py
+++ b/euphonic/cli/intensity_map.py
@@ -1,10 +1,10 @@
-import argparse
-from typing import List
+from argparse import ArgumentParser
+from typing import List, Optional
 
 import numpy as np
 
 import euphonic
-from euphonic import ureg
+from euphonic import ureg, Spectrum2D
 import euphonic.plot
 from euphonic.util import get_qpoint_labels
 from .utils import (_bands_from_force_constants, _calc_modes_kwargs,
@@ -14,7 +14,7 @@
                     matplotlib_save_or_show)
 
 
-def main(params: List[str] = None) -> None:
+def main(params: Optional[List[str]] = None) -> None:
     args = get_args(get_parser(), params)
     calc_modes_kwargs = _calc_modes_kwargs(args)
 
@@ -87,7 +87,7 @@ def main(params: List[str] = None) -> None:
     if x_tick_labels:
         spectrum.x_tick_labels = x_tick_labels
 
-    spectra = spectrum.split(**split_args)
+    spectra = spectrum.split(**split_args)  # type: List[Spectrum2D]
     if len(spectra) > 1:
         print(f"Found {len(spectra)} regions in q-point path")
 
@@ -100,7 +100,7 @@ def main(params: List[str] = None) -> None:
     matplotlib_save_or_show(save_filename=args.save_to)
 
 
-def get_parser() -> argparse.ArgumentParser:
+def get_parser() -> ArgumentParser:
     parser, sections = _get_cli_parser(
         features={'read-fc', 'read-modes', 'q-e', 'map', 'btol', 'ebins',
                   'ins-weighting', 'plotting'})
diff --git a/euphonic/cli/optimise_dipole_parameter.py b/euphonic/cli/optimise_dipole_parameter.py
index 42888be95..cfea4b1e7 100644
--- a/euphonic/cli/optimise_dipole_parameter.py
+++ b/euphonic/cli/optimise_dipole_parameter.py
@@ -7,9 +7,9 @@
 value
 """
 
-import argparse
+from argparse import ArgumentParser
 import time
-from typing import List, Tuple
+from typing import List, Tuple, Optional
 import warnings
 
 import numpy as np
@@ -18,7 +18,7 @@
                                 force_constants_from_file)
 
 
-def main(params: List[str] = None) -> None:
+def main(params: Optional[List[str]] = None) -> None:
     args = get_args(get_parser(), params)
     params = vars(args)
     params.update({'print_to_terminal': True})
@@ -33,7 +33,7 @@ def calculate_optimum_dipole_parameter(
         n: int = 500,
         print_to_terminal: bool = False,
         **calc_modes_kwargs
-        ) -> Tuple[float, float, float, np.array, np.array, np.array]:
+        ) -> Tuple[float, float, float, np.ndarray, np.ndarray, np.ndarray]:
     """
     Calculate the optimum dipole_parameter and other dipole_parameters
     from the filename castep_bin file
@@ -102,7 +102,6 @@ def calculate_optimum_dipole_parameter(
 
         # Time per qpt
         qpts = np.full((n, 3), 0.5)
-        t_total = []
         start = time.time()
         # Need reduce_qpts=False because all q-points are the same,
         # if reduce_qpts=True only one q-point would be calculated
@@ -127,7 +126,7 @@ def calculate_optimum_dipole_parameter(
             dipole_parameters, t_init, t_per_qpt)
 
 
-def get_parser() -> argparse.ArgumentParser:
+def get_parser() -> ArgumentParser:
     parser, sections = _get_cli_parser(
         features={'read-fc', 'dipole-parameter-optimisation'})
     parser.description=(
diff --git a/euphonic/cli/powder_map.py b/euphonic/cli/powder_map.py
index 6bc5190c5..33f78d59c 100755
--- a/euphonic/cli/powder_map.py
+++ b/euphonic/cli/powder_map.py
@@ -1,8 +1,6 @@
-#! /usr/bin/env python3
-
-import argparse
+from argparse import ArgumentParser
 from math import ceil
-from typing import List
+from typing import List, Optional
 
 import numpy as np
 
@@ -29,7 +27,7 @@ def tqdm(sequence):
             return sequence
 
 
-def get_parser() -> 'argparse.ArgumentParser':
+def get_parser() -> ArgumentParser:
 
     parser, sections = _get_cli_parser(
         features={'read-fc', 'pdos-weighting', 'ins-weighting',
@@ -54,7 +52,7 @@ def get_parser() -> 'argparse.ArgumentParser':
     return parser
 
 
-def main(params: List[str] = None):
+def main(params: Optional[List[str]] = None) -> None:
     args = get_args(get_parser(), params)
     calc_modes_kwargs = _calc_modes_kwargs(args)
 
@@ -206,7 +204,3 @@ def update_max(max_val):
         maxbox.on_submit(update_max)
 
     matplotlib_save_or_show(save_filename=args.save_to)
-
-
-if __name__ == '__main__':
-    main()
diff --git a/euphonic/cli/show_sampling.py b/euphonic/cli/show_sampling.py
index 6f8aed8a5..bf7e6fef4 100644
--- a/euphonic/cli/show_sampling.py
+++ b/euphonic/cli/show_sampling.py
@@ -1,6 +1,6 @@
 # -*- coding: UTF-8 -*-
-import argparse
-from typing import List
+from argparse import ArgumentParser
+from typing import List, Optional
 
 import matplotlib.pyplot as plt
 from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
@@ -14,8 +14,8 @@
               'random-sphere'}
 
 
-def get_parser() -> argparse.ArgumentParser:
-    parser = argparse.ArgumentParser()
+def get_parser() -> ArgumentParser:
+    parser = ArgumentParser()
     parser.add_argument('npts', type=int)
     parser.add_argument('sampling', type=str,
                         choices=(choices_2d | choices_3d))
@@ -27,7 +27,7 @@ def get_parser() -> argparse.ArgumentParser:
     return parser
 
 
-def main(params: List[str] = None) -> None:
+def main(params: Optional[List[str]] = None) -> None:
     parser = get_parser()
     args = get_args(parser, params)
 
diff --git a/euphonic/cli/utils.py b/euphonic/cli/utils.py
index ff733eac7..4dc1eeadc 100644
--- a/euphonic/cli/utils.py
+++ b/euphonic/cli/utils.py
@@ -30,11 +30,11 @@ def force_constants_from_file(filename: Union[str, os.PathLike]
 
     Returns
     -------
-    ForceConstants
+    fc
     """
     path = pathlib.Path(filename)
     if path.suffix in ['.hdf5', '.yaml']:
-        phonopy_kwargs = {}
+        phonopy_kwargs: Dict[str, Union[str, os.PathLike]] = {}
         phonopy_kwargs['path'] = path.parent
         if (path.parent / 'BORN').is_file():
             phonopy_kwargs['born_name'] = 'BORN'
@@ -77,7 +77,7 @@ def modes_from_file(filename: Union[str, os.PathLike]
 
     Returns
     -------
-    QpointPhononmodes
+    modes
     """
     path = pathlib.Path(filename)
     if path.suffix == '.phonon':
@@ -113,12 +113,12 @@ def load_data_from_file(filename: Union[str, os.PathLike]
 
     Parameters
     ----------
-    filename : str
+    filename
         The file with a path
 
     Returns
     -------
-    QpointPhononmodes or ForceConstants
+    file_data
     """
     qpoint_phonon_modes_suffixes = ('.phonon')
     force_constants_suffixes = ('.castep_bin', '.check')
@@ -145,22 +145,24 @@ def load_data_from_file(filename: Union[str, os.PathLike]
             " data from Euphonic should have extension '.json'.")
 
 
-def get_args(parser: ArgumentParser, params: List[str] = None):
+def get_args(parser: ArgumentParser, params: Optional[List[str]] = None
+             ) -> Namespace:
     """
     Get the arguments from the parser. params should only be none when
     running from command line.
 
     Parameters
     ----------
-    parser : ArgumentParser
+    parser
         The parser to get the arguments from
-    params : List[str], optional
+    params
         The parameters to get arguments from, if None,
          then parse_args gets them from command line args
 
     Returns
     -------
-    Arguments object for use e.g. args.unit
+    args
+        Arguments object for use e.g. args.unit
     """
     if params is None:
         args = parser.parse_args()
@@ -169,14 +171,14 @@ def get_args(parser: ArgumentParser, params: List[str] = None):
     return args
 
 
-def matplotlib_save_or_show(save_filename: str = None):
+def matplotlib_save_or_show(save_filename: str = None) -> None:
     """
     Save or show the current matplotlib plot.
     Show if save_filename is not None which by default it is.
 
     Parameters
     ----------
-    save_filename : str, optional
+    save_filename
         The file to save the plot in
     """
     import matplotlib.pyplot as plt
@@ -202,7 +204,7 @@ def _get_q_distance(length_unit_string: str, q_distance: float) -> Quantity:
 
 
 def _get_energy_bins(
-        modes: QpointPhononModes,
+        modes: Union[QpointPhononModes, QpointFrequencies],
         n_ebins: int, emin: Optional[float] = None,
         emax: Optional[float] = None,
         headroom: float = 1.05) -> Quantity:
@@ -257,10 +259,8 @@ def _get_break_points(bandpath: dict) -> List[int]:
 
     Returns
     -------
-    list[int]
-
+    break_points
         Indices at which the spectrum should be split into subplots
-
     """
     # Find break points between continuous spectra: wherever there are two
     # adjacent labels
@@ -348,10 +348,10 @@ def _bands_from_force_constants(data: ForceConstants,
 def _grid_spec_from_args(crystal: Crystal,
                            grid: Optional[Sequence[int]] = None,
                            grid_spacing: Quantity = 0.1 * ureg('1/angstrom')
-                           ) -> Sequence[int]:
+                           ) -> Tuple[int, int, int]:
     """Get Monkorst-Pack mesh divisions from user arguments"""
     if grid:
-        grid_spec = grid
+        grid_spec = tuple(grid)
     else:
         grid_spec = crystal.get_mp_grid_spec(spacing=grid_spacing)
     return grid_spec
@@ -375,7 +375,7 @@ def _get_debye_waller(temperature: Quantity,
     return dw_phonons.calculate_debye_waller(temperature)
 
 
-def _get_pdos_weighting(cl_arg_weighting: str) -> str:
+def _get_pdos_weighting(cl_arg_weighting: str) -> Optional[str]:
     """
     Convert CL --weighting to weighting for calculate_pdos
     e.g. --weighting coherent-dos to weighting=coherent
@@ -424,16 +424,22 @@ def _get_cli_parser(features: Collection[str] = {}
                                Dict[str, _ArgumentGroup]]:
     """Instantiate an ArgumentParser with appropriate argument groups
 
-    Args:
-        sections:
-            collection (e.g. set, list) of str for known argument groups.
-            Known keys: read-fc, read-modes, ins-weighting, pdos-weighting,
-            powder, mp-grid, plotting, ebins, adaptive-broadening, q-e,
-            map, btol, dipole-parameter-optimisation
+    Parameters
+    ----------
+    features
+        collection (e.g. set, list) of str for known argument groups.
+        Known keys: read-fc, read-modes, ins-weighting, pdos-weighting,
+        powder, mp-grid, plotting, ebins, adaptive-broadening, q-e,
+        map, btol, dipole-parameter-optimisation
 
-    Returns:
-        Parser and a dict of parser subsections. This allows their help strings
-        to be customised and specialised options to be added.
+    Returns
+    -------
+    parser
+        ArgumentParser with the requested features
+    sections
+        Dictionary of section labels and their argument groups. This
+        allows their help strings to be customised and specialised
+        options to be added
 
     """
     def deprecation_text(deprecated_arg: str, new_arg: str):
diff --git a/euphonic/crystal.py b/euphonic/crystal.py
index 3e48c81af..07406d46a 100644
--- a/euphonic/crystal.py
+++ b/euphonic/crystal.py
@@ -1,6 +1,6 @@
 import inspect
 from math import ceil
-from typing import List, Tuple, TypeVar, Dict, Any
+from typing import List, Tuple, Type, TypeVar, Dict, Any
 from collections import OrderedDict
 
 import numpy as np
@@ -12,10 +12,6 @@
                          _obj_to_dict, _process_dict)
 from euphonic import ureg, Quantity
 
-
-CR = TypeVar('CR', bound='Crystal')
-
-
 class Crystal:
     """
     Stores lattice and atom information
@@ -37,6 +33,7 @@ class Crystal:
         Shape (n_atoms,) float Quantity in mass units. The mass of each
         atom in the unit cell
     """
+    T = TypeVar('T', bound='Crystal')
 
     def __init__(self, cell_vectors: Quantity, atom_r: np.ndarray,
                  atom_type: np.ndarray, atom_mass: Quantity) -> None:
@@ -82,26 +79,26 @@ def __init__(self, cell_vectors: Quantity, atom_r: np.ndarray,
         self.atom_mass_unit = str(atom_mass.units)
                 
     @property
-    def cell_vectors(self):
+    def cell_vectors(self) -> Quantity:
         return self._cell_vectors*ureg('bohr').to(
             self.cell_vectors_unit)
 
     @cell_vectors.setter
-    def cell_vectors(self, value):
+    def cell_vectors(self, value: Quantity) -> None:
         self.cell_vectors_unit = str(value.units)
         self._cell_vectors = value.to('bohr').magnitude
 
     @property
-    def atom_mass(self):
+    def atom_mass(self) -> Quantity:
         return self._atom_mass*ureg('m_e').to(
             self.atom_mass_unit)
 
     @atom_mass.setter
-    def atom_mass(self, value):
+    def atom_mass(self, value: Quantity) -> None:
         self.atom_mass_unit = str(value.units)
         self._atom_mass = value.to('m_e').magnitude
 
-    def __setattr__(self, name, value):
+    def __setattr__(self, name: str, value: Any) -> None:
         _check_unit_conversion(self, name, value,
                                ['cell_vectors_unit', 'atom_mass_unit'])
         super(Crystal, self).__setattr__(name, value)
@@ -137,12 +134,12 @@ def cell_volume(self) -> Quantity:
         Returns
         -------
         volume
-            Scalar float quantity in length**3 units. The cell volume
+            Scalar float Quantity in length**3 units. The cell volume
         """
         vol = self._cell_volume()*ureg.bohr**3
         return vol.to(ureg(self.cell_vectors_unit)**3)
 
-    def _cell_volume(self):
+    def _cell_volume(self) -> float:
         return _cell_vectors_to_volume(self._cell_vectors)
 
     def get_mp_grid_spec(self,
@@ -300,7 +297,7 @@ def to_dict(self) -> Dict[str, Any]:
                                    'atom_type', 'atom_mass'])
         return dout
 
-    def to_json_file(self, filename: str):
+    def to_json_file(self, filename: str) -> None:
         """
         Write to a JSON file. JSON fields are equivalent to
         Crystal.from_dict keys
@@ -313,7 +310,7 @@ def to_json_file(self, filename: str):
         _obj_to_json_file(self, filename)
 
     @classmethod
-    def from_dict(cls: CR, d: Dict[str, Any]) -> CR:
+    def from_dict(cls: Type[T], d: Dict[str, Any]) -> T:
         """
         Convert a dictionary to a Crystal object
 
@@ -338,7 +335,7 @@ def from_dict(cls: CR, d: Dict[str, Any]) -> CR:
                    d['atom_mass'])
 
     @classmethod
-    def from_json_file(cls: CR, filename: str) -> CR:
+    def from_json_file(cls: Type[T], filename: str) -> T:
         """
         Read from a JSON file. See Crystal.from_dict for required fields
 
@@ -354,7 +351,7 @@ def from_json_file(cls: CR, filename: str) -> CR:
         return _obj_from_json_file(cls, filename)
 
     @classmethod
-    def from_cell_vectors(cls: CR, cell_vectors: Quantity) -> CR:
+    def from_cell_vectors(cls: Type[T], cell_vectors: Quantity) -> T:
         """
         Create a Crystal object from just cell vectors, containing no
         detailed structure information (atomic positions, species,
diff --git a/euphonic/debye_waller.py b/euphonic/debye_waller.py
index 0f4a7e9ec..c05dcea89 100644
--- a/euphonic/debye_waller.py
+++ b/euphonic/debye_waller.py
@@ -1,4 +1,5 @@
 import inspect
+from typing import TypeVar, Dict, Any, Type
 
 from euphonic.validate import _check_constructor_inputs, _check_unit_conversion
 from euphonic.io import (_obj_to_json_file, _obj_from_json_file,
@@ -13,28 +14,34 @@ class DebyeWaller:
 
     Attributes
     ----------
-    crystal : Crystal
+    crystal
         Lattice and atom information
-    debye_waller : (n_ions, 3, 3) float Quantity
-        The anisotropic Debye-Waller exponent W_ab for each atom, where
-        the Debye-Waller factor is exp(-W_ab*Q_a*Q_b) where a,b run over
+    debye_waller
+        Shape (n_ions, 3, 3) float Quantity. The anisotropic
+        Debye-Waller exponent W_ab for each atom, where the
+        Debye-Waller factor is exp(-W_ab*Q_a*Q_b) where a,b run over
         the 3 Cartesian directions
-    temperature : float Quantity
-        The temperature the Debye-Waller exponent was calculated at
+    temperature
+        Scalar float Quantity. The temperature the Debye-Waller
+        exponent was calculated at
     """
+    T = TypeVar('T', bound='DebyeWaller')
 
-    def __init__(self, crystal, debye_waller, temperature):
+    def __init__(self, crystal: Crystal, debye_waller: Quantity,
+                 temperature: Quantity) -> None:
         """
         Parameters
         ----------
-        crystal : Crystal
+        crystal
             Lattice and atom information
-        debye_waller : (n_atoms, 3, 3) Quantity
-            The anisotropic Debye-Waller exponent W_ab for each atom,
-            where the Debye-Waller factor is exp(-W_ab*Q_a*Q_b) where
-            a,b run over the 3 Cartesian directions
-        temperature : float Quantity
-            The temperature the Debye-Waller exponent was calculated at
+        debye_waller
+            Shape (n_atoms, 3, 3) Quantity. The anisotropic
+            Debye-Waller exponent W_ab for each atom, where the
+            Debye-Waller factor is exp(-W_ab*Q_a*Q_b) where a,b run
+            over the 3 Cartesian directions
+        temperature
+            Scalar float Quantity. The temperature the Debye-Waller
+            exponent was calculated at
         """
         _check_constructor_inputs([crystal], [Crystal], [()], ['crystal'])
         n_atoms = crystal.n_atoms
@@ -51,31 +58,31 @@ def __init__(self, crystal, debye_waller, temperature):
         self.temperature_unit = str(temperature.units)
 
     @property
-    def debye_waller(self):
+    def debye_waller(self) -> Quantity:
         return self._debye_waller*ureg('bohr**2').to(
             self.debye_waller_unit)
 
     @debye_waller.setter
-    def debye_waller(self, value):
+    def debye_waller(self, value: Quantity) -> None:
         self.debye_waller_unit = str(value.units)
         self._debye_waller = value.to('bohr**2').magnitude
 
     @property
-    def temperature(self):
+    def temperature(self) -> Quantity:
         # See https://pint.readthedocs.io/en/latest/nonmult.html
         return Quantity(self._temperature, ureg('K')).to(self.temperature_unit)
 
     @temperature.setter
-    def temperature(self, value):
+    def temperature(self, value: Quantity) -> None:
         self.temperature_unit = str(value.units)
         self._temperature = value.to('K').magnitude
 
-    def __setattr__(self, name, value):
+    def __setattr__(self, name: str, value: Any) -> None:
         _check_unit_conversion(self, name, value,
                                ['debye_waller_unit', 'temperature_unit'])
         super(DebyeWaller, self).__setattr__(name, value)
 
-    def to_dict(self):
+    def to_dict(self) -> Dict[str, Any]:
         """
         Convert to a dictionary. See DebyeWaller.from_dict for details
         on keys/values
@@ -87,26 +94,26 @@ def to_dict(self):
         dout = _obj_to_dict(self, ['crystal', 'debye_waller', 'temperature'])
         return dout
 
-    def to_json_file(self, filename):
+    def to_json_file(self, filename: str) -> None:
         """
         Write to a JSON file. JSON fields are equivalent to
         DebyeWaller.from_dict keys
 
         Parameters
         ----------
-        filename : str
+        filename
             Name of the JSON file to write to
         """
         _obj_to_json_file(self, filename)
 
     @classmethod
-    def from_dict(cls, d):
+    def from_dict(cls: Type[T], d: Dict[str, Any]) -> T:
         """
         Convert a dictionary to a DebyeWaller object
 
         Parameters
         ----------
-        d : dict
+        d
             A dictionary with the following keys/values:
 
             - 'crystal': dict, see Crystal.from_dict
@@ -125,14 +132,14 @@ def from_dict(cls, d):
         return cls(crystal, d['debye_waller'], d['temperature'])
 
     @classmethod
-    def from_json_file(cls, filename):
+    def from_json_file(cls: Type[T], filename: str) -> T:
         """
         Read from a JSON file. See DebyeWaller.from_dict for required
         fields
 
         Parameters
         ----------
-        filename : str
+        filename
             The file to read from
 
         Returns
diff --git a/euphonic/force_constants.py b/euphonic/force_constants.py
index 90072e4cf..9a8ab052d 100644
--- a/euphonic/force_constants.py
+++ b/euphonic/force_constants.py
@@ -2,7 +2,7 @@
 import os
 import sys
 import warnings
-from typing import Optional, Tuple, Union
+from typing import Optional, Tuple, Union, TypeVar, Sequence, Dict, Any, Type
 from multiprocessing import cpu_count
 
 import numpy as np
@@ -36,39 +36,50 @@ class ForceConstants:
 
     Attributes
     ----------
-    crystal : Crystal
+    crystal
         Lattice and atom information
-    force_constants : (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float Quantity
-        Force constants matrix
-    sc_matrix : (3, 3) int ndarray
-        The supercell matrix
-    n_cells_in_sc : int
+    force_constants
+        Shape (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float Quantity in
+        energy/length**2 units. The force constants matrix
+    sc_matrix
+        Shape (3, 3) int ndarray. The supercell matrix
+    n_cells_in_sc
         Number of cells in the supercell
-    cell_origins : (n_cells_in_sc, 3) int ndarray
-        The locations of the unit cells within the supercell
-    born : (n_atoms, 3, 3) float Quantity or None
+    cell_origins
+        Shape (n_cells_in_sc, 3) int ndarray. The locations of the unit
+        cells within the supercell
+    born
+        Shape (n_atoms, 3, 3) float Quantity in charge units or None.
         The Born charges for each atom
-    dielectric : (3, 3) float Quantity or None
-        The dielectric permittivity tensor
+    dielectric
+        Shape (3, 3) float Quantity in charge**2/(length*energy) units
+        or None. The dielectric permittivity tensor
     """
+    T = TypeVar('T', bound='ForceConstants')
 
-    def __init__(self, crystal, force_constants, sc_matrix, cell_origins,
-                 born=None, dielectric=None):
+    def __init__(self, crystal: Crystal, force_constants: Quantity,
+                 sc_matrix: np.ndarray, cell_origins: np.ndarray,
+                 born: Optional[Quantity] = None,
+                 dielectric: Optional[Quantity] = None) -> None:
         """
         Parameters
         ----------
-        crystal : Crystal
+        crystal
             Lattice and atom information
-        force_constants : (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float Quantity
-            Force constants matrix
-        sc_matrix : (3, 3) int ndarray
-            The supercell matrix
-        cell_origins : (n_cells_in_sc, 3) int ndarray
-            The locations of the unit cells within the supercell
-        born : (n_atoms, 3, 3) float Quantity, optional
-            The Born charges for each atom
-        dielectric : (3, 3) float Quantity, optional
-            The dielectric permittivity tensor
+        force_constants
+            Shape (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float Quantity
+            in energy/length**2 units. The force constants matrix
+        sc_matrix
+            Shape (3, 3) int ndarray. The supercell matrix
+        cell_origins
+            Shape (n_cells_in_sc, 3) int ndarray. The locations of the
+            unit cells within the supercell
+        born
+            Shape (n_atoms, 3, 3) float Quantity in charge units. The
+            Born charges for each atom
+        dielectric
+            Shape (3, 3) float Quantity in charge**2/(length*energy)
+            units. The dielectric permittivity tensor
         """
         # Check independent inputs first
         _check_constructor_inputs(
@@ -105,30 +116,30 @@ def __init__(self, crystal, force_constants, sc_matrix, cell_origins,
                 'e**2/(bohr*hartree)')))
 
     @property
-    def force_constants(self):
+    def force_constants(self) -> Quantity:
         return self._force_constants*ureg(
             'hartree/bohr**2').to(
                 self.force_constants_unit)
 
     @force_constants.setter
-    def force_constants(self, value):
+    def force_constants(self, value: Quantity) -> None:
         self.force_constants_unit = str(value.units)
         self._force_constants = value.to('hartree/bohr**2').magnitude
 
     @property
-    def born(self):
+    def born(self) -> Union[Quantity, None]:
         if self._born is not None:
             return self._born*ureg('e').to(self.born_unit)
         else:
             return None
 
     @born.setter
-    def born(self, value):
+    def born(self, value: Quantity) -> None:
         self.born_unit = str(value.units)
         self._born = value.to('e').magnitude
 
     @property
-    def dielectric(self):
+    def dielectric(self) -> Union[Quantity, None]:
         if self._dielectric is not None:
             return self._dielectric*ureg((
                 'e**2/(bohr*hartree)')).to(self.dielectric_unit)
@@ -136,11 +147,11 @@ def dielectric(self):
             return None
 
     @dielectric.setter
-    def dielectric(self, value):
+    def dielectric(self, value: Quantity) -> None:
         self.dielectric_unit = str(value.units)
         self._dielectric = value.to('e**2/(bohr*hartree)').magnitude
 
-    def __setattr__(self, name, value):
+    def __setattr__(self, name: str, value: Any) -> None:
         _check_unit_conversion(self, name, value,
                                ['force_constants_unit', 'born_unit',
                                 'dielectric_unit'])
@@ -169,10 +180,12 @@ def calculate_qpoint_phonon_modes(
 
         Parameters
         ----------
-        qpts : (n_qpts, 3) float ndarray
-            The q-points to interpolate onto
-        weights : (n_qpts,) float ndarray
-            The weight for each q-point. If not given, equal weights are
+        qpts
+            Shape (n_qpts, 3) float ndarray. The q-points to
+            interpolate onto in reciprocal cell vector units
+        weights
+            Shape (n_qpts,) float ndarray. The weight for each q-point.
+            If not given, equal weights are
             applied
         asr
             One of {'realspace', 'reciprocal'}. Which acoustic sum rule
@@ -368,7 +381,7 @@ def calculate_qpoint_phonon_modes(
 
         .. [3] J. R. Yates, X. Wang, D. Vanderbilt and I. Souza, Phys. Rev. B, 2007, 75, 195121
         """
-        qpts, weights, freqs, evecs, grads = self._calculate_phonons_at_qpts(
+        qpts, freqs, weights, evecs, grads = self._calculate_phonons_at_qpts(
             qpts, weights, asr, dipole, dipole_parameter, eta_scale,
             splitting, insert_gamma, reduce_qpts, use_c, n_threads,
             return_mode_gradients, return_mode_widths,
@@ -402,7 +415,7 @@ def calculate_qpoint_frequencies(
         q-points. See ForceConstants.calculate_qpoint_phonon_modes for
         argument and algorithm details
         """
-        qpts, weights, freqs, _, grads = self._calculate_phonons_at_qpts(
+        qpts, freqs, weights, _, grads = self._calculate_phonons_at_qpts(
             qpts, weights, asr, dipole, dipole_parameter, eta_scale,
             splitting, insert_gamma, reduce_qpts, use_c, n_threads,
             return_mode_gradients, return_mode_widths,
@@ -417,8 +430,8 @@ def calculate_qpoint_frequencies(
     def _calculate_phonons_at_qpts(
             self,
             qpts: np.ndarray,
-            weights: np.ndarray,
-            asr: str,
+            weights: Optional[np.ndarray],
+            asr: Optional[str],
             dipole: bool,
             dipole_parameter: float,
             eta_scale: float,
@@ -429,9 +442,34 @@ def _calculate_phonons_at_qpts(
             n_threads: Optional[int],
             return_mode_gradients: bool,
             return_mode_widths: bool,
-            return_eigenvectors: bool) -> Union[
-                Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray],
-                Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]]:
+            return_eigenvectors: bool) -> Tuple[
+                np.ndarray, Quantity, Optional[np.ndarray],
+                Optional[np.ndarray], Optional[Quantity]]:
+        """
+        Calculates phonon frequencies, and optionally eigenvectors and
+        phonon frequency gradients. See calculate_qpoint_phonon_modes
+        for argument details
+
+        Returns
+        -------
+        qpts
+            Shape (n_qpts, 3) float ndarray. The calculated q-points,
+            note that this may not be the same as the input q-points
+            if splitting and insert_gamma=True and there are gamma
+            points
+        frequencies
+            Shape (n_qpts, 3*n_atoms) float Quantity in energy units,
+            the phonon frequencies. n_qpts may not be the same as the
+            number of input q-points, if splitting and
+            insert_gamma=True and there are gamma points
+        weights
+            Shape (n_qpts,) float ndarray. The weight of each q-point
+        eigenvectors
+            Shape (n_qpts, 3*n_atoms, n_atoms, 3) complex ndarray
+        gradients
+            Shape (n_qpts, 3*n_atoms, 3) float Quantity in
+            energy*length units. The phonon mode gradients
+        """
         if return_mode_widths:
             warnings.warn(
                 'return_mode_widths has been deprecated and will be removed '
@@ -520,21 +558,21 @@ def _calculate_phonons_at_qpts(
         # Get q-directions for non-analytical corrections
         if splitting:
             split_idx = np.where(is_gamma(reduced_qpts))[0]
-            q_dirs = self._get_q_dir(qpts, qpts_i, split_idx)
+            q_dirs = self._get_q_dirs(qpts, qpts_i, split_idx)
         else:
             split_idx = np.array([])
             q_dirs = np.array([])
 
-        lim = 2  # Supercell image limit
-        # Construct list of supercell ion images
+        n_sc_shells = 2  # How many supercells out to search for atom images
+        # Construct list of supercell atom images
         if not hasattr(self, '_sc_image_i'):
-            self._calculate_supercell_images(lim)
+            self._calculate_supercell_images(n_sc_shells)
 
         # Get a list of all the unique supercell image origins and cell
         # origins in x, y, z and how to rebuild them to minimise
         # expensive phase calculations later
         sc_image_r = get_all_origins(
-            np.repeat(lim, 3) + 1, min_xyz=-np.repeat(lim, 3))
+            np.repeat(n_sc_shells, 3) + 1, min_xyz=-np.repeat(n_sc_shells, 3))
         sc_origins = np.einsum('ij,jk->ik', sc_image_r,
                                self.sc_matrix).astype(np.int32)
         unique_sc_origins = [[] for i in range(3)]
@@ -572,17 +610,15 @@ def _calculate_phonons_at_qpts(
                        dipole_parameter != self._dipole_parameter):
             self._dipole_correction_init(dipole_parameter)
 
+        force_constants = self._force_constants
         if asr == 'realspace':
             if not hasattr(self, '_force_constants_asr'):
                 self._force_constants_asr = self._enforce_realspace_asr()
             force_constants = self._force_constants_asr
-        else:
-            force_constants = self._force_constants
         # Precompute fc matrix weighted by number of supercell atom
         # images (for cumulant method)
         n_sc_images_repeat = (
-            self._n_sc_images.repeat(3, axis=2).repeat(3, axis=1)
-        )
+            self._n_sc_images.repeat(3, axis=2).repeat(3, axis=1))
         fc_img_weighted = np.divide(
             force_constants, n_sc_images_repeat, where=n_sc_images_repeat != 0)
 
@@ -597,9 +633,7 @@ def _calculate_phonons_at_qpts(
             if dipole:
                 dyn_mat_gamma += self._calculate_dipole_correction(q_gamma)
             recip_asr_correction = self._enforce_reciprocal_asr(dyn_mat_gamma)
-            if len(recip_asr_correction) == 0:
-                # Finding acoustic modes failed
-                asr = None
+
 
         rfreqs = np.zeros((n_rqpts, 3*n_atoms))
         if return_eigenvectors:
@@ -668,27 +702,27 @@ def _calculate_phonons_at_qpts(
                             '_cells', '_gvec_phases', '_gvecs_cart',
                             '_dipole_q0']
             _ensure_contiguous_attrs(self, attrs, opt_attrs=dipole_attrs)
-            reciprocal_asr = 1 if asr == 'reciprocal' else 0
             with threadpool_limits(limits=1):
                 euphonic_c.calculate_phonons(
                     self, cell_vectors, recip_vectors, reduced_qpts,
                     split_idx, q_dirs, fc_img_weighted, sc_origins,
                     recip_asr_correction, dyn_mat_weighting, dipole,
-                    reciprocal_asr, splitting, rfreqs, reigenvecs,
-                    rmode_gradients, all_origins_cart, n_threads)
+                    splitting, rfreqs, reigenvecs, rmode_gradients,
+                    all_origins_cart, n_threads)
         else:
-            q_independent_args = (
-                reduced_qpts, split_idx, q_dirs, fc_img_weighted,
-                unique_sc_origins, unique_sc_i, unique_cell_origins,
-                unique_cell_i, recip_asr_correction, dyn_mat_weighting,
-                dipole, asr, splitting, all_origins_cart)
-            for q in range(n_rqpts):
-                rfreqs[q], evecs, grads = self._calculate_phonons_at_q(
-                        q, q_independent_args)
+            for qi, qpt in enumerate(reduced_qpts):
+                q_dir = None
+                if splitting and is_gamma(qpt):
+                    q_dir_idx = np.where(split_idx == qi)[0][0]
+                    q_dir = q_dirs[q_dir_idx]
+                rfreqs[qi], evecs, grads = self._calculate_phonons_at_q(
+                    qpt, fc_img_weighted, unique_sc_origins, unique_sc_i,
+                    unique_cell_origins, unique_cell_i, all_origins_cart,
+                    dyn_mat_weighting, recip_asr_correction, dipole, q_dir)
                 if return_eigenvectors:
-                    reigenvecs[q] = evecs
+                    reigenvecs[qi] = evecs
                 if return_mode_gradients:
-                    rmode_gradients[q] = grads
+                    rmode_gradients[qi] = grads
 
         freqs = rfreqs[qpts_i]*ureg('hartree').to('meV')
         if return_eigenvectors:
@@ -719,44 +753,102 @@ def _calculate_phonons_at_qpts(
             mode_gradients = mode_gradients_to_widths(
                 mode_gradients,
                 self.crystal.cell_vectors)
-        return qpts, weights, freqs, eigenvectors, mode_gradients
+        return qpts, freqs, weights, eigenvectors, mode_gradients
 
-    def _calculate_phonons_at_q(self, q, args):
+    def _calculate_phonons_at_q(
+            self,
+            qpt: np.ndarray,
+            fc_img_weighted: np.ndarray,
+            unique_sc_origins: Sequence[Sequence[int]],
+            unique_sc_i: np.ndarray,
+            unique_cell_origins: Sequence[Sequence[int]],
+            unique_cell_i: np.ndarray,
+            all_origins_cart: np.ndarray,
+            dyn_mat_weighting: np.ndarray,
+            recip_asr_correction: np.ndarray,
+            dipole: bool,
+            q_dir: Optional[np.ndarray] = None
+            ) -> Tuple[np.ndarray, np.ndarray, Union[np.ndarray, None]]:
         """
         Given a q-point and some precalculated q-independent values,
         calculate and diagonalise the dynamical matrix and return the
-        frequencies and eigenvalues. Optionally also includes the Ewald
-        dipole sum correction and LO-TO splitting
-        """
-        (reduced_qpts, split_idx, q_dirs, fc_img_weighted, unique_sc_origins,
-         unique_sc_i, unique_cell_origins, unique_cell_i,
-         recip_asr_correction, dyn_mat_weighting, dipole, asr,
-         splitting, all_origins_cart) = args
+        frequencies, eigenvectors and optionally mode gradients
+
+        Parameters
+        ----------
+        qpt
+            Shape (3,) float ndarray. The q-point to calculate
+        fc_img_weighted
+            Shape (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float ndarray.
+            The force constants matrix weighted by the number of
+            supercell atom images for each ij displacement
+        unique_sc_origins
+            A list containing 3 lists of the unique supercell image
+            offsets in each direction. The supercell offset is
+            calculated by multiplying the supercell matrix by the
+            supercell image indices (obtained by _get_all_origins()). A
+            list of lists rather than a Numpy array is used as the 3
+            lists are independent and their size is not known
+            beforehand
+        unique_sc_i
+            Shape ((2*n_sc_shells + 1)**3, 3) int ndarray, where
+            n_sc_shells is how many supercells out to go to search for
+            atom images. The indices needed to reconstruct sc_origins
+            from the unique values in unique_sc_origins
+        unique_cell_origins
+            A list containing 3 lists of the unique cell origins in
+            each direction. A list of lists rather than a Numpy array
+            is used as the 3 lists are independent and their size is
+            not known beforehand
+        unique_cell_i
+            Shape (n_cells_in_sc, 3) int ndarray. The indices needed to
+            reconstruct cell_origins from the unique values in
+            unique_cell_origins
+        all_origins_cart
+            Shape (n_cells_in_sc, n_atoms, n_atoms, n_images, 3) or
+            (0, 3) float ndarray. The Cartesian coordinate of each atom
+            in each supercell image. Is only used if calculating
+            dynamical matrix gradients
+        dyn_mat_weighting
+            Shape (3*n_atoms, 3*n_atoms) float ndarray. The dynamical
+            matrix mass weighting
+        recip_asr_correction
+            Shape (3*n_atoms, 3*n_atoms) or (0,) complex ndarray. Is
+            only used if the reciprocal ASR was requested
+        dipole
+            Whether to apply the dipole correction
+        q_dir
+            Shape (3,) float ndarray, the q-direction to use in LO-TO
+            splitting, if applicable
 
-        qpt = reduced_qpts[q]
+        Returns
+        -------
+        frequencies
+            Shape (3*n_atoms) float ndarray. The phonon frequencies
+        eigenvectors
+            Shape (3*n_atoms, n_atoms, 3) complex ndarray. The
+            eigenvectors
+        mode_gradients
+            Shape (3*n_atoms, 3) complex ndarray. The gradient for
+            each mode
+        """
         n_atoms = self.crystal.n_atoms
 
-        dmat_args = (qpt, fc_img_weighted, unique_sc_origins, unique_sc_i,
-                     unique_cell_origins, unique_cell_i, all_origins_cart)
-        dyn_mat, dmat_grad = self._calculate_dyn_mat(*dmat_args)
+        dyn_mat, dmat_grad = self._calculate_dyn_mat(
+            qpt, fc_img_weighted, unique_sc_origins, unique_sc_i,
+            unique_cell_origins, unique_cell_i, all_origins_cart)
 
         if dipole:
             dipole_corr = self._calculate_dipole_correction(qpt)
             dyn_mat += dipole_corr
 
-        if asr == 'reciprocal':
+        if len(recip_asr_correction) > 0:
             dyn_mat += recip_asr_correction
 
         # Calculate LO-TO splitting by calculating non-analytic
         # correction to dynamical matrix
-        if splitting and is_gamma(qpt):
-            q_dir_idx = np.where(split_idx == q)[0][0]
-            na_corr = self._calculate_gamma_correction(q_dirs[q_dir_idx])
-        else:
-            # Correction is zero if not a gamma point or splitting=False
-            na_corr = np.array([0])
-
-        dyn_mat += na_corr
+        if q_dir is not None:
+            dyn_mat += self._calculate_gamma_correction(q_dir)
 
         # Mass weight dynamical matrix
         dyn_mat *= dyn_mat_weighting
@@ -781,9 +873,16 @@ def _calculate_phonons_at_q(self, q, args):
         else:
             return evals, evecs, None
 
-    def _calculate_dyn_mat(self, q, fc_img_weighted, unique_sc_origins,
-                           unique_sc_i, unique_cell_origins, unique_cell_i,
-                           all_origins_cart):
+    def _calculate_dyn_mat(
+            self,
+            qpt: np.ndarray,
+            fc_img_weighted: np.ndarray,
+            unique_sc_origins: Sequence[Sequence[int]],
+            unique_sc_i: np.ndarray,
+            unique_cell_origins: Sequence[Sequence[int]],
+            unique_cell_i: np.ndarray,
+            all_origins_cart: np.ndarray) -> Tuple[np.ndarray,
+                                                   Optional[np.ndarray]]:
         """
         Calculate the non mass weighted dynamical matrix at a specified
         q-point from the image weighted force constants matrix and the
@@ -792,38 +891,48 @@ def _calculate_dyn_mat(self, q, fc_img_weighted, unique_sc_origins,
 
         Parameters
         ----------
-        q : (3,) float ndarray
-            The q-point to calculate the correction for
-        fc_img_weighted : (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float ndarray
+        qpt
+            Shape (3,) float ndarray. The q-point to calculate
+        fc_img_weighted
+            Shape (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float ndarray.
             The force constants matrix weighted by the number of
             supercell atom images for each ij displacement
-        unique_sc_origins : list of lists of ints
+        unique_sc_origins
             A list containing 3 lists of the unique supercell image
             offsets in each direction. The supercell offset is
             calculated by multiplying the supercell matrix by the
             supercell image indices (obtained by _get_all_origins()). A
             list of lists rather than a Numpy array is used as the 3
-            lists are independent and their size is not known beforehand
-        unique_sc_i : ((2*lim + 1)**3, 3) int ndarray
-            The indices needed to reconstruct sc_origins from the unique
-            values in unique_sc_origins
-        unique_cell_origins : list of lists of ints
-            A list containing 3 lists of the unique cell origins in each
-            direction. A list of lists rather than a Numpy array is used
-            as the 3 lists are independent and their size is not known
+            lists are independent and their size is not known
             beforehand
-        unique_cell_i : (n_cells_in_sc, 3) int ndarray
-            The indices needed to reconstruct cell_origins from the
-            unique values in unique_cell_origins
-        all_origins_cart : (n_cells_in_sc, n_atoms, n_atoms, n_images, 3)
-            or (0, 3) float ndarray
-            The Cartesian coordinate of each atom in each supercell
-            image. Is only used if calculating dynamical matrix gradients
+        unique_sc_i
+            Shape ((2*n_sc_shells + 1)**3, 3) int ndarray, where
+            n_sc_shells is how many supercells out to go to search for
+            atom images. The indices needed to reconstruct sc_origins
+            from the unique values in unique_sc_origins
+        unique_cell_origins
+            A list containing 3 lists of the unique cell origins in
+            each direction. A list of lists rather than a Numpy array
+            is used as the 3 lists are independent and their size is
+            not known beforehand
+        unique_cell_i
+            Shape (n_cells_in_sc, 3) int ndarray. The indices needed to
+            reconstruct cell_origins from the unique values in
+            unique_cell_origins
+        all_origins_cart
+            Shape (n_cells_in_sc, n_atoms, n_atoms, n_images, 3) or
+            (0, 3) float ndarray. The Cartesian coordinate of each atom
+            in each supercell image. Is only used if calculating
+            dynamical matrix gradients
 
         Returns
         -------
-        dyn_mat : (3*n_atoms, 3*n_atoms) complex ndarray
-            The non mass weighted dynamical matrix at q
+        dyn_mat
+            Shape (3*n_atoms, 3*n_atoms) complex ndarray. The non mass
+            weighted dynamical matrix at q
+        dyn_mat_grad
+            Shape (3*n_atoms, 3*n_atoms, 3) complex ndarray. The dynamical
+            matrix gradient in each Cartesian direction
         """
         sc_image_i = self._sc_image_i
 
@@ -837,7 +946,7 @@ def _calculate_dyn_mat(self, q, fc_img_weighted, unique_sc_origins,
         # hence phase of zero can be used
         sc_phases = np.zeros(len(unique_sc_i) + 1, dtype=np.complex128)
         sc_phases[:-1], cell_phases = self._calculate_phases(
-            q, unique_sc_origins, unique_sc_i, unique_cell_origins,
+            qpt, unique_sc_origins, unique_sc_i, unique_cell_origins,
             unique_cell_i)
         sc_phase_sum = np.sum(sc_phases[sc_image_i],
                               axis=3)
@@ -859,7 +968,7 @@ def _calculate_dyn_mat(self, q, fc_img_weighted, unique_sc_origins,
         else:
             return dyn_mat, None
 
-    def _dipole_correction_init(self, dipole_parameter=1.0):
+    def _dipole_correction_init(self, dipole_parameter: float = 1.0) -> None:
         """
         Calculate the q-independent parts of the long range correction
         to the dynamical matrix for efficiency. The method used is based
@@ -868,7 +977,7 @@ def _dipole_correction_init(self, dipole_parameter=1.0):
 
         Parameters
         ----------
-        dipole_parameter : float, optional
+        dipole_parameter
             Changes the cutoff in real/reciprocal space for the dipole
             Ewald sum. A higher value uses more reciprocal terms
         """
@@ -1006,7 +1115,7 @@ def _dipole_correction_init(self, dipole_parameter=1.0):
         self._gvec_phases = gvec_phases
         self._dipole_q0 = dipole_q0
 
-    def _calculate_dipole_correction(self, q):
+    def _calculate_dipole_correction(self, q: np.ndarray) -> np.ndarray:
         """
         Calculate the long range correction to the dynamical matrix
         using the Ewald sum, see eqs 72-74 from Gonze and Lee PRB 55,
@@ -1014,13 +1123,15 @@ def _calculate_dipole_correction(self, q):
 
         Parameters
         ----------
-        q : (3,) float ndarray
-            The q-point to calculate the correction for
+        q
+            Shape (3,) float ndarray. The q-point to calculate the
+            correction for
 
         Returns
         -------
-        corr : (3*n_atoms, 3*n_atoms) complex ndarray
-            The correction to the dynamical matrix
+        corr
+            Shape (3*n_atoms, 3*n_atoms) complex ndarray. The
+            correction to the dynamical matrix
         """
         recip = self.crystal.reciprocal_cell().to('1/bohr').magnitude
         n_atoms = self.crystal.n_atoms
@@ -1089,7 +1200,29 @@ def _calculate_dipole_correction(self, q):
         return np.reshape(np.transpose(dipole, axes=[0, 2, 1, 3]),
                           (3*n_atoms, 3*n_atoms))
 
-    def _get_q_dir(self, qpts, qpts_i, gamma_idx):
+    def _get_q_dirs(self, qpts: np.ndarray, qpts_i: np.ndarray,
+                   gamma_idx: np.ndarray) -> np.ndarray:
+        """
+        Get q-directions at gamma points
+
+        Parameters
+        ----------
+        qpts
+            Shape (n_qpts, 3) float ndarray. The q-points.
+        qpts_i
+            Shape (n_qpts,) int ndarray. Describes how to index
+            the 'reduced' q-points back into the full q-points
+            array
+        gamma_idx
+            Shape (n_gamma_points,) int ndarray. The indices
+            of gamma-points in the full q-points array
+
+        Returns
+        -------
+        q_dirs
+            Shape (n_gamma_points, 3 float ndarray. The q-direction
+            for splitting at each gamma point)
+        """
         q_dirs = np.zeros((len(gamma_idx), 3))
         if len(qpts) > 1:
             for i, idx in enumerate(gamma_idx):
@@ -1109,7 +1242,7 @@ def _get_q_dir(self, qpts, qpts_i, gamma_idx):
                         q_dirs[i] = qpts[idx_in_qpts] - qpts[idx_in_qpts - 1]
         return q_dirs
 
-    def _calculate_gamma_correction(self, q_dir):
+    def _calculate_gamma_correction(self, q_dir: np.ndarray) -> np.ndarray:
         """
         Calculate non-analytic correction to the dynamical matrix at q=0
         for a specified direction of approach. See Eq. 60 of X. Gonze
@@ -1117,14 +1250,15 @@ def _calculate_gamma_correction(self, q_dir):
 
         Parameters
         ----------
-        q_dir : (3,) float ndarray
-            The direction along which q approaches 0, in reciprocal
-            fractional coordinates
+        q_dir
+            Shape (3,) float ndarray. The direction along which q
+            approaches 0, in reciprocal fractional coordinates
 
         Returns
         -------
-        na_corr : (3*n_atoms, 3*n_atoms) complex ndarray
-            The correction to the dynamical matrix
+        na_corr
+            Shape (3*n_atoms, 3*n_atoms) complex ndarray. The
+            correction to the dynamical matrix
         """
         n_atoms = self.crystal.n_atoms
         born = self._born
@@ -1150,21 +1284,21 @@ def _calculate_gamma_correction(self, q_dir):
 
         return na_corr
 
-    def _get_shell_origins(self, n):
+    def _get_shell_origins(self, n: int) -> np.ndarray:
         """
         Given the shell number, compute all the cell origins that lie in
         that shell
 
         Parameters
         ----------
-        n : int
+        n
             The shell number
 
         Returns
         -------
-        origins : ((2*n + 1)**3 - (2*n - 1)**3, 3) int ndarray
-            The cell origins. Note: if n = 0, origins = [[0, 0, 0]]
-
+        origins
+            Shape ((2*n + 1)**3 - (2*n - 1)**3, 3) int ndarray. The
+            cell origins. Note: if n = 0, origins = [[0, 0, 0]]
         """
 
         if n == 0:
@@ -1175,9 +1309,9 @@ def _get_shell_origins(self, n):
         # the edges, the xz plane has (2*n + 1) - 2 rows in z, rather
         # than (2*n + 1). The yz plane also has (2*n + 1) - 2 rows in z
         # and (2*n + 1) - 2 columns in y
-        xy = get_all_origins([n+1, n+1, 1], min_xyz=[-n, -n, 0])
-        xz = get_all_origins([n+1, 1, n], min_xyz=[-n, 0, -n+1])
-        yz = get_all_origins([1, n, n], min_xyz=[0, -n+1, -n+1])
+        xy = get_all_origins((n+1, n+1, 1), min_xyz=(-n, -n, 0))
+        xz = get_all_origins((n+1, 1, n), min_xyz=(-n, 0, -n+1))
+        yz = get_all_origins((1, n, n), min_xyz=(0, -n+1, -n+1))
 
         # Offset each plane by n and -n to get the 6 planes that make up
         # the shell
@@ -1194,7 +1328,7 @@ def _get_shell_origins(self, n):
 
         return origins
 
-    def _enforce_realspace_asr(self):
+    def _enforce_realspace_asr(self) -> np.ndarray:
         """
         Apply a transformation to the force constants matrix so that it
         satisfies the acousic sum rule. Diagonalise, shift the acoustic
@@ -1205,8 +1339,9 @@ def _enforce_realspace_asr(self):
 
         Returns
         -------
-        force_constants : (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float ndarray
-            The corrected force constants matrix
+        force_constants
+            Shape (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float ndarray. The
+            corrected force constants matrix
         """
         n_cells_in_sc = self.n_cells_in_sc
         n_atoms = self.crystal.n_atoms
@@ -1243,7 +1378,8 @@ def _enforce_realspace_asr(self):
 
         return fc
 
-    def _enforce_reciprocal_asr(self, dyn_mat_gamma):
+    def _enforce_reciprocal_asr(self, dyn_mat_gamma: np.ndarray
+            ) -> np.ndarray:
         """
         Calculate the correction to the dynamical matrix that would have
         to be applied to satisfy the acousic sum rule. Diagonalise the
@@ -1254,15 +1390,16 @@ def _enforce_reciprocal_asr(self, dyn_mat_gamma):
 
         Parameters
         ----------
-        dyn_mat_gamma : (3*n_atoms, 3*n_atoms) complex ndarray
-            The non mass-weighted dynamical matrix at gamma
+        dyn_mat_gamma
+            Shape (3*n_atoms, 3*n_atoms) complex ndarray. The non
+            mass-weighted dynamical matrix at gamma
 
         Returns
         -------
-        dyn_mat : (3*n_atoms, 3*n_atoms) complex ndarray or empty array
+        dyn_mat
+            Shape (3*n_atoms, 3*n_atoms) complex ndarray or empty array
             The corrected, non mass-weighted dynamical matrix at q.
-            Returns empty array (np.array([])) if finding the 3 acoustic
-            modes fails
+            Returns empty array if finding the 3 acoustic modes fails
         """
         tol = 5e-15
 
@@ -1284,7 +1421,8 @@ def _enforce_reciprocal_asr(self, dyn_mat_gamma):
 
         return recip_asr_correction
 
-    def _find_acoustic_modes(self, dyn_mat):
+    def _find_acoustic_modes(self, dyn_mat: np.ndarray
+            ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
         """
         Find the acoustic modes from a dynamical matrix, they should
         have the sum of c of m amplitude squared = mass (note: have not
@@ -1292,17 +1430,20 @@ def _find_acoustic_modes(self, dyn_mat):
 
         Parameters
         ----------
-        dyn_mat : (3*n_atoms, 3*n_atoms) complex ndarray
-            A dynamical matrix
+        dyn_mat
+            Shape (3*n_atoms, 3*n_atoms) complex ndarray. A dynamical
+            matrix
 
         Returns
         -------
-        ac_i : (3,) int ndarray
-            The indices of the acoustic modes
-        evals : (3*n_atoms) float ndarray
-            Dynamical matrix eigenvalues
-        evecs : (3*n_atoms, n_atoms, 3) complex ndarray
-            Dynamical matrix eigenvectors
+        ac_i
+            Shape (3,) int ndarray. The indices of the acoustic modes
+        evals
+            Shape (3*n_atoms) float ndarray. Dynamical matrix
+            eigenvalues
+        evecs
+            Shape (3*n_atoms, n_atoms, 3) complex ndarray. Dynamical
+            matrix eigenvectors
         """
         n_branches = dyn_mat.shape[0]
         n_atoms = int(n_branches/3)
@@ -1323,8 +1464,13 @@ def _find_acoustic_modes(self, dyn_mat):
 
         return ac_i, evals, evecs
 
-    def _calculate_phases(self, q, unique_sc_origins, unique_sc_i,
-                          unique_cell_origins, unique_cell_i):
+    def _calculate_phases(
+            self,
+            qpt: np.ndarray,
+            unique_sc_origins: Sequence[Sequence[int]],
+            unique_sc_i: np.ndarray,
+            unique_cell_origins: Sequence[Sequence[int]],
+            unique_cell_i: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
         """
         Calculate the phase factors for the supercell images and cells
         for a single q-point. The unique supercell and cell origins
@@ -1333,36 +1479,39 @@ def _calculate_phases(self, q, unique_sc_origins, unique_sc_i,
 
         Parameters
         ----------
-        q : (3,) float ndarray
-            The q-point to calculate the phase for
-        unique_sc_origins : list of lists of ints
+        qpt
+            Shape (3,) float ndarray. The q-point to calculate
+        unique_sc_origins
             A list containing 3 lists of the unique supercell image
-            offsets in each direction in units of the unit cell vectors.
-            The supercell offset is calculated by multiplying the
-            supercell matrix by the supercell image indices
-            (obtained by _get_all_origins()). A list of lists rather
-            than a Numpy array is used as the 3 lists are independent
-            and their size is not known beforehand
-        unique_sc_i : ((2*lim + 1)**3, 3) int ndarray
-            The indices needed to reconstruct sc_origins from the unique
-            values in unique_sc_origins
-        unique_cell_origins : list of lists of ints
-            A list containing 3 lists of the unique cell origins in each
-            direction in units of the unit cell vectors. A list of
-            lists rather than a Numpy array is used as the 3 lists are
-            independent and their size is not known beforehand
-        unique_cell_i : (cell_origins, 3) int ndarray
-            The indices needed to reconstruct cell_origins from the
-            unique values in unique_cell_origins
+            offsets in each direction. The supercell offset is
+            calculated by multiplying the supercell matrix by the
+            supercell image indices (obtained by _get_all_origins()). A
+            list of lists rather than a Numpy array is used as the 3
+            lists are independent and their size is not known
+            beforehand
+        unique_sc_i
+            Shape ((2*n_sc_shells + 1)**3, 3) int ndarray, where
+            n_sc_shells is how many supercells out to go to search for
+            atom images. The indices needed to reconstruct sc_origins
+            from the unique values in unique_sc_origins
+        unique_cell_origins
+            A list containing 3 lists of the unique cell origins in
+            each direction. A list of lists rather than a Numpy array
+            is used as the 3 lists are independent and their size is
+            not known beforehand
+        unique_cell_i
+            Shape (n_cells_in_sc, 3) int ndarray. The indices needed to
+            reconstruct cell_origins from the unique values in
+            unique_cell_origins
 
         Returns
         -------
-        sc_phases : (unique_sc_i,) float ndarray
-            Phase factors exp(iq.r) for each supercell image coordinate
-            in sc_origins
-        cell_phases : (unique_cell_i,) float ndarray
-            Phase factors exp(iq.r) for each cell coordinate in the
-            supercell
+        sc_phases
+            Shape (unique_sc_i,) float ndarray. Phase factors exp(iq.r)
+            for each supercell image coordinate in sc_origins
+        cell_phases
+            Shape (unique_cell_i,) float ndarray. Phase factors
+            exp(iq.r) for each cell coordinate in the supercell
         """
 
         # Only calculate exp(iq) once, then raise to power to get the
@@ -1370,7 +1519,7 @@ def _calculate_phases(self, q, unique_sc_origins, unique_sc_i,
         # expensive exp calculations
         # exp(iq.r) = exp(iqh.ra)*exp(iqk.rb)*exp(iql.rc)
         #           = (exp(iqh)^ra)*(exp(iqk)^rb)*(exp(iql)^rc)
-        phase = np.exp(2j*math.pi*q)
+        phase = np.exp(2j*math.pi*qpt)
         sc_phases = np.ones(len(unique_sc_i), dtype=np.complex128)
         cell_phases = np.ones(len(unique_cell_i), dtype=np.complex128)
         for i in range(3):
@@ -1382,17 +1531,18 @@ def _calculate_phases(self, q, unique_sc_origins, unique_sc_i,
 
         return sc_phases, cell_phases
 
-    def _calculate_supercell_images(self, lim):
+    def _calculate_supercell_images(self, n_sc_shells: int) -> None:
         """
-        For each displacement of ion i in the unit cell and ion j in the
-        supercell, calculate the number of supercell periodic images
-        there are and which supercells they reside in, and sets the
-        sc_image_i, and n_sc_images ForceConstants attributes
+        For each displacement of ion i in the unit cell and ion j in
+        the supercell, calculate the number of supercell periodic
+        images there are and which supercells they reside in, and sets
+        the _sc_image_i, and _n_sc_images ForceConstants attributes
 
         Parameters
         ----------
-        lim : int
-            The supercell image limit
+        n_sc_shells
+            The number of supercells out to go to search for atom
+            images
         """
 
         n_atoms = self.crystal.n_atoms
@@ -1418,14 +1568,15 @@ def _calculate_supercell_images(self, lim):
 
         # Get Cartesian coords of supercell images and ions in supercell
         sc_image_r = get_all_origins(
-            np.repeat(lim, 3) + 1, min_xyz=-np.repeat(lim, 3))
+            np.repeat(n_sc_shells, 3) + 1, min_xyz=-np.repeat(n_sc_shells, 3))
         sc_image_cart = np.einsum('ij,jk->ik', sc_image_r, sc_vecs)
         sc_atom_cart = np.einsum('ijk,kl->ijl',
                                  cell_origins[:, ax, :] + atom_r[ax, :, :],
                                  cell_vectors)
 
-        sc_image_i = np.full((n_cells_in_sc, n_atoms, n_atoms, (2*lim + 1)**3),
-                             -1, dtype=np.int32)
+        sc_image_i = np.full(
+            (n_cells_in_sc, n_atoms, n_atoms, (2*n_sc_shells + 1)**3),
+            -1, dtype=np.int32)
         n_sc_images = np.zeros((n_cells_in_sc, n_atoms, n_atoms),
                                dtype=np.int32)
 
@@ -1474,7 +1625,7 @@ def _calculate_supercell_images(self, lim):
         # nonexistent images
         self._sc_image_i = sc_image_i[:, :, :, :np.max(n_sc_images)]
 
-    def to_dict(self):
+    def to_dict(self) -> Dict[str, Any]:
         """
         Convert to a dictionary. See ForceConstants.from_dict for
         details on keys/values
@@ -1489,26 +1640,26 @@ def to_dict(self):
                                    'dielectric'])
         return dout
 
-    def to_json_file(self, filename):
+    def to_json_file(self, filename: str) -> None:
         """
         Write to a JSON file. JSON fields are equivalent to
         ForceConstants.from_dict keys
 
         Parameters
         ----------
-        filename : str
+        filename
             Name of the JSON file to write to
         """
         _obj_to_json_file(self, filename)
 
     @classmethod
-    def from_dict(cls, d):
+    def from_dict(cls: Type[T], d: Dict[str, Any]) -> T:
         """
         Convert a dictionary to a ForceConstants object
 
         Parameters
         ----------
-        d : dict
+        d
             A dictionary with the following keys/values:
 
             - 'crystal': dict, see Crystal.from_dict
@@ -1526,7 +1677,7 @@ def from_dict(cls, d):
 
         Returns
         -------
-        ForceConstants
+        forceconstants
         """
         crystal = Crystal.from_dict(d['crystal'])
         d = _process_dict(
@@ -1536,43 +1687,46 @@ def from_dict(cls, d):
                    d['cell_origins'], d['born'], d['dielectric'])
 
     @classmethod
-    def from_json_file(cls, filename):
+    def from_json_file(cls: Type[T], filename: str) -> T:
         """
         Read from a JSON file. See ForceConstants.from_dict for required
         fields
 
         Parameters
         ----------
-        filename : str
+        filename
             The file to read from
 
         Returns
         -------
-        ForceConstants
+        forceconstants
         """
         return _obj_from_json_file(cls, filename)
 
     @classmethod
-    def from_castep(cls, filename):
+    def from_castep(cls: Type[T], filename: str) -> T:
         """
         Reads from a .castep_bin or .check file
 
         Parameters
         ----------
-        filename : str
+        filename
             The path and name of the file to read
 
         Returns
         -------
-        ForceConstants
+        forceconstants
         """
         data = castep.read_interpolation_data(filename)
         return cls.from_dict(data)
 
     @classmethod
-    def from_phonopy(cls, path='.', summary_name='phonopy.yaml',
-                     born_name=None, fc_name='FORCE_CONSTANTS',
-                     fc_format=None):
+    def from_phonopy(cls: Type[T],
+                     path: str = '.',
+                     summary_name: str = 'phonopy.yaml',
+                     born_name: Optional[str] = None,
+                     fc_name: str = 'FORCE_CONSTANTS',
+                     fc_format: Optional[str] = None) -> T:
         """
         Reads data from the phonopy summary file (default phonopy.yaml)
         and optionally born and force constants files. Only attempts to
@@ -1581,27 +1735,28 @@ def from_phonopy(cls, path='.', summary_name='phonopy.yaml',
 
         Parameters
         ----------
-        path : str, optional
+        path
             Path to directory containing the file(s)
-        summary_name : str, optional
+        summary_name
             Filename of phonopy summary file, default phonopy.yaml. By
             default any information (e.g. force constants) read from
             this file takes priority
-        born_name : str, optional
+        born_name
             Name of the Phonopy file containing born charges and
             dielectric tensor (by convention in Phonopy this would be
             called BORN). Is only read if Born charges can't be found in
             the summary_name file
-        fc_name : str, optional
+        fc_name
             Name of file containing force constants. Is only read if
             force constants can't be found in summary_name
-        fc_format : {'phonopy', 'hdf5'} str, optional
-            Format of file containing force constants data.
-            FORCE_CONSTANTS is type 'phonopy'
+        fc_format
+            One of {'phonopy', 'hdf5'}. Format of file containing force
+            constants data, if it can't be determined from the
+            extension. FORCE_CONSTANTS is type 'phonopy'
 
         Returns
         -------
-        ForceConstants
+        forceconstants
         """
         data = phonopy.read_interpolation_data(
             path=path, summary_name=summary_name, born_name=born_name,
diff --git a/euphonic/io.py b/euphonic/io.py
index c5de5ede4..c4d0366fd 100644
--- a/euphonic/io.py
+++ b/euphonic/io.py
@@ -1,13 +1,15 @@
 import copy
 import json
 import os
+from typing import Dict, Type, TypeVar, Sequence, Any
 
 import numpy as np
 
 from euphonic import ureg, Quantity, __version__
 
+T = TypeVar('T')
 
-def _to_json_dict(dictionary):
+def _to_json_dict(dictionary: Dict[str, Any]) -> Dict[str, Any]:
     """
     Convert all keys in an output dictionary to a JSON serialisable
     format
@@ -23,7 +25,8 @@ def _to_json_dict(dictionary):
     return dictionary
 
 
-def _from_json_dict(dictionary, type_dict={}):
+def _from_json_dict(dictionary: Dict[str, Any],
+                    type_dict: Dict[str, Type[Any]] = {}) -> Dict[str, Any]:
     """
     For a dictionary read from a JSON file, convert all list key values
     to that specified in type_dict. If not specified in type_dict, just
@@ -45,7 +48,7 @@ def _from_json_dict(dictionary, type_dict={}):
     return dictionary
 
 
-def _obj_to_dict(obj, attrs):
+def _obj_to_dict(obj: T, attrs: Sequence[str]) -> Dict[str, Any]:
     dout = {}
     for attr in attrs:
         val = getattr(obj, attr)
@@ -73,7 +76,8 @@ def _obj_to_dict(obj, attrs):
     return dout
 
 
-def _process_dict(dictionary, quantities={}, optional={}):
+def _process_dict(dictionary: Dict[str, Any], quantities: Sequence[str] = [],
+                  optional: Sequence[str] = []) -> Dict[str, Any]:
     """
     Process an input dictionary for creating objects. Convert keys in
     'quantities' to Quantity objects, and if any 'optional' keys are
@@ -92,7 +96,7 @@ def _process_dict(dictionary, quantities={}, optional={}):
     return dictionary
 
 
-def _obj_to_json_file(obj, filename):
+def _obj_to_json_file(obj: T, filename: str) -> None:
     """
     Generic function for writing to a JSON file from a Euphonic object
     """
@@ -104,7 +108,8 @@ def _obj_to_json_file(obj, filename):
     print(f'Written to {os.path.realpath(f.name)}')
 
 
-def _obj_from_json_file(cls, filename, type_dict={}):
+def _obj_from_json_file(cls: Type[T], filename: str,
+                        type_dict: Dict[str, Type[Any]] = {}) -> T:
     """
     Generic function for reading from a JSON file to a Euphonic object
     """
diff --git a/euphonic/plot.py b/euphonic/plot.py
index 073fe4cca..dc82320d0 100644
--- a/euphonic/plot.py
+++ b/euphonic/plot.py
@@ -160,7 +160,7 @@ def plot_1d(spectra: Union[Spectrum1D,
 
     """
     if isinstance(spectra, (Spectrum1D, Spectrum1DCollection)):
-        spectra = [spectra]
+        spectra = (spectra,)
     else:
         # Check units are consistent
         for spectrum in spectra[1:]:
@@ -295,7 +295,7 @@ def _get_q_range(data: Quantity) -> float:
 
     intensity_unit = spectra[0].z_data.units
 
-    def _get_minmax_intensity(spectrum: Spectrum2D) -> float:
+    def _get_minmax_intensity(spectrum: Spectrum2D) -> Tuple[float, float]:
         dimensionless_data = spectrum.z_data.to(intensity_unit).magnitude
         assert isinstance(dimensionless_data, np.ndarray)
         return np.min(dimensionless_data), np.max(dimensionless_data)
@@ -322,7 +322,7 @@ def _get_minmax_intensity(spectrum: Spectrum2D) -> float:
     return fig
 
 def _set_x_tick_labels(ax: Axes,
-                       x_tick_labels: List[Tuple[int, str]],
+                       x_tick_labels: Optional[Sequence[Tuple[int, str]]],
                        x_data: Quantity) -> None:
     if x_tick_labels is not None:
         locs, labels = [list(x) for x in zip(*x_tick_labels)]
diff --git a/euphonic/powder.py b/euphonic/powder.py
index c06946c2e..141eeb6d1 100644
--- a/euphonic/powder.py
+++ b/euphonic/powder.py
@@ -204,7 +204,7 @@ def sample_sphere_structure_factor(
     energy_bins: Quantity = None,
     scattering_lengths: Union[str, Dict[str, Quantity]] = 'Sears1992',
     **calc_modes_args
-) -> Spectrum1D:
+    ) -> Spectrum1D:
     """Sample structure factor, averaging over a sphere of constant |q|
 
     (Specifically, this is the one-phonon inelastic-scattering structure
@@ -315,7 +315,7 @@ def sample_sphere_structure_factor(
     return s.calculate_1d_average(energy_bins)
 
 
-def _get_default_bins(phonons: QpointPhononModes,
+def _get_default_bins(phonons: Union[QpointPhononModes, QpointFrequencies],
                       nbins: int = 1000) -> Quantity:
     """Get a default set of energy bin edges for set of phonon modes"""
     max_energy = np.max(phonons.frequencies) * 1.05
diff --git a/euphonic/qpoint_frequencies.py b/euphonic/qpoint_frequencies.py
index 3b793dda1..b25af8c02 100644
--- a/euphonic/qpoint_frequencies.py
+++ b/euphonic/qpoint_frequencies.py
@@ -1,5 +1,5 @@
 import warnings
-from typing import Any, Dict, List, Optional, Tuple, TypeVar
+from typing import Any, Dict, List, Optional, Tuple, TypeVar, Type
 
 import numpy as np
 
@@ -12,28 +12,26 @@
                       Spectrum1DCollection, Spectrum2D)
 
 
-T = TypeVar('T', bound='QpointFrequencies')
-
-
 class QpointFrequencies:
     """
     A class to read and store frequency data at q-points
 
     Attributes
     ----------
-    crystal : Crystal
+    crystal
         Lattice and atom information
-    n_qpts : int
+    n_qpts
         Number of q-points in the object
     qpts
         Shape (n_qpts, 3) float ndarray. Q-point coordinates, in
         fractional coordinates of the reciprocal lattice
     frequencies
-        Shape (n_qpts, n_branches) float Quantity. Frequencies
-        per q-point and mode
+        Shape (n_qpts, n_branches) float Quantity in energy units.
+        Frequencies per q-point and mode
     weights
         Shape (n_qpts,) float ndarray. The weight for each q-point
     """
+    T = TypeVar('T', bound='QpointFrequencies')
 
     def __init__(self, crystal: Crystal, qpts: np.ndarray,
                  frequencies: Quantity,
@@ -46,8 +44,8 @@ def __init__(self, crystal: Crystal, qpts: np.ndarray,
         qpts
             Shape (n_qpts, 3) float ndarray. Q-point coordinates
         frequencies
-            Shape (n_qpts, n_branches) float Quantity. Frequencies
-            per q-point and mode
+            Shape (n_qpts, n_branches) float Quantity in energy units.
+            Frequencies per q-point and mode
         weights
             Shape (n_qpts,) float ndarray. The weight for each q-point.
             If None, equal weights are assumed
@@ -76,15 +74,15 @@ def __init__(self, crystal: Crystal, qpts: np.ndarray,
             self.weights = np.full(self.n_qpts, 1/self.n_qpts)
 
     @property
-    def frequencies(self):
+    def frequencies(self) -> Quantity:
         return self._frequencies*ureg('hartree').to(self.frequencies_unit)
 
     @frequencies.setter
-    def frequencies(self, value):
+    def frequencies(self, value: Quantity) -> None:
         self.frequencies_unit = str(value.units)
         self._frequencies = value.to('hartree').magnitude
 
-    def __setattr__(self, name, value):
+    def __setattr__(self, name: str, value: Any) -> None:
         _check_unit_conversion(self, name, value,
                                ['frequencies_unit'])
         super(QpointFrequencies, self).__setattr__(name, value)
@@ -99,8 +97,8 @@ def calculate_dos(self, dos_bins: Quantity,
         Parameters
         ----------
         dos_bins
-            Shape (n_e_bins + 1,) float Quantity. The energy bin edges
-            to use for calculating the DOS
+            Shape (n_e_bins + 1,) float Quantity in energy units. The
+            energy bin edges to use for calculating the DOS
         mode_widths
             Shape (n_qpts, n_branches) float Quantity in energy units.
             The broadening width for each mode at each q-point, for
@@ -193,6 +191,29 @@ def calculate_dos_map(self, dos_bins: Quantity,
                           mode_widths: Optional[Quantity] = None,
                           mode_widths_min: Quantity = Quantity(0.01, 'meV')
                           ) -> Spectrum2D:
+        """
+        Produces a bandstructure-like plot, using the DOS at each q-point
+
+        Parameters
+        ----------
+        dos_bins
+            Shape (n_e_bins + 1,) float Quantity in energy units. The
+            energy bin edges to use for calculating the DOS
+        mode_widths
+            Shape (n_qpts, n_branches) float Quantity in energy units.
+            The broadening width for each mode at each q-point, for
+            adaptive broadening
+        mode_widths_min
+            Scalar float Quantity in energy units. Sets a lower limit on
+            the mode widths, as mode widths of zero will result in
+            infinitely sharp peaks
+
+        Returns
+        -------
+        dos_map
+            A 2D spectrum containing the q-point bins on the x-axis, energy
+            bins on the y-axis and DOS on the z-axis
+        """
         x_data, x_tick_labels = self._get_qpt_axis_and_labels()
         dos_map = np.zeros((len(self.qpts), len(dos_bins) - 1))
         for i in range(len(self.qpts)):
@@ -228,6 +249,18 @@ def get_dispersion(self) -> Spectrum1DCollection:
 
     def _get_qpt_axis_and_labels(self
             ) -> Tuple[Quantity, List[Tuple[int, str]]]:
+        """
+        Converts the qpts stored in this object an array of
+        scalar distances and applies appropriate symmetry labels
+
+        Returns
+        -------
+        abscissa
+            Shape (n_qpts,) float quantity in 1/length units. The
+            scalar q-point distances
+        x_tick_labels
+            The tick labels
+        """
         abscissa = _calc_abscissa(self.crystal.reciprocal_cell(), self.qpts)
         # Calculate q-space ticks and labels
         x_tick_labels = get_qpoint_labels(
@@ -256,7 +289,7 @@ def to_json_file(self, filename: str) -> None:
         _obj_to_json_file(self, filename)
 
     @classmethod
-    def from_dict(cls: T, d: Dict[str, Any]) -> T:
+    def from_dict(cls: Type[T], d: Dict[str, Any]) -> T:
         """
         Convert a dictionary to a QpointFrequencies object
 
@@ -280,7 +313,7 @@ def from_dict(cls: T, d: Dict[str, Any]) -> T:
                    d['weights'])
 
     @classmethod
-    def from_json_file(cls: T, filename: str) -> T:
+    def from_json_file(cls: Type[T], filename: str) -> T:
         """
         Read from a JSON file. See from_dict for
         required fields
@@ -293,7 +326,7 @@ def from_json_file(cls: T, filename: str) -> T:
         return _obj_from_json_file(cls, filename)
 
     @classmethod
-    def from_castep(cls: T, filename: str) -> T:
+    def from_castep(cls: Type[T], filename: str) -> T:
         """
         Reads precalculated phonon mode data from a CASTEP .phonon file
 
@@ -306,7 +339,7 @@ def from_castep(cls: T, filename: str) -> T:
         return cls.from_dict(data)
 
     @classmethod
-    def from_phonopy(cls: T, path: str = '.',
+    def from_phonopy(cls: Type[T], path: str = '.',
                      phonon_name: str = 'band.yaml',
                      phonon_format: Optional[str] = None,
                      summary_name: str = 'phonopy.yaml') -> T:
diff --git a/euphonic/qpoint_phonon_modes.py b/euphonic/qpoint_phonon_modes.py
index fde40e976..bd0b0cc6b 100644
--- a/euphonic/qpoint_phonon_modes.py
+++ b/euphonic/qpoint_phonon_modes.py
@@ -1,5 +1,5 @@
 import math
-from typing import Dict, Optional, Union, TypeVar, Any
+from typing import Dict, Optional, Union, TypeVar, Any, Type
 from collections.abc import Mapping
 
 import numpy as np
@@ -13,9 +13,6 @@
                       StructureFactor, Spectrum1DCollection)
 
 
-T = TypeVar('T', bound='QpointPhononModes')
-
-
 class QpointPhononModes(QpointFrequencies):
     """
     A class to read and store vibrational data from model (e.g. CASTEP)
@@ -39,10 +36,11 @@ class QpointPhononModes(QpointFrequencies):
         Shape (n_qpts, 3*crystal.n_atoms, crystal.n_atoms, 3) complex
         ndarray. The dynamical matrix eigenvectors
     """
+    T = TypeVar('T', bound='QpointPhononModes')
 
     def __init__(self, crystal: Crystal, qpts: np.ndarray,
-                 frequencies: Quantity, eigenvectors: Quantity,
-                 weights: Optional[np.ndarray] = None):
+                 frequencies: Quantity, eigenvectors: np.ndarray,
+                 weights: Optional[np.ndarray] = None) -> None:
         """
         Parameters
         ----------
@@ -546,7 +544,7 @@ def to_qpoint_frequencies(self) -> QpointFrequencies:
             self.crystal, self.qpts, self.frequencies, self.weights)
 
     @classmethod
-    def from_dict(cls: T, d: Dict[str, Any]) -> T:
+    def from_dict(cls: Type[T], d: Dict[str, Any]) -> T:
         """
         Convert a dictionary to a QpointPhononModes object
 
@@ -571,7 +569,7 @@ def from_dict(cls: T, d: Dict[str, Any]) -> T:
                    d['eigenvectors'], d['weights'])
 
     @classmethod
-    def from_json_file(cls: T, filename: str) -> T:
+    def from_json_file(cls: Type[T], filename: str) -> T:
         """
         Read from a JSON file. See QpointPhononModes.from_dict for
         required fields
@@ -585,7 +583,7 @@ def from_json_file(cls: T, filename: str) -> T:
                                    type_dict={'eigenvectors': np.complex128})
 
     @classmethod
-    def from_castep(cls: T, filename: str) -> T:
+    def from_castep(cls: Type[T], filename: str) -> T:
         """
         Reads precalculated phonon mode data from a CASTEP .phonon file
 
@@ -598,7 +596,7 @@ def from_castep(cls: T, filename: str) -> T:
         return cls.from_dict(data)
 
     @classmethod
-    def from_phonopy(cls: T, path: str = '.',
+    def from_phonopy(cls: Type[T], path: str = '.',
                      phonon_name: str = 'band.yaml',
                      phonon_format: Optional[str] = None,
                      summary_name: str = 'phonopy.yaml') -> T:
diff --git a/euphonic/readers/castep.py b/euphonic/readers/castep.py
index eb0417b26..23680ffb7 100644
--- a/euphonic/readers/castep.py
+++ b/euphonic/readers/castep.py
@@ -1,6 +1,6 @@
 import re
 import struct
-from typing import Dict, Any, TextIO, Tuple, Optional, List, Union
+from typing import Dict, Any, TextIO, BinaryIO, Tuple, Optional, List, Union
 
 import numpy as np
 
@@ -75,8 +75,7 @@ def read_phonon_dos_data(
             dos_data = np.array([[float(elem) for elem in line.split()]
                                       for line in data[:n_bins]])
 
-    data_dict = {}
-
+    data_dict: Dict[str, Any] = {}
     data_dict['crystal'] = {}
     cry_dict = data_dict['crystal']
     cry_dict['n_atoms'] = n_atoms
@@ -179,7 +178,8 @@ def read_phonon_data(
             except EOFError:
                 break
             [f.readline() for x in range(2)]  # Skip 2 label lines
-            evec_lines = [f.readline() for x in range(n_atoms*n_branches)]
+            evec_lines = np.array(
+                [f.readline() for x in range(n_atoms*n_branches)])
             if read_eigenvectors:
                 evec_lines = np.array([x.split()[2:] for x in evec_lines],
                                       dtype=np.float64)
@@ -207,7 +207,7 @@ def read_phonon_data(
                     eigenvecs = np.concatenate((eigenvecs, [qeigenvec]))
             idx += 1
 
-    data_dict = {}
+    data_dict: Dict[str, Any] = {}
     data_dict['crystal'] = {}
     cry_dict = data_dict['crystal']
     cry_dict['n_atoms'] = n_atoms
@@ -231,7 +231,7 @@ def read_phonon_data(
 
 
 def _read_crystal_info(f: TextIO, n_atoms: int
-    ) -> Tuple[int, int, int, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
     """
     Reads the header crystal information from a CASTEP text file, from
     'Unit cell vectors' to 'END header'
@@ -274,7 +274,7 @@ def _read_crystal_info(f: TextIO, n_atoms: int
 
 def _read_frequency_block(
             f: TextIO, n_branches: int, extra_columns: Optional[List] = None
-    ) -> Tuple[np.ndarray, float, np.ndarray, Union[None, np.ndarray]]:
+    ) -> Tuple[np.ndarray, float, np.ndarray, Optional[np.ndarray]]:
     """
     For a single q-point reads the q-point, weight, frequencies
     and optionally any extra columns
@@ -310,23 +310,22 @@ def _read_frequency_block(
         raise EOFError
     float_patt = re.compile('-?\d+\.\d+')
     floats = re.findall(float_patt, qpt_line)
-    qpt = [float(x) for x in floats[:3]]
+    qpt = np.array([float(x) for x in floats[:3]])
     qweight = float(floats[3])
     freq_lines = [f.readline().split()
                   for i in range(n_branches)]
     freq_col = 1
     qfreq = np.array([float(line[freq_col]) for line in freq_lines])
+    extra: Optional[np.ndarray] = None
     if extra_columns is not None:
         if len(qpt_line.split()) > 6:
             # Indicates a split gamma point in .phonon, so there is an
             # additional column
-            extra_columns += 1
+            extra_columns = [x + 1 for x in extra_columns]
         extra = np.zeros((len(extra_columns), n_branches))
         for i, col in enumerate(extra_columns):
             extra[i] = np.array(
                 [float(line[freq_col + col + 1]) for line in freq_lines])
-    else:
-        extra = None
     return qpt, qweight, qfreq, extra
 
 
@@ -401,7 +400,7 @@ def read_interpolation_data(
                 dielectric = np.transpose(np.reshape(
                     _read_entry(f, float_type), (3, 3)))
 
-    data_dict = {}
+    data_dict: Dict[str, Any] = {}
     data_dict['crystal'] = {}
     cry_dict = data_dict['crystal']
     cry_dict['n_atoms'] = n_atoms
@@ -441,34 +440,38 @@ def read_interpolation_data(
     return data_dict
 
 
-def _read_cell(file_obj, int_type, float_type):
+def _read_cell(file_obj: BinaryIO, int_type: str, float_type: str
+               ) -> Tuple[int, np.ndarray, np.ndarray,
+                          np.ndarray, np.ndarray]:
     """
     Read cell data from a .castep_bin or .check file
 
     Parameters
     ----------
-    f : file object
+    file_obj
         File object in read mode for the .castep_bin or .check file
-    int_type : str
+    int_type
         Python struct format string describing the size and endian-ness
         of ints in the file
-    float_type : str
+    float_type
         Python struct format string describing the size and endian-ness
         of floats in the file
 
-        Returns
+    Returns
     -------
-    n_atoms : int
+    n_atoms
         Number of atoms in the unit cell
-    cell_vectors : (3, 3) float ndarray
-        The unit cell vectors in bohr
-    atom_r : (n_atoms, 3) float ndarray
-        The fractional position of each atom within the unit cell
-    atom_mass : (n_atoms,) float ndarray
-        The mass of each atom in the unit cell in units of electron mass
-    atom_type : (n_atoms,) string ndarray
-        The chemical symbols of each atom in the unit cell. Atoms are in
-        the same order as in atom_r
+    cell_vectors
+        Shape (3, 3) float ndarray. The unit cell vectors in bohr
+    atom_r
+        Shape (n_atoms, 3) float ndarray. The fractional position of
+        each atom within the unit cell
+    atom_mass
+        Shape (n_atoms,) float ndarray. The mass of each atom in the
+        unit cell in units of electron mass
+    atom_type
+        Shape (n_atoms,) string ndarray. The chemical symbols of each
+        atom in the unit cell
     """
     header = ''
     while header.strip() != b'END_UNIT_CELL':
@@ -524,22 +527,23 @@ def _read_cell(file_obj, int_type, float_type):
     return n_atoms, cell_vectors, atom_r, atom_mass, atom_type
 
 
-def _read_entry(file_obj, dtype=''):
+def _read_entry(file_obj: BinaryIO, dtype: str = ''
+                ) -> Union[str, int, float, np.ndarray]:
     """
     Read a record from a Fortran binary file, including the beginning
     and end record markers and return the data inbetween
 
     Parameters
     ----------
-    f : file object
+    f
         File object in read mode for the Fortran binary file
-    dtype : str, optional, default ''
+    dtype
         String determining what order and type to unpack the bytes as.
         See 'Format Strings' in Python struct documentation
 
     Returns
     -------
-    data : str, int, float or ndarray
+    data
         Data type returned depends on dtype specified. If dtype is not
         specified, return type is a string. If there is more than one
         element in the record, it is returned as an ndarray of floats or
diff --git a/euphonic/readers/phonopy.py b/euphonic/readers/phonopy.py
index 63c68216f..96d012792 100644
--- a/euphonic/readers/phonopy.py
+++ b/euphonic/readers/phonopy.py
@@ -1,6 +1,6 @@
 import os
 import warnings
-from typing import Optional, Dict, Any
+from typing import Optional, Dict, Any, Union, Tuple, TextIO, Sequence, List
 
 import numpy as np
 
@@ -23,19 +23,19 @@ def __str__(self):
         return self.message
 
 
-def _convert_weights(weights):
+def _convert_weights(weights: np.ndarray) -> np.ndarray:
     """
     Convert q-point weights to normalised convention
 
     Parameters
     ----------
-    weights : (n_qpts,) float ndarray
-        Weights in Phonopy convention
+    weights
+        Shape (n_qpts,) float ndarray. Weights in Phonopy convention
 
     Returns
     ----------
-    norm_weights : (n_qpts,) float ndarray
-        Normalised weights
+    norm_weights
+        Shape (n_qpts,) float ndarray. Normalised weights
     """
     total_weight = weights.sum()
     return weights/total_weight
@@ -43,7 +43,7 @@ def _convert_weights(weights):
 
 def _extract_phonon_data_yaml(filename: str,
                               read_eigenvectors: bool = True
-    ) -> Dict[str, Any]:
+                              ) -> Dict[str, np.ndarray]:
     """
     From a mesh/band/qpoint.yaml file, extract the relevant information
     as a dict of Numpy arrays. No unit conversion is done at this point,
@@ -119,7 +119,7 @@ def _extract_phonon_data_yaml(filename: str,
 
 def _extract_phonon_data_hdf5(filename: str,
                               read_eigenvectors: bool = True
-    ) -> Dict[str, Any]:
+                              ) -> Dict[str, np.ndarray]:
     """
     From a mesh/band/qpoint.hdf5 file, extract the relevant information
     as a dict of Numpy arrays. No unit conversion is done at this point,
@@ -148,8 +148,8 @@ def _extract_phonon_data_hdf5(filename: str,
         raise ImportPhonopyReaderError from e
 
     with h5py.File(filename, 'r') as hdf5_file:
+        data_dict = {}
         if 'qpoint' in hdf5_file.keys():
-            data_dict = {}
             data_dict['qpts'] = hdf5_file['qpoint'][()]
             data_dict['frequencies'] = hdf5_file['frequency'][()]
             if read_eigenvectors:
@@ -167,33 +167,23 @@ def _extract_phonon_data_hdf5(filename: str,
         # Is a band.hdf5 file - q-points are stored in 'path' and need
         # special treatment
         else:
-            data_dict = _extract_band_data_hdf5(hdf5_file)
-
-    return data_dict
-
+            data_dict['qpts'] = hdf5_file['path'][()].reshape(
+                -1, hdf5_file['path'][()].shape[-1])
+            data_dict['frequencies'] = hdf5_file['frequency'][()].reshape(
+                -1, hdf5_file['frequency'][()].shape[-1])
+            try:
+                # The last 2 dimensions of eigenvectors in bands.hdf5 are for
+                # some reason transposed compared to mesh/qpoints.hdf5, so also
+                # transpose to handle this
+                data_dict['eigenvectors'] = hdf5_file[
+                    'eigenvector'][()].reshape(
+                        -1, *hdf5_file['eigenvector'][()].shape[-2:]
+                        ).transpose([0,2,1])
+            except KeyError:
+                pass
 
-def _extract_band_data_hdf5(hdf5_file):
-    """
-    Read data from a Phonopy band.hdf5 file. All information is stored
-    in the shape (n_paths, n_qpoints_in_path, x) rather than (n_qpts, x)
-    so needs special treatment
-    """
-    data_dict = {}
-    data_dict['qpts'] = hdf5_file['path'][()].reshape(
-        -1, hdf5_file['path'][()].shape[-1])
-    data_dict['frequencies'] = hdf5_file['frequency'][()].reshape(
-        -1, hdf5_file['frequency'][()].shape[-1])
-    try:
-        # The last 2 dimensions of eigenvectors in bands.hdf5 are for
-        # some reason transposed compared to mesh/qpoints.hdf5, so also
-        # transpose to handle this
-        data_dict['eigenvectors'] = hdf5_file['eigenvector'][()].reshape(
-            -1, *hdf5_file['eigenvector'][()].shape[-2:]).transpose([0,2,1])
-    except KeyError:
-        pass
     return data_dict
 
-
 def read_phonon_data(
         path: str = '.',
         phonon_name: str = 'band.yaml',
@@ -202,7 +192,8 @@ def read_phonon_data(
         cell_vectors_unit: str = 'angstrom',
         atom_mass_unit: str = 'amu',
         frequencies_unit: str = 'meV',
-        read_eigenvectors: bool = True) -> Dict[str, Any]:
+        read_eigenvectors: bool = True
+        ) -> Dict[str, Union[int, str, np.ndarray]]:
     """
     Reads precalculated phonon mode data from a Phonopy
     mesh/band/qpoints.yaml/hdf5 file and returns it in a dictionary. May
@@ -294,7 +285,7 @@ def read_phonon_data(
                 f'with summary file {summary_pathname}. Please '
                 'check contents'))
 
-    data_dict = {}
+    data_dict: Dict[str, Any] = {}
     data_dict['crystal'] = {}
     cry_dict = data_dict['crystal']
     cry_dict['n_atoms'] = len(phonon_dict['atom_r'])
@@ -320,7 +311,8 @@ def read_phonon_data(
     return data_dict
 
 
-def convert_eigenvector_phases(phonon_dict):
+def convert_eigenvector_phases(phonon_dict: Dict[str, np.ndarray]
+                               ) -> np.ndarray:
     """
     When interpolating the force constants matrix, Euphonic uses a phase
     convention of e^iq.r_a, where r_a is the coordinate of each CELL in
@@ -343,37 +335,44 @@ def convert_eigenvector_phases(phonon_dict):
     return np.reshape(eigvecs, (n_qpts, 3*n_atoms, n_atoms, 3))
 
 
-def _extract_force_constants(fc_pathname, n_atoms, n_cells, summary_name,
-                             cell_origins_map, sc_relative_idx, p2s_map):
+def _extract_force_constants(fc_pathname: str, n_atoms: int, n_cells: int,
+                             summary_name: str, cell_origins_map: np.ndarray,
+                             sc_relative_idx: np.ndarray, p2s_map: np.ndarray,
+                             ) -> np.ndarray:
     """
     Reads force constants from a Phonopy FORCE_CONSTANTS file
 
     Parameters
     ----------
-    fc_pathname : str
+    fc_pathname
         The FORCE_CONSTANTS file to read from
-    n_atoms : int
+    n_atoms
         Number of atoms in the unit cell
-    n_cells : int
+    n_cells
         Number of unit cells in the supercell
-    cell_origins_map : (n_atoms*n_cells, 2) int ndarray
-        In the case of of non-diagonal supercell_matrices, the cell
-        origins are not the same for each atom. This is a map of the
-        equivalent cell origins for each atom, which is required to
-        reorder the force constants matrix so that all atoms in the unit
-        cell share equivalent origins.
-    sc_relative_idx : (n_cells, n_cells) int ndarray
-        The index n of the equivalent vector in cell_origins for each
-        cell_origins[i] -> cell_origins[j] vector in the supercell.
-        See _get_supercell_relative_idx
-    p2s_map : (n_atoms,) int ndarray
-        The index of the primitive atoms in the atoms of the supercell.
-        Used if fc_format = 'full'
+    summary_name
+        Name of input phonopy.yaml file
+    cell_origins_map
+        Shape (n_atoms*n_cells, 2) int ndarray. In the case of
+        non-diagonal supercell_matrices, the cell origins are not the
+        same for each atom. This is a map of the equivalent cell origins
+        for each atom, which is required to reorder the force constants
+        matrix so that all atoms in the unit cell share equivalent
+        origins
+    sc_relative_idx
+        Shape (n_cells, n_cells) int ndarray. The index n of the
+        equivalent vector in cell_origins for each
+        cell_origins[i] -> cell_origins[j] vector in the supercell. See
+        _get_supercell_relative_idx
+    p2s_map
+        Shape (n_atoms,) int ndarray. The index of the primitive atoms
+        in the atoms of the supercell. Used if fc_format = 'full'
 
     Returns
     -------
-    fc : (n_cells, 3*n_atoms, 3*n_atoms) float ndarray
-        The force constants, in Euphonic convention
+    fc
+        Shape (n_cells, 3*n_atoms, 3*n_atoms) float ndarray. The force
+        constants, in Euphonic convention
     """
     with open(fc_pathname, 'r') as f:
         fc_dims =  [int(dim) for dim in f.readline().split()]
@@ -396,8 +395,10 @@ def _extract_force_constants(fc_pathname, n_atoms, n_cells, summary_name,
     return _reshape_fc(fc, n_atoms, n_cells, cell_origins_map, sc_relative_idx)
 
 
-def _extract_force_constants_hdf5(filename, n_atoms, n_cells, summary_name,
-                                  cell_origins_map, sc_relative_idx):
+def _extract_force_constants_hdf5(
+        filename: str, n_atoms: int, n_cells: int, summary_name: str,
+        cell_origins_map: np.ndarray, sc_relative_idx: np.ndarray
+        ) -> np.ndarray:
     try:
         import h5py
     except ModuleNotFoundError as e:
@@ -415,33 +416,39 @@ def _extract_force_constants_hdf5(filename, n_atoms, n_cells, summary_name,
                        sc_relative_idx)
 
 
-def _extract_force_constants_summary(summary_object, cell_origins_map,
-                                     sc_relative_idx, p2s_map):
+def _extract_force_constants_summary(
+        summary_object: Dict[str, Any],
+        cell_origins_map: np.ndarray,
+        sc_relative_idx: np.ndarray,
+        p2s_map: np.ndarray) -> np.ndarray:
     """
     Get force constants from phonopy yaml summary file.
 
     Parameters
     ----------
-    summary_object : dict
+    summary_object
         Dict containing contents of phonopy.yaml
-    cell_origins_map : (n_atoms*n_cells, 2) int ndarray
-        In the case of of non-diagonal supercell_matrices, the cell
-        origins are not the same for each atom. This is a map of the
-        equivalent cell origins for each atom, which is required to
-        reorder the force constants matrix so that all atoms in the unit
-        cell share equivalent origins.
-    sc_relative_idx : (n_cells, n_cells) int ndarray
-        The index n of the equivalent vector in cell_origins for each
+    cell_origins_map
+        Shape (n_atoms*n_cells, 2) int ndarray. In the case of of
+        non-diagonal supercell_matrices, the cell origins are not the
+        same for each atom. This is a map of the equivalent cell
+        origins for each atom, which is required to reorder the force
+        constants matrix so that all atoms in the unit cell share
+        equivalent origins
+    sc_relative_idx
+        Shape (n_cells, n_cells) int ndarray. The index n of the
+        equivalent vector in cell_origins for each
         cell_origins[i] -> cell_origins[j] vector in the supercell.
         See _get_supercell_relative_idx
-    p2s_map : (n_atoms,) int ndarray
-        The index of the primitive atoms in the atoms of the supercell.
-        Used if fc_format = 'full'
+    p2s_map
+        Shape (n_atoms,) int ndarray. The index of the primitive atoms
+        in the atoms of the supercell. Used if fc_format = 'full'
 
     Returns
-    ----------
-    units : dict
-        Dict containing force constants in Euphonic format.
+    -------
+    fc
+        Shape (n_cells, 3*n_atoms, 3*n_atoms) float ndarray. Force
+        constants in Euphonic format
     """
     fc_entry = summary_object['force_constants']
     fc_format = fc_entry['format']
@@ -457,33 +464,39 @@ def _extract_force_constants_summary(summary_object, cell_origins_map,
     fc = _reshape_fc(fc, n_atoms, n_cells, cell_origins_map, sc_relative_idx)
     return fc
 
-def _reshape_fc(fc, n_atoms, n_cells, cell_origins_map, sc_relative_idx):
+def _reshape_fc(fc: np.ndarray, n_atoms: int, n_cells: int,
+                cell_origins_map: np.ndarray, sc_relative_idx: np.ndarray
+                ) -> np.ndarray:
     """
     Reshape force constants from Phonopy convention to Euphonic
     convention
 
     Parameters
     ----------
-    fc : (n_atoms, n_atoms*n_cells, 3, 3) float ndarray
-        Force constants matrix in Phonopy shape
-    n_atoms : int
+    fc
+        Shape (n_atoms, n_atoms*n_cells, 3, 3) float ndarray. Force
+        constants matrix in Phonopy shape
+    n_atoms
         Number of atoms in the unit cell
-    n_cells : int
+    n_cells
         Number of cells in the supercell
-    cell_origins_map : (n_atoms*n_cells, 2) int ndarray
-        This describes how the supercell atoms in Phonopy map onto the
-        cells and atoms in Euphonic. The first column is the index of
-        the equivalent atom in the unit cell, and the second column is
-        the index of the cell in the supercell.
-    sc_relative_idx : (n_cells, n_cells) int ndarray
-        The index n of the equivalent vector in cell_origins for each
+    cell_origins_map
+        Shape (n_atoms*n_cells, 2) int ndarray. This describes how the
+        supercell atoms in Phonopy map onto the cells and atoms in
+        Euphonic. The first column is the index of the equivalent atom
+        in the unit cell, and the second column is the index of the
+        cell in the supercell
+    sc_relative_idx
+        Shape (n_cells, n_cells) int ndarray. The index n of the
+        equivalent vector in cell_origins for each
         cell_origins[i] -> cell_origins[j] vector in the supercell.
         See _get_supercell_relative_idx
 
     Returns
     -------
-    fc : (n_cells, 3*n_atoms, 3*n_atoms) float ndarray
-        Force constants matrix in Euphonic shape
+    fc
+        Shape (n_cells, 3*n_atoms, 3*n_atoms) float ndarray. Force
+        constants matrix in Euphonic shape
     """
     fc_phonopy = np.reshape(fc, (n_atoms, n_atoms*n_cells, 3, 3))
     fc_euphonic = np.full((n_atoms, n_cells, n_atoms, 3, 3), -1.0)
@@ -506,7 +519,9 @@ def _reshape_fc(fc, n_atoms, n_cells, cell_origins_map, sc_relative_idx):
         axes=[1, 0, 3, 2, 4]), (n_cells, 3*n_atoms, 3*n_atoms))
 
 
-def _check_fc_shape(fc_shape, n_atoms, n_cells, fc_filename, summary_filename):
+def _check_fc_shape(fc_shape: Tuple[int, int], n_atoms: int,
+                    n_cells: int, fc_filename: str,
+                    summary_filename: str) -> None:
     """
     Check if force constants has the correct shape
     """
@@ -519,68 +534,67 @@ def _check_fc_shape(fc_shape, n_atoms, n_cells, fc_filename, summary_filename):
             f' and {n_cells} cells in the supercell'))
 
 
-def _extract_born(born_object):
+def _extract_born(born_file_obj: TextIO) -> Dict[str, Union[float, np.ndarray]]:
     """
     Parse and convert dielectric tensor and born effective
     charge from BORN file
 
     Parameters
     ----------
-    born_object : file-like object
-        Object representing contents of BORN
+    born_file_obj
+        Pointer to BORN file
 
     Returns
     ----------
-    born_dict : float ndarray
+    born_dict
         Dict containing dielectric tensor and born effective charge
     """
-    born_lines_str = born_object.readlines()
-    born_lines = [narr.split() for narr in born_lines_str]
-
-    for l_i, line in enumerate(born_lines):
-        try:
-            born_lines[l_i] = [float(i) for i in line]
-            continue
-        except:
+    born_lines_str = born_file_obj.readlines()
+    born_lines: List[List[float]] = []
+    for line in born_lines_str:
+        # Ignore comments
+        if line.startswith('#'):
             pass
+        else:
+            born_lines.append([float(x) for x in line.split()])
 
-    n_lines = len(born_lines)
-    # factor first line
-    factor = born_lines[0][0]
+    born_dict: Dict[str, Union[float, np.ndarray]] = {}
 
-    # dielectric second line
-    # xx, xy, xz, yx, yy, yz, zx, zy, zz.
-    dielectric = np.array(born_lines[1]).reshape([3,3])
+    idx0 = 0
+    if len(born_lines[0]) == 1:
+        # Then this is the NAC conversion factor
+        born_dict['factor'] = born_lines[0][0]
+        idx0 = 1
 
-    # born ec third line onwards
+    # dielectric first line after factor
     # xx, xy, xz, yx, yy, yz, zx, zy, zz.
-    born_lines = born_lines[2:]
-    born = np.array([np.array(bl).reshape([3,3]) for bl in born_lines])
+    born_dict['dielectric'] = np.array(born_lines[idx0]).reshape([3,3])
 
-    born_dict = {}
-    born_dict['factor'] = factor
-    born_dict['dielectric'] = dielectric
-    born_dict['born'] = born
+    # born charges after dielectric
+    # xx, xy, xz, yx, yy, yz, zx, zy, zz.
+    born_dict['born'] = np.array(
+        [np.array(bl).reshape([3,3]) for bl in born_lines[idx0+1:]])
 
     return born_dict
 
 
-def _extract_summary(filename, fc_extract=False):
+def _extract_summary(filename: str, fc_extract: bool = False
+                     ) -> Dict[str, Union[str, int, np.ndarray]]:
     """
     Read phonopy.yaml for summary data produced during the Phonopy
-    post-process.
+    post-process
 
     Parameters
     ----------
-    filename : str
+    filename
         Path and name of the Phonopy summary file (usually phonopy.yaml)
-    fc_extract : bool, optional, default False
+    fc_extract
         Whether to attempt to read force constants and related
         information from summary_object
 
     Returns
     ----------
-    summary_dict : dict
+    summary_dict
         A dict with the following keys: n_atoms, cell_vectors, atom_r,
         atom_type, atom_mass, ulength, umass. Also optionally has the
         following keys: sc_matrix, n_cells_in_sc, cell_origins,
@@ -597,7 +611,6 @@ def _extract_summary(filename, fc_extract=False):
 
     with open(filename, 'r') as summary_file:
         summary_object = yaml.load(summary_file, Loader=SafeLoader)
-
     (cell_vectors, n_atoms, atom_r, atom_mass,
      atom_type, _) = _extract_crystal_data(summary_object['primitive_cell'])
 
@@ -658,7 +671,7 @@ def _extract_summary(filename, fc_extract=False):
         cell_origins_map = np.zeros((n_atoms*n_pcells_in_sc, 2),
                                     dtype=np.int32)
         # Get origins of adjacent supercells in prim cell frac coords
-        sc_origins =  get_all_origins([2,2,2], min_xyz=[-1,-1,-1])
+        sc_origins =  get_all_origins((2,2,2), min_xyz=(-1,-1,-1))
         sc_origins_pcell = np.einsum('ij,jk->ik', sc_origins, p_to_sc_matrix)
         for i in range(n_pcells_in_sc*n_atoms):
             co_idx = np.where(
@@ -709,30 +722,37 @@ def _extract_summary(filename, fc_extract=False):
     return summary_dict
 
 
-def _extract_crystal_data(crystal):
+def _extract_crystal_data(crystal: Dict[str, Any]
+                          ) -> Tuple[np.ndarray, int, np.ndarray,
+                                     np.ndarray, np.ndarray, np.ndarray]:
     """
     Gets relevant data from a section of phonopy.yaml
 
     Parameters
     ----------
-    crystal : dict
+    crystal
         Part of the dict obtained from reading a phonopy.yaml file. e.g.
         summary_dict['unit_cell']
 
     Returns
     -------
-    cell_vectors : (3,3) float ndarray
-        Cell vectors, in same units as in the phonopy.yaml file
-    n_atoms : int
+    cell_vectors
+        Shape (3,3) float ndarray. Cell vectors, in same units as in
+        the phonopy.yaml file
+    n_atoms
         Number of atoms
-    atom_r : (n_atoms, 3) float ndarray
-        Fractional position of each atom
-    atom_mass : (n_atoms,) float ndarray
-        Mass of each atom, in same units as in the phonopy.yaml file
-    atom_type : (n_atoms,) str ndarray
-        String specifying the species of each atom
-    idx_in_pcell : (n_atoms,) int ndarray
-        Maps the atom index on the unit/supercell to the primitive cell
+    atom_r
+        Shape (n_atoms, 3) float ndarray. Fractional position of each
+        atom
+    atom_mass
+        Shape (n_atoms,) float ndarray. Mass of each atom, in same
+        units as in the phonopy.yaml file
+    atom_type
+        Shape (n_atoms,) str ndarray. String specifying the species
+        of each atom
+    idx_in_pcell
+        Shape (n_atoms,) int ndarray. Maps the atom index on the
+        unit/supercell to the primitive cell
     """
     n_atoms = len(crystal['points'])
     cell_vectors = np.array(crystal['lattice'])
@@ -862,7 +882,7 @@ def read_interpolation_data(
     umass = summary_dict['umass']
     ufc = summary_dict['ufc']
 
-    data_dict = {}
+    data_dict: Dict[str, Any] = {}
     data_dict['crystal'] = {}
     cry_dict = data_dict['crystal']
     cry_dict['cell_vectors'] = summary_dict['cell_vectors']*ureg(
diff --git a/euphonic/spectra.py b/euphonic/spectra.py
index bb88d2446..11a5ad4a8 100644
--- a/euphonic/spectra.py
+++ b/euphonic/spectra.py
@@ -5,7 +5,7 @@
 import math
 from numbers import Integral
 from typing import (Any, Dict, List, Optional, overload,
-                    Sequence, Tuple, TypeVar, Union)
+                    Sequence, Tuple, TypeVar, Union, Type, Iterable)
 
 import numpy as np
 from scipy.ndimage import correlate1d, gaussian_filter
@@ -18,13 +18,9 @@
 from euphonic import ureg, Quantity
 
 
-S = TypeVar('S', bound='Spectrum')
-S1D = TypeVar('S1D', bound='Spectrum1D')
-SC = TypeVar('SC', bound='Spectrum1DCollection')
-S2D = TypeVar('S2D', bound='Spectrum2D')
-
-
 class Spectrum(ABC):
+    T = TypeVar('T', bound='Spectrum')
+
     def __setattr__(self, name: str, value: Any) -> None:
         _check_unit_conversion(self, name, value,
                                ['x_data_unit', 'y_data_unit'])
@@ -37,22 +33,22 @@ def _set_data(self, data: Quantity, attr_name: str) -> None:
         setattr(self, f'{attr_name}_data_unit', str(data.units))
 
     @property
-    def x_data(self):
+    def x_data(self) -> Quantity:
         return self._x_data*ureg(self._internal_x_data_unit).to(
             self.x_data_unit)
 
     @x_data.setter
-    def x_data(self, value):
+    def x_data(self, value: Quantity) -> None:
         self.x_data_unit = str(value.units)
         self._x_data = value.to(self._internal_x_data_unit).magnitude
 
     @property
-    def y_data(self):
+    def y_data(self) -> Quantity:
         return self._y_data*ureg(self._internal_y_data_unit).to(
             self.y_data_unit)
 
     @y_data.setter
-    def y_data(self, value):
+    def y_data(self, value: Quantity) -> None:
         self.y_data_unit = str(value.units)
         self._y_data = value.to(self._internal_y_data_unit).magnitude
 
@@ -63,7 +59,7 @@ def to_dict(self) -> Dict[str, Any]:
 
     @classmethod
     @abstractmethod
-    def from_dict(cls: S, d: Dict[str, Any]) -> S:
+    def from_dict(cls: Type[T], d: Dict[str, Any]) -> T:
         """Initialise a Spectrum object from dictionary"""
         ...
 
@@ -80,7 +76,7 @@ def to_json_file(self, filename: str) -> None:
         _obj_to_json_file(self, filename)
 
     @classmethod
-    def from_json_file(cls: S, filename: str) -> S:
+    def from_json_file(cls: Type[T], filename: str) -> T:
         """
         Read from a JSON file. See from_dict for required fields
 
@@ -93,12 +89,12 @@ def from_json_file(cls: S, filename: str) -> S:
         return _obj_from_json_file(cls, filename, type_dict)
 
     @abstractmethod
-    def _split_by_indices(self: S, indices: Union[Sequence[int], np.ndarray]
-                          ) -> List[S]:
+    def _split_by_indices(self: T, indices: Union[Sequence[int], np.ndarray]
+                          ) -> List[T]:
         """Split data along x axis at given indices"""
         ...
 
-    def _split_by_tol(self: S, btol: float = 10.0) -> List[S]:
+    def _split_by_tol(self: T, btol: float = 10.0) -> List[T]:
         """Split data along x-axis at detected breakpoints"""
         diff = np.diff(self.x_data.magnitude)
         median = np.median(diff)
@@ -107,7 +103,7 @@ def _split_by_tol(self: S, btol: float = 10.0) -> List[S]:
 
     @staticmethod
     def _ranges_from_indices(indices: Union[Sequence[int], np.ndarray]
-                             ) -> List[Tuple[int, int]]:
+                             ) -> List[Tuple[int, Optional[int]]]:
         """Convert a series of breakpoints to a series of slice ranges"""
         if len(indices) == 0:
             ranges = [(0, None)]
@@ -120,9 +116,10 @@ def _ranges_from_indices(indices: Union[Sequence[int], np.ndarray]
         return ranges
 
     @staticmethod
-    def _cut_x_ticks(x_tick_labels: Union[List[Tuple[int, str]], None],
+    def _cut_x_ticks(x_tick_labels: Union[Sequence[Tuple[int, str]], None],
                      x0: int,
-                     x1: Union[int, None]) -> List[Tuple[int, str]]:
+                     x1: Union[int, None]) -> Union[List[Tuple[int, str]],
+                                                    None]:
         """Crop and shift x labels to new x range"""
         if x_tick_labels is None:
             return None
@@ -130,8 +127,8 @@ def _cut_x_ticks(x_tick_labels: Union[List[Tuple[int, str]], None],
             return [(int(x - x0), label) for (x, label) in x_tick_labels
                     if (x >= x0) and ((x1 is None) or (x < x1))]
 
-    def split(self: S, indices: Union[Sequence[int], np.ndarray] = None,
-              btol: float = None) -> List[S]:
+    def split(self: T, indices: Union[Sequence[int], np.ndarray] = None,
+              btol: float = None) -> List[T]:
         """Split to multiple spectra
 
         Data may be split by index. Alternatively, x-axis data may be
@@ -141,9 +138,9 @@ def split(self: S, indices: Union[Sequence[int], np.ndarray] = None,
 
         Parameters
         ----------
-        indices:
+        indices
             positions in data of breakpoints
-        btol:
+        btol
             parameter used to identify breakpoints. This is a ratio
             between the gap in values and the median gap between
             neighbouring x-values. If neither indices nor btol is
@@ -151,10 +148,10 @@ def split(self: S, indices: Union[Sequence[int], np.ndarray] = None,
 
         Returns
         -------
+        split_spectra
             Separated spectrum regions. If passed to the appropriate
             functions in euphonic.plot this would be interpreted as a
             series of subplots.
-
         """
 
         if indices is None:
@@ -291,15 +288,19 @@ class Spectrum1D(Spectrum):
         == (n_x_data + 1,))
     y_data
         Shape (n_x_data,) float Quantity. The plot data in y
-    x_tick_labels : Sequence[Tuple[int, str]] or None
-        Special tick labels e.g. for high-symmetry points. The int
-        refers to the index in x_data the label should be applied to
-    metadata : Dict[str, Union[int, str]]
-        Contains metadata about the spectrum. Keys should be strings
-        and values should be strings or integers
+    x_tick_labels
+        Sequence[Tuple[int, str]] or None. Special tick labels e.g. for
+        high-symmetry points. The int refers to the index in x_data the
+        label should be applied to
+    metadata
+        Dict[str, Union[int, str]]. Contains metadata about the
+        spectrum. Keys should be strings and values should be strings
+        or integers
         There are some functional keys:
           - 'label' : str. This is used label lines on a 1D plot
     """
+    T = TypeVar('T', bound='Spectrum1D')
+
     def __init__(self, x_data: Quantity, y_data: Quantity,
                  x_tick_labels: Optional[Sequence[Tuple[int, str]]] = None,
                  metadata: Optional[Dict[str, Union[int, str]]] = None
@@ -336,7 +337,7 @@ def __init__(self, x_data: Quantity, y_data: Quantity,
         self.x_tick_labels = x_tick_labels
         self.metadata = {} if metadata is None else metadata
 
-    def __add__(self, other: S1D):
+    def __add__(self, other: 'Spectrum1D') -> 'Spectrum1D':
         """
         Sums the y_data of two Spectrum1D objects together,
         their x_data axes must be equal, and their y_data must
@@ -349,9 +350,9 @@ def __add__(self, other: S1D):
         spec_col = Spectrum1DCollection.from_spectra([self, other])
         return spec_col.sum()
 
-    def _split_by_indices(self: S1D,
+    def _split_by_indices(self: T,
                           indices: Union[Sequence[int], np.ndarray]
-                          ) -> List[S1D]:
+                          ) -> List[T]:
         """Split data along x-axis at given indices"""
         ranges = self._ranges_from_indices(indices)
 
@@ -374,7 +375,7 @@ def to_dict(self) -> Dict[str, Any]:
                                    'metadata'])
 
     @classmethod
-    def from_dict(cls: S1D, d: Dict[str, Any]) -> S1D:
+    def from_dict(cls: Type[T], d: Dict[str, Any]) -> T:
         """
         Convert a dictionary to a Spectrum1D object
 
@@ -403,8 +404,8 @@ def from_dict(cls: S1D, d: Dict[str, Any]) -> S1D:
                    metadata=d['metadata'])
 
     @classmethod
-    def from_castep_phonon_dos(cls: S1D, filename: str,
-                               element: Optional[str] = None) -> S1D:
+    def from_castep_phonon_dos(cls: Type[T], filename: str,
+                               element: Optional[str] = None) -> T:
         """
         Reads DOS from a CASTEP .phonon_dos file
 
@@ -423,11 +424,11 @@ def from_castep_phonon_dos(cls: S1D, filename: str,
         else:
             metadata['species'] = element
         metadata['label'] = element
-        return Spectrum1D(data['dos_bins']*ureg(data['dos_bins_unit']),
-                          data['dos'][element]*ureg(data['dos_unit']),
-                          metadata=metadata)
+        return cls(data['dos_bins']*ureg(data['dos_bins_unit']),
+                   data['dos'][element]*ureg(data['dos_unit']),
+                   metadata=metadata)
 
-    def broaden(self: S1D, x_width: Quantity, shape: str = 'gauss') -> S1D:
+    def broaden(self: T, x_width: Quantity, shape: str = 'gauss') -> T:
         """
         Broaden y_data and return a new broadened spectrum object
 
@@ -453,7 +454,7 @@ def broaden(self: S1D, x_width: Quantity, shape: str = 'gauss') -> S1D:
             [self.get_bin_centres().magnitude],
             [x_width.to(self.x_data_unit).magnitude],
             shape=shape)
-        return Spectrum1D(
+        return type(self)(
             np.copy(self.x_data.magnitude)*ureg(self.x_data_unit),
             y_broadened*ureg(self.y_data_unit),
             copy.copy((self.x_tick_labels)),
@@ -479,19 +480,23 @@ class Spectrum1DCollection(collections.abc.Sequence, Spectrum):
     y_data
         Shape (n_entries, n_x_data) float Quantity. The plot data in y,
         in rows corresponding to separate 1D spectra
-    x_tick_labels : Sequence[Tuple[int, str]] or None
-        Special tick labels e.g. for high-symmetry points. The int
-        refers to the index in x_data the label should be applied to
-    metadata : Dict[str, Union[int, str, LineData]] or None
-        Contains metadata about the spectra. Keys should be strings
-        and values should be strings or integers. There are some
-        functional keys:
+    x_tick_labels
+        Sequence[Tuple[int, str]] or None. Special tick labels e.g. for
+        high-symmetry points. The int refers to the index in x_data the
+        label should be applied to
+    metadata
+        Dict[str, Union[int, str, LineData]] or None. Contains metadata
+        about the spectra. Keys should be strings and values should be
+        strings or integers.
+        There are some functional keys:
           - 'line_data' : LineData
                           This is a Sequence[Dict[str, Union[int, str]],
                           it contains metadata for each spectrum in
                           the collection, and must be of length
                           n_entries
     """
+    T = TypeVar('T', bound='Spectrum1DCollection')
+
     def __init__(
             self, x_data: Quantity, y_data: Quantity,
             x_tick_labels: Optional[Sequence[Tuple[int, str]]] = None,
@@ -542,7 +547,7 @@ def __init__(
                     f'{len(metadata["line_data"])} entries')
         self.metadata = {} if metadata is None else metadata
 
-    def __add__(self, other: SC):
+    def __add__(self: T, other: T) -> T:
         """
         Appends the y_data of 2 Spectrum1DCollection objects,
         creating a single Spectrum1DCollection that contains
@@ -555,11 +560,11 @@ def __add__(self, other: SC):
         spectra are retained in the top level dictionary, any
         others are put in the individual 'line_data' entries
         """
-        return Spectrum1DCollection.from_spectra([*self, *other])
+        return type(self).from_spectra([*self, *other])
 
-    def _split_by_indices(self: SC,
+    def _split_by_indices(self,
                           indices: Union[Sequence[int], np.ndarray]
-                          ) -> List[SC]:
+                          ) -> List[T]:
         """Split data along x-axis at given indices"""
 
         ranges = self._ranges_from_indices(indices)
@@ -578,14 +583,14 @@ def __getitem__(self, item: int) -> Spectrum1D:
         ...
 
     @overload  # noqa: F811
-    def __getitem__(self, item: slice) -> SC:
+    def __getitem__(self, item: slice) -> T:
         ...
 
     @overload  # noqa: F811
-    def __getitem__(self, item: Union[Sequence[int], np.ndarray]) -> SC:
+    def __getitem__(self, item: Union[Sequence[int], np.ndarray]) -> T:
         ...
 
-    def __getitem__(self, item: Union[int, slice, Sequence[int]]):  # noqa: F811
+    def __getitem__(self, item: Union[int, slice, Sequence[int], np.ndarray]):  # noqa: F811
         new_metadata = copy.deepcopy(self.metadata)
         line_metadata = new_metadata.pop('line_data',
                                          [{} for _ in self._y_data])
@@ -616,7 +621,7 @@ def __getitem__(self, item: Union[int, slice, Sequence[int]]):  # noqa: F811
                           metadata=new_metadata)
 
     @classmethod
-    def from_spectra(cls: SC, spectra: Sequence[Spectrum1D]) -> SC:
+    def from_spectra(cls: Type[T], spectra: Sequence[Spectrum1D]) -> T:
         if len(spectra) < 1:
             raise IndexError("At least one spectrum is needed for collection")
 
@@ -648,7 +653,7 @@ def _type_check(spectrum):
 
     @staticmethod
     def _combine_metadata(all_metadata: Sequence[Dict[str, Union[int, str]]]
-                          ) -> Dict[str, Any]:
+                          ) -> Dict[str, Union[int, str, LineData]]:
         """
         From a sequence of metadata dictionaries, combines all common
         key/value pairs into the top level of a metadata dictionary,
@@ -729,7 +734,7 @@ def to_dict(self) -> Dict[str, Any]:
                                    'metadata'])
 
     @classmethod
-    def from_dict(cls: SC, d) -> SC:
+    def from_dict(cls: Type[T], d) -> T:
         """
         Convert a dictionary to a Spectrum1DCollection object
 
@@ -758,7 +763,7 @@ def from_dict(cls: SC, d) -> SC:
                    metadata=d['metadata'])
 
     @classmethod
-    def from_castep_phonon_dos(cls: SC, filename: str) -> SC:
+    def from_castep_phonon_dos(cls: Type[T], filename: str) -> T:
         """
         Reads total DOS and per-element PDOS from a CASTEP
         .phonon_dos file
@@ -783,7 +788,7 @@ def from_castep_phonon_dos(cls: SC, filename: str) -> SC:
             y_data*ureg(data['dos_unit']),
             metadata=metadata)
 
-    def broaden(self, x_width: Quantity, shape: str = 'gauss') -> SC:
+    def broaden(self, x_width: Quantity, shape: str = 'gauss') -> T:
         """
         Individually broaden each line in y_data, returning a new
         Spectrum1DCollection
@@ -817,7 +822,7 @@ def broaden(self, x_width: Quantity, shape: str = 'gauss') -> SC:
             copy.copy((self.x_tick_labels)),
             copy.deepcopy(self.metadata))
 
-    def group_by(self, *line_data_keys: str) -> SC:
+    def group_by(self, *line_data_keys: str) -> T:
         """
         Group and sum y_data for each spectrum according to the values
         mapped to the specified keys in metadata['line_data']
@@ -878,7 +883,7 @@ def sum(self) -> Spectrum1D:
                           metadata=metadata)
 
     def select(self, **select_key_values: Union[
-            str, int, Sequence[str], Sequence[int]]) -> SC:
+            str, int, Sequence[str], Sequence[int]]) -> T:
         """
         Select spectra by their keys and values in metadata['line_data']
 
@@ -941,20 +946,24 @@ class Spectrum2D(Spectrum):
         == (n_y_data + 1,))
     z_data
         Shape (n_x_data, n_y_data) float Quantity. The plot data in z
-    x_tick_labels : Sequence[Tuple[int, str]] or None
-        Special tick labels e.g. for high-symmetry points. The int
-        refers to the index in x_data the label should be applied to
-    metadata : Dict[str, Union[int, str]]
-        Contains metadata about the spectrum. Keys should be strings
-        and values should be strings or integers
+    x_tick_labels
+        Sequence[Tuple[int, str]] or None. Special tick labels e.g. for
+        high-symmetry points. The int refers to the index in x_data the
+        label should be applied to
+    metadata
+        Dict[str, Union[int, str]]. Contains metadata about the
+        spectrum. Keys should be strings and values should be strings
+        or integers
     """
+    T = TypeVar('T', bound='Spectrum2D')
+
     def __init__(self, x_data: Quantity, y_data: Quantity,
                  z_data: Quantity,
                  x_tick_labels: Optional[Sequence[Tuple[int, str]]] = None,
                  metadata: Optional[Dict[str, Union[int, str]]] = None
                  ) -> None:
         """
-        Attributes
+        Parameters
         ----------
         x_data
             Shape (n_x_data,) or (n_x_data + 1,) float Quantity. The
@@ -992,12 +1001,12 @@ def __init__(self, x_data: Quantity, y_data: Quantity,
         self.metadata = {} if metadata is None else metadata
 
     @property
-    def z_data(self):
+    def z_data(self) -> Quantity:
         return self._z_data*ureg(self._internal_z_data_unit).to(
             self.z_data_unit)
 
     @z_data.setter
-    def z_data(self, value):
+    def z_data(self, value: Quantity) -> None:
         self.z_data_unit = str(value.units)
         self._z_data = value.to(self._internal_z_data_unit).magnitude
 
@@ -1006,9 +1015,9 @@ def __setattr__(self, name: str, value: Any) -> None:
                                ['z_data_unit'])
         super(Spectrum2D, self).__setattr__(name, value)
 
-    def _split_by_indices(self: S2D,
+    def _split_by_indices(self,
                           indices: Union[Sequence[int], np.ndarray]
-                          ) -> List[S2D]:
+                          ) -> List[T]:
         """Split data along x-axis at given indices"""
         ranges = self._ranges_from_indices(indices)
         return [type(self)(self.x_data[x0:x1], self.y_data,
@@ -1018,9 +1027,9 @@ def _split_by_indices(self: S2D,
                            metadata=self.metadata)
                 for x0, x1 in ranges]
 
-    def broaden(self: S2D, x_width: Optional[Quantity] = None,
+    def broaden(self, x_width: Optional[Quantity] = None,
                 y_width: Optional[Quantity] = None, shape: str = 'gauss'
-                ) -> S2D:
+                ) -> T:
         """
         Broaden z_data and return a new broadened Spectrum2D object
 
@@ -1120,7 +1129,7 @@ def to_dict(self) -> Dict[str, Any]:
                                    'x_tick_labels', 'metadata'])
 
     @classmethod
-    def from_dict(cls: S2D, d: Dict[str, Any]) -> S2D:
+    def from_dict(cls: Type[T], d: Dict[str, Any]) -> T:
         """
         Convert a dictionary to a Spectrum2D object
 
diff --git a/euphonic/structure_factor.py b/euphonic/structure_factor.py
index 08230c280..64d244ebe 100644
--- a/euphonic/structure_factor.py
+++ b/euphonic/structure_factor.py
@@ -1,4 +1,4 @@
-from typing import Optional, TypeVar, Dict, Any
+from typing import Optional, TypeVar, Dict, Any, Type
 import warnings
 
 import numpy as np
@@ -9,9 +9,6 @@
                       Spectrum2D)
 
 
-T = TypeVar('T', bound='StructureFactor')
-
-
 class NoTemperatureError(Exception):
     pass
 
@@ -42,6 +39,8 @@ class StructureFactor(QpointFrequencies):
         population factor). None if no temperature-dependent effects
         have been applied
     """
+    T = TypeVar('T', bound='StructureFactor')
+
     def __init__(self, crystal: Crystal, qpts: np.ndarray,
                  frequencies: Quantity, structure_factors: Quantity,
                  weights: Optional[np.ndarray] = None,
@@ -91,17 +90,17 @@ def __init__(self, crystal: Crystal, qpts: np.ndarray,
             self.temperature_unit = str(ureg.K)
 
     @property
-    def structure_factors(self):
+    def structure_factors(self) -> Quantity:
         return self._structure_factors*ureg('bohr**2').to(
             self.structure_factors_unit)
 
     @structure_factors.setter
-    def structure_factors(self, value):
+    def structure_factors(self, value: Quantity) -> None:
         self.structure_factors_unit = str(value.units)
         self._structure_factors = value.to('bohr**2').magnitude
 
     @property
-    def temperature(self):
+    def temperature(self) -> Quantity:
         if self._temperature is not None:
             # See https://pint.readthedocs.io/en/latest/nonmult.html
             return Quantity(self._temperature,
@@ -110,11 +109,11 @@ def temperature(self):
             return None
 
     @temperature.setter
-    def temperature(self, value):
+    def temperature(self, value: Quantity) -> None:
         self.temperature_unit = str(value.units)
         self._temperature = value.to('K').magnitude
 
-    def __setattr__(self, name, value):
+    def __setattr__(self, name: str, value: Any) -> None:
         _check_unit_conversion(self, name, value,
                                ['frequencies_unit', 'structure_factors_unit',
                                 'temperature_unit'])
@@ -290,7 +289,8 @@ def _bose_corrected_structure_factor(self, e_bins: Quantity,
 
         return sqw_map*sf_conv/e_conv
 
-    def _bose_factor(self, temperature: Optional[Quantity] = None):
+    def _bose_factor(self, temperature: Optional[Quantity] = None
+                     ) -> np.ndarray:
         """
         Calculate the Bose factor for the frequencies stored in
         StructureFactor
@@ -358,7 +358,7 @@ def to_qpoint_frequencies(self) -> QpointFrequencies:
             self.crystal, self.qpts, self.frequencies, self.weights)
 
     @classmethod
-    def from_dict(cls: T, d: Dict[str, Any]) -> T:
+    def from_dict(cls: Type[T], d: Dict[str, Any]) -> T:
         """
         Convert a dictionary to a StructureFactor object
 
@@ -389,11 +389,11 @@ def from_dict(cls: T, d: Dict[str, Any]) -> T:
                    d['temperature'])
 
     @classmethod
-    def from_castep(cls):
+    def from_castep(cls: Type[T]) -> None:  # type: ignore
         ''
         raise AttributeError
 
     @classmethod
-    def from_phonopy(cls):
+    def from_phonopy(cls: Type[T]) -> None:  # type: ignore
         ''
         raise AttributeError
diff --git a/euphonic/util.py b/euphonic/util.py
index 8c9c53697..09cb37419 100644
--- a/euphonic/util.py
+++ b/euphonic/util.py
@@ -4,7 +4,7 @@
 import sys
 import warnings
 import os.path
-from typing import Dict, Sequence, Union, Tuple
+from typing import Dict, Sequence, Union, Tuple, Optional, List
 from collections import OrderedDict
 
 import numpy as np
@@ -17,11 +17,26 @@
 import euphonic.data
 
 
-def direction_changed(qpts, tolerance=5e-6):
+def direction_changed(qpts: np.ndarray, tolerance: float = 5e-6
+                      ) -> np.ndarray:
     """
-    Takes a N length list of q-points and returns an N - 2 length list of
-    booleans indicating whether the direction has changed between each pair
-    of q-points
+    Determines whether the q-direction has changed between each pair of
+    q-points
+
+    Parameters
+    ----------
+    qpts
+        Shape (n_qpts, 3) float ndarray. The q-points to use
+    tolerance
+        The minimum difference between the dot product and magnitude of
+        both vectors multiplied together for each pair of direction
+        vectors required for the direction to have changed
+
+    Returns
+    -------
+    has_direction_changed
+        Shape (n_qpts - 2,) bool ndarray. Indicates whether the
+        direction vector has changed between each pair of q-points
     """
 
     # Get vectors between q-points
@@ -38,42 +53,42 @@ def direction_changed(qpts, tolerance=5e-6):
     return direction_changed
 
 
-def is_gamma(qpt):
+def is_gamma(qpt: np.ndarray) -> Union[bool, np.ndarray]:
     """
     Determines whether the given point(s) are gamma points
 
     Parameters
     ----------
-    qpts: (3,) or (N, 3) float ndarray
-        The q-point or q-points
+    qpts
+        Shape (3,) or (N, 3) float ndarray. The q-point or q-points
 
     Returns
     -------
-    isgamma: bool or (N,) bool ndarray
-        Whether the input q-points(s) are gamma points. Returns a scalar
-        if only 1 q-point is provided
+    isgamma
+        bool or shape (N,) bool ndarray. Whether the input q-points(s)
+        are gamma points. Returns a scalar if only 1 q-point is
+        provided
     """
     tol = 1e-15
     isgamma = np.sum(np.absolute(qpt - np.rint(qpt)), axis=-1) < tol
     return isgamma
 
 
-def mp_grid(grid):
+def mp_grid(grid: Tuple[int, int, int]) -> np.ndarray:
     """
     Returns the q-points on a MxNxL Monkhorst-Pack grid specified by
     grid
 
     Parameters
     ----------
-    grid : (3,) int ndarray
-        Length 3 array specifying the number of points in each direction
+    grid
+        The number of points in each direction
 
     Returns
     -------
-    qgrid : (M*N*L, 3) float ndarray
-        Q-points on an MP grid
+    qgrid
+        Shape (M*N*L, 3) float ndarray. Q-points on an MP grid
     """
-
     # Monkhorst-Pack grid: ur = (2r-qr-1)/2qr where r=1,2..,qr
     qh = np.true_divide(
         2*(np.arange(grid[0]) + 1) - grid[0] - 1, 2*grid[0])
@@ -87,24 +102,27 @@ def mp_grid(grid):
     return np.column_stack((qh, qk, ql))
 
 
-def get_all_origins(max_xyz, min_xyz=[0, 0, 0], step=1):
+def get_all_origins(max_xyz: Tuple[int, int, int],
+                    min_xyz: Tuple[int, int, int] = (0, 0, 0),
+                    step: int = 1) -> np.ndarray:
     """
     Given the max/min number of cells in each direction, get a list of
     all possible cell origins
 
     Parameters
     ----------
-    max_xyz : (3,) int ndarray
+    max_xyz
         The number of cells to count to in each direction
-    min_xyz : (3,) int ndarray, optional
+    min_xyz
         The cell number to count from in each direction
-    step : integer, optional
+    step
         The step between cells
 
     Returns
     -------
-    origins : (prod(max_xyz - min_xyz)/step, 3) int ndarray
-        The cell origins
+    origins
+        Shape (prod(max_xyz - min_xyz)/step, 3) int ndarray. The cell
+        origins
     """
     diff = np.absolute(np.subtract(max_xyz, min_xyz))
     nx = np.repeat(range(min_xyz[0], max_xyz[0], step), diff[1]*diff[2])
@@ -115,23 +133,27 @@ def get_all_origins(max_xyz, min_xyz=[0, 0, 0], step=1):
     return np.column_stack((nx, ny, nz))
 
 
-def get_qpoint_labels(qpts, cell=None):
+def get_qpoint_labels(qpts: np.ndarray,
+                      cell: Optional[Tuple[List[List[float]],
+                                           List[List[float]],
+                                           List[int]]] = None
+                      ) -> List[Tuple[int, str]]:
     """
     Gets q-point labels (e.g. GAMMA, X, L) for the q-points at which the
     path through reciprocal space changes direction
 
     Parameters
     ----------
-    qpts : (n_qpts, 3) float ndarray
-        The q-points to get labels for
-    cell : (list, list, list), optional
+    qpts
+        Shape (n_qpts, 3) float ndarray. The q-points to get labels for
+    cell
         The cell structure as defined by spglib. Can be obtained by
         Crystal.to_spglib_cell. If not provided, the labels will be
         generic e.g. '1/3 1/2 0' rather than high-symmetry point labels
 
     Returns
     -------
-    x_tick_labels : list (int, string) tuples or None
+    x_tick_labels
         Tick labels and the q-point indices that they apply to
     """
     xlabels, qpts_with_labels = _recip_space_labels(qpts, cell=cell)
@@ -283,8 +305,8 @@ def _cell_vectors_to_volume(cell_vectors: np.ndarray) -> float:
 
 
 def _get_unique_elems_and_idx(
-        all_elems: Sequence[Tuple[Union[int, str]]]
-        ) -> 'OrderedDict[Union[str, int, tuple], np.ndarray]':
+        all_elems: Sequence[Tuple[Union[int, str], ...]]
+        ) -> 'OrderedDict[Tuple[Union[int, str], ...], np.ndarray]':
     """
     Returns an ordered dictionary mapping the unique sequences of
     elements to their indices
@@ -299,22 +321,26 @@ def _get_unique_elems_and_idx(
         ) for elem in unique_elems)
 
 
-def _calc_abscissa(reciprocal_cell, qpts):
+def _calc_abscissa(reciprocal_cell: Quantity, qpts: np.ndarray
+                       ) -> Quantity:
     """
     Calculates the distance between q-points (e.g. to use as a plot
     x-coordinate)
 
     Parameters
     ----------
-    reciprocal_cell : (3, 3) float Quantity
-        The reciprocal cell, can be calculated with
-        Crystal.reciprocal_cell
-    qpts : (n_qpts, 3) float ndarray
-        The q-points to get the distance between, in reciprocal lattice
-        units
-
-    abscissa : (n_qpts) float Quantity
-        The distance between q-points in 1/crystal.cell_vectors_unit
+    reciprocal_cell
+        Shape (3, 3) float Quantity. The reciprocal cell, can be
+        calculated with Crystal.reciprocal_cell
+    qpts
+        Shape (n_qpts, 3) float ndarray. The q-points to get the
+        distance between, in reciprocal lattice units
+
+    Returns
+    -------
+    abscissa
+        Shape (n_qpts) float Quantity. The distance between q-points
+        in 1/crystal.cell_vectors_unit
     """
     recip = reciprocal_cell.to('1/bohr').magnitude
     # Get distance between q-points in each dimension
@@ -353,28 +379,34 @@ def _calc_abscissa(reciprocal_cell, qpts):
     return abscissa*ureg('1/bohr').to(reciprocal_cell.units)
 
 
-def _recip_space_labels(qpts, cell=None):
+def _recip_space_labels(qpts: np.ndarray,
+                        cell: Optional[Tuple[List[List[float]],
+                                             List[List[float]],
+                                             List[int]]]
+                        ) -> Tuple[np.ndarray, np.ndarray]:
     """
     Gets q-points point labels (e.g. GAMMA, X, L) for the q-points at
     which the path through reciprocal space changes direction
 
     Parameters
     ----------
-    qpts : (n_qpts, 3) float ndarray
-        The q-points to get labels for
-    cell : (list, list, list), optional
+    qpts
+        Shape (n_qpts, 3) float ndarray. The q-points to get labels for
+    cell
         The cell structure as defined by spglib. Can be obtained by
         Crystal.to_spglib_cell. If not provided, the labels will be
         generic e.g. '1/3 1/2 0' rather than high-symmetry point labels
 
     Returns
     -------
-    labels : (n_qpts_direction_changed,) string ndarray
-        List of the labels for each q-point at which the path through
-        reciprocal space changes direction
-    qpts_with_labels : (n_qpts_direction_changed,) int ndarray
-        List of the indices of the q-points at which the path through
-        reciprocal space changes direction
+    labels
+        Shape (n_qpts_direction_changed,) string ndarray. The labels
+        for each q-point at which the path through reciprocal space
+        changes direction
+    qpts_with_labels
+        Shape (n_qpts_direction_changed,) int ndarray. The indices of
+        the q-points at which the path through reciprocal space changes
+        direction
     """
 
     # First and last q-points should always be labelled
@@ -405,7 +437,7 @@ def _recip_space_labels(qpts, cell=None):
     return labels, qpts_with_labels
 
 
-def _generic_qpt_labels():
+def _generic_qpt_labels() -> Dict[str, Tuple[float, float, float]]:
     """
     Returns a dictionary relating fractional q-point label strings to
     their coordinates e.g. '1/4 1/2 1/4' = [0.25, 0.5, 0.25]. Used for
@@ -419,12 +451,14 @@ def _generic_qpt_labels():
         for j, s2 in enumerate(label_strings):
             for k, s3 in enumerate(label_strings):
                 key = s1 + ' ' + s2 + ' ' + s3
-                value = [label_coords[i], label_coords[j], label_coords[k]]
+                value = (label_coords[i], label_coords[j], label_coords[k])
                 generic_labels[key] = value
     return generic_labels
 
 
-def _get_qpt_label(qpt, point_labels):
+def _get_qpt_label(qpt: np.ndarray,
+                   point_labels: Dict[str, Tuple[float, float, float]]
+                   ) -> str:
     """
     Gets a label for a particular q-point, based on the high symmetry
     points of a particular space group. Used for labelling the
@@ -432,9 +466,9 @@ def _get_qpt_label(qpt, point_labels):
 
     Parameters
     ----------
-    qpt : (3,) float ndarray
-        3 dimensional coordinates of a q-point
-    point_labels : dictionary
+    qpt
+        Shape (3,) float ndarray. Single q-point coordinates
+    point_labels
         A dictionary with N entries, relating high symmetry point labels
         (e.g. 'GAMMA', 'X'), to their 3-dimensional coordinates (e.g.
         [0.0, 0.0, 0.0]) where N = number of high symmetry points for a
@@ -442,7 +476,7 @@ def _get_qpt_label(qpt, point_labels):
 
     Returns
     -------
-    label : string
+    label
         The label for this q-point. If the q-point isn't a high symmetry
         point label is just an empty string
     """
@@ -476,7 +510,8 @@ def _get_qpt_label(qpt, point_labels):
     return label
 
 
-def _get_supercell_relative_idx(cell_origins, sc_matrix):
+def _get_supercell_relative_idx(cell_origins: np.ndarray,
+                                sc_matrix: np.ndarray) -> np.ndarray:
     """"
     For each cell_origins[i] -> cell_origins[j] vector in the supercell,
     gets the index n of the equivalent cell_origins[n] vector, where
@@ -485,16 +520,18 @@ def _get_supercell_relative_idx(cell_origins, sc_matrix):
 
     Parameters
     ----------
-    cell_origins : (n_cells, 3) int ndarray
-        The vector to the origin of each cell in the supercell, in unit
-        cell fractional coordinates
-    sc_matrix : (3, 3) int ndarray
-        The matrix for converting from the unit cell to the supercell
+    cell_origins
+        Shape (n_cells, 3) int ndarray. The vector to the origin of
+        each cell in the supercell, in unit cell fractional coordinates
+    sc_matrix
+        Shape (3, 3) int ndarray. The matrix for converting from the
+        unit cell to the supercell
 
     Returns
     -------
-    supercell_relative_idx : (n_cells, n_cells) in ndarray
-        The index n of the equivalent vector in cell_origins
+    supercell_relative_idx
+        Shape (n_cells, n_cells) int ndarray. The index n of the
+        equivalent vector in cell_origins
     """
     n_cells = len(cell_origins)
     ax = np.newaxis
diff --git a/euphonic/validate.py b/euphonic/validate.py
index 64699fba5..75c123870 100644
--- a/euphonic/validate.py
+++ b/euphonic/validate.py
@@ -1,4 +1,4 @@
-from typing import List, Any
+from typing import List, Any, Union, Type, Tuple, Optional
 
 import numpy as np
 from pint import DimensionalityError
@@ -6,28 +6,32 @@
 from euphonic import ureg
 
 
-def _check_constructor_inputs(objs, types, shapes, names):
+def _check_constructor_inputs(
+        objs: List[object],
+        types: List[Union[Type, List[Type]]],
+        shapes: List[Union[Tuple[int, ...], List[Tuple[int, ...]]]],
+        names: List[str]) -> None:
     """
     Make sure all the inputs are all the expected type, and if they are
     an array, the correct shape
 
     Parameters
     ----------
-    objs : list of objects
+    objs
         The objects to check
-    types : list of types or lists of types
+    types
         The expected class of each input. If multiple types are
         accepted, the expected class can be a list of types. e.g.
         types=[[list, np.ndarray], int]
-    shapes : list of tuples
+    shapes
         The expected shape of each object (if the object has a shape
         attribute). If the shape of some dimensions don't matter,
         provide -1 for those dimensions, or if none of the dimensions
         matter, provide an empty tuple (). If multiple shapes are
         accepted, the expected shapes can be a list of tuples. e.g.
-        shapes=[[(n, 3), (n + 1, 3)], 3]
-    names : list of strings
-        The name of each array
+        shapes=[[(n, 3), (n + 1, 3)], (3,)]
+    names
+        The name of each input variable
 
     Raises
     ------
@@ -85,18 +89,20 @@ def _check_unit_conversion(obj: object, attr_name: str, attr_value: Any,
                     f'unit for "{attr_name}"'))
 
 
-def _replace_dim(expected_shape, obj_shape):
+def _replace_dim(expected_shape: Tuple[int, ...],
+                 obj_shape: Tuple[int, ...]) -> Tuple[int, ...]:
     # Allow -1 dimensions to be any size
     idx = np.where(np.array(expected_shape) == -1)[0]
     if len(idx) == 0 or len(expected_shape) != len(obj_shape):
         return expected_shape
     else:
-        expected_shape = np.array(expected_shape)
-        expected_shape[idx] = np.array(obj_shape)[idx]
-        return tuple(expected_shape)
+        expected_shape_replaced = np.array(expected_shape)
+        expected_shape_replaced[idx] = np.array(obj_shape)[idx]
+        return tuple(expected_shape_replaced)
 
 
-def _ensure_contiguous_attrs(obj, required_attrs, opt_attrs=[]):
+def _ensure_contiguous_attrs(obj: object, required_attrs: List[str],
+                             opt_attrs: List[str] = []) -> None:
     """
     Make sure all listed attributes of obj are C Contiguous and of the
     correct type (int32, float64, complex128). This should only be used
@@ -105,12 +111,12 @@ def _ensure_contiguous_attrs(obj, required_attrs, opt_attrs=[]):
 
     Parameters
     ----------
-    obj : Object
+    obj
         The object that will have it's attributes checked
-    required_attrs : list of strings
+    required_attrs
         The attributes of obj to be checked. They should all be Numpy
         arrays
-    opt_attrs : list of strings, optional
+    opt_attrs
         The attributes of obj to be checked, but if they don't exist
         will not throw an error. e.g. Depending on the material
         ForceConstants objects may or may not have 'born' defined
@@ -129,7 +135,7 @@ def _ensure_contiguous_attrs(obj, required_attrs, opt_attrs=[]):
             pass
 
 
-def _ensure_contiguous_args(*args):
+def _ensure_contiguous_args(*args: np.ndarray) -> Tuple[np.ndarray, ...]:
     """
     Make sure all arguments are C Contiguous and of the correct type
     (int32, float64, complex128). This should only be used internally,
@@ -138,33 +144,30 @@ def _ensure_contiguous_args(*args):
 
     Parameters
     ----------
-    *args : any number of ndarrays
+    *args
         The Numpy arrays to be checked
 
     Returns
     -------
-    args_contiguous : the same number of ndarrays as args
+    contiguous_args
         The same as the provided args, but all contiguous.
     """
-    args = list(args)
-    for i in range(len(args)):
-        args[i] = args[i].astype(_get_dtype(args[i]), order='C', copy=False)
-
-    return args
+    return tuple(arg.astype(_get_dtype(arg), order='C', copy=False)
+                 for arg in args)
 
 
-def _get_dtype(arr):
+def _get_dtype(arr: np.ndarray) -> Optional[Type]:
     """
     Get the Numpy dtype that should be used for the input array
 
     Parameters
     ----------
-    arr : ndarray
+    arr
         The Numpy array to get the type of
 
     Returns
     -------
-    dtype : Numpy dtype
+    dtype
         The type the array should be
     """
     if np.issubdtype(arr.dtype, np.integer):


Test patch:
diff --git a/tests_and_analysis/test/data/phonopy_files/NaCl/BORN_nacl_nofactor b/tests_and_analysis/test/data/phonopy_files/NaCl/BORN_nacl_nofactor
new file mode 100644
index 000000000..a5062e44f
--- /dev/null
+++ b/tests_and_analysis/test/data/phonopy_files/NaCl/BORN_nacl_nofactor
@@ -0,0 +1,10 @@
+# Testing behaviour with comments and no factor
+2.43533967 0 0 0 2.43533967 0 0 0 2.43533967
+1.086875 0 0 0 1.086875 0 0 0 1.086875
+1.086875 0 0 0 1.086875 0 0 0 1.086875
+1.086875 0 0 0 1.086875 0 0 0 1.086875
+1.086875 0 0 0 1.086875 0 0 0 1.086875
+-1.086875 0 0 0 -1.086875 0 0 0 -1.086875
+-1.086875 0 0 0 -1.086875 0 0 0 -1.086875
+-1.086875 0 0 0 -1.086875 0 0 0 -1.086875
+-1.086875 0 0 0 -1.086875 0 0 0 -1.086875
diff --git a/tests_and_analysis/test/euphonic_test/test_force_constants.py b/tests_and_analysis/test/euphonic_test/test_force_constants.py
index 3fe9c9794..ae5717268 100644
--- a/tests_and_analysis/test/euphonic_test/test_force_constants.py
+++ b/tests_and_analysis/test/euphonic_test/test_force_constants.py
@@ -224,6 +224,9 @@ def create_from_dict(self, request):
         ('NaCl', {'summary_name': 'phonopy_nofc_noborn.yaml',
                   'fc_name': 'FORCE_CONSTANTS_nacl',
                   'born_name': 'BORN_nacl'}),
+        ('NaCl', {'summary_name': 'phonopy_nofc_noborn.yaml',
+                  'fc_name': 'FORCE_CONSTANTS_nacl',
+                  'born_name': 'BORN_nacl_nofactor'}),
         # Explicitly test the default behaviour (if fcs aren't found
         # in phonopy.yaml they should be read from FORCE_CONSTANTS. BORN
         # is not read by default as it is not required, so it is difficult
diff --git a/tests_and_analysis/test/euphonic_test/test_force_constants_calculate_qpoint_phonon_modes.py b/tests_and_analysis/test/euphonic_test/test_force_constants_calculate_qpoint_phonon_modes.py
index 9c8c16f2a..998248fe4 100644
--- a/tests_and_analysis/test/euphonic_test/test_force_constants_calculate_qpoint_phonon_modes.py
+++ b/tests_and_analysis/test/euphonic_test/test_force_constants_calculate_qpoint_phonon_modes.py
@@ -265,9 +265,9 @@ def test_calc_qpt_ph_modes_with_large_complex_mode_gradients_raises_warning(
         # Mock force constants to return nonsense gradients with high
         # imaginary terms
         class MockFC(ForceConstants):
-            def _calculate_phonons_at_q(self, q, args):
+            def _calculate_phonons_at_q(self, *args):
                 freqs, evecs, grads = ForceConstants._calculate_phonons_at_q(
-                    self, q, args)
+                    self, *args)
                 mocked_grads = np.ones(grads.shape) + np.ones(grads.shape)*1j
                 return freqs, evecs, mocked_grads
         mocker.patch('euphonic.ForceConstants', MockFC)