Linear-scaling Density Functional Theory Research Articles

Strain-dependent one-dimensional confinement channels in twisted bilayer 1T′−WTe2

The low-symmetry and anistropic lattice of 1T′WTe2 is responsible for the existence of parallel one-dimensional channels in the moiré patterns of twisted bilayers. This gives the opportunity to explore moiré physics of a different nature to that widely observed in twisted bilayers of materials with hexagonal symmetries. Here, we combine plane-wave and linear-scaling density functional theory calculations to describe the electronic properties of twisted bilayer 1T′WTe2. We investigate the strain dependence of both aligned bilayer configurations across a range of stacking geometries, and of the twisted bilayer itself. For a small change in the lattice parameters of the constituent 1T′WTe2 monolayers, we find a very substantial increase in the moiré-induced striped electrostatic potential landscape in the twisted bilayer. We measure a peak-to-trough magnitude >200 meV between parallel linear channels, sufficient to induce Luttinger liquid behavior. Published by the American Physical Society 2024

Noise reduction of stochastic density functional theory for metals.

Density Functional Theory (DFT) has become a cornerstone in the modeling of metals. However, accurately simulating metals, particularly under extreme conditions, presents two significant challenges. First, simulating complex metallic systems at low electron temperatures is difficult due to their highly delocalized density matrix. Second, modeling metallic warm-dense materials at very high electron temperatures is challenging because it requires the computation of a large number of partially occupied orbitals. This study demonstrates that both challenges can be effectively addressed using the latest advances in linear-scaling stochastic DFT methodologies. Despite the inherent introduction of noise into all computed properties by stochastic DFT, this research evaluates the efficacy of various noise reduction techniques under different thermal conditions. Our observations indicate that the effectiveness of noise reduction strategies varies significantly with the electron temperature. Furthermore, we provide evidence that the computational cost of stochastic DFT methods scales linearly with system size for metal systems, regardless of the electron temperature regime.

Linear-scaling density functional theory (DFT) simulations of point, Frenkel and Schottky defects in CeO2

CeO2 (ceria) is a material of significant industrial and technological importance, used in solid oxide fuel cells and catalysis. Here, we explore the usage of linear-scaling density functional theory as implemented in the ONETEP code, which allows to use larger simulation cells. By using DFT+U calculations we revise the defect chemistry of ceria, including point defects, Frenkel and Schottky defects.We found that the ground state of an oxygen vacancy is associated to two neighbouring reduced cerium sites. A cerium vacancy is the least favourable point defect, where holes localise on neighbouring oxygen sites. It is more favourable to displace an oxygen interstitial defect away from the octahedral interstitial site, with the formation of a stable peroxide species. Our simulations show that a cerium interstitial is best accommodated in the octahedral interstitial site, as this minimises the distortion of the lattice.Placing a vacancy and an interstitial defect at a separation of 5.18 Å for the OF¡110¿ and 4.77 Å for the CeF¡111¿, stable Frenkel defects can be formed. We also studied the effect of different supercell size on the energetic ordering of Schottky defects, where the S¡111¿ is more favourable than the S¡110¿ for a given simulation cells containing 324 or more atoms.

A variational formulation of the Harris functional as a correction to approximate Kohn-Sham density functional theory.

Accurate descriptions of intermolecular interactions are of great importance in simulations of molecular liquids. We present an electronic structure method that combines the accuracy of the Harris functional approach with the computational efficiency of approximately linear-scaling density functional theory (DFT). The non-variational nature of the Harris functional has been addressed by constructing a Lagrangian energy functional, which restores the variational condition by imposing stationarity with respect to the reference density. The associated linear response equations may be treated with linear-scaling efficiency in an atomic orbital based scheme. Key ingredients to describe the structural and dynamical properties of molecular systems are the forces acting on the atoms and the stress tensor. These first-order derivatives of the Harris Lagrangian have been derived and implemented in consistence with the energy correction. The proposed method allows for simulations with accuracies close to the Kohn-Sham DFT reference. Embedded in the CP2K program package, the method is designed to enable ab initio molecular dynamics simulations of molecular solutions for system sizes of several thousand atoms. Available subsystem DFT methods may be used to provide the reference density required for the energy correction at near linear-scaling efficiency. As an example of production applications, we applied the method to molecular dynamics simulations of the binary mixtures cyclohexane-methanol and toluene-methanol, performed within the isobaric-isothermal ensemble, to investigate the hydrogen bonding network in these non-ideal mixtures.

Development of Nanoscale Graphene Oxide Models for the Adsorption of Biological Molecules.

Graphene oxide (GO), a nanomaterial with promising applications that range from water purification to enzyme immobilization, is actively present in scientific research since its discovery. GO studies with computational methodologies such as molecular dynamics are frequently reported in the literature; however, the models used often rely on approximations, such as randomly placing functional groups and the use of generalized force fields. Therefore, it is important to develop new MD models that provide a more accurate description of GO structures and their interaction with an aqueous solvent and other adsorbate molecules. In this paper, we derived new force field non-bonded parameters from linear-scaling density functional theory calculations of nanoscale GO sheets with more than 10,000 atoms through an atoms-in-molecules (AIM) partitioning scheme. The resulting GAFF2-AIM force field, derived from the bonded terms of GAFF2 parameterization, reproduces the solvent structure reported in ab initio MD simulations better than the force field nowadays widely used in the literature. Additionally, we analyzed the effect of the ionic strength of the medium and of the C/O ratio on the distribution of charges surrounding the GO sheets. Finally, we simulated the adsorption of natural amino acid molecules to a GO sheet and estimated their free energy of binding, which compared very favorably to their respective experimental values, validating the force field presented in this work.

First-Principles Study on the Effect of Lithiation in Spinel LixMn2O4 (0 ≤ x ≤ 1) Structure: Calibration of CASTEP and ONETEP Simulation Codes

Lithium–manganese–oxide (Li-Mn-O) spinel is among the promising and economically viable, high-energy density cathode materials for enhancing the performance of lithium-ion batteries. However, its commercialization is hindered by its poor cyclic performance. In computational modelling, pivotal in-depth understanding of material behaviour and properties is sizably propelled by advancements in computational methods. Hence, the current work compares traditional DFT (CASTEP) and linear-scaling DFT (ONETEP) in a LiMn2O4 electronic property study to pave way for large-scale DFT calculations in a quest to improve its electrochemical properties. The metallic behaviour of LixMn2O4 (0.25 ≤ x ≤ 1) and Li2Mn2O4 was correctly determined by both CASTEP and ONETEP code in line with experiments. Furthermore, OCV during the discharge cycle deduced by both codes is in good accordance and is between 5 V and 2.5 V in the composition range of 0 ≤ x ≤ 1. Moreover, the scaling of the ONETEP code was performed at South Africa’s CHPC to provide guidelines on more productive large-scale ONETEP runs. Substantial total computing time can be saved by systematically adding the number of processors with the growing structure size. The study also substantiates that true linear scaling of the ONETEP code is achieved by a systematic truncation of the density kernel.

Linear Weak Scalability of Density Functional Theory Calculations without Imposing Electron Localization.

Linear scaling density functional theory (DFT) approaches to the electronic structure of materials are often based on the tendency of electrons to localize in large atomic and molecular systems. However, in many cases of actual interest, such as semiconductor nanocrystals, system sizes can reach a substantial extension before significant electron localization sets in, causing a considerable deviation from linear scaling. Herein, we address this class of systems by developing a massively parallel DFT approach which does not rely on electron localization and is formally quadratic scaling yet enables highly efficient linear wall-time complexity in the weak scalability regime. The method extends from the stochastic DFT approach described in Fabian et al. (WIRES: Comp. Mol. Sci.2019, e1412) but is entirely deterministic. It uses standard quantum chemical atom-centered Gaussian basis sets to represent the electronic wave functions combined with Cartesian real-space grids for some operators and enables a fast solver for the Poisson equation. Our main conclusion is that when a processor-abundant high-performance computing (HPC) infrastructure is available, this type of approach has the potential to allow the study of large systems in regimes where quantum confinement or electron delocalization prevents linear scaling.

Efficiently Computing Excitations of Complex Systems: Linear-Scaling Time-Dependent Embedded Mean-Field Theory in Implicit Solvent.

Quantum embedding schemes have the potential to significantly reduce the computational cost of first-principles calculations while maintaining accuracy, particularly for calculations of electronic excitations in complex systems. In this work, I combine time-dependent embedded mean field theory (TD-EMFT) with linear-scaling density functional theory and implicit solvation models, extending previous work within the ONETEP code. This provides a way to perform multilevel calculations of electronic excitations on very large systems, where long-range environmental effects, both quantum and classical in nature, are important. I demonstrate the power of this method by performing simulations on a variety of systems, including a molecular dimer, a chromophore in solution, and a doped molecular crystal. This work paves the way for high accuracy calculations to be performed on large-scale systems that were previously beyond the reach of quantum embedding schemes.

Forces from Stochastic Density Functional Theory under Nonorthogonal Atom-Centered Basis Sets.

We develop a formalism for calculating forces on the nuclei within the linear-scaling stochastic density functional theory (sDFT) in a nonorthogonal atom-centered basis set representation (Fabian et al. Wiley Interdiscip. Rev.: Comput. Mol. Sci.2019, 9, e1412, 10.1002/wcms.1412) and apply it to the Tryptophan Zipper 2 (Trp-zip2) peptide solvated in water. We use an embedded-fragment approach to reduce the statistical errors (fluctuation and systematic bias), where the entire peptide is the main fragment and the remaining 425 water molecules are grouped into small fragments. We analyze the magnitude of the statistical errors in the forces and find that the systematic bias is of the order of 0.065 eV/Å (∼1.2 × 10–3Eh/a0) when 120 stochastic orbitals are used, independently of system size. This magnitude of bias is sufficiently small to ensure that the bond lengths estimated by stochastic DFT (within a Langevin molecular dynamics simulation) will deviate by less than 1% from those predicted by a deterministic calculation.

Ab initio study of lithium intercalation into a graphite nanoparticle

Ab initio workflow for prediction of Li intercalation, with minimal calculations, in anode-like graphite nanoparticles using linear-scaling DFT. The workflow is able to reproduce key experimental data including staging, charge transfer, and OCVs.

Unravelling the binding affinity and selectivity of molybdenum(II) phenanthroline complexes with DNA G-quadruplexes by using linear-scaling DFT studies. The important role of ancillary ligands.

We have used near linear-scaling density functional theory (LS-DFT) methods including dispersion, for the first time, to study the interaction of two isomers, equatorial (Eq) and axial (Ax), of the [Mo(η3-C3H5)Br(CO)2(phen)] metal complex with the DNA G-quadruplexes (GQ) to gain insight into its cytotoxicity. The LMKLL/DZDP level of calculation, which includes van der Waals contributions, with the SIESTA software was used to treat by means of first-principles computations the whole biological studied model system with ∼1000 atoms. Computed formation energies point to systems containing the Ax isomer as the most stable although the nearest system in energy containing the Eq isomer is only 7.5 kcal mol-1 above. On the other hand, the energy decomposition analysis (EDA) favours interaction energies for the systems containing the Eq isomer. However, when solvent effects are taken into account the systems containing the Ax isomer are again the most stable. This Ax isomer was found interacting by means of end-stacking with the GQ and surprisingly totally inside the non-canonical secondary structure, where all the ligands of the metal complex produce several weak interactions with the DNA structure. On the other hand, the Eq isomer prefers to interact from outside by means of intercalation in which the ancillary ligands also have some role in the interaction. Such features and comparison with the results regarding the interaction of the [Mo(η3-C3H5)Br(CO)2(phen)] metal complex with duplex DNA suggest that the [Mo(η3-C3H5)Br(CO)2(phen)] would have a higher affinity and eventual selectivity for non-canonical DNA GQ structures.

Semi-empirical and linear-scaling DFT methods to characterize duplex DNA and G-quadruplexes in the presence of interacting small molecules

SIESTA LS-DFT method at LMKLL/DZDP level and PM6-DH2 semi-empirical method are suitable to obtain correct geometries for DNA structures interacting with small molecules. However, only the SIESTA method reproduces the interaction energies when comparing to DLPNO-CCSD(T) and data bases.

Electrochemistry from first-principles in the grand canonical ensemble.

Progress in electrochemical technologies, such as automotive batteries, supercapacitors, and fuel cells, depends greatly on developing improved charged interfaces between electrodes and electrolytes. The rational development of such interfaces can benefit from the atomistic understanding of the materials involved by first-principles quantum mechanical simulations with Density Functional Theory (DFT). However, such simulations are typically performed on the electrode surface in the absence of its electrolyte environment and at constant charge. We have developed a new hybrid computational method combining DFT and the Poisson-Boltzmann equation (P-BE) capable of simulating experimental electrochemistry under potential control in the presence of a solvent and an electrolyte. The charged electrode is represented quantum-mechanically via linear-scaling DFT, which can model nanoscale systems with thousands of atoms and is neutralized by a counter electrolyte charge via the solution of a modified P-BE. Our approach works with the total free energy of the combined multiscale system in a grand canonical ensemble of electrons subject to a constant electrochemical potential. It is calibrated with respect to the reduction potential of common reference electrodes, such as the standard hydrogen electrode and the Li metal electrode, which is used as a reference electrode in Li-ion batteries. Our new method can be used to predict electrochemical properties under constant potential, and we demonstrate this in exemplar simulations of the differential capacitance of few-layer graphene electrodes and the charging of a graphene electrode coupled to a Li metal electrode at different voltages.

Intercalation voltages for spinel LixMn2O4 (0≤x≤2) cathode materials: Calibration of calculations with the ONETEP linear-scaling DFT code

Density functional theory (DFT) has become an engine for driving ab-initio quantum mechanical simulations spanning a vast range of applications. However, conventional DFT has limitations of the accessible system size due to computational expense. Recent progress on linear scaling DFT methods has enabled us to investigate larger systems. In this paper, we investigate the numerical agreement between conventional DFT codes and ONETEP, a linear-scaling DFT code, for important materials properties particularly calculating the intercalation potential of spinel electrode materials. Modulating materials with high energy density is an important aspect that contributes to the significant gap in our knowledge of the factors. We provide typical simulation results based on calculated intercalation potentials for discharging the spinel LixMn2O4, which plays a key role in developing high energy density lithium-ion batteries. The structural properties obtained after geometry optimisation with CASTEP yielded materials with volume within 3 % and lattice parameters within 1 % relative error with experimental values. The average intercalation potentials calculated with the CASTEP and ONETEP codes are within 3 % agreement with each other.

Atomic and electronic structure of two-dimensional Mo(1− x )W x S2 alloys

Alloying enables engineering of the electronic structure of semiconductors for optoelectronic applications. Due to their similar lattice parameters, the two-dimensional semiconducting transition metal dichalcogenides of the MoWSeS group (MX2 where M = Mo or W and X = S or Se) can be grown as high-quality materials with low defect concentrations. Here we investigate the atomic and electronic structure of Mo(1−x)W x S2 alloys using a combination of high-resolution experimental techniques and simulations. Analysis of the Mo and W atomic positions in these alloys, grown by chemical vapour transport, shows that they are randomly distributed, consistent with Monte Carlo simulations that use interaction energies determined from first-principles calculations. Electronic structure parameters are directly determined from angle resolved photoemission spectroscopy measurements. These show that the spin–orbit splitting at the valence band edge increases linearly with W content from MoS2 to WS2, in agreement with linear-scaling density functional theory predictions. The spin–orbit splitting at the conduction band edge is predicted to reduce to zero at intermediate compositions. Despite this, polarisation-resolved photoluminescence spectra on monolayer Mo0.5W0.5S2 show significant circular dichroism, indicating that spin-valley locking is retained. These results demonstrate that alloying is an important tool for controlling the electronic structure of MX2 for spintronic and valleytronic applications.

Discovery of SARS-CoV-2 Mpro peptide inhibitors from modelling substrate and ligand binding.

The main protease (Mpro) of SARS-CoV-2 is central to viral maturation and is a promising drug target, but little is known about structural aspects of how it binds to its 11 natural cleavage sites. We used biophysical and crystallographic data and an array of biomolecular simulation techniques, including automated docking, molecular dynamics (MD) and interactive MD in virtual reality, QM/MM, and linear-scaling DFT, to investigate the molecular features underlying recognition of the natural Mpro substrates. We extensively analysed the subsite interactions of modelled 11-residue cleavage site peptides, crystallographic ligands, and docked COVID Moonshot-designed covalent inhibitors. Our modelling studies reveal remarkable consistency in the hydrogen bonding patterns of the natural Mpro substrates, particularly on the N-terminal side of the scissile bond. They highlight the critical role of interactions beyond the immediate active site in recognition and catalysis, in particular plasticity at the S2 site. Building on our initial Mpro-substrate models, we used predictive saturation variation scanning (PreSaVS) to design peptides with improved affinity. Non-denaturing mass spectrometry and other biophysical analyses confirm these new and effective ‘peptibitors’ inhibit Mpro competitively. Our combined results provide new insights and highlight opportunities for the development of Mpro inhibitors as anti-COVID-19 drugs.

Parallel Implementation of Large-Scale Linear Scaling Density Functional Theory Calculations With Numerical Atomic Orbitals in HONPAS.

Linear-scaling density functional theory (DFT) is an efficient method to describe the electronic structures of molecules, semiconductors, and insulators to avoid the high cubic-scaling cost in conventional DFT calculations. Here, we present a parallel implementation of linear-scaling density matrix trace correcting (TC) purification algorithm to solve the Kohn–Sham (KS) equations with the numerical atomic orbitals in the HONPAS package. Such a linear-scaling density matrix purification algorithm is based on the Kohn's nearsightedness principle, resulting in a sparse Hamiltonian matrix with localized basis sets in the DFT calculations. Therefore, sparse matrix multiplication is the most time-consuming step in the density matrix purification algorithm for linear-scaling DFT calculations. We propose to use the MPI_Allgather function for parallel programming to deal with the sparse matrix multiplication within the compressed sparse row (CSR) format, which can scale up to hundreds of processing cores on modern heterogeneous supercomputers. We demonstrate the computational accuracy and efficiency of this parallel density matrix purification algorithm by performing large-scale DFT calculations on boron nitrogen nanotubes containing tens of thousands of atoms.

Probing the Operating Voltage and Structural Transitions in Li-Mn-O Layered-Spinel Composite Nanoarchitectures during the Discharge Process.

Li-Mn-O layered-spinel (LS) composite is one of the utmost promising cathode materials for advancing lithium-ion batteries, which are lacking in comparison to technological advancements. The currently used cathode materials are one of the limiting factors in realising the envisaged affordable, durable, and safe rechargeable lithium-ion batteries. Their limitations include capacity fading transpiring during cycling, inherent low operational capacities, safety issues due to thermal instabilities, and high cost which inhibit their full utilisation in large scale systems. The LS composite demonstrates enhanced electrochemical properties owing to the synergism of the layered and spinel components. In particular, the 0.7Li2MnO3.0.3Li4Mn5O12 composite, which has demonstrated a higher specific capacity, >250 mAh/g. Recent studies have focused mainly on enhancing the specific capacity of LS composites. However, there is limited knowledge of the structural changes that take place during the discharge process and how they affect the working voltage of the LS composite cathode material. Li-Mn-O LS composites indicating different stages of the discharge process have been generated using the molecular dynamics simulated amorphisation and recrystallisation method, which utilises empirical interatomic potentials. Different structural features of the LS composites in all significant stages of the discharge process were visualised and captured. X-ray diffraction patterns (XRDs) and radial distribution functions were used to characterise the formed layered and spinel components of the LS blend nanoparticles. Moreover, their average voltages were also computed using the linear scaling density functional theory code.

Large-scale first-principles quantum transport simulations using plane wave basis set on high performance computing platforms

As the characteristic lengths of advanced electronic devices are approaching the atomic scale, ab initio simulation method, with full consideration of quantum mechanical effects, becomes essential to study the quantum transport phenomenon in them. The widely used non-equilibrium Green’s function (NEGF) combined with the density functional theory (DFT) approach prefers a localized basis set. As many states of the art DFT calculations for solid state systems are carried out in plane waves, it is thus worth to investigate the feasibility of using the plane wave basis set for large-scale quantum transport calculations. Here we present a plane wave method for large-scale transport calculations based on a previously developed scattering state calculation approach (Wang, 2005). We address the unique computational challenges of applying that approach for large-scale systems where it is too expensive to calculate all the occupied eigenstates of the system as in conventional DFT calculations. By applying several high-efficiency parallel algorithms, including linear-scaling DFT algorithm, folded spectrum method, and Chebyshev filter technique, we demonstrate that it is possible to use this approach to simulate a system with several thousand atoms on high performance computing platforms. This method is not only used to study nanowire interconnects, showing how the shape and point defect affects their transport properties, but also used to study nanoscale Si transistor. Such quantum transport simulation method will be useful for investigating and designing nanoscale devices.

Electronic structure calculations in electrolyte solutions: Methods for neutralization of extended charged interfaces.

Density functional theory (DFT) is often used for simulating extended materials such as infinite crystals or surfaces, under periodic boundary conditions (PBCs). In such calculations, when the simulation cell has non-zero charge, electrical neutrality has to be imposed, and this is often done via a uniform background charge of opposite sign ("jellium"). This artificial neutralization does not occur in reality, where a different mechanism is followed as in the example of a charged electrode in electrolyte solution, where the surrounding electrolyte screens the local charge at the interface. The neutralizing effect of the surrounding electrolyte can be incorporated within a hybrid quantum-continuum model based on a modified Poisson-Boltzmann equation, where the concentrations of electrolyte ions are modified to achieve electroneutrality. Among the infinite possible ways of modifying the electrolyte charge, we propose here a physically optimal solution, which minimizes the deviation of concentrations of electrolyte ions from those in open boundary conditions (OBCs). This principle of correspondence of PBCs with OBCs leads to the correct concentration profiles of electrolyte ions, and electroneutrality within the simulation cell and in the bulk electrolyte is maintained simultaneously, as observed in experiments. This approach, which we call the Neutralization by Electrolyte Concentration Shift (NECS), is implemented in our electrolyte model in the Order-N Electronic Total Energy Package (ONETEP) linear-scaling DFT code, which makes use of a bespoke highly parallel Poisson-Boltzmann solver, DL_MG. We further propose another neutralization scheme ("accessible jellium"), which is a simplification of NECS. We demonstrate and compare the different neutralization schemes on several examples.