School of Mathematics


Software created by members of the Applied and Computational Mathematics group

Ben Goddard


2DChebClass is a Matlab implementation of a novel, efficient pseudo-spectral collocation scheme for the solution of stiff, non-local, non-linear PDEs in one and two spatial dimensions. In particular, we compute the non-local terms in real space with the help of a specialised Gauss quadrature. Due to the exponential accuracy of the quadrature and a convenient choice of collocation points near interfaces, we can use grids with a significantly lower number of nodes than most other reported methods. We have demonstrated the capabilities of our numerical methodology by applying it to (Dynamical) Density Functional Theory problems, studying equilibrium and dynamic two-dimensional test cases with single-and multispecies hard-sphere and hard-disc particles modelled with fundamental measure theory, with and without van der Waals attractive forces, in bounded and unbounded physical domains. We have shown that our results satisfy statistical mechanical sum rules. See the associated paper for details.

2DChebClass was jointly developed with Andreas Nold (Max Plank Institute for Brain Research) and Serafim Kalliadasis (Imperial College London, Chem. Eng.). It is freely available through Edinburgh DataShare and GitHub.

Superadiabatic quantum molecular dynamics

Matlab code for computation of transitions through avoided crossings using superadiabatic representations. See e.g. 'Wave packet dynamics in the optimal superadiabatic approximation' V. Betz, B. D. Goddard, and U. Manthe, The Journal of Chemical Physics 144, 224109 (2016). Open access pdf.

In the above paper, we explain the concept of superadiabatic representations and show how in the context of electronically non-adiabatic transitions they lead to an explicit formula that can be used to predict transitions at avoided crossings. Based on this formula, we present a simple method for computing wave packet dynamics across avoided crossings. Only knowledge of the adiabatic potential energy surfaces near the avoided crossing is required for the computation. In particular, this means that no diabatization procedure is necessary, the adiabatic electronic energies can be computed on the fly, and they only need to be computed to higher accuracy when an avoided crossing is detected. We test the quality of our method on the paradigmatic example of photo-dissociation of NaI, finding very good agreement with results of exact wave packet calculations.  Both the code and data for NaI are available.

Ben Leimkuhler

Extasy Project

Computational simulations are a powerful approach to addressing one of the "Grand Challenges in the Chemical Sciences", namely the ability to understand and control the properties of complex macromolecular systems. In turn, one of the great challenges for computational simulations is to sample the conformational space of complex macromolecular systems sufficiently to enable reliable predictions to be made. The ExTASY project involves the development of a set of software tools that, by coupling large ensemble Molecular Dynamics calculations with novel analysis methods and efficient integrators, can provide orders-of-magnitude improvements in sampling in comparison to standard approaches. ExTASY tools are designed to be flexible (able to be applied to many different types of simulation problem) and extensible (able to be integrated with many other simulation software tools)

MIST (Molecular Integration Simulation Toolkit) is a C++ library which provides two things - a simple API for integrator developers which hides the complexities of real MD codes, enabling rapid development of new integration algorithms at a higher level of abstraction; and a plug-in interface, enabling the new integrators to be deployed directly in several MD codes, taking advantage of their highly efficient force evaluation code, file format support, and supporting tool chains.


Source files, accompanying "Efficient molecular dynamics using geodesic integration and solvent-solute splitting" (Proceedings of the Royal Society, 2016) are available

MD.M (Molecular Dynamics in Matlab)

MD.M (pronounced em dee dot em) is an extremely easy to use MATLAB toolkit for experimenting with molecular dynamics integrators and simple molecular models.

Although MATLAB is an interpreter-based programming system, recent improvements in the software mean that M-files run reasonably fast and it is possible to do some elementary studies involving a handful of atoms (and certainly it is enough for understanding the main differences among timestepping algorithms). For larger systems (with, say, more than a few hundred atoms, or if long simulation times are needed) it is recommended to either take advantage of the CMEX interface (see submenu at left) or to use an advanced molecular modelling software package such as NAMD/VMD, Gromacs, Amber, LAMMPS or similar which will offer much greater computational efficiency as well as having a force field interface, a more varied selection of algorithms, visualization tools, etc. MD.M has been specifically formulated to accompany the Molecular Dynamics book of B. Leimkuhler and C. Matthews. 

James Maddison


Dr. Maddison was the PI of the ARCHER embedded CSE eCSE03-8, "Parallel supermeshing for multimesh modelling", for which the libsupermesh parallel supermeshing library was developed. Details regarding the library can be found in the manual and in the report.

The latest version can be found in the libsupermesh Bitbucket repository.

timestepping library

Dr. Maddison is the developer of the "timestepping" extension to the FEniCS system, which is included with the dolfin-adjoint library. This is used to optimise time dependent finite element problems, with automated adjoining and similar optimisation of the time-dependent adjoint.

The principles used by the library are described in Maddison and Farrell, CMAME 276, pp. 95‒121, 2014. Further details can be found in the manual.

The latest version, compatible with FEniCS 1.5 and 1.6, can be found in the dolfin-adjoint Bitbucket repository.  

John Pearson


Code for computing the confluent and Gauss hypergeometric functions numerically can be found on Edinburgh DataShare. This accompanies the paper "Numerical methods for the computation of the confluent and Gauss hypergeometric functions" (Numerical Algorithms, 2017).


Code for implementing a deferred correction scheme for time-dependent PDE-constrained optimization problems can be found here. This accompanies the paper "A rational deferred correction approach to parabolic optimal control problems" (IMA Journal of Numerical Analysis, 2017).


Code applying an iterative solver for the Ohta-Kawasaki equation can be found on Bitbucket. It relies on FEniCS (version 1.6) and PETSc/petsc4py (version 3.6). This accompanies the paper "A preconditioner for the Ohta-Kawasaki equation" (SIAM Journal on Matrix Analysis and Applications, 2017).


Code applying refined saddle-point preconditioners for the Stokes equations can be found here. This accompanies the paper "Refined saddle-point preconditioners for discretized Stokes problems" (Numerische Mathematik, 2017).


Code implementing tensor train solvers for optimization problems with fractional differential equation constraints can be found here. This accompanies the paper "Fast tensor product solvers for optimization problems with fractional differential equations as constraints" (Applied Mathematics and Computation, 2016).


Dr Pearson wrote the codes for solving PDE-constrained optimization problems within the IFISS software system, which can be found here.