An electric field-based approach for quantifying effective volumes and radii of chemically affected space.

Chemical shape and size play a critical role in chemistry. The van der Waals (vdW) radius, a familiar manifold used to quantify size by assuming overlapping spheres, provides rapid estimates of size in atoms, molecules, and materials. However, the vdW method may be too rigid to describe highly polarized systems and chemical species that stray from spherical atomistic environments. To deal with these exotic chemistries, numerous alternate methods based on electron density have been presented. While each boasts inherent generality, all define the size of a chemical system, in one way or another, by its electron density. Herein, we revisit the longstanding problem of assessing sizes of atoms and molecules, instead through examination of the local electric field produced by them. While conceptually different than nuclei-centered methods like that of van der Waals, the field assesses chemically affected volumes. This approach implicitly accounts for long-range fields in highly polar systems and predicts that cations should affect more space than neutral counterparts. This journal is © The Royal Society of Chemistry.

Introduction

Chemical shape, size, and accompanying surface area and volume are fundamental properties that govern a wealth of intermolecular phenomena in atoms, molecules, and materials.[1-7] Determination of chemical volumes and surface area is, however, obscured by the definition of the atomic “surface”.[8,9] Today, there are numerous approaches to quantify chemical size and while the van der Waals method is certainly the most prevalent, alternatives have been developed through a synergy of experiment and theory.[9] Understanding the limitations and applicability of these alternatives is important because numerous advanced measurements implicitly rely on size in some form (e.g., specific surface area of gas molecules used in surface area and volumetric measurements of porous materials within the Brunauer–Emmett–Teller (BET) formalism,[10] quantifying void-space docking sites in enzymes,[11,12] and so forth). Hence, there remains intrinsic value in revisiting this age-old problem: how big are molecules?

Initial experimental measurements of atomic size were performed by Meyer in 1870,[13] where he identified a relationship between material density and atomic size, and obtained a periodic trend in atomic volumes. These values were later refined by Bragg[14] and Pauling,[15] who developed methods for assessing atomic radii through X-ray scattering. In two separate works, Bondi and Batsanov revisited the radii presented by Pauling and Bragg, and it is these works that are synonymous with the “van der Waals (vdW) radii” of atoms and ions in the solid state.[16-18] An alternate approach was taken by Alvarez,[19,20] Biswas, and Ghosh,[21] who extracted atomic radii using statistical analyses of online databases. The generality of this approach is limited however because of the uncertainty for chemical environments not represented in the those data,[9] for systems featuring regions of high polarity,[22-24] molecules with elongated bonds,[25-27] atoms under high pressure,[28,29] and other exotic chemistries. Further, it is not immediately evident that crystal-derived atomic size[30] is necessarily applicable to gas and solution phase atoms and ions.[31]

First principle simulations have provided another alternative, boasting inherent generality to recover volumetric data of both known and undiscovered molecules from their computed electronic structures while also in principle providing some insight into gas phase radii.[32] Early examples of these calculations were presented by Slater who employed the maximum radial density of outermost single particle wavefunctions to define atomic radii,[33-35] and several other related methods have also been reported.[36-41] Among these advances, Bader computed elemental size at the Hartree–Fock level of theory.[42] There, the surface of a chemical system was defined by the deduced electron density at a cutoff = 0.002 e bohr−3. This cutoff was refined by Boyd[43] to be 0.001 e bohr−3 within the Density Function Theory (DFT) construct,[44,45] with the justification that any smaller value of electron density would result in a negligible change in calculated radii. Later, Rahm and Hoffmann furthered Boyd's work, applying the method to atomic ions using electron densities obtained from hybrid-GGA DFT[9] (PBE0,[46-48] with a large basis set).

Yet the use of electron density alone poses problems for modeling cations, which certainly interact with their surroundings beyond their electron cloud; they create a large electric field.[49] With this in mind, we thought to revisit the size quantification problem through examination of the electric field, a value computed from the derivative of the electrostatic potential. While other approaches concern the space that a chemical system occupies, an electric field description captures the volume that a chemical system affects. In this regard, both cations and anions should be larger than their charge neutral counterparts (i.e. they affect more space) because their coulombic charge increases, polar bonds should produce larger fields than non-polar analogues, and size/shape should be affected by external fields. Additionally, the effective size of the atom should still depend on quantum chemical properties, such as orbital filling and effective core charge (i.e. the charge experienced by valence electrons as described by Slater's rule), and so forth. Thus, an electric field metric should provide a conceptually different description of chemical size, and a unique approach for defining the edge of a chemical system. Herein, we explore the generality and implications of quantifying chemical size through examination of atomic and polyatomic electric fields, and provide some examples of the utility of this approach.

Results and discussion

Method description

To sample the electric field, we have developed a post-electronic structure processing software called STREUSEL (Structure Topology REcovery Using Sampling of the ELectric field),[50] which computes chemically “affected” volumes and surface areas of atoms, molecules, and materials. Since the electric field is defined as the negative gradient of the electrostatic potential it is highly sensitive to subtle changes in polarization. Here, the edge of the chemical system is defined as the point in space where there is near-zero variance in the electric field magnitude (variance is computed by considering magnitude of the field in neighboring volumetric pixels, voxels). Conventional and ab initio calculations return reliable electrostatic potential values on the order of 10−6 eV per electron (i.e. ∼2.3 × 10−5 kcal mol−1), thus we consider a change of less than 10−5 eV per electron (∼2.3 × 10−4 kcal mol−1) to be conservative. Conceptually, as cutoff decreases more vacuum space becomes included in the regions associated with the molecule, and while it is undoubtedly a variable, we justify our 10−5 eV potential cutoff thermodynamically; the energy cutoff is on the order of kB, around 10−5 eV K−1. For reference, a typical van der Waals interaction is on the order of 0.956–1.912 kcal mol−1.[51,52] Our cutoff accounts for fluctuations in the third decimal place, in principle providing sensitive, but still experimentally measurable information.

To employ this cutoff, one must first sample the electrostatic potential generated from a DFT or ab initio calculation, with the density computed at a discrete number of voxels. Like the field cutoff, there is a dependence of the voxel size on the computed molecular size coming at a trade-off between time-to-solution and voxel resolution, see Fig. S1.† From these data, a voxel of volume 0.008 Å3 yields a desirable balance between computation time and volumetric resolution, while being sufficiently high fidelity to describe rapid changes in field across conventional chemical bonds.

To illustrate the conceptual difference between our approach and other size metrics, we present the Mg0 and Mg2+ vdW radii alongside the electric field radii computed using CCSD-full[53-55]/aug-cc-pVTZ,[56]Fig. 1. Conventional chemical tenets suggested by Batsanov would indicate that Mg2+, and cations in general, are smaller than their neutral counterparts (see Table S1, and Fig. S3† for a Mg-specific example). Yet, from an electrical field perspective, Mg2+ should be significantly larger than Mg0, because the electric field ultimately depends on the ratio of the number of protons-to-electrons. In other words, the area affected by a cation should be large, while the density of electrons should be small. Hence, the Mg2+ example serves as an illustration that the field-defined size is not “atomic size” in the conventional sense and provides different insights than other existing atomic size models.

Fig. 1

Comparison of atomic radii of Mg and Mg2+ between the Batsanov metal compressibility method and the electric field-based method presented in this work (isosurface drawn from our cutoff using data from CCSD-full/aug-cc-pVTZ data). A coordination number of 8 is used for the Mg2+ from Batsanov. See Fig. S3 and Table S1† for a comparison of other size comparisons of Mg0 and Mg2+.

The computed volume, however, should heavily depend on familiar quantum chemical variables such as functional, basis set, density grid size, etc. Size should also depend on quantum chemical aspects of orbital filling. For cations, we expect that the radius should be inversely proportional to the electronic screening of the nuclear charge,[57] whereas for anions we hypothesize that the radius will be governed by both the magnitude of charge as well as some consideration for which orbitals are occupied (e.g. one might expect that anions with 4s valence to be larger than an isoelectronic 3d valence, see Fig. S9†). We will explore these and other parameters throughout this paper.

Atoms

Most sizing methods delineate between charged and charge neutral chemical systems, and free atoms and those in bonded environments. For example, the Pyykkö family of atomic radii[58] should only be used for covalently bonded systems, the Shannon and Prewitt[59,60] formalism is reserved for oxidized metals, and the vdW radii presented by Bondi are useful for charge neutral atoms. To contrast our method to these and other reported approaches, we compute the periodic table of elements using CCSD-full/def2-TZVP as implemented in Gaussian09 (ref. 61) (see Fig. S4† for a periodic table of elements) alongside some comparable and widely used methods, Fig. 2. Although other calculations throughout this paper invoke the aug-cc-pVTZ basis, the def2-TZVP afforded up to Rn. Fortunately, differences in these triple zeta bases yielded negligible difference in computed size (see Table S2†).

Fig. 2

A comparison of atomic radii recovered using various size metrics. (a) Sizes computed from STREUSEL are like other electronic structure-based methods, and (b) comparable to the vdW (Batsanov) sizes. Alternative size metrics (singly-bound data are presented for Pyykkö, and Boyd's unscaled approach) define the upper and lower limits. Inset: the first row transition metals show a contraction between V and Cr, attributed to the transition from 4s2 to 4s1 orbital filling.

Generally, the atomic radius computed from the electric field is comparable with other electronic structure-based approaches (those presented by Alvarez,[19,20] Boyd[43] and Rahm,[9]Fig. 2a). For the most part, our approach seems to predict similar size trends for free neutral atoms. Fig. 2b reveals larger deviations between compared methods; Boyd's unscaled DFT-derived sizes and Pyykkö radii for singly-bonded atoms bracket the radii of neutral elements. Pyykkö's singly-bound radii serves as an important reminder that the local environment plays a significant role in determining the interatomic radii. We would expect that the electric field produced by atoms in interacting environments would also deviate from the neutral sizes, Fig. 2. The other striking conclusion is that the vdW radii seem to align well with electric field-derived sizing for the first 20 atoms, with subtle deviations in heavier main group elements. Indeed, the general shape of the predicted atomic sizes using STREUSEL seems to obey a trend common to other electronic structure-based methods (per Fig. 2a).

Ions

Most atoms exist in some formal, non-zero oxidation state. Previous work by Sen and Politzer studied the radii of anions by examination of their electrostatic potential.[62] In that study, the authors highlighted that anions have a minimum radius in which the nuclear potential is cancelled by the radial electron charge. Moving beyond this minimum, the potential progresses towards zero. In contrast, cationic potential decays steadily towards zero and governs the potential produced by both neutral and cationic species. By examining the magnitude of the field, however, we can make comparisons between anions and cations directly. To complement Sen and Politzer's study, feedback provided on the preprint of this paper noted that if classical electrostatics (i.e. Coulombs law) were to dominate the field magnitude, then one could expect the radius of C2+ to be half that of C4+. To assess this, we computed the radius of the range of carbon-based anions, C4−⋯4+, Fig. 3. These data point to two outcomes; (i) the radius of C2+ is much greater than half that of C4+ indicating that the effect is not purely classical, and (ii) anionic and cationic C of equal charge magnitude have roughly the same radii. However, the anions are consistently smaller than the cations, with the most striking difference being the singly charged systems C− and C+. We interpret this result to be determined by the extent of 2s–2p electron mixing: C− should maximize mixing between 2s and 2p states (as each p orbital is singly occupied), resulting in field contraction. Since both C− and C+ carry the same magnitude of charge, but the valence electron in C+ has reduced s–p mixing, the valence electron should have elevated energy hence creating a marginally larger field.

Fig. 3

The radius of carbon as a function of oxidation state. Cations and anions of equal charge magnitude have approximately similar radii, but differences can be attributed to the mixing of s–p orbitals. Radii are presented in italics, effective core charge is shown for C− and C+.

Additionally, we computed ionic radii for some transition metals in the periodic table, Fig. 4, and compared our computed size to yet another conventional volumetric approach presented by Shannon and Prewitt.[59,60] The Shannon–Prewitt method shows a clear dependence on the oxidation state of the ion, but convolutes a direct comparison through the inclusion of a number of inner-sphere ligands (with a general relationship between increasing number of inner-sphere ligands and reducing ionic radii). Given our approach computes size in the absence of ligands, we expected that STREUSEL radii should be larger than those from Shannon–Prewitt, Fig. 4b. From our method, the key observation is that radii are proportional to atomic charge (e.g. the extent of polarization), not necessarily identity of the atom. As an example, this contrasts with the Shannon–Prewitt finding that both Ti2+ and Zn2+ are larger than Fe3+. There is also a dependence on orbital filling, with larger radii predicted for Sc, Ti, V, and Cr when computed in non-ground state electronic configurations (see Table S4 and Fig. S9†). Another conclusion from these data is that the oxidation states appear in bands. For example, Fe2+/Fe3+ and Co2+/Co3+ have similar differences in radii. These values will of course dramatically reduce in the presence of external fields, highlighting the Shannon–Prewitt-type dependence on local coordination.

Fig. 4

A comparison of ionic radii in the d-block. (a) The Shannon–Prewitt sizes depend on oxidation state and show a reduced size for mid-block elements (computed for 8-inner sphere ligands). (b) The electric field sizes are computed using the free ion and show a general relationship in radius to oxidation state.

Bonded atoms

Like the Shannon–Prewitt method, an electric field-derived approach should show a change in chemical size depending on proximity to other electric fields (e.g. those produced by ligands, atoms, molecules, surfaces of materials etc.). For example, the field produced by a free Li+ should be larger than Li+ in proximity of Cl− simply because they interact with one another. We can demonstrate that the size of ions depends on atomic proximity through progressive increase in bond lengths of some simple diatomics (LiF, LiCl, LiBr), Fig. 5a.

Fig. 5

(a) STREUSEL volumes for lithium halides at arbitrarily expanded bond lengths. At infinite separation the volumes are equivalent to Li+ and the corresponding halide anion. Even at 4 Å, separation the ions are affected by one another, as evidenced by reduced field size, and atomic volume. (b) Comparison of molecular volumes for a series of geometrically equilibrated diatomic molecules with varying dipole moments. Divergence between electron density-based size metrics and STREUSEL depends on the system polarity and atomic electronegativity. (c) The difference in size between Boyd and STREUSEL for LiCl.

Here, the volumes of LiF, LiCl, and LiBr increase with increasing bond length. They inevitably converge to the size of the sum of Li+ and X−, computed either separately or very far apart within the same calculation (>10 Å, to exceed the long-range cation interactions[63,64]). While the trend of increasing size with increasing interatomic separation is common to both the vdW and STREUSEL approaches, Fig. S5,† the volumetric trend of the dissociated lithium halides reveals an opposing trend in atomic sizes based on electron density. STREUSEL predicts that the radius of F− > Cl− > Br−. This can be readily understood from Slater's rules by examination of the effective screening of core charge. The valence electrons experience an effective nuclear charge (Zeff) of 4.85, 5.75, and 7.25 for F−, Cl− and Br− respectively. Hence, a significantly lighter ion of similar charge should create a larger field. Noting that this is at odds with the generally accepted metric that Br− is the largest of the three, this serves as a prime example of the unique perspective provided by the electric field. It is an inherent measure of affected size, and seemingly follows the general trend offered by Lewis' hard/soft acid/base theory (i.e. soft ions are indicative of small fields).

We can further highlight this relationship by examining the predicted volumes of HF, HCl, HBr, Fig. 4b. These data reveal two key observations; (i) the acidic halides have near-zero dipole moment, resulting in volumes predicted by STREUSEL to align more closely with electron density-based methods (like data shown in Fig. 2a for charge neutral atoms), and (ii) large dipole moment results in significant increase in volume. The latter also results in a significant shift in the regions of the molecule producing this volume, Fig. 5c and S6.† In summary, ions of comparable mass with identical oxidation state will have similar electric fields, and hence similar radii. As the atom becomes heavy, the effective core charge increases, resulting in smaller fields for identically charged ions. Finally, as two chemical system approach, the fields polarize one another. To investigate the latter we now turn more complex molecular systems.

DFT and applications to complex systems

The utility of the presented approach is somewhat limited by the dependence on using CCSD-full and a triple zeta basis. The following section seeks to identify alternate, less computationally intensive DFT-based methods that provide reliable size comparisons which can be used for large molecules and materials. Such an approach would boast the benefit of increased time-to-solution, and the ability to be applied to more diverse chemical systems.

As with all electronic structure methods, the size and shape of molecules depends on both the functional and basis set used for geometric equilibration. For this purpose, our basis remains constant (aug-cc-pVTZ), and the electronic structure method is explored. In a recent assessment of 128 DFT functionals, the authors canvassed the method performance for recovery of total energy and, separately, electron density.[65] That paper highlighted that approximations in the DFT formalism[66] leads to significant energetic dependency on functional (having an indirect effect on the shape of the molecule).[67] The latter is true because atomic position is determined by electron density, and electron density is determined by atomic position (i.e. the self-consistent field and geometry optimization routine). In this context we are interested in molecular volumes, surface areas, and shape, with energetics playing a secondary role (although there are certainly obvious future studies that may harness the energetics of interactions from electric field overlap).

To arrive at an ideal method to recover size using STREUSEL, we examined 15 neutral molecules (Ne, H2, N2, F2, Cl2, Br2, H2O, H2S, NH3, CO, CO2, CH4, C2H2, CNCl, and SO2) using 49 DFT and ab initio methods (Table S3†). The selection of small molecules spans a range of polarizations and oxidation states. While the methods in Fig. 6 only canvas the upper rungs of Jacob's ladder,[67] eight electronic structure classes are included (ab initio, generalized gradient approximation (GGA), generalized gradient exchange (GGE), hybrid-GGA (H-GGA), hybrid-meta-GGA (HM-GGA), local-density approximation (LDA), meta-GGA (M-GGA), and range separated functionals). Like before, CCSD-full is used as our geometric and electronic reference for the exact solution, and each molecule was geometrically equilibrated using the stated functional. Its volume, radius, surface area, and other topological properties were then computed from the electrostatic potential.

Fig. 6

The calculated sizes for the small molecule systems shown, compared to the volumes computed using the structure and volume from CCSD-full. The x-axis presents the natural logarithm of the volume fraction, y-axis groups functionals based on their electronic structure method.

As we ascend Jacob's ladder, higher level DFT functionals do not immediately appear to outperform lower-level ones, excluding ab initio methods which are highly accurate. Generally, GGE, HM-GGA, and M-GGA functionals appear to systematically underestimate molecular volumes (revealed in the mean volume deviations, Fig. S7†), and F2 and Ne are anomalously overestimated independent of DFT method. It should be noted that the axis in Fig. 6 is deliberately presented as logarithmic to show differences between functionals, and perhaps a better comparison is achieved by examining average mean deviations between functional classes, Fig. S7.† It becomes apparent that there is no clear preferred functional, at least not predictably so. It also highlights one plausible reason why the vdW approach has been widely adopted—there is very little dependence on minor fluctuations in bond length in predicted molecule volume (see Table S4†). In addition the mean absolute error for molecules optimized with B3LYP, a widely used functional for small organic molecules, appears to reliably approximate that of CCSD-full.

Instead of naming a clear champion, we acknowledge that it is difficult to compare sizes of molecular systems computed using different levels of theory. This, however, is no surprise—one would not compare the energies of a chemical system computed with two different functionals, nor should one compare their size. The utility of our approach is rather self-contained; the user should simply select a functional from a tractable rung of the ladder. One could select a functional using game theory,[68] or perhaps the old-fashioned way of simply balancing both accuracy and time-to-solution. With respect to Fig. 6, we are guided to favor functionals and ab initio methods displaying lowest mean absolute deviation within their functional class: MP3 and MP4 (ab initio); CAM-B3LYP (range separated); BMK (HM-GGA): BHandH, BHandHLYP, X3LYP, B1LYP, and B3LYP (H-GGA); TPSS and VSXC (M-GGA), BPL, G96LYP, and TPSSLYP1W (GGA); BVWN5 (GGE); Xalpha (LDA). However, one can use any functional within a given study. By way of example, the following case study invokes M06L, one of the worst performing functionals to obtain sizes comparable to that produced from CCSD-full.

Case study

In the case of ionic liquids, molar volume is a critical parameter which is thought to govern the physical properties of the bulk, including density, viscosity, and so forth.[69] The density and molar volume should depend on the size of the ions, which itself may depend on whether they are computed together or separately. It may be useful to predict the packing volumes from small molecule calculations, ideally of the free ions themselves, to overcome sampling of the various geometric configurations. By way of example, we consider two simple chemical systems, [BMIM][VCl4] and [BMIM]2[CoCl4] (BMIM = 1-butyl-3-methylimidazolium, Table 1).[70] Using the geometries computed using M06L, several volume metrics are presented. Both STREUSEL and vdW appear to predict similar sizes of the free ions, while electron density methods significantly under-estimate the sizes of the molecular ions.

Boyd, Bader, Batsanov and STREUSEL-derived molecular volumes (Å3) for the [BMIM]2[CoCl4] ionic liquid, BMIM = 1-butyl-3-methylimidazolium. The sum of single molecule volumes for the individual constituents ([BMIM]+, [CoCl4]2−) are presented, as well an equilibrium ion pair. The percent difference (%) is presented for each ionic liquid model

		Boyd		Bader		STREUSEL		Batsanov		Dipole
Single molecules	[BMIM]+	73.98		66.62		158.8		161.3		5.32
	[VCl4]−	64.37		58.24		87.94		107.3		0.055
	[BMIM]+ + [VCl4]−	138.3	19%	124.9	17%	246.7	0.16%	268.6	21%
Ion pair	[BMIM][VCl4]	114.8		105.1		247.1		217.6
Single molecules	[BMIM]+	73.98		66.62		158.8		161.3		5.32
	[CoCl4]2−	68.35		61.95		86.63		107.6		0.403
	2[BMIM]+ + [CoCl4]2−	216.3	22%	195.2	20%	404.2	3.6%	430.2	38%
Ion pair	[BMIM]2[CoCl4]	173.1		159.6		390.1		293.6

From Fig. 5, one may initially assume that the free molecular ions would always be larger than their paired analogues. Yet, within polynuclear systems there are three competing phenomena determining the size of the ions;

(i) Anionic and cationic electric fields interfering with one another. As electric fields of opposing magnitude interact, they should screen one another leading to a decrease in net volume.

(ii) A change in molecular shape due to intermolecular interactions. As molecules (charged or neutral) interact, they progress towards new equilibrium geometries. As a result, intermolecular interactions lead to an unpredictable change in shape.

(iii) The possibility of the external field polarizing the neighboring molecule. Depending on the polarity of the field some bond may increase in length, other decrease. Some may become more polarized than others. Together, increasing polarity should lead to an increase in volume.

To explore these nuances, we further computed the same ionic liquid components as ion pairs in an equilibrium geometry, Table 1 (entry 2). There, the vdW, Boyd, and Bader methods show a dramatic decrease in volume due to intermolecular overlap between the species. However, STREUSEL appears to recover similar volumes for the molecules computed separately or together (<4% change). Unlike those data presented in Fig. 5, the lack of change is likely because the atoms carrying the formal charge are heavily buried within the molecule, thereby being shielded from external fields. We can imagine a multitude of follow-up studies that isolate the screening effects of bound substituents in diverse chemical systems.

Undoubtedly, the utility of this approach may prove to be more general than a single example we imagined for this paper. This finding may prove to be a useful design principle for assessing volumes, densities, and other 3D chemical properties a priori. These data may feed into computer accelerated materials discovery[71-73] by rapidly screening molecular libraries for ideal multimolecular aggregates, forgoing arduous sampling of intermolecular interacting geometries.

Conclusions

The calculation of atomic, molecular, and material size has always been defined by an arbitrary cutoff where one chemical system ends and another begins. In prior approaches, the determination of atomic radii, volume, and associated surface area, has depended on electron density (either experimental or from simulation). Yet the limitations of these models have been highlighted by comparison of highly polarized systems, with the most striking being a difference in affected volume of free cations. Indeed, the contrast between the work presented herein, and the van der Waals metric is the description of the atomic surface being defined by a field smaller than the thermal energy provided at near absolute 0 K. The value of the latter is multifold: it allows for an assessment of how large an area is affected by the presence of the ion (an effect that may prove to be significant in applied electrochemistry,[74] and other complex systems[75,76]). It also allows for direct estimates of dispersion interactions (by redefining the cutoff to account for temperature effects such that materials can interact via their fields more strongly than kBT).

Beyond atomic cations, there are clear pathways to future fundamental studies examining how electric fields interact with, and are shielded by, chemical moieties (e.g. examining the field dissipation in a series of related cations, (CH3)4N+, (C2H5)4N+etc.). We can imagine studies that broaden our understanding of how fields permeate both through space and through bond, elucidating limitations of electric field-based descriptions of atomic, molecular, and material effective size. While our method undoubtedly does not replace the vdW radius or other size definitions, it does provide an alternative framework to contemplate the boundaries of chemistry.

Data availability

STREUSEL is available on github (see ESI†). Example input and output files are also hosted on the same site.

Author contributions

The computations were performed by AMM, AMD, JS and JD. The project was conceived by AMM and CHH. All authors partook in the drafting of the manuscript.

Conflicts of interest

There are no conflicts to declare.