diff --git a/doc/source/conf.py b/doc/source/conf.py
index 17ed96341288345a545cbc235a3ad625b04473d0..d28e97fce63ee03eee1686228c3677ef40272b05 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -16,6 +16,8 @@
 import sys
 import os
 
+import sphinx_rtd_theme
+
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
@@ -113,7 +115,7 @@ todo_include_todos = False
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
-html_theme = 'alabaster'
+html_theme = 'sphinx_rtd_theme'
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
@@ -121,7 +123,7 @@ html_theme = 'alabaster'
 #html_theme_options = {}
 
 # Add any paths that contain custom themes here, relative to this directory.
-#html_theme_path = []
+html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 
 # The name for this set of Sphinx documents.  If None, it defaults to
 # "<project> v<release> documentation".
diff --git a/pycgtool.py b/pycgtool.py
index 52cfcdc1b479bd3a1b919d0ae71509fb5b660f9f..d2676b00d4ddc3f4d681b904555ec3f6bc38752c 100755
--- a/pycgtool.py
+++ b/pycgtool.py
@@ -9,6 +9,13 @@ from pycgtool.forcefield import ForceField
 
 
 def main(args):
+    """
+    Main function of the program PyCGTOOL.
+
+    Performs the complete AA->CG mapping and outputs a GROMACS forcefield directory.
+
+    :param args: Arguments from argparse
+    """
     frame = Frame(gro=args.gro, xtc=args.xtc)
 
     if args.bnd is not None:
diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py
index 4a3461581f7ac26bfdeaf86b8ef3eedff80f6298..d45912b714cb0094e0655954df1376e18834ece3 100644
--- a/pycgtool/bondset.py
+++ b/pycgtool/bondset.py
@@ -1,3 +1,9 @@
+"""
+Module containing classes to calculate bonded properties from a Frame.
+
+BondSet contains a dictionary of lists of Bonds.  Each list corresponds to a single molecule.
+"""
+
 import numpy as np
 
 from .util import stat_moments, sliding
@@ -7,9 +13,21 @@ np.seterr(all="raise")
 
 
 class Bond:
+    """
+    Class holding the properties of a single bonded term.
+
+    Bond lengths, angles and dihedrals are all equivalent, distinguised by the number of atoms present.
+    """
     __slots__ = ["atoms", "atom_numbers", "values", "eqm", "fconst"]
 
     def __init__(self, atoms=None, atom_numbers=None):
+        """
+        Create a single bond definition.
+
+        :param atoms: List of atom names defining the bond
+        :param atom_numbers: List of atom numbers defining the bond
+        :return: Instance of Bond
+        """
         self.atoms = atoms
         self.atom_numbers = atom_numbers
         self.values = []
@@ -18,6 +36,11 @@ class Bond:
         return len(self.atoms)
 
     def boltzmann_invert(self, temp=310):
+        """
+        Perform Boltzmann Inversion using measured values of this bond to calculate equilibrium value and force constant.
+
+        :param temp: Temperature at which the simulation was performed
+        """
         mean, var = stat_moments(self.values)
 
         rt = 8.314 * temp / 1000.
@@ -42,6 +65,14 @@ class Bond:
 
 
 def angle(a, b, c=None):
+    """
+    Calculate the signed angle between two/three vectors.
+
+    :param a: First vector
+    :param b: Second vector
+    :param c: Optional third orientation vector
+    :return: Signed angle in radians
+    """
     if c is None:
         c = np.cross(a, b)
     dot = np.dot(a, b)
@@ -50,7 +81,18 @@ def angle(a, b, c=None):
 
 
 class BondSet:
+    """
+    Class used to perform bond measurements in a Frame.
+
+    BondSet contains a dictionary of lists of Bonds.  Each list corresponds to a single molecule.
+    """
     def __init__(self, filename=None):
+        """
+        Read in bonds from a file.
+
+        :param filename: File to read
+        :return: Instance of BondSet
+        """
         self._molecules = {}
 
         if filename is not None:
@@ -62,6 +104,13 @@ class BondSet:
                         molbnds.append(Bond(atoms=atomlist))
 
     def _populate_atom_numbers(self, mapping):
+        """
+        Add atom numbers to all bonds.
+
+        Uses previously defined atom names.
+
+        :param mapping: AA->CG mapping used to collect atom/bead numbers
+        """
         for mol in self._molecules:
             molmap = mapping[mol]
             index = [bead.name for bead in molmap]
@@ -69,6 +118,12 @@ class BondSet:
                 bond.atom_numbers = [index.index(atom.lstrip("+-")) for atom in bond.atoms]
 
     def write_itp(self, filename, mapping):
+        """
+        Output a GROMACS .itp file containing atoms/beads and bonded terms.
+
+        :param filename: Name of output file
+        :param mapping: AA->CG Mapping from which to collect bead properties
+        """
         self._populate_atom_numbers(mapping)
 
         with open(filename, "w") as itp:
@@ -120,6 +175,11 @@ class BondSet:
                     ), file=itp)
 
     def apply(self, frame):
+        """
+        Calculate bond lengths/angles for a given Frame and store into Bonds.
+
+        :param frame: Frame from which to calculate values
+        """
         def calc_length(atoms):
             vec = atoms[1].coords - atoms[0].coords
             return np.sqrt(np.sum(vec*vec))
@@ -162,6 +222,9 @@ class BondSet:
                     pass
 
     def boltzmann_invert(self):
+        """
+        Perform Boltzmann Inversion of all bonds to calculate equilibrium value and force constant.
+        """
         for mol in self._molecules:
             for bond in self._molecules[mol]:
                 bond.boltzmann_invert()
diff --git a/pycgtool/forcefield.py b/pycgtool/forcefield.py
index 6e0f39988e352823bda4ab6b25b13a47e0650754..d552cf9053b6c1acfc6617b1a5f443e0f1a2d4f3 100644
--- a/pycgtool/forcefield.py
+++ b/pycgtool/forcefield.py
@@ -1,3 +1,7 @@
+"""
+This module contains a single class ForceField used to output a GROMACS .ff forcefield.
+"""
+
 import os
 
 # Python 3.2 doesn't have FileExistsError
@@ -11,7 +15,15 @@ except FileExistsError:
 
 
 class ForceField:
+    """
+    Class used to output a GROMACS .ff forcefield
+    """
     def __init__(self, name):
+        """
+        Open a forcefield directory at name.  If it does not exist it is created.
+
+        :param name: Directory name to open/create
+        """
         try:
             os.makedirs(name)
         except FileExistsError as e:
@@ -33,6 +45,13 @@ class ForceField:
             print("PyCGTOOL produced force field", file=doc)
 
     def write_atp(self, mapping):
+        """
+        Write a GROMACS forcefield .atp file.
+
+        This file lists the atomtypes used in the forcefield.
+
+        :param mapping: AA->CG mapping from which to collect atomtypes
+        """
         with open(os.path.join(self.dirname, "atomtypes.atp"), "w") as atp:
             types = set()
             for mol in mapping:
@@ -42,6 +61,15 @@ class ForceField:
                         types.add(bead.type)
 
     def write_rtp(self, name, mapping, bonds):
+        """
+        Write a GROMACS .rtp file.
+
+        This file defines the residues present in the forcefield and allows pdb2gmx to be used.
+
+        :param name: Name of the .rtp file to create
+        :param mapping: AA->CG mapping from which to collect molecules
+        :param bonds: BondSet from which to collect bonds
+        """
         # TODO print everything to file
         with open(os.path.join(self.dirname, name), "w") as rtp:
             print("[ bondedtypes ]", file=rtp)
@@ -82,4 +110,3 @@ class ForceField:
                         bond.atoms[0], bond.atoms[1], bond.atoms[2], bond.atoms[3],
                         bond.eqm, bond.fconst, 1
                     ), file=rtp)
-
diff --git a/pycgtool/frame.py b/pycgtool/frame.py
index 8efba67adf5ec97d9429b04ffb40c55a377f3d27..364e015a36fac5b2e3a366d4513fc192e66a89e4 100644
--- a/pycgtool/frame.py
+++ b/pycgtool/frame.py
@@ -1,3 +1,10 @@
+"""
+This module contains classes for storing atomic data.
+
+The Frame class may contain multiple Residues which may each contain multiple Atoms.
+Both Frame and Residue are iterable. Residue is indexable with either atom numbers or names.
+"""
+
 import numpy as np
 
 from simpletraj import trajectory
diff --git a/pycgtool/mapping.py b/pycgtool/mapping.py
index 628054f1e753fce763e83a0f5dd3ac4b2355522a..9a9532fbb170bef558c60b88674aa414f183e5bb 100644
--- a/pycgtool/mapping.py
+++ b/pycgtool/mapping.py
@@ -1,3 +1,10 @@
+"""
+This module contains classes required to perform an atomistic to coarse-grain mapping.
+
+The Mapping class contains a dictionary of lists of BeadMaps.
+Each list corresponds to a single molecule.
+"""
+
 import numpy as np
 
 from .frame import Atom, Residue, Frame
@@ -7,9 +14,22 @@ np.seterr(all="raise")
 
 
 class BeadMap:
+    """
+    POD class holding values relating to the AA->CG transformation for a single bead.
+    """
     __slots__ = ["name", "typ", "type", "atoms", "charge", "mass"]
 
     def __init__(self, name=None, type=None, atoms=None, charge=0, mass=0):
+        """
+        Create a single bead mapping
+
+        :param name: The name of the bead
+        :param type: The bead type
+        :param atoms: The atom names from which the bead is made up
+        :param charge: The net charge on the bead
+        :param mass: The total bead mass
+        :return: Instance of BeadMap
+        """
         self.name = name
         self.type = type
         self.atoms = atoms
@@ -17,6 +37,11 @@ class BeadMap:
         self.mass = mass
 
     def __iter__(self):
+        """
+        Iterate through the atom names from which the bead is made up.
+
+        :return: Iterator over atoms
+        """
         return iter(self.atoms)
 
     def __len__(self):
@@ -24,11 +49,25 @@ class BeadMap:
 
 
 class EmptyBeadError(Exception):
+    """
+    Exception used to indicate that none of the required atoms are present.
+    """
     pass
 
 
 class Mapping:
+    """
+    Class used to perform the AA->CG mapping.
+
+    Contains a dictionary of lists of BeadMaps.  Each list corresponds to a single molecule.
+    """
     def __init__(self, filename=None):
+        """
+        Read in the AA->CG mapping from a file.
+
+        :param filename: File from which to read mapping
+        :return: Instance of Mapping
+        """
         self._mappings = {}
 
         if filename is not None:
@@ -52,7 +91,14 @@ class Mapping:
     def __iter__(self):
         return iter(self._mappings)
 
-    def apply(self, frame, exclude=set()):
+    def apply(self, frame, exclude=None):
+        """
+        Apply the AA->CG mapping to an atomistic Frame.
+
+        :param frame: Frame to which mapping will be applied
+        :param exclude: Set of molecule names to exclude from mapping - e.g. solvent
+        :return: A new Frame instance containing the CG frame
+        """
         cgframe = Frame()
         cgframe.name = frame.name
         cgframe.box = frame.box
@@ -60,7 +106,7 @@ class Mapping:
         cgframe.residues = []
 
         for aares in frame:
-            if aares.name in exclude:
+            if exclude is not None and aares.name in exclude:
                 continue
             try:
                 molmap = self._mappings[aares.name]
diff --git a/pycgtool/parsers/cfg.py b/pycgtool/parsers/cfg.py
index d287ae6ba5e5474d97da5fd794b55adb63b4f240..32fa4f31bda67f6b92a7ffd799398a32ebd6710d 100644
--- a/pycgtool/parsers/cfg.py
+++ b/pycgtool/parsers/cfg.py
@@ -1,10 +1,23 @@
+"""
+Module containing classes used to parse custom CFG file format.
+
+Format is based upon GROMACS .itp files but does not support nesting of sections.
+"""
 import os
 
 
 class Section:
+    """
+    Class representing a single section of the config file.
+    """
     __slots__ = ["name", "_lines"]
 
     def __init__(self, name=None):
+        """
+        Create a section and storage for the lines it contains.
+
+        :param name: Name of section
+        """
         self.name = name
         self._lines = []
 
@@ -19,14 +32,28 @@ class Section:
 
 
 class DuplicateSectionError(Exception):
+    """
+    Exception used to indicate that a section has appeared twice in a file.
+    """
     def __repr__(self):
         return "Section {0} appears twice in file {1}.".format(*self.args)
 
 
 class CFG:
+    """
+    Class representing a CFG file.
+
+    Contains a dictionary of Sections.
+    """
     __slots__ = ["filename", "_sections", "_section_names", "_iter_section"]
 
     def __init__(self, filename=None):
+        """
+        Parse a config file and extract Sections.
+
+        :param filename: Name of file to read
+        :return: Instance of CFG
+        """
         self.filename = filename
         self._sections = {}
         self._section_names = []
@@ -76,10 +103,6 @@ class CFG:
         sec = self._section_names[self._iter_section]
         return self._sections[sec]
 
-    # For Python 2
-    def next(self):
-        return self.__next__()
-
     def __contains__(self, item):
         return item in self._section_names
 
diff --git a/pycgtool/util.py b/pycgtool/util.py
index 0ac391bdbf167554b8c5b645fdc475797e45d433..a1733d8775fb21ed3313ee5f7f6fa3fc51285366 100644
--- a/pycgtool/util.py
+++ b/pycgtool/util.py
@@ -1,3 +1,7 @@
+"""
+This module contains some general purpose utility functions used in PyCGTOOL.
+"""
+
 import numpy as np
 np.seterr(all="raise")
 
@@ -11,7 +15,7 @@ def stat_moments(vals, ignore_nan=True):
     :return: Numpy array of moments - population mean and variance
     """
     if ignore_nan:
-        vals_tmp = [val for val in vals if not np.isnan(val)]
+        vals_tmp = [val for val in vals if np.isfinite(val)]
     else:
         vals_tmp = vals
 
diff --git a/test/test_forcefield.py b/test/test_forcefield.py
index 8c847eb7609b8c2638cd5877b2953e6dcae75c7b..306ff610cce6daf86ba340b576f89883ec96f6fd 100644
--- a/test/test_forcefield.py
+++ b/test/test_forcefield.py
@@ -6,7 +6,7 @@ from pycgtool.forcefield import ForceField
 
 class ForceFieldTest(unittest.TestCase):
     def test_create(self):
-        name = "test.ff"
+        name = "fftest.ff"
         ff = ForceField(name)
         self.assertTrue(os.path.exists(name))
         self.assertTrue(os.path.isdir(name))