Predicting interaction properties of a many-body system

Abstract

A system and method for predicting electronic structure properties of an atomic system of arbitrary definition. The method includes receiving parameters of the atomic system and determining a total electronic free energy of the atomic system, by inter alia constructing a matrix of coefficients of plane-wave expansions of electronic orbitals corresponding to the atomic system. A plurality of energies are then calculated, including Hartree potential energy. The total electronic free energy and derivatives of the total electronic free energy with respect to the matrix of coefficients, are found by summing the energies. If the derivatives do not reflect a stable local minimum, then the above process is repeated until convergence (i.e., stability), by applying gradient descent to the matrix using the derivatives. If a stable local minimum is found, a prediction of interaction properties of the atomic system with other atomic systems is generated, by reference to the total electronic free energy, and derivative.

Description

PREDICTING INTERACTION PROPERTIES OF A MANY-BODY SYSTEMTECHNIC AL FIELD

[0001] The present invention provides a system and method for determining minimum free energy of an atomic or many-body system. The invention can be used in, but is not limited to, predicting interaction properties of such many-body systems, for use in development of materials, drugs, chemicals and other products BACKGROUND

[0002] The reference in this specification to any prior publication (or information derived from it), or to any matter which is known, is not, and should not be taken as an acknowledgment or admission or any form of suggestion that that prior publication (or information derived from it) or known matter forms part of the common general knowledge in the field of endeavour to which this specification relates.

[0003] Modem development processes have increasingly involved methods for reducing work, time to development and other factors. In addition, in racing to develop new products, it is critical to understand how those products will behave. For new materials, drugs and chemicals, it has become critical to avoid delays, liability and cost associated with investigating unsuitable materials, compounds without the intended functional effect or with undesirable functional effects, and chemicals with unknown or unexpected interactions.|0004| Electronic structure theory has emerged as a candidate for reducing wasted investigation, directed at defining or predicting the behavior of atomic structures based on their electronic structure. Relatedly, the explosion of interest in machine learning and artificial intelligence has driven dramatic advances in numerical optimization techniques, including through electronic structure elucidation. Naturally, these techniques are increasingly being applied to the solution of complex problems in materials science. A key problem in this context is the determination of the ground state (equilibrium) properties of electrons in crystalline materials.|0005| Historically, self-consistent solutions have been employed, generally involving multiple computational iterations to derive the electronic structure, properties of electrons in crystalline materials. The most commonly used approach is the self-consistent field (SCF) Kohn- Sham density functional theory (KS-DFT) in which the ground state properties are extracted from an auxiliary system of non-interacting electrons subjected to an effective periodic potential thatdepends on the electronic density (see - and The Kohn-Sham equations for the occupied and potentially occupied single-particle orbitals (Bloch waves in the case of periodic systems) are solved to self-consistency yielding properties like the equilibrium density, the energy and most important, the Kohn-Sham eigenvalues, which are generally considered reliable proxies for the single-particle excitation energies, i.e., the band structure of the material.

[0006] A most attractive feature of the KS-DFT method is the ease with which the orthogonality of the occupied and potentially occupied orbitals is guaranteed, thus ensuring the satisfaction of the Pauli exclusion principle. On the other hand, SCF techniques often struggle with convergence issues due to their susceptibility to encountering saddle points on the energy landscape, which can result in erroneous outcomes or even complete convergence failure as when, for example, the solution keeps oscillating back and forth between two densities, neither of which is self-consistent. Furthermore the process of computing orbitals and eigenvalues by direct diagonalisation of the Kohn-Sham Hamiltonian is very time-consuming (scaling as N3) when the number of electrons per unit cell, N, is large.

[0007] More recently, direct optimization has been used as an alternative to self-consistent solutions, such as the Kohn-Sham equations. This approach requires a minimization of the free energy with respect to molecular orbitals and an occupation matrix, subject to the constraint of orthonormality for the orbitals and the Pauli exclusion principle for the occupation matrix. It is of interest primarily because of its promise to yield more efficient and more reliable convergence than the conventional SCF method. In fact, the KS approach is just a special case of the broader problem of finding the global minimum of the energy functional with respect to orthonormal orbitals and their occupations. Simply put, the KS equation is the condition for the existence of a stationary point (not a minimum but, in general, a saddle point) of the energy functional, whereas direct optimization targets the absolute minimum of the energy functional.

[0008] A crucial aspect of direct optimization is the need to optimize simultaneously the potentially occupied orbitals and their fractional occupations. This is especially important in metals due to ambiguous occupation of energy levels across k-space. It turns out that by a slight extension of the original objective function, i.e., replacing the energy by the free energy through the additionof a suitable entropy term, it is possible to guarantee that the optimal orbitals are simultaneous eigenstates of the effective KS hamiltonian and of the occupation matrix, as they should. Tn this manner the KS eigenvalues emerge “spontaneously” from the optimization, without requiring the explicit solution of the time-consuming eigenvalue problem. This would not be the case in the plain energy minimization formalism.

[0009] To perform the minimization, some technologies use a two-loop optimization scheme in which the free energy is first minimized with respect to the occupations matrix (at fixed orbitals) and only in a second stage with respect to the orbitals. Despite its ingenuity, the two-loop method requires one to work with a non-diagonal occupation matrix during the intermediate stages of the calculation, which therefore is computationally very demanding.

[0010] Accordingly, there exists a need for a simpler approach to direct optimization, that avoids or reduces the requirement for iterative improvement, and reduces computational resource demands, particularly where such improvement requires time-consuming calls to external code libraries.SUMMARY

[0011] In view of the abovementioned drawbacks of the prior art, a simpler approach is proposed herein whereby, in some embodiments, the orbitals and the occupations are generated by a QR decomposition algorithm. The QR decomposition algorithm automatically enforces the physical constraints and replaces the two-loop minimization with a single direct minimization with respect to two independent sets of rectangular matrices.[0012| The present technology provides a method for use in predicting electronic structure properties of an atomic system of arbitrary definition, comprising the steps of: a) receiving parameters of the atomic system; b) determining a total electronic free energy of the atomic system, by: i) constructing a matrix of coefficients of plane-wave expansions of electronic orbitals corresponding to the atomic system, ii) calculating a plurality of energies, the plurality of energies comprising Hartree potential energy;iii) obtaining the total electronic free energy and derivative of the total electronic free energy with respect to the matrix of coefficients, based on summing the energies; iv) determine if the derivatives reflect a stable local minimum, and v) repeating steps b)i) to b)iv) if the derivatives do not reflect a stable local minimum, by applying gradient descent to the matrix using the derivatives; and c) generating a prediction of interaction properties of the atomic system with other atomic systems, by reference to the total electronic free energy, and derivative.

[0013] Also disclosed is a system for use in predicting interaction properties of an atomic system, comprising: memory; at least one processor; a receiver for receiving parameters of the atomic system; a total electronic free energy estimation module for determining a total electronic free energy of the atomic system, by: constructing a matrix of coefficients of plane- wave expansions of electronic orbitals corresponding to the atomic system; calculating a plurality of energies, the plurality of energies comprising Hartree potential energy; and obtaining the total electronic free energy and derivatives of the total electronic free energy with respect to the matrix of coefficients, based on summing the energies; a derivative module for determining if the derivatives reflect a stable local minimum; and a predictor for predicting interaction properties of the atomic system with other atomic systems, by reference to the total electronic free energy, and derivatives, wherein the memory stores instructions that, when executed by the at least one processor, cause the total electronic free energy estimation module and derivative module to repetitively determine the total electronic free energy of the atomic system and determineif corresponding derivatives reflect a stable local minimum, where, for each iteration, gradient descent is applied to the matrix using the derivatives of the total electronic free energy

[0014] Throughout this specification and the claims which follow, unless the context requires otherwise, the word “comprise”, and variations such as “comprises” and “comprising”, will be understood to imply the inclusion of a stated integer or step or group of integers or steps but not the exclusion of any other integer or step or group of integers or steps.BRIEF DESCRIPTION OF DRAWINGS

[0015] Some embodiments of the present invention will not be described, by way of nonlimiting example only, with reference to the accompanying drawings in which:

[0016] FIG. 1 is a flowchart of an algorithm for predicting interaction properties of a manybody system.

[0017] FIG. 2 illustrates a computational graph for total electronic free energy minimization.10018| FIG. 3 illustrates the fine-tuning process for band structure calculations for diamond.

[0019] FIGs. 4 and 5 show comparative analyses for band structures ascertained using various methods.

[0020] FIG. 6 illustrates geometry optimization.

[0021] Fig. 7 shows convergence curves of various optimizers on lithium and carbon (diamond) crystals - the y-axis represents the logarithmic scale of the shifted total electronic free energy.[ 0022| FIG. 8 illustrates a system for performing the method of FIG. 1.DETAILED DESCRIPTION

[0023] Embodiments of the methods described below, and used for the elucidation of electronic structure for the purposes of, for example, compound, material and drug discovery, leverage the invariance of free energy under unitary transformations to work with a diagonal occupation matrix at all stages of the calculation. Two physical constraints (1) orthonormality of the orbitals and (2) fractional occupation numbers comprised between 0 and 1, are “hard-wired” in the calculation where orbitals and occupation numbers are generated through an operation known as “QR decomposition”. QR decomposition generates a rectangular matrix of orthonormal columns out of an arbitrary rectangular matrix of the same shape.

[0024] Using QR decomposition and the resultant rectangular matrices, the conventional two- loop minimization can be replaced by a single direct minimization (a single-loop) with respect to two independent sets of rectangular matrices. While the occupation matrix remains diagonal throughout the process, the effective KS hamiltonian turns into a diagonal matrix “spontaneously” as the minimum of the free energy is reached - i.e., self-diagonalises. Consequently, the optimal orbitals are automatically eigenfunctions of the KS equation without warranting solving the KS equation itself.

[0025] Some embodiments are implemented in JAX - a machine learning framework for transforming numerical functions - to afford end-to-end differentiability As a result, the present framework can be used, for example, to compute interatomic forces as derivatives of the energy with respect to atomic coordinates.

[0026] The single loop process maintains efficacy, in particular, the ability to extract heuristic band structures without explicit diagonalisation of the Kohn-Sham Hamiltonian, and also the flexibility to handle nonlocal energy functionals, e.g., hybrid functionals, without a significant increase in computational effort. Moreover, the complexity of the calculation remains formally unchanged while the actual computational time can be greatly reduced using the QR decomposition algorithm: in other words, the prefactor of N3can be reduced by as much as an order of magnitude.

[0027] An embodiment of the method is shown in FIG. 1. The method 100 is used in predicting interaction properties of a many-body system, by producing information that facilitates material property analysis. As used herein, a "material" comprises an atomic system, including many-body systems, such as chemical compounds, proteins, DNA, drugs and corresponding receptor sites, catalytic surfaces, semi-conductors and others.

[0028] The method 100 takes input variables through a receiver system - Step 1. The receiver system may take the input variables from an external database or memory, user input or other source. The input variables include the geometry and identity of the atomic system - i.e., many- body system for which the electronic structure is being elucidated. Other input variables can include the desired electronic structure basis, and the choice of gradient descent algorithm - the rate of convergence depends on the gradient descent algorithm such as Adam, Yogi or Stochastic Gradient Descent, and can in some embodiments be fixed (i.e., not chosen by the user). The method100 employs direct gradient descent, returning the electron density corresponding to the total electronic free energy and all derivatives associated with any input variables - e.g., if the desired derivatives (such as the derivative of energy with respect to nuclear coordinates, which facilitates geometry optimisation - see Equation (37)) have been specific in the input variables. The processes performed by the method 100 are implemented in a deep learning framework, to facilitate back- propagation of gradients, all functions are differentiable so that the total electronic free energy can be optimized by gradient descent, and does not require self-consistent field iterations, even though the convergence result - i.e., the stable local minimum of the total electronic free energy, and its related derivatives - is self-consistent.

[0029] After receiving inputs at Step 1, Steps 2 to 8 involve constructing a matrix of coefficients of plane-wave expansions of electronic orbitals corresponding to the atomic system (C). Various energies, including Hartree potential energy, are then calculated at Steps 9 to 14. A total electronic free energy [ A(n)] and derivative (with respect to C) for the atomic system are then obtained at Step 15, by summing the energies calculated at Steps 9 to 14. The method 100 necessarily drives towards identifying a minimum free energy in the electronic structure of the atomic structure. Therefore, after Step 15, Step 16 determines whether the derivatives of the total electronic free energy reflect a stable local minimum - e.g., the derivatives are each "0" and represent the base of a saddle, or are otherwise sufficiently close to zero that further refinement will achieve negligible change (e.g., below a predetermined threshold such as 10-6Hartree). If the derivatives do not reflect a stable local minimum then, per Step 17, the derivatives of the total electronic free energy (i.e., the gradient of a function representing the total electronic free energy) is used to adapt the matrix Q and thus modify C constructed in Steps 4 and 8, respectively. If the derivatives reflect a stable local minimum, then the electronic structure has been found and the total electronic free energy and derivatives can be used to predict interaction properties of the atomic system - e.g., to output a material property analysis. Material properties can include, but are not limited to, phonon spectra, vibrational spectra, nuclear magnetic resonance shielding tensors, dipole moment, quadrupole moments, octupole moments, excitation energies, charge transfer properties, binding energies, electron densities, electron spin-densities, others. These properties can be determined by querying a library of electronic structures and corresponding material properties, to generate a report comprising the electronic structure, or portions thereof,and corresponding material properties as extracted from the library. The desired derivatives inputted with the input variables received at Step 1 will depend on the material properties sought to be discovered - e g.. Nuclear Magnetic Resonance Shielding Tensors require the derivative of the magnetic field, and the nuclear magnetic spin momentum (which would need to be programmed as inputs), phonon and vibrational spectra both use variations of the derivatives with respect to nuclear coordinates that geometry optimization uses, and polarisability and hyperpolarisability can be found as derivatives of the energy with respect to an electric field (the electric field will have to be programmed as an input too).

[0030] With regard to identifying local minima, the derivative of energy with respect to C is a multidimensional object. The derivative in every direction is 0 (or effectively 0) at a stable local minimum. Non-stable local minima will have at least one or more directions that have negative gradient. Negative derivatives mean it is possible to decrease the energy further by moving in that direction. Positive derivatives are also possible, and mean the energy will increase if that direction is followed. By consistently following negative derivative directions one will necessarily find a stable local minimum, eventually.

[0031] Step 16 and 17 seek to minimize the energy functional, being the sum of the energies calculated at Steps 9 to 14. While method 100 calculates the kinetic energy, exchange correlation energy, external energy and Hartree energy, some embodiments do not require all four energies to be calculated, and other embodiments may calculate different energies instead of one or more of those referenced in Fig. 1.

[0032] Step 3 involves initializing a random matrix (X). Per Step 4, QR decomposition is then performed on matrix X, in the first iteration, and is updated using derivatives of the total electronic free energy, for each subsequent iteration To apply QR decomposition, to consider the Hohenberg-Kohn density functional theory (DFT) as generalized by Mermin to thermal ensembles. According to this the equilibrium density of the many-electron system is obtained by minimizing the free energy A[n], regarded as a functional of the electronic density n(r):(1)

[0033] Tn Equation (1), Tsis the non-interacting electron kinetic energy of the non-interacting electronic system, EHis the Hartree energy, Excis the exchange correlation energy, Vextis the potential energy from the external field due to positively charge nuclei, S and r being the entropy of the system and r being the specified atomic coordinates and T being the temperature as input parameter

[0034] The equilibrium density is assumed to be obtainable as the density of an auxiliary noninteracting system, the Kohn-Sham system, described by a single particle Hamiltonian.

[0035] For periodic systems, let Ψi(k, r) denote a sufficiently large set of orthonormal eigenfunctions (collectively denoted by {Ψ}) characterized by a band index i = 1, ..., I (I larger than the number of electrons, N, per unit cell) and a Bloch wave vector k, with K wave vectors in each band, where K is the number of unit cells across which Born-vonKarman periodic boundary conditions are applied. This set of eigenfunctions will represent orbitals referred to herein as “potentially occupied orbitals” referred to herein as the actual occupations are specified by an hermitian occupation matrix f with matrix elements Γij(k) = Γij(k). such that density takes the form: (2)(sum over k in the first Brillouin zone) where: (3)

[0036] The occupation matrix must have the following properties: (i) its eigenvalues {fi} are all comprised between 0 and 1, in compliance with Pauli’s exclusion principle: (4)(ii) its trace (summed over k) equals the total number of electrons:(5)(iii) It must be diagonalisable together with Hs, i.e. it must be possible to find eigenstates of Hs, with eigenvalues εi(k) , which are also eigenstates of T , with eigenvalues fi(k) comprised between 0 and 1 and adding up to NK when summed over all i and k. As shown below, this condition is essential for the existence of stationary points of the functional Ω[n] .

[0037] The free energy is now represented as a functional, referred to as the "objective function" that is intended to be optimised:where E [Γ , {Ψ }] is the energy and the density is given by Equation (2).

[0038] The objective function is invariant with respect to a unitary transformation of the potentially occupied orbitals among themselves, combined with the corresponding transformation of the occupation matrix leaves the density, the energy, and in general all the physical properties of the system unchanged. This implies that the solution of the minimization problem cannot be unique, since a unitary transformation can be applied that does not change the physical properties.

[0039] This enables the Kohn-Sham Hamiltonian and hermitian matrix to be brought to diagonal form using a unitary transformation, such that εij(k) -> εi(k)δijThus, a set of optimal orbitals {Ψ } is unitarily equivalent to the set of eigenfunctions of the Kohn-Sham Hamiltonian, with eigenvalues εi(k).

[0040] In addεtion, by constraining the objective function such that it is stationary with respect to infinitesimal variation of Γij(k) yields, by differentiation of equation (6):where μ (chemical potential) is the Lagrange multiplier that enforces the trace constraint, and we have relied on Janak’s theorem (https: / / journals.aps.org / prb / abstract / 10.1103 / PhysRevB.18.7165) to obtain

[0041] A unitary transformation that diagonalises the right hand side of Equation (7) must also diagonalise the left hand side, which is possible only if Γ is disagonalised by the same transformation.

[0042] The free energy can ultimately be written in a way that the interaction contribution of the Kohn-Sham eigenvalues is removed and the interaction energy is restored. As a result, every stationary point of the free energy (even if it is not the true minimum) is associated with a self- consistent solution of the Kohn-Sham equation in which the Kohn-Sham Hamiltonian and the occupation matrix are simultaneously diagonal.

[0043] When directly optimizing the objective function, Equation (1 ), there is no loss of generality in assuming that the occupation matrix is diagonal in the basis of the potentially occupied orbitals. Thus searching for the minimumA will be conducted over diagonal occupation matrices, with diagonal matrix elements fi(k), while the physical constraints (positivity, Pauli exclusion principle, orthonormality of the orbitals) will be enforced in the parameterizations of the occupations and the potentially occupied orbitals. A single steepest descent search can then be used, with a single objective. In particular, the optimization procedure commences by expanding the potentially occupied Bloch wave functions Ψi(k, r) in plane waves across the material, presently assumed to be a periodic crystal lattice:where the sum runs over reciprocal lattice vectors G and Ω is the volume of the unit cell. The coefficients CG i(k) are indexed by a reciprocal lattice vector G (row index), and a band index i, and are functions of k in the first Brillouin zone. These coefficients are sub) ect to an orthonormality constraint for a given k:(9)The coefficients are then arranged in an elongated rectangular array C(k) of dimensions Gx I, where G is the number of basis states (plane waves) and I is the number of potentially occupied bands. In this notation the orthonormality constraint takes the form: (10)where C^(k) is a rectangular I x G matrix, IIis the I x I identity matrix and the dot denotes the usual matrix multiplication

[0044] In the standard KS-DFT approach the orthogonality of the Kohn-Sham orbitals is guaranteed by the hermiticity of the Kohn-Sham equation. Moving away from KS-DFT a crucial consideration is whether the orthogonality constraint can be incorporated for CG(k). To address the challenges of exponential translation methods only being applicable to square matrices - resulting in the number of basis functions greatly exceeding the number of bands, which is inefficient where only the eigenstates yielding lowest eigenvalues are required - and the requirement to perform multiple eigendecompositions and matrix multiplications, a reparametrisation of the orthogonality constraint for the matrix of the coefficients C(k) is performed. This is achieved by applying QR decomposition at Step 4, to transform the problem to a more computationally efficient form. Starting with an arbitrary Gxl matrix X(k) that, in the first iteration is randomly initialized at Step 3, the QR decomposition implemented at Step 4 expresses X(k) as the product of a G x I matrix with orthonormal columns (Q) and an I x I upper triangular matrix (R). Discarding R results in an orthogonal matrix C. Thus, the reparameterization of the orthonormal coefficients can be expressed as:(11)where Xkdenotes the set of variational parameters that are subject to optimization and QR[Xk] denotes the orthonormal part of the QR decomposition of the matrix Xk. This form ensures that the orthonormality of the wave function coefficients is maintained during the optimization process.

[0045] By virtue of the invariance of the free energy with respect to unitary transformations, a search for minimum free energy can be restricted to diagonal occupation matrices; (12)

[0046] These occupations can be arranged in a square diagonal matrix of dimension IK x IK, denoted by F (i.e. F is the practical implementation of Γ), in:

[0047] Diagonal elements of F are between 0 and 1 and add up to NK, wherein N< I. To satisfy this constraint, the following reparametrization is used: (14)matrix of orthonormal columns which is generated by applying the QR decompositions to an arbitrary rectangular matrix Y, being an overlap matrix initialized at Step 5, of dimensions IK x NK and discarding the upper triangular part:(15)

[0048] The elements of the matrix Y are variational parameters subject to optimization. Since, by construction: (16)where INKis the NK x NK identity matrix we immediately see that: (17)

[0049] In addition, the IK x IK square matrix F = V • Vtis idempotent: (18)

[0050] This implies its diagonal elements (which coincide with the diagonal elements of F, hence with the occupation numbers fi(k)) are all comprised between 0 and 1 Indeed:from which the desired inequality: (20)follows.

[0051] Thus, reparametrization transforms the constrained search problem into an unconstrained minimization problem with objective function: (21)with variational parameters Y and Xk[here we have set fik= fi(k)]. This approach incorporates the constraints of finite-temperature KS-DFT. Moreover, the entire process is differentiable, enabling straightforward application of gradient-based algorithms.

[0052] The unconstrained minimization problem set out at Equation (21) has an optimal solution where the Hamiltonian matrix is diagonalised. The trace of the product of the occupation matrix f and the Hamiltonian matrix H has a variational lower bound:(22)where εi(k) are the eigenvalues of the Hamiltonian matrix listed in ascending order: (23)and fi(k) are the eigenvalues of the occupation matrix Γ, also listed in ascending order: (24)

[0053] The left-hand-side trace will achieve the minimum only when Γ and H share the same eigenvectors, i.e., they are simultaneously diagonalisable. By enforcing the occupation matrix to be diagonal, and ensuring it satisfies the conditions outlined in Equations (4) and (5), the optimization of the trace on the left-hand side will inherently lead to the diagonalisation of the Hamiltonian matrix. This is because the inequality becomes an equality only when the eigenvalues of the Hamiltonian matrix are aligned appropriately; otherwise, the left-hand side is always greater than the sum of the products of the eigenvalues of Γ and H . This methodology provides a variational lower bound for the energy in the direct optimization approach.

[0054] Using the foregoing, the method 100 seeks to minimize the total free electronic energy subject to the orthonormal constraint, defined by:using the Kohn-Sham formalism, employing a single slater determinant formed by orthonormal single-particle wave functions as ansatz. The electronic density is then given byand the ground state density can then be obtained by Equation (25), whereEkin, Eext, Ehar, Exc, Esdenote the kinetic energy calculated at Step 12, external potential energy calculated at Step 14, Hartree energy also calculated at Step 14, and the exchange- correlation energy calculated at Step 13, and other energies as applicable such as entropic energy calculatedat Step 11 This constraint optimization problem can be solved via the Lagrange multiplier method, and the first order condition of the Lagrangian gives the Kohn-Sham equation: (26)where HKS[n] is the KS Hamiltonian that depends on the electronic density n(r), andare the KS eigenvalues. The KS equation can be solved via SCF.

[0055] Applying the above to a periodic crystalline system input at Step 1, the system is simulated by solving Schrodinger's equation in a Bravais lattice in Steps 2 to 8, by capturing the periodic structure of the material by tiling a parallelepiped (unit cell) defined by lattice constantsin space. The lattice is subjected to periodic boundary condition (PBC), which means the lattice itself is infinitely tiled in space. We denote the number of unit cell in the lattice as M = x M2x M3, and each cell in the lattice can be indexed by m: = (mlrm2, m3) where in, — 0,1,— 1 for iThe external potential is periodic in the lattice, so byBloch’s theorem, the electronic wavefunction only differs by a phase shift between different cells, i.e.:is a function periodic in the lattice, bi are the reciprocal lattice vector. The wavevectors kmtakes discrete values due to the PBC.

[0056] Similarly, if KS formalism is applied to a periodic system, the electronic density will be distributed among different wavevector km:

[0057] While Equations (25) and (26) still hold, the periodic system produces I X M orbitals Equation (27) then simplifies the orthonormal constraint sinceorbitals with different wavevector are automatically orthogonal, and the KS equation decouples into Al equations: (29)which can be solved separately for each km. The KS eigenvalues £i mhere determines the band structure of the material.

[0058] is the parametrized - e.g., via a Fourier basis. The mesh of the Fourier basis isdefined with size N =within the confines of the unit cell, then at each lattice point whcre n: =, the value of can berepresented by a linear combination of plane waves;n£>. |les onthe reciprocal FFT-mesh, andare the Fourier coefficients of With this parameterization, direct optimization of Equation (25) can be performed. Specifically, to enforce the constraint for each mwhere QR is a function that returns the orthogonal part from a QR decomposition, isthe learnable parameter. The integrals in the energy functional has analytical formulaein Fourier space.

[0059] Per Step 12, the kinetic energy can be directly calculated from the wavefunction coefficients in C. At Step 9, the inverse Fourier transform of C is then calculated, yielding the wavefunction in real space, and is then squared to determine the electronic orbitals and electronic density function in real space. From the electronic density function, the exchange correlation energy can be calculated per Step 13. Step 10 applies a fast Fourier transform to the electronic orbitals to find the density function in reciprocal space. The density function in reciprocal space isused to calculated the Hartree energy and external potential energy per Step 14. In the present embodiment, entropy (S) is also calculated from the occupation matrix F These energies are summed at Step 15 to produce the total electronic free energy, and derivatives of the total electronic free energy are calculated with respect to matrix C

[0060] Fig. 2 shows the computational graph for total electronic free energy minimization, with each operation being differentiable. In particular, after calculating learnable parameter X (at 200), QR decomposition (at 202) is applied to generate C (204), from which the kinetic energy Eki is determined (at 206). Inverse fast Fourier transform is then applied (at 208) to produce the wavefunction in real space (at 210). That wavefunction is squared (at 212) to produce the electronic orbitals from which the density function is determined (at 214), that enables the exchange correlation energy to be calculated (at 216). A fast Fourier transform is applied (at 218) to theelectronic orbitals to generate the density function in reciprocal space (at 220) from which the Hartree energy,and external potential energy,are calculated (at 222 and 224, respectively). The energies are the summed (at 226). This process can be performed in multiple iterations, where each iteration, other than the first iteration, starts with X updated using the gradients (derivatives) of the total electronic free energy calculated from the previous iteration.

[0061] In particular, after the total electronic free energy and its derivatives are calculated, Step 16 involves determining if the derivatives reflect a stable local minimum. If not, per Step 17, the derivatives are used to modify QR decompositions at Steps 4 and 6. The updated matrices Q and V (used to construct occupation matrix F) are then used in the next iteration of energy calculations. If the derivatives reflect a stable local minimum, then the total electronic free energy and derivatives are used for material property analysis. If a stable local minimum is found, per Step 18, the method reports the total electronic free energy (TFE), along with the derivatives of TFE (e.g., with respect to X and Y), and results of any material property analyses - e g. consulting an external reference library for comparable structures (or substructures - parts of the elucidated structure), and extracting material properties from the reference library.

[0062] As mentioned above, with regard to Equation (29), methods herein can be used to determine the electronic band structure of a material. The electronic band structure is critical in the domain of solid-state physics, instrumental in delineating the electronic behavior within crystalline solids. Serving as a holistic representation, the band structure imparts pivotalinformation regarding a material’s intrinsic properties. It describes whether a material exhibits conductive, insulating, or semiconducting behavior. Additionally, it provides insights into the material’s optical properties, magnetic tendencies, and thermal conductivity attributes, thereby offering a comprehensive understanding vital for diverse technological applications.

[0063] The abovementioned processes can be adapted specifically for the computation of electronic band structure. In particular, for a fixed point the KS Hamiltonian represented inthe Fourier basis is a matrixused to help find the electronic band structure, which differs from the F used previous to now), where the matrix element is given by:

[0064] Using the Fourier coefficient the KS equation can be written in matrixform:which require a N x N matrix to be diagonalised for each sampled However, to minimizeenergy, on the lowest I eigenvalues of this matrix are of interest. To extract these, the sum of eigenvalues is minimized for the smaller I x I matrix i.e., calculate:where the coefficients Cmare obtained from the QR decomposition of XmThere are two main advantages of this approach: (1) The orthogonality constraint is seamlessly converted to an unconstrained form using the QR reparameterization trick, which makes the band structure optimization problem a direct optimization process that aligns well with the deep learning paradigm; (2) Typically, I « N, so by only storing the I x I matrix the memory requirements are drastically reduced.

[0065] Thus, the band structure elucidation, the method comprises: selecting k-mesh size M1X M2X M3, then running the direct optimization process to obtain the ground state electron density n*; for each kmin the A-path, use the ground state density n* to construct the matrix representation of the KS Hamiltonian, then obtain coefficients Cm; and perform a one-step diagonalisation oftoobtain KS eigenvalues which represent the electronicband structure.

[0066] Although in the above algorithm, the calculation for each kmcan be parallelized, when running on hardware with limited memory each kmmust be processed sequentially. Instead of randomly initializing the parameters Xmfor each km, the parameter Xm' obtained from the previous step as an initial guess and fine-tune it to determine the parameters for the subsequent k- point. Given the continuity of the band structure, only a few fine-tuning iterations are needed, as shown in Fig. 3.

[0067] Fig. 3 illustrates the fine-tuning process for band structure calculation on the diamond structure using LDA functional. In the top image, the band structure is optimized via the present direct optimization method The valence and conduction bands are represented by red and blue lines, respectively. A distinct band gap is highlighted in the grey area. The bottom image presents the training curve corresponding to the minimization process. For optimization, the Adam optimizer is used, employing 4000 iterations for the first k point and 100 for fine-tuning the subsequent k points.

[0068] This method significantly cuts down computation time without demanding additional computational resources.

[0069] Differentiability is especially useful for jointly optimizing other parameters in the system as in the task of geometry optimization. The problem of geometry optimization is to find the equilibrium geometry, which is atom configurationthat minimizes the total electronic and ionic energy of the system where C is the basis coefficientsof the electronic wavefunction and is subjected to the orthonormal constraint, and is thenuclear repulsion energy:where are atomic positions, are the corresponding charge distributions, indexed by = 1, 2,. . . , L, where L represents the total number of atoms within the unit-cell,allows consideration of interactions beyond the primitive cell where m — (m1m2, m3) represents an integer translation of the lattice vectors, we exclude the terms whereGeometry optimization is an important problem since equilibrium geometry provides a fundamental understanding of a material. An accurate atomic structure, obtained through optimization, is essential to precisely compute electronic structure, and other properties of the material like mechanical, thermal, and vibrational properties.

[0070] From the problem definition, the most natural way to solve geometry optimization is to perform constraint minimization of the total electronic free energyby varying the wavefunction C and geometry R at the same time, i.e.

[0071] However, when using the SCF approach, to ensure convergence the update in C is usually dampened, i.e. history-dependent, and interleaving geometry update will disrupt the convergence of SCF. Therefore when using SCF, the force method is usually used, which is usually given by the Hellman-Feynman theorem (note FHFis a new matrix of Hellman-Feynmann forces, which differs from the F matrices used previous to this).

[0072] The differences between the algorithms of the two approaches are set out below.

[0073] The present approach is fully differentiable, allowing direct joint optimization of Equation (36), utilizing gradient descent. It is important to note that the nuclear repulsion indicated in Equation (35) is not straightforwardly evaluable in real space due to its slow convergence. Consequently, the employment of a technique known as Ewald summation is useful.

[0074] EXPERIMENTS

[0075] Experiments were conducted on real-world crystals and materials, focusing primarily on three subtasks: electronic band structure prediction, geometry optimization and convergence analysis.

[0076] The method was implemented in deep learning framework JAX. All the experiments were conducted on an NVIDIA Al 00 GPU with 40GB memory. All remaining experiments with other implementations were conducted on a server powered by an Intel Xeon CPU @ 2.10GHz and furnished with 64GB of memory.

[0077] Firstly dealing with electron band structure prediction, the accuracy and efficacy of the present direct optimization method was compared to existing implementations in the context of electron band structure prediction, which is a crucial step for uncovering the electronic properties of materials. Four distinct crystal structures were tested: lithium, aluminum, with crystal structures obtained from the Open Crystallographic Database (OCD), and carbon (diamond), and silicon, with crystal structures obtained from the critic2 library.|0078| To generate baselines, comparisons were conducted with several packages, including Quantum ESPRESSO (QE), Fritz Haber Institute ab initio materials simulation (FHIaims) package, and GPAW QE conducts DFT calculations using plane-wave basis sets and pseudopotentials to represent electron-ion interactions. Since the present methods use a planewave only basis, in order to make the closest comparison possible with an established method, in some cases the number of electrons contained in the pseudo-potential is reduced to zero toapproximate an all-electron plane-wave basis. Where unmodified pseudopotentials are used, they are generated from the pslibrary. FHI-aims implements numerical atomic orbital basis sets and provides the option to use periodic boundary conditions. All FHIaims calculations in this work use tight basis settings. GPAW utilizes a plane-wave basis with the pseudopotentials for which the recommended defaults were used. All calculations presented in this work use the LDA functional.

[0079] The results of the lithium and carbon band structures are illustrated in Fig. 4 and aluminium and silicon band structures are illustrated in Fig. 5. In Fig. 4, insets (a) and (c) show a comparative analysis of the valence / conduction bands with QE, FHIaims, and GPAW, implemented on various crystals. Insets (b) and (d) of Fig. 4 exhibit the comprehensive band structures, encompassing the core bands. In Fig. 5, insets (a) and (c) show a comparative analysis of the valence / conduction bands with QE, FHIaims, and GPAW, implemented on various crystals. Insets (b) and (d) of Fig. 5 show an exhibition of the comprehensive band structures, encompassing the core bands. In these calculations a 1 x 1 * 1k -point mesh is used. A cut-off energy of 800 Hartree the present method and QE, GPAW uses a cut-off energy of 400 Hartree. Adam optimizer was adopted. QE uses an empty pseudopotential in all cases, and cold smearing with a smearing temperature of 0 01 Ry for aluminium and lithium, and fixed occupations for diamond and silicon.

[0080] The band structure indicates whether a crystal is an insulator, semiconductor, or conductor. A substantial, clear band gap suggests that the crystal primarily exhibits insulating properties. The results of the present method, aligning consistently with existing implementations, reveal no band gap for lithium and aluminum, suggesting metallic characteristics. In contrast, a significant band gap is observable for carbon, indicative of insulating properties for diamond, while a small band gap for silicon implies semiconductor behavior. This large band gap of diamond means that, at room temperature, there are no available energy states for electrons in the conduction band, and thus, no electrical current can flow through, rendering diamond an excellent insulator, whereas for metals, the valence band and the conduction band overlap, meaning there are always available energy states for electrons to move to. The conclusions drawn from our study confirm the conventional understanding of these materials and align with the theoretical predictions provided by other established implementations.

[0081] As an all-electron plane-wave method, the present approach illustrates the complete band energy level spectrum of crystals, encompassing the core band. The examination of corebands is pivotal for achieving complete representation of a material’s electronic structure, which, in turn, facilitates the understanding of various associated physical and chemical properties of the materials. The illustrative results of this analysis are depicted in Figs. 4 and 5, insets (b) and (d), with the bands names annotated based on the atomic orbital of origin on the right side for reference. The comparison reveals that the present method exhibits a high degree of alignment with QE, yielding a congruent band spectrum. This congruence underscores the reliability and accuracy of the present method in analyzing and representing the intricate features of band structures, inclusive of the core bands.

[0082] Secondly, dealing with geometry optimization: the structure of a primitive unit cell from a diamond-structured carbon crystal was examined, consisting of two carbon atoms. In the experiments, one carbon atom is anchored at the origin, and the method seeks to ascertain the position of the second carbon atom. This atom is moving along the plane highlighted in red, as depicted in Fig. 6, inset (d). The corresponding potential energy surface (PES) is computed and illustrated in Fig. 6, inset (a). The optimum is denoted by the star point, accurately representing the experimentally known diamond structure. Observations reveal the presence of multiple local optima and saddle points interspersed across the PES, underscoring the complexity of the energy landscapes typically found in atomic systems. In the experiment, a cut-off energy of 40 Hartree is employed along with a 24 x 24 x 24 FFT-mesh grid, and a 2 x 2 x 2 k-poiiil mesh. Three distinct optimizers were used: SGD, Adam, and Yogi, with a learning rate at 1 x 10'4.

[0083] In Fig. 6, the PES shown in insets (a) to (c) is computed by moving the atom in the centre of the unit-cell in inset (d) of Fig. 6 around the highlighted plane in a regular 50 x 50 grid.

[0084] yogi’, ‘adam’, and ‘SGD’ highlight optimization paths by the respective optimization method from their starting points as shown in insets (b) and (c) of Fig. 6. The global minimum is highlighted with a black star, and the geometry of the crystal at point A is given in inset (e) of Fig. 6, and consistent with the known structure of diamond.

[0085] The geometry of the local minima at point B to which the SGD method converges to is presented in inset (f) of Fig. 6, which can be seen to approach the known structure of graphite. The different band structures of A and B are shown in inset (g) of Fig. 6. Insets (d) to (f) of Fig. 6 were generated using xcrysden.

[0086] The results of optimizations with varied initializations are illustrated in Fig. 6, insets (b) and (c). Tn inset (b), both the Yogi and Adam algorithms successfully reach the global optimum, while the SGD optimizer lands at a local minimum. Inset (c) demonstrates situations where both Yogi and SGD become stuck at the saddle point, whereas Adam locates the global optima in the adjacent unit cell, reflected back in the plot. This local minimum is appears to be approaching a graphite-like atomic structure, showcased in inset (f). To emphasise the types of material analysis that geometry optimization can allow, the band structure of inset (e)-diamond and inset (f)- graphite-like are shown in inset (g). In particular, there is a clear band gap in inset (g), Part A, indicating insulator-like properties, while the disappearance of a band gap in inset (g), Part B, implies metallic properties.

[0087] Thirdly, convergence of different optimisers was analysed. In this task, attention is divided between two sub-tasks: firstly, analyzing the performance of different optimizers, and secondly, contrasting direct optimization with the SCF approach.

[0088] A variety of optimizers were tested within the present optimization framework, focusing on lithium and carbon crystal structures. All the optimizers employed the same initialization and a learning rate of 1 x 10'5. Each optimizer was run ten times, with the average and 95% confidence interval displayed in Fig. 7. Note that the energy is shifted and is depicted in a logarithmic scale.

[0089] The results illustrate that Yogi consistently outperforms the other optimizers in terms of achieving lower energy, even though it does not exhibit the most rapid decrease. It is also observable that some optimizers, particularly those incorporating a momentum mechanism, are prone to experiencing bounces, manifesting as an initial faster decrease in energy, followed by an increase, before resuming the decrease. Additionally, fluctuations are more pronounced in the convergence patterns of Adam, AdaBelief, and Adamax, showcasing a higher degree of variability as they approach convergence.

[0090] When considering convergence to a stable solution using SCF, it can sometimes be challenging using traditional approaches, which is not incurred using direct optimization as presently described. Tests were conducted to attempt to reproduce the PES generated in Fig. 7 using computational settings in QE and GPAW. As above a cut-off energy of 40 Hartree is used with a 48 x 48 x 48 FFT-mesh, 2 x 2 x 2 A-point mesh. QE begins all computations with an‘atomic+random’ initial guess, allowing a maximum of 300 SCF steps, and using a two-electron n-kjpaw pseudopotential, all other settings remaining default. Where not mentioned above, GP AW calculations use program default settings.

[0091] Where 100% convergence of all points on the PES is reporting using the present method, this test reveals 88.3% convergence for QE and 99.1% convergence for GPAW.

[0092] Fig. 8 shows a system 800 for implementing the method 100, in its various embodiments. The system 800 may be implemented in physical circuit components, or in software code, and all such implementations and combinations thereof are envisaged to be within the understanding of the skilled person in view of present teachings.

[0093] The system 800 comprises memory 802 that stores program code 804, execution of which by the processor or processors 806, cause the system 800 to operate. Operation involves receiving, at receiver 808, inputs such as those defined with reference to Step 1 of method 100, in particular the geometry and identity of the atomic system - as used herein, "identity" refers to the number of protons associated with an atom, by which traditional chemists can identify the element corresponding to the atom. Notably, in the present context, the method can be applied despite the number of protons of each atom not needing to be an integer, and indeed the derivative of the energy with respect to the proton number can nevertheless be computed to aid material discovery. The total electronic free energy estimation module 810 then constructs the matrix of coefficients corresponding to the atomic system, calculates the various energies of interest, including Hartree potential energy, and obtains the total electronic free energy and derivatives of the total electronic free energy with respect to the matrix of coefficients, by summing the energies. The system 800 also includes a derivative module 812 that takes the total electronic free energy and derivatives calculated by the total electronic free energy estimation module 810, and determines whether the derivatives reflect a stable local minimum. If not, the derivatives are used in a gradient descent approach to update the matrix for recalculation of the energies and thus the total electronic free energy and its derivatives. Once the derivatives reflect a stable local minimum, the derivative module passes the total electronic free energy and derivatives to the predictor 814. The predictor 814 predicts the electron band structure or optimizes the structure of the atomic system, or performs other operations, based on the total electronic free energy and derivatives.

[0094] The application of such a method are expansive, allowing for calculation of binding between molecules (as in drugs to cell receptors, or metal coordination inside a battery electrode), calculation of force constants (as in how materials respond to pressure, for applications such as testing new armors at the atomic scale), calculation of nuclear magnetic resonance shielding tensors (for aiding characterisation of unknown, but proposed materials undergoing analysis), or exploration of unknown material space through dynamic adjustment of nuclear identity as defining a target property through differentiable optimization. Adjustment of nuclear identity in the present context is significant, since non-integer nuclear values cannot occur in the real world, and hence cannot be argued as being pre-existing properties of original systems. While there is thus no original system to reference against, such non-integer analyses can help find new real systems nonetheless - i.e., systems comprising integer nuclear values.

[0095] It will be appreciated that many further modifications and permutations of various aspects of the described embodiments are possible. Accordingly, the described aspects are intended to embrace all such alterations, modifications, and variations that fall within the spirit and scope of the appended claims.

Claims

CLAIMS1. A method for use in predicting electronic structure properties of an atomic system of arbitrary definition, comprising the steps of: a) receiving parameters of the atomic system; b) determining a total electronic free energy of the atomic system, by: i) constructing a matrix of coefficients of plane-wave expansions of electronic orbitals corresponding to the atomic system; ii) calculating a plurality of energies, the plurality of energies comprising Hartree potential energy; iii) obtaining the total electronic free energy and derivative of the total electronic free energy with respect to the matrix of coefficients, based on summing the energies; iv) determine if the derivatives reflect a stable local minimum; and v) repeating steps b)i) to b)iv) if the derivatives do not reflect a stable local minimum, by applying gradient descent to the matrix using the derivatives; and c) generating a prediction of interaction properties of the atomic system with other atomic systems, by reference to the total electronic free energy, and derivative.

2. The method of claim 1, wherein the atomic system parameters comprise a geometry and an identity of the atomic system.

3. The method of claim 1 or 2, wherein calculating the plurality of energies comprises calculating one or more of a kinetic energy, exchange correlation energy, and external potential energy.

4. The method of claim 3, wherein calculating the plurality of energies comprises calculating the kinetic energy from the matrix of coefficients.

5. The method of claim 3 or 4, wherein calculating the plurality of energies comprises calculating an inverse discrete Fourier Transform and then a square of the matrix of coefficients, to determine the electronic orbitals and a density, and calculating the exchange correlation energy from the electronic orbitals and the density.

6. The method of any one of claims 3 to 5, wherein calculating the plurality of energies comprises performing a fast Fourier Transform of the electronic orbitals to find a density function in reciprocal space, and calculating the external potential energy from the density function in reciprocal space7. A system for use in predicting interaction properties of an atomic system, comprising: memory; at least one processor; a receiver for receiving parameters of the atomic system; a total electronic free energy estimation module for determining a total electronic free energy of the atomic system, by: constructing a matrix of coefficients of plane- wave expansions of electronic orbitals corresponding to the atomic system; calculating a plurality of energies, the plurality of energies comprising Hartree potential energy; and obtaining the total electronic free energy and derivatives of the total electronic free energy with respect to the matrix of coefficients, based on summing the energies; a derivative module for determining if the derivatives reflect a stable local minimum; and a predictor for predicting interaction properties of the atomic system with other atomic systems, by reference to the total electronic free energy, and derivatives, wherein the memory stores instructions that, when executed by the at least one processor, cause the total electronic free energy estimation module and derivative module to repetitively determine the total electronic free energy of the atomicsystem and determine if corresponding derivatives reflect a stable local minimum, where, for each iteration, gradient descent is applied to the matrix using the derivatives of the total electronic free energy.

8. The system of claim 7, wherein the atomic system parameters comprise a geometry and an identity of the atomic system.

9. The system of claim 7 or 8, wherein calculating the plurality of energies comprises calculating one or more of a kinetic energy, exchange correlation energy, and external potential energy.

10. The system method of claim 9, wherein calculating the plurality of energies comprises calculating the kinetic energy from the matrix of coefficients.

11. The system of claim 9 or 10, wherein calculating the plurality of energies comprises calculating an inverse discrete Fourier Transform and then a square of the matrix of coefficients, to determine the electronic orbitals and a density, and calculating the exchange correlation energy from the electronic orbitals and the density.

12. The system of any one of claims 9 to 11, wherein calculating the plurality of energies comprises performing a fast Fourier Transform of the electronic orbitals to find a density function in reciprocal space, and calculating the external potential energy from the density function in reciprocal space

Images