Search tips
Search criteria 


Logo of wtpaEurope PMCEurope PMC Funders GroupSubmit a Manuscript
J Phys Chem Lett. Author manuscript; available in PMC 2017 March 21.
Published in final edited form as:
PMCID: PMC5360233

Chemisorption of Hydroxide on 2D Materials From DFT Calculations: Graphene Versus Hexagonal Boron Nitride


Recent nanofluidic measurements revealed strongly different surface charge measurements for boron-nitride and graphitic nanotubes when in contact with saline and alkaline water. 1,2 These observations contrast with the similar reactivity of a graphene layer and its boron nitride counterpart, using Density Functional Theory (DFT) framework, for intact and dissociative adsorption of gaseous water molecules. Here, we investigate, by DFT in implicit water, single and multiple adsorption of anionic hydroxide on single layers. A differential adsorption strength is found in vacuum for the first ionic adsorption on the two materials – chemisorbed on BN while physisorbed on graphene. The effect of implicit solvation reduces all adsorption values resulting in a favorable (non-favorable) adsorption on BN (graphene). We also calculate a pKa [similar, equals] 6 for BN in water, in good agreement with experiments. Comparatively, the unfavorable results for graphene in water echoes the weaker surface charge measurements, but points to an alternative scenario.

The successful isolation of single-layer graphene sheet3 has led to tremendous progress in the discovery of new 2D materials including boron nitride, silicene and transition metal dichalcogenides to cite a few. In particular hexagonal boron nitride monolayers share the planar honeycomb structure of graphene, with boron and nitrogen atoms alternating in the vortices of the honeycomb structure. Despite of their similar crystallographic structure, the electronic structure of BN and graphene nano sheets are drastically different : semimetallic for graphene and insulating for BN. Until recently both pristine layered materials were assumed to present inherently low chemical reactivity. However recent experimental works showed that covalent chemical fonctionalization could be achieved under drastic reactive conditions, either via reduction4,5 or reaction with oxidative reagents (H2O2,6 oxygen radical7) or via fluorination8 Furthermore recent nanofluidic experiments have altered the long-standing picture of the chemical inertness of BN nanotubes by showing a large negative surface charge when the BN materials is in contact with aqueous saline solutions. In contrast, similar experimental investigations based on nanofluidic investigations showed that carbon nanotubes2 and planar graphitic surfaces9 exhibit smaller and even minute surface charge.

The magnitude of the maximum charging - ranging from 0.1 C / m2 for C to 1 C/m2 for BN depending on the solution pH - rules out point defects as the sole origin of the measured surface charging. These experimental results therefore call for two challenging questions: (i) what are the potential chemical species present in basic and saline solutions that could adhere on the surface and strongly charge the pristine nanotubes ? (ii) how can we explain the different charging behavior of graphene and BN ideal nanotubes ?

A recent DFT study has investigated the dissociation barrier of a water molecule in contact with pristine graphene and BN nano sheets and concluded that the process is highly endothermic with a larger barrier on graphene than on BN.10 They also proposed facile ways to enhance the dissociation, either via including substitutional defects or either via pre-adsorption of H, which again amounts to the presence of defects.

Here in this DFT work we consider pristine nano sheets and take into account the presence of hydroxide ions in water. We demonstrate that adding one electron (via the charged adsorbate) to these 2D materials permits to tune the adsorption strength of one hydroxide : a strong ionic adsorption on BN and a weak ionic adsorption on graphene. Our atomistic results corroborate the acidity (pKa) of the BN sheet, pointing to the hydroxide anion as the best candidate to explain the charging of BN material immersed in a saline electrolyte.

The total energies of the systems have been calculated using the periodic DFT code Vienna Ab initio Simulation Package VASP (see Methods). The graphene and h-BN single layers were modeled using the same supercell corresponding to an (3×5) orthorhombic unit cell containing 60 atoms in Periodic Boundary Conditions (PBC) (see figure 1). The slab monolayer is separated by a vacuum space of 15 Å to avoid the interaction between the periodic images. The gamma point has been used in the k-point sampling of the first Brillouin zone of the supercell. The graphene and h-BN bond distances are evaluated to 1.43 Å and 1.45 Å with the PBE functional in agreement with previous calculations with the same method.11 We study the adsorption of anionic hydroxide molecule onto one B (or C) atom of the sp2 BN (graphene) layer respectively. The resulting adsorption state is therefore negatively charged and an homogeneous background charge is automatically added in the periodic supercell to ensure its neutrality. As a result the total energy of a charged solid system includes an additional electrostatic interaction with the charged background. To compute the first adsorption energy of OH called Eads1, we consider the adsorbed OH on a negatively charged surface (by adding an extra electron to the total valence electrons) as the final state and the non-interacting species taken in the same supercell as the initial state as shown in Figure 2A-B. By doing so, the additional electrostatic interaction with the charged vacuum cancels out in the resulting adsorption energy.


where the indices a and f stand for adsorbed and desorbed OH respectively (see figure caption 2). The second adsorption energy proceeds similarly as illustrated in Figure 2 C-D.

Figure 1
(A) BN hexagonal primitive unit cell. (B) BN Orthorhombic primitive unit cell. (C) BN orthorhombic supercell adopted in the DFT calculations. Distances, expressed in Å, are written in red (black) for the case of BN (graphene). Boron, Carbon, Nitrogen ...
Figure 2
(top) Simulation cells for the BN monolayer (same on graphene not shown) used to determine the energy of the adsorption of the first anionic hydroxide. (A) Sub + OHf : substrate and one desorbed anionic hydroxide in gaseous or aqueous phase. (B) OHa: ...

Table 1 summarizes the energetics computed within periodic DFT for the first OH adsorption on pristine BN and graphene monolayers subjected to different environments. In vacuum the adsorption is strongly favorable with - 2.3 eV on BN layer, while it is less favored on the graphene layer (-0.5 eV). The obtained values are strikingly different for both materials. The difference can be assigned to the charged ligand. Comparatively the neutral adsorption computed with spin-polarized calculations leads to similar moderate values (around -0.6 eV) in full agreement with previous studies.10,12 Interestingly the adsorption values do not depend on the inclusion, or not, of van der waals (vdw) corrections (see Eads1(vdw) in Table 1 ). Indeed in the final adsorption state the species are short-range bonded precluding vdw interactions and in the initial desorbed state the species are positioned to avoid any spurious vdw interaction.

Table 1
Adsorption energies (in eV) of the first OH species onto BN and Graphene layers in vacuum and in implicit water. Negative (positive) values mean favorable (unfavorable) cases. For the adsorption in vacuum, values in parentheses are spin-polarized ...

In contrast the adsorption values in implicit water are greatly reduced by around 1 eV maintaining a favorable adsorption value for BN while turning to an unfavorable adsorption value for graphene. Notably including the solvent in the DFT modelling permits to discriminate both materials regarding the adsorption of one anionic hydroxyl : BN will spontaneously attach, while graphene should remain inert at least at low temperatures. We note in passing an alternative way to model the adsorption of a charged ligand with zero net charge in the supercell : the insertion of a counterion such as one potassium cation in the vacuum space at long distance from the adsorbate and the surface. The resulting energetics including the potassium cation are very similar to the charged procedure used in this study (see Supporting Information). Hence the discriminative energetics between the two materials are robust against different modeling schemes of the charging effect.

The detailed DFT geometries of the adsorbed structures are displayed in Figure 3. Side views permit to visualize the sp3 rehybridization of the bonding atom inducing an out-of-plane buckling somewhat larger for BN than for graphene. The bonding distances between the O atom and the anchoring site on the substrates are 1.5 Å for BN and 1.54 Å for graphene, typical for single C-O bonds. A pyramidalization angle θp 13 can be determined and amounts to 18° for BN and 15° for graphene corresponding to 93 % and 77 % of sp3 hybridization respectively. The geometrical differences are not as striking as the energy differences but it confirms a weaker adsorption of OH on graphene than on BN.

Figure 3
DFT adsorption structures. Top views of BN (A) and graphene (B). Insets with side views of the bonding site for each material with key distances and angles expressed in Å and degrees, respectively.

We have investigated the spatial distribution of the additional electron using the Mulliken charge analysis (see Methods). It strongly differs between BN and graphene nano sheets (see supplementary Table). The grafted C atom on graphene gains only 0.06 electron and the rest of C atoms remains practically non-charged. On graphene the additional charge is totally delocalized over the nano sheet and the possible charge-activating role to bind hydroxyl is therefore masked. In contrast, on BN the grafted B gains 0.22 electron and becomes less positively charged than the rest of electron-deficient B atoms. The additional charge on BN is therefore substantially localized on the ligand and the grafted atom and the local reduction of the binding B atom can explain its improved reactivity.

We then studied the adsorption of a second OH ion onto BN and graphene monolayers. Figure 4 displays the second adsorption energy in vacuum Eads2 (dashed lines) and in implicit solvant Eads2sol (solid lines) as a function of the distance L from the first adsorption site. Circles mark the successive next nearest neighbors. Because only every two atoms are available for OH adsorption on h-BN, there are twice more of such sites on graphene as compared to BN.

Figure 4
A) Circles display possible sites (labelled with letters) for the second adsorption of OH at a given distance from the first adsorption site labelled A. B) Same for graphene. C) Adsorption energy of a second hydroxide versus the distance between ...

As shown on Fig. 4, the addition of a second OH ion on the nearest neighbor located in the same hexagon destabilizes the adsorption energy by typically ~1 eV as compared to the single ion adsorption. It remains favorable for BN (point B), but becomes readily unfavorable for graphene (point A).

The second adsorption energy is then found to decrease as a function of the distance L from the first ion. The value reaches a plateau in energy for large values of L, slightly above the first adsorption value. Typically the plateau is reached roughly beyond L = 3 Å for graphene and beyond L = 5 Å for BN. Now considering the effect of the solvant, one may observe that it leads to a global destabilizing shift for the second adsorption energy, typically ~ 0.3 − 0.5 eV. The difference in energy between the first adsorption and the plateau is however reduced as compared to vacuum. One may attribute this effect to the screening of charge.

In the following we compute the thermodynamics of adsorption of OH on the BN or graphene layers following successful approaches adapted to solid interfaces1416 (see Supporting Information). The correction terms that link Eads1sol to ΔGads were evaluated at T = 293.15 K and considering a solute standard concentration c0 = 1 mol L−1. The thermodynamic corrections are summarized in Table 2 while more details are available in Supporting Information. As a result we find a favorable ΔGads1 = - 0.46 eV for BN and a non-favorable ΔGads1 = + 0.79 eV for graphene. Thermodynamic corrections also worsen the scenario for the graphene material. This result rules out this charging mechanism on bare graphene but permits to envision it for BN.

Table 2
Thermodynamic corrections and the Gibbs free energy of the first adsorption evaluated in eV at 297,15 K and 1 bar. ϵb designates the dielectric permittivity of the solvent, G(OHa) is the Gibbs free energy of the system with OH adsorbed ...

We further quantified how the predicted adsorption energies ΔGads1 compare with the surface charge data extracted from the BN measurements in.1 From the value of ΔGads1, one can define the surface pKa as


As an overall result, Eq. (15) predicts pKa = 6 for BN and pKa = 27.7 for graphene. The value for BN agrees well with the post-treatment of experimental data,1 in which the electrostatic potential at the surface was evaluated using the Poisson-Boltzmann model. The value proposed for graphene is non-realistic precluding the same mechanism on graphene.

Our results support the strong reactivity of BN monolayer in the presence of hydroxyl anions contained in saline water. The contrasted result on graphene confirms the differential chemical behavior between the two materials. This echoes the experimental results using nanofluidic tools where much smaller adsorption is found on carbon nanotubes and planar graphitic surface.2,9 Still, the presence of a surface charge on carbon tubes, may suggest that other species than OH could adsorb onto the graphene surface and charge it.

Another interesting parameter is the maximum coverage for the ions on the surface. Such a quantity enters the global thermodynamics when considering for example the charge-pH relationship in the experiments. The adsorption equilibrium


is usually modeled by a Langmuir adsorption model, which introduces a maximal fraction θmax < 1 of B or C atoms (for h-BN or graphene respectively) that constitute possible adsorption sites. This value for θmax fixes the maximal charge expected on the surface, which is an interesting quantity per se. Denoting θ the fraction of bonded B or C atoms to an adsorbed OH, the chemical potential of OH at the surface is then modelled as


where ϕ is the electrostatic potential at the surface. The chemical potential of OH− in the bulk at the concentration [OH] expressed in mol.L−1 is fixed at equilibrium by its bulk value.


so that


The equilibrium condition then becomes


which relates the surface coverage to pH of the solution and the electrostatic potential at the surface.

A rough estimate for θmax can be done for BN on the basis of our DFT calculations. To this aim, we performed preliminary investigations of further adsorption of hydroxyl anions keeping the same minimum favorable distance revealed in Fig. 4. First results are presented on Fig. 5. Considering a third adsorption on the same hexagon leads to an energetically unfavorable configuration (plot B). The adsorption of three adsorbates in a raw is energetically favorable, although the zig zag configuration (plot C) is less favorable than the straight segment (plot D). Finally completing the lines with a fourth adsorbate becomes non-favorable (plots E and F). Of course these are first configurations and the stability of many other configurations should be evaluated. This goes beyond the present work and will be the object of a forthcoming study. However, on the basis of these results, one may estimate that 3 hydroxides can be grafted at most in a rectangle of (L+l3OHL, which amounts to a surface of 50 Å2. This corresponds to a site density of 0.06 per Å2, and the associated full coverage site density is ~ 0.18 per Å2(30 B atoms per supercell). This leads to θmax [similar, equals] 0.33 and Σmax [similar, equals]0.96 Cm−2. These values are in good agreement with the values reported for the maximum charge mesured on BN surfaces.1

Figure 5
Beyond the second adsorption on BN monolayer : from two adsorbates in the same honeycomb at the minimum favorable distance (A) exploration of three adsorbates (B-C-D) and four adsorbates (E-D) keeping the same favorable distance. The indicated energy ...

In conclusion, we have investigated comparatively within ab initio DFT framework a monolayer of BN and graphene embedded in implicit water and in contact with one hydroxyl anion. We predict for the first time a contrasted chemical reactivity for graphene and BN : a chemisorption of OH on BN while this mechanism is not highlighted on graphene. This striking difference between the two pristine materials is linked to their electronic structure and the lack of screening for graphene surface as compared to the insulating BN. Furthermore, on the basis of atomistic thermodynamics including all the vibrational contributions to energy and entropy, we have derived a free energy of adsorption of - 0.46 eV on BN which corresponds to a pKa [similar, equals] 6, in quantitative agreement with experimental measurements.1 Hence this charging scenario is viable on BN and our study leans strong support to the recently-revealed non-inertness of the BN material in contact with aqueous saline phase. On the other hand, the case of graphene requires further investigations. Indeed planar graphitic surfaces show minute surface charge (as deduced from surface conductance measurements)9 while the experiments on carbon nanotubes still exhibit a surface charge,2 although much weaker as compared to BN tubes.1 This suggests the possibility of adsorption of other charged species on the carbon surfaces, not considered in the present study. Future work will involve the exploration at finite T of the charging mechanism using explicit water solvent and the inclusion of counterions using Ab Initio Molecular Dynamics AIMD, in order to evaluate a more realistic pKa value. Also taking the full stacking of these 2D materials into account, their curvature effect (since the measurements are made on multi-wall nanotubes of different lengths and diameters) are other ingredients which will be of high interest in understanding the molecular mechanisms at stake in the chemical reactivity of these 2D materials.


Ab initio calculations

In the VASP code,1719 the wave functions of the system contained in the supercell have been expanded in a plane wave basis set with an energy cutoff of 800 eV. The common cutoff for hybrid interfaces is 400 eV in vacuum and taking twice the value is requested for sake of precision when dealing with the implicit solvent scheme (see below).The electronic cores are described by the projector augmented wave method.20 The exchange-correlation energies have been evaluated within the generalized gradient approximation (GGA).We employ GGA functionals as suggested by Perdew, Burke, and Ernzerhof (PBE).21,22 However we have also verified that the key results (adsorption energies) obtained here are not particularly sensitive to the choice of exchange-correlation functional, in particular using the vdW-inclusive optB86b-vdW functional23,24 which include vdw interactions. We use the smearing technique using a Gaussian broadening of 0.2 eV to achieve electronic convergence. To check the possible influence of the calculation settings, we have modified sequentially the smearing value to 0.1 eV, the k-mesh to 3x3x1 and the vacuum space to 20 Å. The adsorption values varie between 0.1 to 0.2 eV but the striking difference between graphene and BN adsorption energies holds true.

In a second step the solvation energies in implicit water (Esol) are evaluated within the joint density functional theory framework as implemented into VASP by Mathew and Hennig25 and successfully used recently.2629 In this technique the dielectric permittivity of the medium is defined as a functional of the electronic density. Indeed dielectric permittivity approaches asymptotically the bulk water value at 1 bar, ϵb = 78.4, in the region where the electronic density is lower than ρcut = 0.0025 Å3 (the default cutoff charge density), while ϵb =1 in the region when the electronic density is high (inside the slab treated as a solute). The cavitation energies which account for the solvant contributions are calculated with a surface tension parameter of 0.525 meV/Å2. The systems are fully optimized in vacuum with atomic forces lower than 0.02 eV per Å. Then the optimized vacuo geometries are subjected to the implicit solvant interaction and the total energies named Esol are corrected after an energy and geometry optimization. For the second adsorption Eads2 however we use the frozen geometry from vacuum for sake of computational time. As a result the adsorption energy in implicit water transforms into :


To evaluate the vibrational entropic corrections to this adsorption energy, we have run a vibrational analysis of state A and B (see figure 2) in order to determine the vibrational modes and frequencies. The classical procedure in VASP is to calculate and diagonalize the Hessian matrix using the finite-difference method.

In addition we derived Mulliken atomic charges from a static Dmol3 calculation3032 based of the VASP optimized structures in order to understand the spatial location of the excess charge.

Supplementary Material

Supporting Information Available

Additional figures and information.

This material is available free of charge via the Internet at

Supporting Information


B.G., C.P. and M.-L.B. thank the national computational center IDRIS, CINES and TGCC (Grant 2015- [x2015087364]) for the computer grant. L.B. and A.S. thank ANR Blue Energy for support. A.S. and B.G. acknowledge funding from the European Union’s H2020 Framework Programme / ERC Starting Grant agreement number 637748 - NanoSOFT.


(1) Siria A, Poncharal P, Biance A-L, Fulcrand R, Blase X, Purcell ST, Bocquet L. Giant Osmotic Energy Conversion Measured in a Single Transmembrane Boron Nitride Nanotube. Nature. 2013;494:455–458. [PubMed]
(2) Secchi E, Niguès A, Jubin L, Siria A, Bocquet L. Scaling Behavior for Ionic Transport and its Fluctuations in Individual Carbon Nanotubes. Phys Rev Lett. 2016;116:154501. [PMC free article] [PubMed]
(3) Novoselov KS, Geim AK, Morozov SV, Jiang D, Zhang Y, Dubonos SV, Grigorieva IV, Firsov AA. Electric Field Effect in Atomically Thin Carbon Films. Science. 2004;306:666–669. [PubMed]
(4) Sun Q, Li Z, Searles DJ, Chen Y, Lu GM, Du A. Charge-Controlled Switchable CO2 Capture on Boron Nitride Nanomaterials. J Am Chem Soc. 2013;135:8246–8253. [PubMed]
(5) Shin H, Guan J, Zgierski MZ, Kim KS, Kingston CT, Simard B. Covalent Functionalization of Boron Nitride Nanotubes via Reduction Chemistry. ACS Nano. 2015;9:12573–12582. [PubMed]
(6) Zhi C, Bando Y, Terao T, Tang C, Kuwahara H, Golberg D. Chemically Activated Boron Nitride Nanotubes. Chem Asian J. 2009;4:1536–1540. [PubMed]
(7) Sainsbury T, Satti A, May P, Wang Z, McGovern I, Gunko YK, Coleman J. Oxygen Radical Functionalization of Boron Nitride Nanosheets. J Am Chem Soc. 2012;134:18758–18771. [PubMed]
(8) Tang C, Bando Y, Huang Y, Yue S, Gu C, Xu F, Golberg D. Fluorination and Electrical Conductivity of BN Nanotubes. J Am Chem Soc. 2005;127:6552–6553. [PubMed]
(9) Radha B, Esfandiar A, Wang FC, Rooney AP, Gopinadhan K, Keerthi A, Mishchenko A, Janardanan A, Blake P, Fumagalli L, et al. Molecular Transport Through Capillaries Made With Atomic-Scale Precision. Nature. 2016;538:222–225. [PubMed]
(10) Al-Hamdani YS, Alfè D, von Lilienfeld OA, Michaelides A. Tuning Dissociation Using Isoelectronically Doped Graphene and Hexagonal Boron Nitride: Water and Other Small Molecules. J Chem Phys. 2016;144:154706. [PubMed]
(11) Graziano G, Klimeš J, Fernandez-Alonso F, Michaelides A. Improved Description of Soft Layered Materials with Van der Waals Density Functional Theory. J Phys: Condens Matter. 2012;24:424216. [PubMed]
(12) Ghaderi N, Peressi M. First-Principle Study of Hydroxyl Functional Groups on Pristine, Defected Graphene, and Graphene Epoxide. J Phys Chem C. 2010;114:21625–21630.
(13) Niyogi S, Hamon MA, Hu H, Zhao B, Bhowmik P, Sen R, Itkis ME, Haddon RC. Chemistry of Single-Walled Carbon Nanotubes. Acc Chem Res. 2002;35:1105–1113. [PubMed]
(14) Reuter K, Scheffler M. Composition and Structure of the RuO2(110) Surface in an O2 and CO Environment: Implications for the Catalytic Formation of CO2. Phys Rev B. 2003;68:045407.
(15) Loffreda D. Theoretical Insight of Adsorption Thermodynamics of Multifunctional Molecules on Metal Surfaces. Surf Sci. 2006;600:2103–2112.
(16) Hill T. An Introduction to Statistical Thermodynamics. Addison-Wesley series in chemistry; Dover Publications; 1960.
(17) Kresse G, Hafner J. Ab Initio Molecular Dynamics for Liquid Metals. Phys Rev B. 1993;47:558–561. [PubMed]
(18) Kresse G, Furthmüller J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput Mater Sci. 1996;6:15–50.
(19) Kresse G, Furthmüller J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys Rev B. 1996;54:11169–11186. [PubMed]
(20) Kresse G. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys Rev B. 1999;59:1758–1775.
(21) Perdew JP, Burke K, Ernzerhof M. Generalized Gradient Approximation Made Simple. Phys Rev Lett. 1996;77:3865–3868. [PubMed]
(22) Ziesche P, Kurth S, Perdew JP. Density Functionals from LDA to GGA. Comput Mater Sci. 1998;11:122–127.
(23) Klimeš J, Bowler DR, Michaelides A. Van der Waals Density Functionals Applied to Solids. Phys Rev B. 2011;83:195131.
(24) Klimeš J, Michaelides A. Perspective: Advances and Challenges in Treating Van der Waals Dispersion Forces in Density Functional Theory. J Chem Phys. 2012;137:120901. [PubMed]
(25) Mathew K, Sundararaman R, Letchworth-Weaver K, Arias TA, Hennig RG. Implicit Solvation Model for Density-Functional Study of Nanocrystal Surfaces and Reaction Pathways. J Chem Phys. 2014;140:084106. [PubMed]
(26) Lespes N, Filhol J-S. Using Implicit Solvent in Ab Initio Electrochemical Modeling: Investigating Li+/Li Electrochemistry at a Li/Solvent Interface. J Chem Theory Comput. 2015;11:3375–3382. [PubMed]
(27) Sakong S, Naderian M, Mathew K, Hennig RG, Gro A. Density Functional Theory Study of the Electrochemical Interface Between a Pt Electrode and an Aqueous Electrolyte using an Implicit Solvent Method. J Chem Phys. 2015;142:234107. [PubMed]
(28) Steinmann SN, Michel C, Schwiedernoch R, Filhol J-S, Sautet P. Modeling the HCOOH/CO2 Electrocatalytic Reaction: When Details Are Key. ChemPhysChem. 2015;16:2307–2311. [PubMed]
(29) Garcia-Ratés M, López N. Multigrid-Based Methodology for Implicit Solvation Models in Periodic DFT. J Chem Theory Comput. 2016;12:1331–1341. [PubMed]
(30) Delley B. An All-Electron Numerical Method for Solving the Local Density Functional for Polyatomic Molecules. J Chem Phys. 1990;92:508.
(31) Delley B. From Molecules to Solids with the DMol3 Approach. J Chem Phys. 2000;113:7756.
(32) Delley B. Fast Calculation of Electrostatics in Crystals and Large Molecules. J Phys Chem. 1996;100:6107–6110.