{'fraction_non_alphanumeric': 0.01960635017322115, 'fraction_numerical': 0.010754939658126166, 'mean_word_length': 4.6856060606060606, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Machine Learning Approach for Increased Throughput of Density Functional Theory Substitutional Alloy Studies

In this study, a machine learning based technique is developed to reduce the computational cost required to explore large design spaces of substitional alloys. The first advancement is based on a neural network approach to predict the initial position of ions for both minority and majority ions prior to ion relaxation. The second advancement is to allow the neural network to predict the total energy for every possibility minority ion position and select the most stable configuration in the absence of relaxing each trial position. This study a bismuth oxide materials system, (BixLayYbz)2 MoO6, is used as an model system to demonstrate the developed method and potential computational speedup. Comparing a brute force method that requires calculation of every possible minority concentration location and subsequent relaxation there is a 1.3x speedup if the NN is allowed to predict the initial position prior to relaxation. This speedup is a result in an average decrease of 4 hour reduction in supercell relaxation wall time for all trials. Implementation of the second advancement allowed the NN to predict the total energy for all possible trials prior to relaxation resulting in a speed up of approximately 37x. Validation was done by comparing both position and energy between the NN to DFT calculation. A maximum vector mean squared error (MSE) of 1.6x10−2 and a maximum energy MSE of 2.3x10−7 was predicted. This method demonstrates a significant computational that even more impressive for larger design spaces where the size of the design space is a function of a factorial number of minority components.


Introduction
Advancements in first-principle material modeling techniques is finding applications for modern day materials development that is required effective use of resources for new devices and structures. Modern materials require precise control of properties such as rapid phase response to external stimuli like pressure, light, magnetic field, so that meaningful uses are possible in modern day or can be expected in the near future. These modern materials (magnetic, ferroelectric, superconducting) are often multicomponent systems such as but not limited to high temperature superconductors, magnetic tunnel-junctions, and perovskite materials with complex magnetic structures 3 . The underlying principle cohered with for these first-principle methods is that the parameters of the formulated theory are fixed by the basic assumptions and equations of quantum mechanics.
During the past two decades first-principle calculations based on density-functional theory (DFT) 19 in the generalized gradient approximation (GGA) 4,5,15 or the local density approximation (LDA) unfolded as a successful approach to solve the electronic structure of matter 3 . DFT is a widely used electronic structure method implement to assist in understanding a a Department of Mechanical  wide range of material properties. The theory is able to reduce the many body Schrödinger equation to an effective single electron problem by relying on Hohenberg-Kohn theorem 12 and Kohn-Sham method 19 , thus making material property predictions computationally feasible 16 . The profound success of DFT for describing ground-state properties for vast classes of materials such as semiconductors, insulators, half metals, semimetals, transition metals, etc., at the nanostructure scale makes it one of the most used method for modern electronic structure analyses 3 . Its noted that the goal of these calculations is to gain insight on a well defined model so that studies can find and predict trends that can better assist in developing different levels of understanding for any system in question 3 .
Due to the extreme computational costs of most theoretical studies, limitations can and do arise when using approximation methods because accuracy is compromised in exchange for speed up time 24 . One of the most challenging aspects in modern theoretical calculations is to develop and apply an approximation method that expedite first-principle calculations speed up time without the loss of accuracy. Methods such as fragmenting the system 18 , construction of empirical potentials, corrections through statistical methods 20 , linear scaling 21 , or semi-emperical (SE) methods have been applied to attain accurate first-principle calculations that are still computationally effective 20,24 . In the case of SE Tight Binding approximations to DFT, the time expense is reduced by treating the Hamiltonian elements as parameters adjusted for the desired properties of the system such as band gap, effective mass, etc, also SE Tight Binding to DFT simplifies the Hamiltonian to nearest neighbor interactions among atoms 11 . Introduction of the exchange and correlation functionals allowed DFT to improve cost even further when compared to high level first principle methods like MP2 and CI, with very similar level of accuracy 20 . However, even the fastest DFT techniques, such as O(N) approach that is based on Local Orbitals, use up most of the computational time to iteratively formulate the Hamiltonian and solve self consistently for the ground state electron eigenstates 11 . Thus relies heavily on the initial state of the system, the closer the final state of the system is to the initial state specified as the input, the less iterations spent formulating the Hamiltonian and solving self consistently the ground state.
The reader should note that each molecule is unique and thus to explore different configurations, tailor or substitute various atoms into a structure, will only increase the computational cost. This is not to say one can avoid computational cost but in practicality greatly reduce it if the initial state of the system being explored is not just referenced based on some configuration but slightly refined based on the referenced configurations and based on what is being tailored or substituted. Thus exploring new configurations can be more computationally feasible, especially in a time sensitive world where industry applications relies heavily on a number of material systems such as geometries, boundary conditions need, and so on to be evaluated quickly and effectively to meet with the growing demand for production and application. Machine learning approaches have demonstrated viable solutions in attaining various forms of interacting and non interacting atomistic potential by utilizing regression algorithms in resent years 24 . Variety of applications such as chemistry and physics have successfully applied machine learning methods to predict reaction pathways 14 , formation energies 6 , excited energy states 10 , atomic forces and resonance chemical shifts, etc., to assist in searching and classifying material systems 24 . Also, computational material science has had several advancements made in applying various machine learning techniques such as predictions of DFT functionals 25 , mapping of spacial atomic data for predicting total energies 2 , and computation of electronic properties 23 in recent years. Thus machine learning express potential to predict molecular interactions with accuracy and reduction of computational cost.
To address the issue of computational expense, this study proposes a machine learning based method to predict positions of atoms (vectors that represent each atom) and the total ground state energy for each structure (unit cell). The predicted atomic positions that make up the unit cell and energy is then compared to that of the predicted DFT calculations. Note that this is predicting the fractional coordinates Fig. 1 Illustration of the Bi 2 MoO 6 in orthorhombic configuration that is used as the example structure system to test the machine learning approach for this study. The orthorhombic unit cell is made up of 24 O atoms 8 transition metal atoms (Bi) and 4 Mo atoms making a total of 36 atoms in a unit cell. La, and Yb are subbed into the 8 metal atom positions (Bi) (atom position numbers 3,16,17,18,19,20,21,22). These is the initial state of all DFT calculated structures and thus the numbers that correspond to each atom are the same for any configuration. All calculated DFT data is referencing the same atom position number to keep consistency, thus all attained calculations of position are for the same atom position number.
for each atom that makes up one primitive unit cell and the respected lattice constants. The machine learning model uses training sets based on DFT simulations done for various structures. This paper focuses on bismuth-based photocatalysts (Bi 2 MoO 6 ) in the orthorhombic configuration ( Figure 1) and predictions made are validated with DFT calculations for configurations of (Bi x M y ) 2 MoO 6 where (M = La, Yb). However, this method is independent of material system and depends only on referenced DFT simulations for a specific structure. The aim of this study is to predict fractional coordinates and the total ground state energy for a bismuth-based photocatalysts (Bi 2 MoO 6 ) by taking already attained DFT calculations for (Bi x La y ) 2 MoO 6 configurations as the training set to a neural network and then predict the ground state energy and the fractional coordinate vectors with respect to the DFT calculations for (Bi x Yb y ) 2 MoO 6 in the orthorhombic configuration. By using the training set, a much better approximation of the initial state of the system can be attained which reduces the effective iterative steps taken in the DFT calculations to solve 1-?? | 2 self consistently the ground state.

Methodology
Applying a neural network to any problem at hand requires input features (specified by user) to be mapped to some target output in some non-trivial way. In this paper, a supervised neural network (NN) is applied to a data set consisting of inputs (features that describe atomic classification in a unit cell) to be mapped on desired output (expressing atomic positions in a unit cell and total ground state energy).  2 Illustration A is the vector representation of input (charge number Z) and output (positions along the lattice vectors (a, b, c) in each unit cell). Illustration B is the artificial neural network architecture implemented for this paper. It's noted that the numbers indicated in this figure (such as a 1 or Z 36 ) are referencing the same numbers initially started for each unit cell in Figure 1 to keep referencing atomic number positions in the NN consistent.

Model Design Space
This study focused on the substitution of La and Yb into the bismuth-based structure (Bi 2 MoO 6 ) as shown in Figure 1. The main structure modification that was looked at is the substitution of La and Yb atoms for the Bi atoms in various combinations of the 8 possible positions occupied in the unit cell expressed in Figure 1 for atom position numbers 3, 16,17,18,19,20,21,22. Thus all combination of Bi and Yb, Bi and La are explored in an attempt to express a method for machine learning that can reduce the computational cost of exploring the substitution for a structure system. A usual method implemented for substitution is to calculate all possible combinations that the substitution can occur in the structure, which is a very computationally heave method. To put this in perspective, for the structure system explored in this study, more then 300 unique configurations are evaluated, this is a classical permutation problem. For example, if to say we are looking at 1 Yb atom and all 7 Bi atoms implemented into the structure, Yb could potentially occupy any of the 8 positions in the unit cell, thus 8 unique combinations can be expressed for the Yb 1/7 Bi ratio. Thus this study aims to express a viable solution to the permutation problem so that not all structures are needed to be explored in a system, and if structures are explored, to express a better production of the initial state (initial configuration of the unit cell) of the system to reduce computational cost. The way proposed in this study is to use machine learning technique to predict the fractional coordinates for all atoms that make up a unit cell (Figure 1), and predict the total ground state energy for structure system which expresses a more viable solution to screen large design spaces. The model for the machine learning approach is DFT simulations of atomic structures, and also the validated results are compared to the DFT simulations.

Input and Output Data Sets
The data set consists of inputs in a set of atomic charges (atomic numbers) that make up a structure while the output of the network consists of fractional coordinates along lattice vectors (a, b, c) of each atom position that make up one unit cell and the corresponding total energy. The reference/training [(Bi x La y ) 2 MoO 6 ] and validation [(Bi x Yb y ) 2 MoO 6 ] data is attained by consistent DFT calculations. It is noted that all calculated DFT structures initially start in the same configuration, Figure 1, which allows referencing input and output features consistent with the number that correspond to each atom position for all unit cells as expressed by numbers in Figure 1. The way the input and output features are expressed for the data set is by decomposing the input features (charges of atoms in the unit cell) and output features (atomic positions in the unit cell) as the vector representation expressed in Figure 2A. It is noted, that operations such as rotation, translation, and permutation to atoms in the data set (input and output) would ultimately change the ordering and positions of referenced atoms that make up the unit cell. Thus it is very important that the input and output features formed for the neural network reflect the same corresponding number in the sequence of vector input and output given to the neural network for each unit cell. This is why the order in which the numbers expressed in Figure 1 is maintained when classifying input and output vectors for each unit cell represented as the vectors in Figure 2A. Thus allowing the network to progressively see how individual output feature changes as a consequence of input vectors.
To reduce the complexity of solving the atomic positions and gain insight into fundamental and meaningful solutions to the problem given. This study imposes dimensionless quantities given to the neural network for data set interpretation, thus fractional coordinates are used for atomic positions. Understanding the energy as substitution occurs in the structure plays quite a central role in studies of chemical and biological systems. Thus the ground state energy of the relaxed structure became another very important feature the neural network was trained to calculate. Keeping the same idea of dimensionless quantity in mind, this study proposes to express the energy as a ratio of the calculated sum of the individual energy that make up the structure constituent parts divided by the total ground state energy, which is referred to as the energy ratio in this study.

Energy Ratio = ∑ Individual Energy o f Constituent Parts
Total Energy o f Structure (1) Thus the ground state energy is the DFT calculated energy and the individual energy of atoms that make up the system is just solving the self-consistent calculation for each individual atom. From this point foreword energy is expressed as energy ratio. Unlike the output vector expressed for each input vector for fractional coordinates, the energy ratio is just a scalier representation of energy for each input vector. Thus the total ground state energy for each unit cell can be calculated simply by using the sum of energies that make up the individual part of the unit cell divided by the energy ratio.

Artifical Neural Network (ANN) Model
This study aims to train a neural network with all configurations of (Bi x La y ) 2 MoO 6 , inputs being the atomic number (charge) of atoms that make up the unit cell structure and outputs being the fractional coordinates for each atom in the unit cell, then to validate this with configurations of (Bi x Yb y ) 2 MoO 6 . General ANN model is a simple information processing unit with multiple inputs and output. A simple neuron within the architecture of the network is attaining inputs from other neurons or from the exterior through path modeling. The artificial neuron output is computed as the weighted sum of all inputs modified by an activation function. These weights are adjusted through the learning process of the ANN. The idea is to train the network with a given set of non linear input to output data set, in order to express patterns in the data set. The method implemented in this study is that of Radial Basis Function (RBF) ANN, where the network is comprised of three layers (Figure 2B), the input layer (IL), the hidden layer (HL) and the output layer (OL) 22 .
Each layer in the RBF-NN has a different task 13,27 , the general architecture of the RBF is expressed in Figure 2B. From the IL to the HL of RBF network the distance between the network input and hidden layer centers is calculated. From the HL to OL the weighted sum is computed for each neuron. Each neuron of the HL has a vector parameter called center, and the general expression of the network is given as 13,27 : where, N is the number of neurons in the HL (i ∈ {1, 2, ..., N}), w i j are the weight of the i th neuron and j th output,ȳ j is the neural network's response to the j th output, and φ i is the radial basis function for the i th neuron. The radical basis function is refereed to as the activation function which is taken as the Gaussian function defined as, where, σ i is the spread parameter of the i th neuron, x is the input data vector, c i is the center vector of the i th neuron, and x − c i 2 is the Euclidean distance. Figure 2B expresses the architecture of the used RBF network. Note that for the structure system this study is exploring the input vector (x) is a 1x36 vector (Z 1 ,...,Z 36 ), the vector is passed to the HL comprised of N neurons. At the HL each neuron calculates the Euclidean distance between its center vector and the input vector, then calculates the Gaussian function (activation function) with the spread parameter specified (Equation 3). Then the output layer calculates the weighted sum for each input multiplied by the neurons activation response (Equation 2) to that input vector to get the output of the RBF network (ȳ 1 , ...,ȳ j ), where M denotes the number of outputs ( j ∈ {1, 2, ..., M}), being in this study the fractional coordinates.
The training part of the NN involves determining the number of neurons in the HL, and to attain the desired output for the network, the w, σ , and c parameters can be adjusted and attained. The most common error reference response of the NN typically used is mean square error and sum square error to train the NN. Thus for this study the error-based expression used (supervised learning) is that defined as: where, y j indicates the desired output, andȳ j is the network output. Thus ultimately the training steps for the network entails the minimization of the error function. The training algorithm for the NN utilizes three approaches. The first, is to use k-means clustering 9 for initially attaining the centers (c) for each neuron based on the input vectors. The second, is to have the weights (w) updated based on the activation function for the training set by using pseudo-inverse 26 of the activation function matrix used for all training sets. Lastly, using a gradient decent algorithm (GD) 17 to progressively update the spread for each neuron (σ ), thus ultimately minimizing the error function.

Training Algorithm used for the Model
To attain the desired result for the neural network output to the training parameters, the training algorithm follows three steps. The reader should note that this section is only referencing the training set in order to teach the NN model, and not for validation set. The first step used in the training algorithm is a k-means clustering 9 for input vectors of all training set to find the centers. The point of k-means clustering is to partition all the input vectors of training set into N clusters, N being the number of neurons specified, in which each input vector belongs to the nearest mean cluster, cluster is the center specified for the i th neuron. However if the user wants to increases or decreases the number of neurons for the model a new k-means clustering has to be done to attain the new center for each neuron.
Before the second step the spread (σ ) is initialized. Since the output of the HL multiplied by the weights is supposed to approximate the output training data as given in Equation 2 or represented as follows, where, O is the number of training data set's, N being the number of neurons, c being the center that corespondents to each neuron, j represents the output feature approximating for (such as a 1 in Figure 2B). For exampleȳ (1) O is the a 1 value for the O vector input (x O ) where asȳ (2) O is the b 1 for the O vector input, reference Figure 2B. A much simpler way to wright Equation 5 is as follows, For the second step the weights are computed by using the inverse of Φ ( j) if it is a square matrix, being the number of neurons is the same as the number of the tanning set (N=O), note that in this case the centers will just be the input vectors and the ANN will provide exact estimations for the training data, this is not recommended because the estimation will tend to over fit all desired inputs to the training data. However, if less neurons are used then will need to compute the pseudo-inverse 26 of Φ ( j) , which is done for most cases, thus the weights are expressed as: Its noted that Y ( j) is the desired output for the training data and thus the weights are updated this way for each updated spread of the Gaussian function (σ ). The final step in the training algorithm is updating the spread for each neuron (σ ) for the next iteration of the training algorithm to further minimize the error in the NN output to the desired output of the training set. The method used is gradient descent algorithm (GD) 17 which is a first order derivative based optimization algorithm used for finding local minimums for a function.
The method takes small steps proportional to the negative of the gradient for the error function at current iteration to update the spread to the next iteration. The error function is defined as follows, where, Y is the desired output of the training set,Ȳ is the NN output, and j is the output feature of the tanning set, same j represented in Equation 5-9. Thus the GD 17 algorithm is used to minimize the error and optimize the adjusting spread for each neuron by iteratively computing the partial derivative and updating spreads (σ ) in parallel and follows this update, where j is the output feature of the tanning set, i is the current iteration step, η is a small step size (referred to as the learning rate). Thus a simple summery of training is as follows: 1) Select the number of neurons and initialize the spread parameter for each neuron.
2) For the training set do k-means 9 algorithm to attain the center for each neuron.
3) Compute the activation for each neuron and get the weights for each neuron by pseudo-inverse approach 26 .

4)
Check to see convergence, if not converged do a GD 17 algorithm to update the spread for next iteration and repeat step 3.

1-?? | 5
The convergence criteria specified for this study is that the training epoch are iterated until the NN training output data set stops improving in accuracy for 100 epochs. The optimization process implemented, is to carry out the training method 5 times using an order of magnitude smaller learning rate each time. Its noted that w, σ , and c parameters very much relies on the number of neurons used for the network. Ultimately having more neurons than the needed amount causes the model to over fit the output to the training data set and increases the complexity of the network. Therefore the number of neurons used directly affects the performance of the network and has to be investigated based on desired result.
Once training is done approximating any set becomes relatively simple. First, the input for the desired approximation wanting to be evaluated is taken through each neuron by computing the Euclidean distance to each center of each neuron specified in the training set. Second the neuron computes the activation function (Equation 3) which describes the relationship of input to the center of that neuron. Then lastly the approximated output for the desired input is computed as the sum of all neurons activation function multiplied by the corresponding weights attained in the training set, thus solving Equation 6 for desired input and follows Figure 2B.

DFT Computational Details
Ground state properties and total ground state energies are approximated for relaxed configurations of a bismuth oxide structure system (see Figure 1), which assisted this study to implement a neural network that can analyze the relevant trends in data set by means of density functional theory (DFT) approach 7 . The DFT calculations used in this study implemented pseudo-wave function functional representation based on Perdew-Burke-Ernzerhof (PBE) exchange-correlation function. The PBE exchange-correlation function implemented at potentials with a cut-off energy wave function of 1496 eV (110 Ry) that demonstrated quite accurate and stable results for the studied unit cells of this paper. The various configurational combinations that the bismuth structure had and this study explored, benefited greatly in the reduction of computational expense due to the implementation of pseudized wave function. It is noted that all the DFT simulations for all configurations had the initial system to be that of bismuth based structure expressed in Figure 1. For all expressed DFT simulations, a Monkhorst-Pack with a k-point mesh sampled at 2x2x2 grid with 1/2,1/2,1/2 offset was implemented. To account for Van der Waals interaction, a Van der Waals correction term 1,8 is implemented in DFT simulations, however this correction did introduce some empiricism into the calculations.
For the dispersion interaction, a cut-off radius of 12 times that of cut-off wave function is implemented, 1320Å, with the scaling parameter specified as 0.7 for DFT simulations.
The unit cell for each configuration geometry were relaxed to a relative total energy less than 1x10 −10 with an overall unit cell pressure less than 0.5 kBar by computationally solving the electronic density self-consistently. It is noted that DFT predictions of energies, band bap, etc. are known to under predict results because of the exchange-correlation terms that are needed to calculate and due to the over-analyticity of the used functionals. Thus the calculated configurations and energies are used not as absolute but to express a relationship between input and output features that the network model will be trained on to predict the trends for various configurations.

Results and Discussion
Understanding the characteristics that describe chemical and biological system behavior and attributes plays quite a role in guiding and predicting new material systems. The modeling tools used to simulate these quantum systems relies on firstprinciple calculations of electronic structures to predict interactions and attributes. This study proposes to use a NN to predict and refine the initial state of structure configurations and ground state energy to reduce the computational expense needed for first-principle calculations based on DFT simulations. As mentioned in previous sections, we express a NN to a training set consisting of all configurations of (Bi x La y ) 2 MoO 6 , and test (validate) results with configurations of (Bi x Yb y ) 2 MoO 6 . The input to the NN is the vector input that describe the charges that make up a structure and output being the desired feature. The first feature the NN is trained to approximate are the fractional coordinates (atomic positions) that make up a unit cell based on the bismuth structure ( Figure 1). The second feature the NN predicts is the energy ratio that describes the total ground state energy of a unit cell.

Approximating Atomic Positions
Predicting atomic positions that make up a unit cell in a sense determines the computational time required to achieve convergence for that structure. In DFT, the total energy of any given system of interacting atoms and electrons is a function of the atomic positions that make up the structure and the electron density 3 . The external potential used in DFT explicitly depends on the atomic position, which is changed by a small step to find optimized atomic structure. Thus the Hamiltonian and the wavefunctions used in DFT are also dependent on the atomic positions. The initial routine of DFT code is to solve the charge density self-consistently, which is solved iteratively by computing the potential terms in the Hamiltonian by initial guess of the input density and comparing that to the output density attained by using Kohn-Sham approach. This iterative step is considered converged when the self-consistent energy is within a specified accuracy. After this self-consistent calculation is done the atomic positions are moved by a small step, then re-evaluation of the density and solving the selfconsistent calculation is done till the problem is solved within some accuracy specified. Thus its noted that the iterations requires to solve this problem relies heavily on the atomic positions that make up the structure. It is quite simple to deduce that if a structure system has atomic positions initially close to that of the final atomic positions then the iterations requires to solve DFT calculations would be reduced. In practice the initial atomic positions are specified by the structure system wanting to be tailored and or substituted, as in this study all DFT structure systems had initial atomic positions specified by Figure 1 which is the bismuth based structure. Figure 3 illustrates the response of the neural network er-ror for vector positions as Yb is substituted in the 8 atomic positions specified in Figure 1. Its noted that for each unit cell there is a total of 36 vectors that correspond to the 36 atomic positions specified in Figure 1. Thus the maximum vector error is plotted for all cases as Yb is substituted, Figure 3A, and the configuration that yielded the maximum overall error, Figure 3B. Initial inspection of Figure 3A shows a random error of the NN response. However, that is not the case when looking carefully, there seems to be a trend, as more Yb substitution occurs the error progressively builds up. Which is to be expected, as lower concentration of Yb tends to be dominated by Bi, thus the training set which has seen high concentration of Bi will yield good results for higher Bi combination with Yb. Where as higher concentration of Yb will result in error build up due to the training set not consisting of any Yb substitution. Thus vectors of low concentration Yb will be evaluated more accurately, and based on the results the concentration of 1 Yb 7 Bi resulted in the lowest maximum overall vector MSE as 1.0x10 −3 Figure 3A. However, for the maximum error response overall of the NN as expressed by '*' in Figure 3A and Figure 3B, the highest maximum vector MSE is 1.6x10 −2 of 3 Yb 5 Bi. Its noted for both the maximum and minimum MSE for all Yb substitution, the vectors are compared to that of DFT, thus the error is reasonable compared to the computational expense DFT takes to attain the final structure. The only downfall is that the NN needs a good interpretation for training set in order to accurately expresses input and output features.
Looking at the worst cause, the overall error in the NN response, as in Figure 3B, atomic vector 35 is contributing to the maximum overall error for the (Bi 5/8 Yb 3/8 ) 2 MoO 6 structure. Further inspection of the atomic position that corresponds to vector 35 in Figure 1 revels the location in the unit cell it occupies. This position is located at the boundary of the unit cell, at the boundary the unit cell tends to be periodic in nature. Thus from a NN perspective the training set had periodicity in the unit cell structures but was never told how this periodicity existed for these structure. The NN only saw how the atomic position would change from one input vector to the next, and so it was difficult for the NN to approximate the behavior of this structure at the boundary condition, a solution to this can be incorporating periodicity in the NN algorithm. However, for the intent of this study, it was more interesting to see if the NN could potentially predict the behavior at the boundary, and not by explicitly telling the network that behavior from user input. After further inspection of the tested set, most of the error was expressed at the boundary conditions, but was reasonably low. A visual interpretation of the error expressed for the vectors in Figure 3B of (Bi 5/8 Yb 3/8 ) 2 MoO 6 structure is expressed in Figure 4.
It is very interesting to see how the NN tries to predict the structure configurations of Yb substitution. Comparing the A B Fig. 4 Considering the worst case, structure A represents the DFT attained optimized structure and B is the NN approximated structure for (Bi 5/8 Yb 3/8 ) 2 MoO 6 . This configuration had the overall maximum error based on the NN approximations. Note that the individual vector MSE expressed in Figure 3B is comparing structure A (DFT) and B (NN). optimized DFT structure, Figure 4A, and the NN predicted structure, Figure 4B, we see that the NN tends to distribute the atoms more evenly in the unit cell. This goes back to the training set, where configurations containing La expressed very similar interactions as configurations containing Bi, were structures expressed more even distribution of atoms in the unit cell. It is noted that all configurations of La, including pure La in the 8 positions proved stable from a DFT convergence perspective. However, not all Yb configurations were stable, and from DFT simulations anything more the 50% would prove to not converge at all. Thus, this study tested (validate) Yb configurations to 50% substitution in the 8 possible positions in the unit cell from DFT simulations. We also see when comparing A and B of Figure 4 that the NN had progressively done a better job at predicting the atomic positions that are within the unit cell. Its noted that atoms closer to the boundary tended to be more evenly unformed for the NN, such as the 4 oxygen sites located at the top of the unit cell in Figure 4 A and B, this feature was expressed in configurations of the training set (La substitution) thus the network though it to be true for most cases. However, with this error in mind, this study found it to be quite difficult for the NN to attaining exactly the atomic positions, this is quite an arduous task that requires very precise calculations of interactions in a structure system. Thus the NN predicted structure is used as a refinement of the (Bi 5/8 Yb 3/8 ) 2 MoO 6 Initial Positions Specified Computational Time Figure 1 ∼30 hr Figure 4B ∼26 hr Table 1 The table expresses the computational time required for DFT simulation to find optimized structure, for the worst case with maximum error, for (Bi 5/8 Yb 3/8 ) 2 MoO 6 configuration starting from initial atomic positions specified by Figure 1 and NN positions specified by Figure 4B. All parameters for DFT calculations are expressed by DFT Computational Details section, the only thing that is changed are the initial atomic positions.
initial atomic positions in the aim of reducing computational time spent for DFT calculations which has proven from time to time to have very accurate results. The structure with the maximum overall error (worst case), (Bi 5/8 Yb 3/8 ) 2 MoO 6 , was re-evaluated using DFT simulation with the initial positions specified by the NN, Figure 4B, and evaluation of the computational time spent revels an overall improvement of 4 hours in wall time compared to the initial positions specified by Figure 1 for DFT simulation (Table 1). This improvement may seem trivial for a single structure, however when looking at the whole design space, this is quite an improvement. Overall the tested set consisted of 100 configurations and each of the configurations had required an average of 30 hours of continuous running time for each simulation. Thus the 4 hours would account for over 400 hour in reduction time in the overall study, if to say each structure would be improved on average of 4 hours. However, this reduction is based on the time spend for the worst case with maximum error, and in practicality most structures had much lower error, thus the improvement would be more then 400 hours overall. For instance, when testing the NN specified initial positions for the best case with the minimum error, being (Bi 7/8 Yb 1/8 ) 2 MoO 6 , the overall reduction in computational time is more then 9 hours, thus overall this study would account for way over 400 hours in computational wall time reduced. It is noted that most design studies of structure systems account for hundreds and even thousands of configurations, thus every possible method of improvement is explored, such as the expressed refinement of the initial system positions. Ultimately tailoring a structure system is a means to find a beater structure that expresses very specific behaviors. In structure simulations, first principle calculations can be used to predict reaction steps that proceed chemical reactions such as the hydrogen evolution reaction (HER) or oxygen evolution reaction (OER). These reactions steps relies on the computation of the total energy for the structure system, thus effectively simulating energy becomes a priority in this study.

Approximating Energy Ratio
Total energy of the unit cell is ultimately effected by the atomic positions that make up the unit cell. Thus it is possible to calculate the energy based on the predicted structure attained from the NN. However, this study intended for the NN to find relevant trends between input (charge) and the output features.
With the error found based on the atomic vectors, it was reasoned that ultimately back calculating the energy would ultimately have similar error. So ultimately energy became a feature for the NN to calculate and not attained by back calculating it relative to the NN predicted positions. Its noted that energy is given as a energy ratio to the NN, thus dimensionless and more manageable for this study (reference Input and Output Data Sets section). Having the NN to predict the energy and not back calculating it proved quite ideal as expressed by the error for Yb substitution in Figure 5. Being able to predict the energy without the computational expense really solves the permutation problem mentioned in the Model Design Space section. This is because not all configurations will require DFT simulation to evaluate the energy. In some cases its more ideal to take the structure system with the lowest total energy, such as the example stated in the Model Design Space section. For instance, finding the best ratio (configuration with lowest energy) of Yb 1/7 Bi ratio requires DFT simulations of 8 unique configurations. Thus by using the NN to evaluate the energy one can get a seance of the relevant trend as Yb is configured in different positions in a unit cell. Figure 5A expresses the overall energy ratio MSE comparing the NN and DFT calculated energy for Yb substitution in the 8 positions. Its very interesting to see that the error expressed in Figure 5A behaves like a step increase from low to high concentration of Yb. This is due to the evaluation of energy by the NN for similar configuration ratios like the 1/7 Yb to Bi. Thus the energy for 1/7 Yb to Bi in one configuration is quite similar to a different configuration with the same ratio of 1/7 Yb to Bi. However, the error expressed for energy is of great accuracy, minimum error being 1.1x10 −8 MSE for low concentration of Yb (12.5%) and maximum error of 2.3x10 −7 MSE for high concentration of Yb in the 8 positions (50%). Figure 5B expresses the quality of fit for the evaluated energy ratio expressed by the error in Figure 5A, where comparison of the reference (DFT) versus predictions (NN) of energy for the training and tested set. The ideal line (x=y) in Figure 5B indicates the quality for expressed predictions, thus predictions matching reference data (training set) would lie on the ideal line. It is quite reasonable to see that the energy attained from the NN is way within the energy computed by DFT simulations (ideal line in Figure 5B), without the computational expense required to evaluate each structure configuration. In Figure 5B, the red points indicate how well the NN responded to the training set, and expresses a visual interpretation of the NN performance for the training algorithm steps taken. The points indicated in blue represent the tested set, being all the configurations of Yb up to 50% substitution in the 8 positions.
Ideally the network should be within the line indicated as Ideal in Figure 5B, this means that the energy calculated from DFT is exactly the same as the predicted NN energy (E DFT =E NN ). Thus having the network express energy based on training set proved quite ideal, and more so having the energy expressed as a ratio between the sum of the individual energy that makes up the structure constituent parts divided by the total ground state energy even more ideal. Thus attain energy is feasible as long as the individual energy of the constituent parts is evaluated.

Computational Efficiency of the Proposed NN Approach
The reader should note that most structure based studies encompass huge design spaces, on the order of hundreds and even thousands of unique configurations. Thus to relatively evaluate the effective efficiency of the proposed NN approach, a control model must be compared and evaluated for the overall wall time required. The control model for this study is all the configurations that make up the training and tested set. The training set is comprised of 230 unique configurations of (Bi x La y ) 2 MoO 6 and the tested set are 100 configurations of (Bi x Yb y ) 2 MoO 6 . Note that on average the wall time required for relaxation of each configuration took 30 hours of consistent DFT calculations. That is roughly 6900 hours of continuous calculations for DFT simulation for the tested set and 3000 hours for the training set. The total wall time required for the control model is mostly used to find the stable energy ratio configurations as expressed by Equation 1, and in practicality there are vast amounts of structure analyses needed to be done for most structure based studies. However, to demonstrate the potential of the NN approach proposed by this study, the overall wall time is reduced by refining the initial atomic positions and predicting the energies prior to DFT simulation for structure system. For the first method of reducing the overall wall time by refining the initial atomic positions as expressed on this chapter in Section 3.1. The overall wall time was shown to reduce by 4 hours for the structure with the maximum error (worst case) and 9 hours for the structure with the minimum error (best case). Thus since most structures expressed error way lower then the worst case as shown in Figure 3A for the tested set. The overall improvement in the wall time is way over 400 hours for the tested set, because this is taking the worst case improvement. The overall wall time reduction for the tested set was calculated to averaged around 7 hours for each trail. That is calculating the average wall time for all the trails of the tested set, which ranged from 4 hours (worst case) to 9 hours (best case) in wall time reduction. Based on the calculated average wall time reduced, an overall 700 hours improvement for the tested set is shown. Thus reducing the 3000 hours of overall wall time to 2300 hours, which is a 1.3x speed up in computation for the tested set.
For the second method, the overall wall time is reduced by predicting the energies for all the energy ratio trials for the tested set as expressed by Section 3.2. This way only the potential trial that represent a specific ratio can be evaluated as expressed in Section 2.1. Note that there is 100 trials for the tested set and the combination of Bi and Yb goes up to 50% substitution. That is to say only rations of Yb occupying 1/8, 2/8, 3/8, 4/8 for the 8 possible positions in the structure system as expressed in Figure 1. Thus by being able to predict the energies for the 100 trials prior to DFT simulation, vast amounts of trials can be eliminated and only trials that represent the specific trial ratio will need to be evaluated. That is only 4 trials of the 100 for the tested set because there is only 4 ratios being evaluated for in the Yb case (1/8, 2/8, 3/8, 4/8). Eliminating most of these trials greatly reduces the overall wall time of 3000 hours to 120 hours, which is a 25x speed up in computation for the tested set. Reference Table 2 for a summery of the overall wall time reduced by the NN approach.

Overall Feasibility
Results in previous sections demonstrate not only the NN feasibility in learning the expressed system, but can provide accurate and meaningful features that greatly reduced DFT simulations overall and DFT computational expense. It would be ideal to demonstrate the proposed method to other DFT simulations of structure systems, for understanding the limitations and feasibility of the proposed method. The two features that the network was able to predict is the atomic positions and the total energy of the unit cell. The first feature was that of fractional coordinates for atomic positions. Where the network  Table 2 The table expresses the computational wall time required for DFT simulation to find optimized structure for the overall configurations of the tested set by implementing the NN proposed approach. The first method is to refine the atomic positions prior to DFT simulation. The second method is to predict the energy ratio for all configurations of the tested set thus eliminating majority of the required trials for the tested set. Note that there are 100 unique configurations (trials) for the tested set and on average each trial required continuous wall time of 30 hours for DFT simulation.
saw various input vectors and the corresponding fractional coordinates as the output feature as described by the training algorithm. Outside of these features no additional assumptions are given to the network. The computed output features relied on the DFT optimized structure, where all assumptions such as boundary conditions are taken care of, thus the network saw pure input and output characteristics. Further investigations regarding the best set of input features that could express the uniqueness of the structure may need to be looked at in future studies. However the method expressed for this study, which are the inputs to the NN as charges that make up individual atoms in the unit cell, proves to express reasonable results, especially for the second feature the network was trained to predict (total energy). The second feature the NN was able to predict is the total energy, no modification to the NN architecture is done, except that the energy was the only output to each input vector. Advantages of the expressed method versus alternate NN methods is the simplicity and ease of the algorithm used. When it comes down to the algorithm, the most arduous part was attaining the data input and output feature. Thus the referenced and training sets were obtained by time consuming selfconsistent calculations using DFT. The method expressed in this study links charge classification directly to structural information and total energy. While this approach is quite powerful and efficient for specific case studies, it does require having access to accurate DFT simulations to explore design space. Thus limitations are based on the data set used and having to train the NN on multiple configurations in order for a relevant trend to be presented for the network, and thus ideal for large design studies. For instance the training algorithm consisted of approximately 230 training configurations of the (Bi x La y ) 2 MoO 6 and evaluated 100 configurations of (Bi x Yb y ) 2 MoO 6 for the tested set. It is also noted that introducing any Yb configuration in the training part would prove to reduce the overall error because the network would see configurations similar to the tested set. However, this study focused on the feasibility of predicting and refining initial positions on configurations that are not introduced to the NN. Thus this method is independent of material system and only depends on referencing DFT simulations for specific class of structure system like the bismuth bases structure (Bi 2 MoO 6 ) explored in this study. The expressed method is thus ideal for large design space based studies where structure systems could potently take multiple substantiational attributes, thus making such studies reliant on computationally effective methods.

Conclusion
This study proposed a neural network technique based on machine learning for evaluating large design spaces to reduce computational cost for DFT simulations. Predicting the optimized atomic structure and the total energy is ultimately the bases for most first principle calculations. Thus this study looks at bismuth bases system (Bi 2 MoO 6 ) as the example structure system and the design space consists of configurations in (Bi x M y ) 2 MoO 6 where (M = La, Yb). We have demonstrated that this method can provide accurate predictions that will ultimately reduce the computational cost required to explore the structure unique configurations by applying one of the two methods. The first method to reduce the computational cost is by having the neural network better predict and refine the initial structure system based on a given set of training data. This way the DFT simulations will have a more refined method to initialize the unit cell being evaluated, which will allow less iterative steps to be taken in evaluating the optimized configuration for that specific structure. The second method proposed, exposes the network to a training data set for a given structure system, in order for the total energy of the optimized configuration to be predicted. This way the total energy of any configuration that makes up a given structure system can be predicted accurately based on already attained training data set for that structure system.
The training set used for this study are configurations of (Bi x La y ) 2 MoO 6 , and the tested set used are of (Bi x Yb y ) 2 MoO 6 configurations. Both the training and the tested set are evaluated with DFT simulations in order to see the error response from the DFT simulations to the neural network. Ultimately for the bismuth based system explored for this study, an overall maximum vector mean squared error of 1.6x10 −2 was expressed by the neural network for structure configuration of (Bi 5/8 Yb 3/8 ) 2 MoO 6 being the worst case. Then the predicted structure for this configuration was used as the refined structure for initializing the unit cell for DFT simulation, and allowed an overall 4 hour reduction in computational time. For the tested set, this overall reduction would account for 400 hours in reduction if all structures had on average 4 hour improvement, however most structures demonstrated over 4 hours reduction when implementing the refined positions from the NN. Thus over 400 hours in reduction is expressed for the design study of the tested set. Also, the network was trained to predict the total energy for the same tested configurations, and yielded a maximum mean squared error of 2.3x10 −7 . The described method is independent of material system and depends only on referenced DFT simulations for a specific structure system. Thus ultimately the machine learning technique expressed for this study provides a promising starting point for high-throughput electronic structure predictions for DFT simulations in large design space based studies.