Comparison between Density Functional Theory and Density Functional Tight Binding approaches for finding the muon stopping site in organic molecular crystals

TL;DR

This study compares DFT and DFTB+ methods for locating muon stopping sites in organic crystals, demonstrating DFTB+'s efficiency and proposing a Python toolkit for standardized site-finding in muon spectroscopy.

Contribution

It introduces a comparative analysis of DFT and DFTB+ for muon site prediction and presents a Python software suite for routine and standardized calculations.

Findings

DFTB+ offers a faster alternative to DFT for muon site prediction.

Both methods show comparable accuracy with certain limitations.

The Python toolkit facilitates standardized and accessible computations.

Abstract

Finding the possible stopping sites for muons inside a crystalline sample is a key problem of muon spectroscopy. In a previous work, we suggested a computational approach to this problem, using Density Functional Theory software in combination with a random structure searching approach using a Poisson sphere distribution. In this work we test this methodology further by applying it to three organic molecular crystals model systems: durene, bithiophene, and tetracyanoquinodimethane (TCNQ). Using the same sets of random structures we compare the performance of Density Functional Theory software CASTEP and the much faster lower level approximation of Density Functional Tight Binding provided by DFTB+, combined with the use of the 3ob-3-1 parameter set. We show the benefits and limitations of such an approach and we propose the use of DFTB+ as a viable alternative to more cumbersome simulations for routine site-finding in organic materials. Finally, we introduce the Muon Spectroscopy Computational Project software suite, a library of Python tools meant to make these methods standardized and easy to use.

Figures

Tables

Table 1. Table 1: Structural characteristics of bithiophene, durene and TCNQ.

Compound	Formula	Mol. per cell	Space Group	Lattice Parameters [$\AA$]	Angles [Degrees]
bithiophene	$C_8{H}_6{S}_2$	2	P 21/c	(7.734, 5.729, 8.933)	(90.00, 106.72, 90.00)
durene	$C_{10}{H}_{14}$	2	P 21/a	(11.590, 5.740, 7.040)	(90.00, 112.80, 90.00)
TCNQ	$C_{12}{H}_4{N}_4$	4	C 2/c	(8.896, 6.913, 16.439)	(90.00, 98.29, 90.00)

Table 2. Table 2: Van der Waals radii for elements appearing in the tested molecules.

Element	H	C	N	S
$r_{vdw}$ (Å)	1.2	1.95	1.55	2

Table 3. Table 3: Clusters obtained from CASTEP and DFTB + + calculations on bithiophene, TCNQ and durene muoniated structures, split by molecular site. Each cluster represents a crystallographically inequivalent version of the same molecular site.

CASTEP
System	Cluster Type	N∘ of Clusters	Structures per Cluster
	S1	3	82,35,8
	C2	2	14,19
Bithiophene	C3	3	34,68,1
	C4	3	3,17,16
	C5	2	49,35
	float	3	9,1,1
	N1	4	57,60,49,50
	C2	4	11,10,19,23
TCNQ	C3	2	31,26
	C4	2	14,10
	C5	4	30,29,32,36
	float	2	32,1
	C1	3	2,1,1
	C2	4	24,27,19,16
Durene	C3	2	80,124
	float	4	3,83,65,1
DFTB$+$
System	Cluster Type	N∘ of Clusters	Structures per Cluster
	S1	4	31,32,1,1
	C2	3	13,17,1
	C3	3	50,69,1
	C4	3	2,11,15
Bithiophene	C5	2	43,70
	S(1,7)	1	1
	S(5,6)	2	4,1
	S(1,3)	3	3,1,1
	float	2	2,15
	N1	4	57,45,28,46
	C2	3	7,6,11
	C3	3	28,39,1
TCNQ	C4	3	7,5,1
	C5	6	28,4,34,29,7,24
	float	2	7,1
	other	1	5
	C1	4	4,2,1,1
	C2	5	29,26,17,23,2
Durene	C3	2	97,186
	C-CC	5	2,3,1,1,1
	float	4	13,1,7,18

Equations

$$r_{min}^{(ion)}=s\frac{r_{vdw}^{(Mu)}+r_{vdw}^{(ion)}}{2}$$

$${\mathbf{f}}_i={\mathbf{R}}_i\mathbf{f}+{\mathbf{t}}_i$$

$$M_{ij}=\frac{|C_i\cap C_j^{\prime }|}{\sqrt{n_i{n}_j^{\prime }}}$$

Full text

In this work we test this methodology further by applying it to three organic molecular crystals model systems: durene, bithiophene, and tetracyanoquinodimethane (TCNQ). Using the same sets of random structures we compare the performance of Density Functional Theory software CASTEP and the much faster lower level approximation of Density Functional Tight Binding provided by DFTB+, combined with the use of the 3ob-3-1 parameter set. We show the benefits and limitations of such an approach and we propose the use of DFTB+ as a viable alternative to more cumbersome simulations for routine site-finding in organic materials. Finally, we introduce the Muon Spectroscopy Computational Project software suite, a library of Python tools meant to make these methods standardized and easy to use.

if keyword

1 Introduction

Muon spin spectroscopy is a technique in which a beam of polarized positive muons is implanted in a solid sample, often with the additional application of an external magnetic field, and the pattern of positron emission caused by their decay is then observed and used to measure the properties of the material. Since the positron emission depends on the direction of the muon’s spin at the time of its decay, observation of its time evolution can allow one to infer information about the magnetic structure of the sample thanks to the hyperfine interaction between electrons and the muon.

A key step of understanding muon spectra is to predict where in a crystal would an incoming muon effectively come to rest. This stopping site problem is well known and has been tackled in a number of different ways 1, including by use of simulations of various levels of accuracy. In our previous work 2 we proposed a new method to approach this problem, by combining the Ab-Initio Random Structure Search methodology, or AIRSS 3, 4 with a new random generation technique making use of a Poisson sphere distribution 5 and the clustering techniques implemented in the Python library Soprano 6.

The main performance cost in the proposed approach consisted in the use of Density Functional Theory (DFT) software to perform the necessary geometry optimization calculations. In this paper, we compare the results obtained with that approach to those obtained with a higher level approximation, Density Functional Tight Binding (DFTB) for a specific class of materials, namely, organic molecular crystals. Our proposal is that DFTB can be used as a computationally cheaper substitute for DFT, thus making calculations for systems under one hundred atoms easily feasible in a few hours on a single desktop machine, and opening the way to the study of larger systems like proteins and polymers. We also introduce new clustering techniques more suitable to treating systems with a large amount of possible stopping sites.

2 Background

While muon spectroscopy is often used to study inorganic systems with peculiar magnetic properties such as high temperature superconductors, there is an interest in using them in organic materials as well 7, 8. When a positive muon is implanted in a material it can either stop as is or capture an electron and form a bound state that is known as ”muonium”. Muonium acts effectively like a lighter hydrogen atom, meaning it possesses the same chemical properties and tends to react in the same ways. The main differences in behavior concern its dynamical properties and the entity of quantum delocalization effects it undergoes - both things dependent on its mass. In this paper, we focus on equilibrium stopping sites of muonium treated classically, and thus both these differences are irrelevant, and muonium can be treated as if it was a hydrogen atom.

In our previous work, we exploited this similarity to treat crystalline systems with an added muon with the DFT software package CASTEP. The muon, in these calculations, was simply represented by a labeled hydrogen atom. The core idea was to use the AIRSS methodology - create a large number of copies of the same host structure in each of which the muon occupied a different starting position and then carrying out a geometry optimization to find the closest energy minimum. This method is justified by the observation that in general the energy landscape for atomistic systems seems to be structured in such a way that the attraction basins for each energy minimum have a size that scales with the depth of that minimum 9, 10, 11. Therefore, by starting with random configurations, one has a high likelihood to find the lowest energy minima.

If we consider the nature of this problem, it is not unreasonable to think that it could be acceptable to replace DFT with a cheaper computational method, as long as some key properties were guaranteed. While a new method could describe an overall different energy landscape for the same system, one only needs it to satisfy two conditions for it to be a viable replacement for DFT in this scenario:

that the resulting landscape has at least as many or more energy minima as the DFT landscape; 2. 2.

that for each minimum in the DFT landscape there exists at least one minimum in the new landscape that falls close to it, or at least within the same attraction basin as defined in the DFT landscape.

Given these two conditions, it is easy to see how one can gain much by replacing DFT with the new, cheaper method. If, for example, 100 random structures are generated, one could find that after optimization they only result in, say, 10 minima. These 10 minima can then be further optimized with DFT, leading to the same results that we would have obtained from a pure AIRSS+DFT calculation. If the new method is cheap enough that its computational cost is trivial with respect to that of a DFT calculation, this results effectively in an almost tenfold increase in speed compared with optimizing all the original 100 structures with DFT.

In this paper, we propose that, at least for organic molecular crystals, Density Functional Tight Binding (DFTB) 12, 13 might be such a method. To validate this idea we compare the muon stopping sites found with both DFT and DFTB starting from the same random generated configurations for three organic molecular crystals.

DFTB is an electronic structure method that, much like DFT, makes use of the Kohn-Sham approximation to solve the quantum many body problem for electrons. It is however approximate in that it only represents these electrons as bound to the ions, and does not describe the full spatial wavefunction. This allows it to be computationally cheaper than a DFT calculation. It is also not ab initio, because it makes use of parametrizations computed from pure DFT calculations to describe interactions between chemical species 14. These parametrizations are stored in sets of so-called Slater-Koster files which are computed to cover only specific groups of elements, and thus can not treat any other species. The choice of organic molecular crystals was in fact prompted also by the availability of a well-documented Slater-Koster parameter set for organic compounds in the 3ob-3-1 set 15, 16, that covers among the other elements carbon, hydrogen, oxygen, nitrogen and sulphur. Success with this specific subset of chemistry could encourage the creation of more system-specific parametrizations in the future to treat other classes of compounds.

3 Calculations

3.1 Choice of molecules

In this paper, we focused on three exemplar, crystalline, molecular-organic systems: tetracyanoquinodimethane (TCNQ), bithiophene and durene. These systems represent the class of organic crystals that are most studied using muon techniques while, at the same time, being of enough individual practical interest.

TCNQ is a cyanocarbon with formula $\mathrm{(NC)_{2}CC_{6}H_{4}C(CN)_{2}}$ (Figure 2). It crystallises around a temperature of 566 K in a structure with four molecules per unit cell (Z=4) and space group C2/c (Figure 3). It can form charge-transfer salts that have relatively large electrical conductivity. For instance, when TCNQ is combined with the electron donor tetrathiafulvene (TTF), it forms the TTF-TCNQ complex, in which TCNQ is the electron acceptor 17. Muon experiments were performed on this compound 18 to investigate its electronic structure; a key step of interpreting them is then to determine the relationship between the muon radical states and the electronic states that are formed in the muoniated molecule.

2-2’-bithiophene is the dimer of thiophene, an heterocyclic compound whose polymers are of great interest due to their conductive properties. Bithiophene has formula $\mathrm{(C_{4}H_{3}S)_{2}}$ (Figure 4) and crystallizes at 306 K in a structure with two molecules per unit cell (Z=2) and space group P2_1/c (Figure 5). Using the electrochemical polymerization method, bithiophene is used as a precursor for polythiophenes (PTs) more often than thiophene itself19, and muons have been used to study the microscopic charge transport processes in polythiopene-based polymers. In particular, muons have been successful in explaining the 1-dimensional intra-chain charge transport mechanisms and 3-dimensional inter-chain charge hopping effects that are related to the electrical and optical properties of these polythiopene-based polymers20.

Finally, durene (1,2,4,5 tetramethyl-benzene) has formula $\mathrm{C_{6}H_{2}(CH_{3})_{4}}$ and crystallizes at 352.3 K in a structure with two molecules per unit cell (Z=2) and space group P2_1/c (Figure 7). It was one of the first substituted aromatic compounds to be studied with muons, and two muon stopping sites were suggested, located at the C-H sites in the ring, from above and below molecular plane 21. Aromatic compounds with different substituting groups such as durene are a useful model to study classical and quantum dynamics of muons in molecular solids. Combined experimental and theoretical studies on crystalline benzene showed the effect of classical dynamics, but no final evidence of quantum tunneling was observed 22. Comparison to similar results in a family of selectively substituted benzene-based crystalline solids, with modified dynamics due to mass and steric effects, can help elucidate the issue.

In all these materials, the main effects of muon implantation are the geometrical distortion of the surroundings of the muon stopping site and the presence of an unpaired electron in the crystalline structure. The unpaired electron in these muoniated compounds goes to the highest un-occupied molecular state. If the extent of the distortion is small, the newly occupied molecular state in the muoniated compounds can be used to study the occupied electronic states that result from charge transfer. Hence, knowledge about the muon stopping sites in these materials, as well as about the parameters of the hyperfine interactions between the muon and the molecular electronic states, can be very helpful. The main structural characteristics of TCNQ, bithiophene and durene are summarized in Table 1.

3.2 Geometry Optimization

Random starting configurations have been generated using the Poisson sphere algorithm originally described in a previous work 2. This algorithm guarantees that all starting positions for the muon are at least at a distance $r_{min}$, here chosen to be 0.8 Å. Additionally, a constraint was enforced to keep the generated muons at a distance from the ions equal to:

[TABLE]

where $r_{vdw}$ are the Van der Waals radii used for the given elements (in this context, muonium is treated as hydrogen) and $s$ is a scaling factor that in this case was set to $s=0.5$. The values of the used Van der Waals radii for various elements are reported in Table 2.

The generation process continues until the entire free available space is filled and any further attempt to add another muon fails to meet the constraints. The process generated 396 different starting structures for bithiophene, 448 for durene, and 529 for TCNQ.

Geometry optimization on these structures was then performed with a BFGS algorithm, with fixed unit cell parameters, and using both CASTEP 8.0 and DFTB+ 17.1, to a tolerance of 0.05 eV/Å for the forces. For CASTEP 23, the PBE exchange-correlation functional was used, in combination with a standard Tkatchenko-Scheffler scheme 24, 25 for dispersion forces and auto-generated ultrasoft pseudopotentials.

DFTB+26 is a code that implements the Density Functional Tight Binding algorithm plus additional features, such as dispersion forces corrections. For these calculations, we made use of the self-consistent charges scheme 27 and 3rd order corrections 28. As mentioned already, the Slater-Koster files were those of the 3ob-3-1 parametrization 15, 16. We also experimented with the use of dispersion correction 29, but it caused convergence to become extremely slow, possibly because of the small size of the unit cell used. Since our main purpose was to test an alternative to CASTEP that would be as quick and cheap as possible -while still producing similar results- we decided not to use dispersion corrections in DFTB+ calculations.

The plane wave cutoff for CASTEP as well as the size of the k-point grid for both CASTEP and DFTB+ were chosen by converging energy and forces to a tolerance inferior to that required by the calculation for each individual system. Values of 700 eV for the CASTEP plane wave cutoff and 2x2x2 for either k-point grid were found satisfying and used for all systems.

3.3 Clustering

In a previous work 2 we already showed how random structure searching could be used fruitfully to identify muon sites in some simple semiconductor crystals. In that work, the optimized structures were analyzed using the software library Soprano 6 and applying basic clustering techniques.

In this paper we expand on that methodology by establishing a work flow of successive classification processes based on discrete variables or the clustering of continuous ones. The main reason for this is that the number of predicted sites for the molecular crystals is much higher than anything seen in the simple crystals treated before. While a single clustering step can indeed be used, meaningful interpretation of the data with that approach becomes much harder. It was therefore deemed more appropriate to use this multi-step process in order to highlight the physically meaningful differences between groups of clusters. The full workflow is illustrated in Figure 8. The key steps are as following:

the number of individual organic molecules in the optimized structure with added muonium, $Z^{\prime}$, is computed by making use of Van der Waals radii (methods included in Soprano) and compared to the number of molecules in the structure without defects, $Z$. If $Z$ and $Z^{\prime}$ are different, this is a sign that:

a) the muon has not bonded with any molecules (in which case $Z^{\prime}=Z+1$) and it is sitting in an intermolecular equilibrium configuration that we termed a ’floating site’; or

b) that the muon has acted as a bridge joining two molecules (in which case $Z^{\prime}=Z-1$), or that some other process by which, for example, molecular fragments have detached has taken place. Of these processes, the ’floating site’ one is the only one observed consistently through both CASTEP and DFTB+ simulations in significant numbers of structures. The three groups ($Z^{\prime}=Z$, $Z^{\prime}=Z+1$, and every other case) are separated, and the latter two are subjected directly to step 3. 2. 2.

in structures in which $Z^{\prime}=Z$, the muon has bonded with one of the existing molecules. For these, the molecule containing the muon is considered and the bonding site is identified automatically analysing the bonding structure, using a method suitable to each system and the library NetworkX for network analysis 30. For TCNQ and bithiophene, the network distance from the muon and the two nitrogen and sulphur atoms respectively was considered. For durene, it was the bonding structure of the nearest carbon. The details of these methods and how they identify the main sites are explained in the Supplementary Information (S1). Figures (2), (4) and (6) show the numbered molecular sites where the muon can bond to the different molecules while maintaining the structural integrity of the molecules. There are, however, particular cases where the muon bonds and breaks the structural integrity of the molecule, i.e.: breaks an aromatic ring. In these cases, the identification string produced by the method was left as is. 3. 3.

finally, a number of sites can be reduced to their ‘asymmetric unit coordinates’ in order to identify equivalent sites. This is done by taking the symmetry operations for the ‘pure’ structure as a series of rotation matrices and translation vectors $\{\mathbf{R}_{i},\mathbf{t}_{i}\}$. Then for a muon site expressed in fractional coordinates $\mathbf{f}$ we can compute all the periodic images:

[TABLE]

and retain only one that is closest to a chosen reference point. In this way, all apparent differences due to symmetry operations are removed, and effectively equivalent defects will be squashed in close proximity one to each other, making clustering operations based on position easier. This functionality is provided in Soprano, by the function soprano.utils.compute_asymmetric_distmat and makes use of the symmetry analysis library Spglib 31.

Clustering was carried out with the k-means method, using a previous hierarchical step as a way to estimate the number of expected clusters, similarly to what was done in our previous work 2. After classifying the sites with clustering, the clusterings obtained from CASTEP and DFTB+ results need to be compared. This is done in two ways, using both a distance matrix and a normalized confusion matrix 32. Both matrices have the purpose of comparing each possible pair of clusters (one from the CASTEP, the other from the DFTB+ results) with some metric describing their similarity. For the distance matrix we first compute an average muon fractional position for each cluster in both sets, and then compute the distance between all pairs of these centers. This is done while carefully accounting for all symmetries of the crystal to make sure that all such distances are computed between the shortest equivalent sites, using the same asymmetric unit coordinates approach described above. Since these are only fractional coordinate distances, they don’t really have a physical meaning, but they correspond to the same criterion with which the clusters were formed, and thus we expect two clusters describing the same stopping site to come out very close with this metric. The closer to zero the distance, the more similar the site represented by the two clusters.

On the other hand, the normalized confusion matrix has the purpose of evaluating how many of the structures, starting from the same point, end up converging to the same cluster through the optimization process. For two clusterings $\mathcal{C}$ and $\mathcal{C}^{\prime}$, the normalized confusion matrix element $M_{ij}$ is defined as:

[TABLE]

where $|C_{i}\cap C_{j}^{\prime}|$ represents the amount of structures shared by clusters $C_{i}$ and $C^{\prime}_{j}$, which belong to the first and second clusterings respectively, and $n_{i}$ and $n^{\prime}_{j}$ are their respective number of elements. If clusters $C_{i}$ and $C^{\prime}_{j}$ are identical, then $n_{i}=n^{\prime}_{j}=|C_{i}\cap C_{j}^{\prime}|$ and $M_{ij}=1$. Conversely, if no structures are shared between $C_{i}$ and $C^{\prime}_{j}$, $M_{ij}=0$. Therefore, the confusion matrix for two identical clusterings would be square, have only one element set to one per row and column, and every other element would be zero.

4 Results and Analysis

Table 3 shows the main features of all the clusters obtained using CASTEP and DFTB$+$ calculations.

The clusters are grouped by ‘type’, meaning whether they’re floating sites (float) bonded to some other atom, or anything else (‘other’). Bonded sites are labelled by element and site number for those that are recognizable (the numbers are the same as in figures 4, 2 and 6 respectively for bithiophene, TCNQ and durene) (see Suplementary Information S2 for examples), and by the signature of the function used to identify them when it did not correspond to any recognizable site. The first important thing to notice from table 3 is that all the bonded and float-type clusters predicted by CASTEP were also predicted by DFTB$+$, i.e. there are no false negatives. In the case of bithiophene and durene, DFTB$+$ predicted some extra bonded structures. For bithiophene these are indicated by the distances between the muon and both sulfur atoms in each molecule. These distances are indicated between brackets as the number of consecutive bonds present between the muon and each sulphur atom in each molecule. The extra bonded structures in bithiophene then are S(1,7), S(5,6) and S(1,3). These linking patterns would not be possible unless some of the molecular bonds were broken - some specific examples are offered in S1 in the Supplementary Information. In general, it seems like DFTB+ has a higher tendency than CASTEP to break bonds in unphysical ways. However, it’s easy to dismiss such sites on account of their tiny numbers compared to the other clusters. The C-CC cluster in durene represents a muon bonded to a carbon that is, in turn, bonded only to two more carbons, probably the result of the muon replacing one hydrogen in the aromatic ring, another unlikely reaction which relies on bond breaking. Similarly, the ’other’ cluster identified for TCNQ, corresponds to a muon site that is neither bonded nor floating. This cluster is formed by a small number of structures where the muon has, for instance, connected two TCNQ molecules. All of these examples can be considered as false positives when using DFTB+ and, while they clutter the space of solutions, they don’t cause us to lose any real important information.

It becomes important then to compare the CASTEP and DFTB+ clusterings to verify how well they map onto each other. As explained in Section 3.3, we do this in two ways: a distance matrix and a normalised confusion matrix. These two tools are complementary in covering different ways in which the clusters can be considered similar. The distance matrix describes how close two clusterings are to describing the same muon stopping sites; the confusion matrix on the other hand gives us a sense of how much two different clusterings overlap in terms of which structures end up in any given class. Since the structures here are labelled based on their initial configuration (identical for both CASTEP and DFTB+ calculations), the latter gives us an idea of how similar the potential landscape is for the two calculation methods. Fundamentally, if the distance matrix shows a good correspondence, that means the two random searches end up converging to similar minima; but if the confusion matrix shows the same correspondence, that also means that the same starting points end up in the same minima, which is a far stronger sign of overlap between the two methods.

Figure 9 shows the distance and normalised confusion matrices for the clusterings obtained with DFTB+ and CASTEP. The correspondence between clusters obtained with the two methods, for the molecular and floating sites, can be seen clearly in the form of black squares (representing low distance or high overlap, respectively, for distance and confusion matrices). The sites that tend to lack an assignment are generally the ones that appeared only in DFTB+ (such as the unconventional $S(1,3)$ in bithiophene or the C-CC in durene). The C1 site in durene does the same, and while in this case it appears in CASTEP too, it still requires the unbinding of a hydrogen from a $CH_{3}$ group to happen, a process that seems unlikely to happen on the time scales of a real muon experiment. Correspondence between distance and confusion matrices, for the molecular and floating sites, is also clearly evident, an excellent indicator of the fact that the potential landscape between CASTEP and DFTB+ is similar. In other words, not only the minima are in similar places, but in general, if the muonium atom is placed in a similar starting spot, it tends to end up in the same one too after geometry optimization.

Figure 10 shows the clusters as dots positioned by average and standard deviation of the energies of the structures they contain. Labels are added to show the average energies for structures classified by site for CASTEP and DFTB+ results. Low standard deviations indicate a more ‘consistent’ cluster, with less dispersed energies, thus better convergence to the minimum. Energies are naturally different from different sites, though they tend to vary very little for crystallographically inequivalent realizations of the same molecular site. More importantly, CASTEP and DFTB+ results tend to match remarkably well in energy in a lot of cases. Low energy sites tend to match especially well (e.g. C5 and C3 in bithiophene, C3 and C2 in durene, and C3 in TCNQ). Some of the higher energy sites (S1 in bithiophene, C4 and C2 in TCNQ, and the floating sites in general) have growing discrepancies as the energies get higher. The biggest error appears for the N1 site in TCNQ, which is the lowest in energy in CASTEP but is approximately 0.6 eV higher in DFTB+.

5 Conclusions

In this work, the muon stopping sites predicted for bithiophene, TCNQ and durene, by using a DFT-based methodology, were compared with those predicted using a higher level approximation that uses Density Functional Tight Binding (as implemented in the DFTB$+$ code).

The potential muon stopping sites found by the DFTB$+$-based methodology agree very well with those found by the CASTEP-based methodology. Moreover, DFTB$+$ calculations are computationally much cheaper than DFT calculations and could potentially be used for treating very large organic systems such as polymers or proteins. This work offers a strong case that it would be reasonable to do so, while also highlighting what are the limitations that a user of the tight binding approach should be watching out for.

6 Acknowledgements

The authors are grateful for the computational support provided by (a) STFC Scientific Computing Department’s SCARF cluster; and (b) the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC (No. EP/P020194/1). Funding for this work was provided by the STFC-ISIS muon source and by the CCP for NMR Crystallography, funded by EPSRC Grant Nos. EP/J010510/1 and EP/M022501/1.