Rapid and Automated Alloy Design with Graph Neural Network-Powered LLM-Driven Multi-Agent Systems

Abstract

A multi-agent AI model is used to automate the discovery of new metallic alloys, integrating multimodal data and external knowledge including insights from physics via atomistic simulations. Our multi-agent system features three key components: (a) a suite of LLMs responsible for tasks such as reasoning and planning, (b) a group of AI agents with distinct roles and expertise that dynamically collaborate, and (c) a newly developed graph neural network (GNN) model for rapid retrieval of key physical properties. A set of LLM-driven AI agents collaborate to automate the exploration of the vast design space of MPEAs, guided by predictions from the GNN. We focus on the NbMoTa family of body-centered cubic (bcc) alloys, modeled using an ML-based interatomic potential, and target two key properties: the Peierls barrier and solute/screw dislocation interaction energy. Our GNN model accurately predicts these atomic-scale properties, providing a faster alternative to costly brute-force calculations and reducing the computational burden on multi-agent systems for physics retrieval. This AI system revolutionizes materials discovery by reducing reliance on human expertise and overcoming the limitations of direct all-atom simulations. By synergizing the predictive power of GNNs with the dynamic collaboration of LLM-based agents, the system autonomously navigates vast alloy design spaces, identifying trends in atomic-scale material properties and predicting macro-scale mechanical strength, as demonstrated by several computational experiments. This approach accelerates the discovery of advanced alloys and holds promise for broader applications in other complex systems, marking a significant step forward in automated materials design.

Keywords Multi-agent system, Large language model  $\cdot$ Deep learning  $\cdot$ Graph neural networks  $\cdot$ complex alloys  $\cdot$ Materials design  $\cdot$ Scientific machine learning

1 Introduction

Multi-principal element alloys (MPEAs) represent a relatively new and innovative class of materials comprising three or more elements and showing remarkable mechanical properties superior to their pure or dilute counterparts, such as mechanical strength, fracture toughness, ductility, and resistance to hydrogen embrittlement. [1, 2, 3, 4]. The BCC refractiry MPEAs of the family Cr-Mo-W-Nb-V-Ta-Ti-Zr-Hf have recently gained particular interest due to their exceptional high-temperaute strength retention, exceeding the capabilities of current superalloys. [5, 6, 7, 8, 9, 10, 11, 12] With millions of possible compositions within just one single crystal, these complex systems offer significant potential to achieve tailored properties for specialized applications. However, navigating the immense multicomponent HEA compositional space to guide alloy discovery with optimized or desired properties poses significant challenges. Developing mechanistic multi-scale theories that connect atomisitc-level phenomenon to microscopic material properties, such as temperature dependent yield sterss, have been instrumental in exploring this design space. [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]

While theoretical models hold significant potential for predicting the macroscopic properties of MPEAs and navigating their vast design space, a key limitation lies in the computational cost of obtaining the necessary input parameters, which typically rely on atomistic simulations. This challenge is particularly acute for body-centered cubic (BCC) materials, where plasticity is controlled by the motion of screw dislocations [24, 21]. In contrast to edge dislocation theories, which offer a simplified approach by using misfit volume as a primary parameter [13, 14, 25], no such simplification exists for screw dislocation theories. As a result, these models require the computation of critical parameters via expensive atomistic simulations. Two key quantities in screw dislocation theories are the Peierls barrier and the solute-screw interaction energy [24, 21]. The Peierls barrier represents the intrinsic lattice resistance to dislocation motion, while the solute-screw interaction energy describes the influence of solute atoms on the dislocation’s behavior. Both factors serve as energy barriers to dislocation motion, and their calculation is further complicated by random solute fluctuations around the dislocation. Atomistic simulations are commonly used to determine these parameters in alloys [26], but applying them to multi-component alloys presents two significant challenges: (a) the vast design space of these alloys and (b) the need for numerous realizations to obtain a statistically accurate average due to the random solute environments. To address these limitations, machine learning (ML) and deep learning (DL) models offer a promising alternative by streamlining the computation of these input parameters, potentially reducing the reliance on computationally expensive brute-force methods.

The advent of ML and DL methods have revolutionized the ways of materials design, physical modeling, and properties measurement [27, 28, 29, 30, 31, 32]. These methods can discover hidden patterns in the training data and thus have been merged into all disciplines. The use of DL methods within the field of crystalline materials is wide, from machine learning interatomic potential development [33, 34, 35, 36] to intrinsic properties calculation [38] to dynamic crack path prediction [39]. Among all the DL architectures, GNN models have been developed to deal with graph structures which model a series of objects (nodes) and the relationships between them (edges), making them an ideal solution for crystalline materials where the atoms are nodes and the edges represent the intermediate bond [40, 41].

While machine learning (ML) methods have accelerated the exploration of multi-principal element alloys (MPEAs), they often target specific material properties in isolation. This narrow focus can limit their ability to incorporate the broader, interdisciplinary knowledge that is critical for true breakthroughs in alloy design and scientific discovery. To overcome these limitations, multi-agent systems have emerged as a transformative approach, facilitating the integration of multi-modal data and external knowledge—such as new developments in physics and materials science including theoretical models—into the design process in a more holistic and adaptive way [42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52]. In multi-agent systems, a group of AI agents—each powered by Large Language Models (LLMs)—work together dynamically to solve complex, multi-faceted problems. By assigning distinct roles to each agent via targeted prompts, the system can divide and conquer different aspects of the alloy design challenge.

In this paper, we present an LLM-based multi-agent system that leverages the power of: (a) a suite of advanced Large Language Models (LLMs), responsible for fundamental tasks such as planning, reasoning, and decision-making; (b) a group of specialized agents, each playing distinct roles within the system; and (c) a set of external tools for various tasks, including a newly developed Graph Neural Network (GNN) model to predict atomic-level material properties. This multi-agent model builds upon AtomAgents, a multimodal multi-agent system with advanced simulation capabilities for alloy design and discovery. AtomAgents was capable of retrieving new physics from atomistic simulations, solving complex multi-faceted alloy discovery problems. However, as the complexity of the alloy systems increased, the computational cost of atomistic simulations became prohibitive.

To address this, we developed a GNN model that bypasses the limitations of brute-force calculations by providing rapid physics predictions. Although trained on a relatively small portion of the compositional space of the ternary Nb-Mo-Ta family of body-centered cubic (BCC) alloys, the GNN model accurately predicts fundamental quantities, including the Peierls barrier and solute/screw interaction energy. Several examples are provided to demonstrate the proficiency of this multi-agent approach in not only exploring the design space but also tackling complex alloy design challenges, such as predicting macroscopic yield stress.

The plan of the paper is as follows. We start by providing details about our GNN model to predict atomic-level material properties in Section 2. We then discuss our proposed LLM-driven multi-agent system for automated mterials design, guided by the predictions from our GNN model. We outline the main components of the multi-agent approach and provide several examples to illustrate the proficiency of our model in solving complex alloy design problems. We then present the key findings in Section 3 and discuss the implications of our multi-agent system for future research in materials discovery.

2 Results and Discussion

The main tool implemented in our multi-agent system is the GNN model that accurately predicts the fundamental material properties in multi-component alloys, i.e. Peierls barrier and potential energy change due to the screw dislocation movement. We first delve into the basics of development of the GNN model and then discuss our proposed LLM-based multi-agent system to design alloys solicitating physics data from the GNN model.

2.1 GNN model

The workflow of the GNN approach for predicting the Peierls barrier and energy changes due to dislocation movement is schematically presented in Figure 1. Initially, the dislocated structures for 200 samples per composition within the training compositional space are generated, as shown in Figure S1 in the SI. After minimizing these structures, the potential energy change is calculated by subtracting the final potential energy from the initial one. Additionally, the Peierls barrier for each sample is extracted from the minimum energy path connecting the initial and final states, as derived from NEB simulations. The potential energy changes and the Peierls barriers serve as the ground truth or labels used to train the deep learning (DL) model. In this model, each random configuration is represented as a graph, where nodes correspond to atoms and edges to chemical bonds. Node features consist of chemical (solute types) and configurational features (screw dislocation displacement). The solute types at each site are one-hot encoded, with vectors $[1,0,0]$, $[0,1,0]$, and $[0,0,1]$ representing Nb, Mo, and Ta, respectively. Additionally, the $z$-component of the screw dislocation displacement, $\delta_{z}$, is calculated for pure Mo and included as a node feature. Notably, only the chemical features vary across compositions, while the $\delta_{z}$ feature remains constant, thus eliminating the need for atomic relaxation in new random configurations during inference. Furthermore, edge features are constructed based on the solute types of the neighboring atoms. After constructing the graph input, the GNN model is trained to predict either the potential energy change or the Peierls barrier, which serves as the label for the entire graph. Once trained, the model can predict the potential energy change and Peierls barrier for new structures, drastically reducing the need for time-consuming atomistic simulations, particularly NEB calculations. This enables us to explore the full compositional space of NbMoTa ternary alloys to identify candidates with enhanced properties, such as a higher Peierls barrier.

Peierls barrier and potential energy change predictions are graph regression tasks. In this work, we adopt a GNN with the Principal Neighborhood Aggregation graph convolution operator (PNAConv),[53] which has outperformed many popular GNN models in the literature, such as GCN,[54] GAT,[55], and GIN [56] on benchmark graph regression tasks. The PNA model combines multiple aggregators, which affect how messages between nodes are passed, with a degree-scaler that generalizes the sum aggregator. Specifically, the architecture employs four different aggregators: mean, maximum, minimum, and standard deviation. The improved performance of the GNN is attributed to this combination, where the degree-based scalers amplify or attenuate signals in the network based on node degree.

More details about the construction of the dataset, the architecture of the GNN model, and the training procedure are provided in the Methods section. In the following sections, we evaluate the model’s performance in predicting material properties on the test dataset and assess its ability to generalize to unseen configurations.

Figure 1: (a) Overview of the workflow used in this work to train GNN models for an end-to-end prediction of the Peierls barrier and potential energy changes. The process begins with the generation of initial and final dislocation structures in random BCC multi-component alloys, which are then minimized to compute the initial and final potential energies. Next, NEB simulations are performed to determine the minimum energy path and calculate the Peierls barrier. The input to the machine learning model is the graph representation of the alloys, with nodes encoding spatial and chemical information, and edges representing bond types. Two supervised GNN models are trained to predict graph-level labels: the Peierls barrier and potential energy change. (b) GNN architecture. The graph input is first passed into an input block to enlarge the dimension of node features. The node embedding, edge features and connectivity of input graphs are then input to the message passing block where the information from the neighbors of each node in the graph is aggregated to update the hidden features of the nodes. The output of the message passing block is then input into a global pooling later which outputs a graph-level embedding by adding node embeddings across the node, and it is connected to a multilayer perceptron (MLP) that returns a predicted graph label, i.e. Peierls barrier or potential energy change.

2.1.1 Peierls barrier

Random fluctuations in the solute environment surrounding the screw dislocation play a crucial role in determining the Peierls barrier, introducing a degree of complexity that must be captured for accurate predictions. In this section, we assess how well our GNN model predicts the Peierls barrier for various random compositions within the test set. To ensure the model’s generalizability, the compositions in the test set are strategically distributed to cover a broad range of the compositional space, including the equimolar composition $\text{Nb}_{33}\text{Mo}_{33}\text{Ta}_{33}$ and binary systems. This comprehensive distribution allows us to verify that the model can accurately capture both the complexity of ternary alloys and the simpler binary ones, which is critical for extending the model to new, unseen compositions. Moreover, each composition consists of 50 realizations of random solute configurations, ensuring a robust evaluation of the model’s predictive capabilities across different atomic arrangements.

We evaluate the model’s performance by comparing the predicted Peierls barrier values with the ground truth data, as illustrated in Figure 2(a). The model achieves a low mean absolute error of 37 meV, indicating that it has effectively captured the fundamental physical mechanisms driving dislocation motion. This accuracy, despite the inherent randomness in the solute distribution, demonstrates the model’s ability to reliably predict the Peierls barrier across varying atomic environments. This result highlights the model’s capacity to grasp complex atomic interactions that influence the Peierls barrier.

From these predictions, we can compute the mean energy barrier for each composition, a key parameter in solute strengthening theories, which rely on these barriers to predict material strength. Figure 2(b) compares the mean Peierls barrier values predicted by the DL model against the ground truth values, showing a strong agreement across both ternary and binary alloys. This alignment highlights the model’s ability to predict individual configurations while also capturing broader trends in solute-strengthening behavior, offering valuable insights into how solute distributions influence material properties across various compositions.

In summary, the model’s low error and consistent performance across diverse alloy compositions affirm its potential as a powerful tool for predicting material behavior, significantly reducing the need for computationally expensive atomistic simulations such as NEB calculations.

Figure 2: Evaluation of the GNN model for Peierls barrier prediction on the test set. Test set contains new compositions that never appeared in the training set. (a) Comparison of machine learning (ML) predictions and NEB results for the Peierls energy barrier. (b) Comparison of the ML results and NEB results for the mean energy barriers of binary and ternary compositions in the test set. The correlation coefficient $R^{2}$ between the predictions and the ground truth is 0.97.

2.1.2 Potential energy change predictions

Another important quantity is the potential energy change as the dislocation moves from one energy minimum to the next. After training the GNN using the potential energy changes as graph labels, we compared the model’s predictions with the ground truth values for the graphs in the test set (see Fig. S1), as shown in Figure 3(a). The relatively low mean absolute error of 60 meV demonstrates the strong performance of our DL model in accurately predicting potential energy changes.

The random fluctuations in the solute environment lead to a distribution of energy changes. In random alloys, the average of this distribution is expected to be zero, and according to solute strengthening theories, the energy barrier against dislocation motion scales with the standard deviation of this distribution [26, 24]. This can be expressed as:

$$\Delta\tilde{E}_{p}(a)=\left(\frac{b}{\zeta}\right)^{\frac{1}{2}}\sigma_{\Delta U}$$ (1)

where $\zeta$ is the dislocation segment length (with $\zeta=5b$ in our case), and $b$ is the screw dislocation Burgers vector.

Figure 3(b) shows the values of the solute/screw interaction energy parameter computed by our DL model compared to the ground truth values, demonstrating good agreement for both ternary and binary alloys. This enables the rapid prediction of $\Delta\tilde{E}_{p}(a)$, a fundamental input parameter in solute-strengthening theories, providing critical insights for predicting material strength in multi-component alloys.

Figure 3: Evaluation of the GNN model for the potential energy change prediction on the test set. Test set contains new compositions that never appeared in the training set. (a) Comparison of ML predictions and NEB results for the potential energy change. (b) Comparison of the ML results and NEB results for the solute/screw interaction energy parameter, Eq. 1 of binary and ternary compositions in the test set. The correlation coefficient $R^{2}$ between the predictions and the ground truth is 0.89.

2.2 Multi-modal multi-agent model

Figure 4 outlines our multi-agent model for rapid multi-component alloy design and analysis. Our model builds upon AtomAgents [52], a multi-modal multi-agent system designed to extract new physical insights from atomistic simulations. In this version, AtomAgents has been enhanced with our new graph neural network model that predicts fundamental physics data. This advancement offers an efficient way to explore the massive compositional space of multi-component alloys compared to costly atomistic simulations, allowing for solving more challenging problems and rapid design of alloys with enhanced or desired material properties.

Figure 4: Overview of the GNN-powered, LLM-based multi-agent system developed here. The system consists of a human and an AI assistant agent at its core, where the human poses queries, and the AI assistant provides responses, seamlessly steering the problem-solving process with the help of integrated tools. These tools are responsible for various tasks, including planning, coding, and multi-modal analysis, and each incorporates a set of AI agents that dynamically collaborate to solve complex tasks. A key component is the physics tool, which includes newly developed graph neural network (GNN) models to retrieve essential physical parameters (such as the Peierls barrier and potential energy change) as well as physics-based theories (like solute-strengthening theory). The GNN models enable the rapid prediction of fundamental material properties, bypassing the need for costly atomistic simulations. The iterative collaboration between agents within the tools and the seamless interaction between the human and AI assistant allows for efficient resolution of complex materials design challenges.

In the current multi-agent framework, the core agents are the "User" and the "AI Assistant," with the User posing the queries and the AI Assistant responding by utilizing various tools. These tools are designed to handle a wide range of tasks, from planning and coding to multi-modal data analysis. Additionally, these tools engage a set of autonomous agents that collaborate to respond effectively to the queries provided by the Assistant. A crucial component of this toolset is the physics tools, which empower the model by leveraging both newly developed Graph Neural Network (GNN) models and physics-based theoretical frameworks, such as solute-strengthening theories. These theories are essential for connecting atomistic-level features to macroscopic properties like yield stress. Each agent in our multi-agent system serves a specialized role, defined by a unique profile, and is powered by advanced general-purpose language models from the GPT family, accessed via the OpenAI API. A detailed description of each agent’s profile and the implementation details of the AI agents and tools are provided in Materials and Methods.

The tools integrated within our multi-agent system enable the resolution of complex challenges in multi-component alloy design, allowing for automated exploration of the vast compositional design space to identify candidates with superior or targeted performance. These tools are automatically activated by the AI Assistant to perform various tasks throughout the alloy design and analysis process. A general workflow of our model is demonstrated in Figure 5. When the User submits a query, the Assistant calls the planning tool to generate a comprehensive plan, detailing the necessary steps, identifying the required tools, and specifying their input parameters. This planning process is driven by two key agents: the planner and the reviewer. The planner drafts the initial version of the plan, while the reviewer critiques and refines it through continuous iterations. This collaborative process ensures that the final plan is precise and thoroughly detailed, which is crucial as the entire problem-solving process depends on this plan.

Once the plan is finalized, the AI Assistant executes it, typically involving tasks such as predicting material properties using the physics tools, collecting and organizing data, and passing them to the coding tool for visualization and plotting. Finally, the Assistant engages the multi-modal agent to conduct a comprehensive analysis of the results.

Figure 5: Overview of a typical workflow executed by our multi-agent system for accelerated and automated multi-component alloy design and analysis. Upon receiving a query from the user, the process begins with the assistant agent, which calls the planning tool to generate a detailed plan. The assistant agent follows this plan by invoking the relevant tools and providing the necessary input functions. These tools are integrated into our system to equip the agent with advanced capabilities, including material property predictions through physics-based retrieval tools (GNN models), data visualization via code execution (through a code-executor pair of agents), and multi-modal analysis (through a multi-modal agent-user pair) for enhanced visualization and reasoning. The system incorporates a human-in-the-loop design, fostering dynamic collaboration between humans and AI. This enables continuous refinement of ideas and queries, resulting in a more adaptive and efficient problem-solving process.

The advantage of incorporating our newly developed GNN model into the multi-agent system is twofold: (a) it enables the rapid exploration of vast design spaces, providing deeper insights into material behavior across a wide range of compositions, and (b) it significantly accelerates the materials design process by allowing the fast computation of critical input parameters for physics-based theoretical models, such as solute-strengthening theories. To demonstrate the effectiveness of our multi-agent system in these key areas, we conducted two experiments, each detailed below.

Variation of the Peierls barrier in binary and ternary alloys

In this experiment, the multi-agent system was tasked with exploring the Peierls barrier across the ternary compositional space, as illustrated in Figure 6(a). Upon receiving the query, the assistant agent activates the planning tool to devise a detailed plan for the task. The embedded planner and reviewer agents collaborate, with the planner drafting the initial plan and the reviewer providing feedback. This feedback loop iterates twice to ensure the accuracy and robustness of the plan. The resulting comprehensive plan outlines the following key steps: (a) Composition Generation: creating a list of all possible binary and ternary compositions in the Nb-Mo-Ta system at 5% intervals (231 total compositions); (b) Computation of Peierls Barrier: using the GNN model to compute the Peierls barrier for all generated compositions; (c) Plotting Results: utilizing the coding tool to plot the computed Peierls barrier values on a ternary diagram representing the Nb-Mo-Ta compositional space; and (d) Data Analysis: applying the analysis tool to examine the plotted data for trends, correlations, and significant insights. The plan also includes a detailed list of functions and input parameters to ensure consistency and precision during execution.

The Assistant agent then executes the plan, generating the possible compositions in the ternary space and activating the relevant tools, as shown in Figure 6(b). We observe that the GNN model enables rapid predictions of the Peierls barrier across the entire compositional space (231 compositions at 5% intervals). Notably, calculating the mean Peierls barrier for a single composition typically requires computationally expensive NEB simulations over numerous random configurations to capture the statistical nature of the problem. The GNN-based approach, however, completes these calculations within seconds per composition, demonstrating significant efficiency. After computing the Peierls barriers, the coding tool is called to write a Python function for plotting the values on a ternary diagram. Despite the complexity of ternary plots, the code generator, powered by the o1-mini LLM, produces the plot, as shown in Figure 6(c) , depicting the variation of the Peierls barrier across Nb-Mo-Ta ternary and binary compositions.

The results of the multi-agent system provide valuable insights into the variation of the Peierls barrier across the entire compositional space. Regions with high Peierls barrier values are identified around compositions with intermediate levels of Nb and Mo and low levels of Ta. As the final step of the workflow, the multi-modal agent performs detailed analysis, identifying the highest Peierls barrier in compositions with a high percentage of Mo and the lowest in those with a high percentage of Ta. Additionally, the agent conducts a comparative analysis of the influence of Nb, Mo, and Ta concentrations on the Peierls barrier.

However, some inconsistencies are observed in the agent’s analysis. For example, the model suggests that the Peierls barrier decreases with increasing Nb content, which is inaccurate based on the predictions. This highlights challenges in the multi-agent system’s ability to detect patterns in complex ternary plots, requiring further experimentation and refinement. The complete results of the ternary plot analysis are presented in Figure S5 in the Supplementary Material.

The full spectrum of Peierls barrier values across the ternary compositional space offers many opportunities for further insights. For instance, it allows exploration of the barrier’s variation within binary systems or specific ternary systems. The human-in-the-loop capability of the multi-agent system makes it possible to issue follow-up queries. Here, we instruct the model to plot the variation of the Peierls barrier with solute concentration for specific binary and ternary systems, including the Nb-Mo and Nb-Ta binary systems, as well as the (NbMo)_2xTa_1-2x and (NbTa)_2xMo_1-2x ternary systems, as shown in Figure 6(e).

The assistant agent identifies the relevant compositions from prior results and passes their Peierls barrier values to the coding agent for plotting, generating the plots shown in Figure 6(e). These plots reveal a non-linear variation of the Peierls barrier with solute concentration. Key insights include: (a) In the Nb-Mo system, the Peierls barrier increases with Nb concentration, peaking around 50%, before decreasing as Nb content approaches 100%. (b) In the Nb-Ta system, the Peierls barrier also increases with Nb concentration, but at a slower rate. (c) In both ternary systems, the Peierls barrier shows an increasing trend as $x$ rises from 0 to 0.5.

The entire workflow is repeated for the solute/screw interaction energy parameter, $\Delta\tilde{E}{p}$, a key material property derived from potential energy changes. The problem-solving approach remains identical, but the model now uses the potential energy GNN to compute $\Delta\tilde{E}{p}$. The results are visualized in the ternary plot of Figure 6(d), showing non-linear variations of the interaction energy parameter with solute concentrations.

Further insights are gained by examining the variation of $\Delta\tilde{E}{p}$ for selected binary and ternary systems, as shown in Figure 6(g). The key findings include: (a) In the Nb-Mo system, $\Delta\tilde{E}{p}$ increases with Nb concentration, peaking around 30%, and then gradually decreases as Nb content continues to rise. This suggests that the interaction energy parameter is optimized at intermediate Nb concentrations. (b) In the Nb-Ta system, $\Delta\tilde{E}{p}$ is higher at high Nb concentrations compared to the Nb-Mo system, indicating stronger interactions at these concentrations. (c) In the (NbTa)_2xMo_1-2x system, $\Delta\tilde{E}{p}$ increases steadily with $x$ and decreases as Mo content becomes minimal. (d) In the (NbMo)_2xTa_1-2x system, $\Delta\tilde{E}_{p}$ rises with increasing $x$, eventually stabilizing, suggesting a saturation point in interaction energy as Ta content diminishes.

These results highlight the complex, non-linear interactions between solute concentrations and the screw dislocation, which influence both the Peierls barrier and the solute/screw interaction energy parameter. The GNN-powered multi-agent system offers an effective framework for efficiently exploring the vast compositional space of multi-component alloys. It enables the identification of key trends and intricate behaviors, significantly enhancing the materials design process through faster and more intelligent exploration.

Figure 6: Overview of the multi-agent collaboration to explore the variation of the Peierls barrier and solute/screw interaction energy parameter in binary and ternary alloys. (a) User query requesting exploration of the Peierls barrier across a wide range of compositions, with a similar query repeated for the solute/screw interaction energy parameter. (b) Workflow detailing the computations and analyses performed by the multi-agent system. (c) and (d) Ternary plots displaying the Peierls barrier and the solute/screw interaction energy parameter, respectively, across the compositional space generated by the multi-agent system. (e) Follow-up tasks requested by the user, involving additional data plotting and analysis. (f) and (g) Plots illustrating the variation of the Peierls barrier and $\Delta\tilde{E_{p}}$, respectively, with solute concentration for binary and ternary alloys, generated by the multi-agent model in response to the follow-up tasks.

Yields stress in multi-component BCC alloys

Solute-strengthening theories focus on the atomic-level mechanisms governing dislocation motion and their interactions with solutes, aiming to develop mechanistic models that predict temperature-dependent mechanical strength, specifically the yield stress, of materials. These models typically rely on a number of input parameters derived from first-principles calculations such as Density Functional Theory (DFT) or atomistic simulations based on empirical potentials. A prominent theory for screw dislocation strengthening in BCC alloys is the Maresca-Curtin theory [24], which requires inputs like the Peierls potential, the solute/screw interaction energy parameter, the kink formation energy, and vacancy/interstitial formation energies. However, due to the complex energy landscapes in multi-component systems, accurately computing these parameters poses a significant challenge. The newly developed GNN model addresses this by enabling the rapid and accurate computation of two critical parameters: the Peierls potential and the solute-dislocation interaction energy. This capability allows us to evaluate these parameters in mere seconds for each composition in the NbMoTa alloy space, offering a major breakthrough in efficiently predicting yield stress in multi-component alloys as a function of temperature.

The objective of this experiment is to demonstrate the potential of our multi-agent system in automating and significantly accelerating the materials design process by seamlessly integrating physics-based theoretical models with advanced DL-driven material prediction tools. In this task, the multi-agent system is assigned the calculation of yield stress for a series of alloy compositions over a broad range of temperatures, as shown in Figure 7(a). Upon receiving the task from the User, the dynamic collaboration between agents is initiated to accomplish the goal, as illustrated in Figure 7(b).

The process begins with the planning tool, which generates a detailed step-by-step plan outlining the required functions and their input parameters. The plan includes the use of GNN-powered tools to predict the Peierls barrier and the solute-dislocation interaction energy—key parameters for yield stress prediction. Additionally, for other input parameters such as the lattice constant, kink formation energy, and vacancy/interstitial formation energies, the system computes their values by averaging over the pure elements. Once all necessary input parameters are obtained, they are fed into the yield stress prediction tool for final calculations. The results are then passed to the coding tool, which generates Python code to plot the predictions. Finally, the multi-modal agent, powered by the GPT-4o LLM, analyzes the resulting plots to evaluate trends across compositions.

Figure 7: Overview of the multi-agent collaboration to predict the yield stress in binary and ternary alloy BCC systems. (a) The input task involving computing the yield stress for a set of compositions. Some material properties for the pure metals are provided as well as the strain rates. Moreover, the experimental data for binary alloys are provided to compare with the predictions. (b) The workflow of computations performed by our multi-agent system featuring key tasks like planning, tool calling for material property predictions from the GNN model, and finally plotting and analyzing the results. (c) Predictions of our model for the ternary alloys and (d) the predictions of our model for binary alloys along with the experimetal results.

Despite the complexity of the task, which requires numerous calculations and input parameters, we observe that all sub-tasks are executed seamlessly. The primary outcomes of the multi-agent collaboration are presented in Figures 7(c) and (d), which show the variation of yield stress with temperature for binary and ternary systems, respectively. The binary plot also includes experimental data from the input task for comparison.

Following this, the multi-modal agent conducts a comprehensive analysis of the results, as detailed in Figures S6 and S7 in the Supplementary Information (SI), covering binary and ternary alloys, respectively. Upon completing the analysis, the assistant provides a summary of key insights, which is presented in Figure 8.

The analysis indicates a typical decrease in alloy strength with increasing temperature, attributed to enhanced atomic mobility at higher temperatures. Additionally, the AI agent highlights that the Maresca-Curtin model predictions align with experimental data for some alloys, such as Nb81Mo19, but diverge for others, like Nb95Mo5 and Nb75Mo25. The agent suggests that these discrepancies may result from variations in alloy composition, microstructure, temperature effects, or experimental variability. To address these issues, the agent proposes refining the model to incorporate more detailed microstructural and compositional factors for improved accuracy. Furthermore, the AI agent emphasizes the role of solid solution strengthening, with lattice distortions impeding dislocation motion as a key strengthening mechanism.

This multi-agent collaboration exemplifies the rigorous integration of multiple modalities, combining AI-driven physics predictions with established theoretical frameworks, along with plotting and visualization reasoning. This multifaceted approach not only enhances the precision and depth of the materials design process but also represents a transformative pathway that could revolutionize how materials design is approached and executed.

Figure 8: A summary of the analysis performed by the multi-modal agent of the yield stress versus temperature predictions as shown in Figures 7c and d.

3 Summary and future perspective

First, we developed a GNN model that directly links dislocated crystalline structures, caused by screw dislocations, to material properties in multi-component BCC alloy systems. Despite being trained on a small subset of the massive compositional space in ternary alloys, the model achieves high accuracy in predicting key properties. These include the potential energy change between consecutive screw dislocations and the Peierls barrier, which represents the highest energy at the transition state. One of the main applications of this energy change is the calculation of the solute/screw interaction energy parameter, a fundamental quantity in solute strengthening theories, which our model captures with high precision. The model’s predictions are generated in a matter of seconds, representing a significant improvement over conventional atomistic simulations, which can take days or even months when applied to the large design space of multi-component alloys. This speed, combined with minimal loss of accuracy, establishes our approach as a feasible alternative to costly simulations for investigating the effects of chemical fluctuations around defects on crystalline material properties. Furthermore, this method can be easily extended to high-entropy alloys and other crystalline structures, such as FCC and HCP systems, as well as to different types of defects, including edge dislocations. Our GNN model enables rapid exploration of the design space, paving the way for the development of alloys with enhanced mechanical properties, particularly those with higher strength.

Next, we constructed a multi-modal artificial intelligence (AI) system that integrates three core components: (a) Large Language Models (LLMs), which excel at a wide range of tasks, including multi-modal reasoning, strategic planning, rational thinking, and even coding; (b) AI agents, each powered by LLMs and external tools, designed with specialized roles and expertise, working collaboratively in a dynamic environment and autonomously tackling multi-faceted problems; and (c) a Graph Neural Network (GNN) model that provides rapid predictions of fundamental material properties. This multi-agent system is an enhanced version of our previous model, AtomAgents [52], which relied on direct atomistic simulations for physics-based insights. In the modified version, the GNN model is integrated to mitigate the high computational cost associated with atomistic simulations, thereby enabling faster exploration of the vast compositional design space inherent to multi-component alloys. This new framework offers a comprehensive method for addressing more complex problems in the design and analysis of multi-component alloys. Its capabilities are demonstrated through a series of experiments, including the exploration of the Peierls barrier and solute/screw interaction energy across the entire ternary composition space, as well as predictions of yield stress in BCC alloys, which are then validated against experimental data.

Our automated multi-agent system incorporates several key features that make it highly effective for complex materials science challenges. First, it is powered by advanced large language models (LLMs) like o1-preview and o1-mini, which provide exceptional capabilities in reasoning, strategic planning, and coding. These models excel at handling multi-modal data, evaluating emerging trends, and integrating information from diverse sources, making them essential for solving multi-scale problems. Second, the system supports the seamless integration of various data modalities, incorporating theoretical models, sophisticated deep learning algorithms, and even generative tools into the design and analysis process. This enables our system to stay at the cutting-edge of scientific and technological advancements, provideing a more holistic perspective on materials design. Third, the system is designed to support continuous evaluation and refinement, allowing for iterative problem-solving where errors are identified and corrected through mutual feedback between AI agents, ensuring that the system remains robust and adaptive throughout the problem-solving process. Finally, the human-in-the-loop functionality embedded within the system further enhances its flexibility and usability. This feature allows users to interact with the system directly, providing feedback, refining model parameters, or introducing additional requests. This AI-human collaboration is particularly valuable for complex, multi-faceted problems, ensuring that the system remains adaptable and responsive to specific user needs. Together, these advancements position the multi-agant modeling of materials as a powerful and innovative tool for tackling some of the most challenging problems in materials science and engineering, traditionally reserved for human experts.

Furthermore, the multi-agent system’s high degree of adaptability makes it suitable for application across a broad range of domains and disciplines, extending its potential impact beyond materials science. One of the primary aspects of this adaptability lies in the underlying LLM, which orchestrates the system’s internal operations and interactions. The performance of the entire system is heavily influenced by the capabilities of this LLM, and as LLM technology continues to evolve, we anticipate even greater improvements in system performance. For example, the transition from earlier versions to o1-preview resulted in significant gains in efficiency and accuracy. This ongoing advancement in LLMs suggests that future iterations of the system could become even more capable, especially as LLMs are further customized to specific tasks and integrated into different agents. Additionally, the system’s flexibility extends to its ability to incorporate a wide range of deep learning and generative models, each developed for various disciplines. This capability allows the system to address domain-specific challenges while maintaining a unified problem-solving framework. Another important feature of this system is its capacity to integrate with other multi-agent systems, such as SciAgents, [42, 57] which specialize in generating advanced research hypotheses. This creates a cohesive scientific discovery process, where one system is responsible for exploring the hidden space of ideas, generating novel hypotheses, while another system is equipped to test, validate, or refine these ideas using advanced tools. Together, this interconnected approach provides a comprehensive and rigorous pathway for scientific discovery, making the multi-agent system an adaptable and powerful framework for tackling interdisciplinary problems.

4 Materials and Methods

Dataset generation

Random compositions of ternary NbMoTa refractory alloys with screw dislcoation are generated to compute the potential energy and Peierls barrier as follows. First pristine crystall strcuctures are generated by creating a rectangular simulation cell using the corresponding lattice constant and then randomly assigning solutes based on the desired concentrations. The simulation cell is oriented with dislocation glide direction x||[1,1,$\bar{2}$], glide plane normal direction, y||[$\bar{1}$10], and dislocation line direction, z||[111] with periodic boundary conditions along x and z and free surface along y. Atomic positions are relaxed by using a combination of the FIRE algorithm [97] and relaxation of the cell dimensions until the convergence is achieved-the norm of the force vector fell below 10e-6 eV/$\AA$ and stresses $\sigma_{xx}$ , $\sigma_{xz}$ , and $\sigma_{zz}$ fell below 0.1 MPa. All the atomistic simulations are performed using LAMMPS [58] at zero temperature. A machine learning MTP potential is used to describe the interatomic interactions in the NbMoTa alloys [59, 34].

We then introduce a screw dislocation in the center of each relaxed simulation cell using the PAD method [16] and minimize the energy along with the relaxation of the pressures. This serves as our initial dislocated configuration. To generate the final dislocated configurations, we use the same initial pristine random structure and insert the dislocation at the adjacent position (next Peierls valley) relative to the initial dislocation at the distance of a along the glide direction where $a$ is the Peierls valley distance. The potential energy changes is computed by subtracting the total potential energies of final and initial configurations. We then perform nudged elastic band (NEB) computations [60, 61, 62, 63] as implemented in LAMMPS [58] on these initial and final structures to compute the minimum energy path among these two energy states. The maximum value along this curve is stored as the corresponding Peierls barrier. The entire process is repeated for all the compositions in the training set (200 realizations for each composition) and test set (50 realizations for each composition) as depicted in Figure S1.

Graph representation

The results from atomistic simulations of atomic configurations and alloy properties are represented as graphs, which are used as inputs and labels to train the GNN. The graph structure is constructed from the pristine pure Mo configuration, where each Mo atom is a node in the graph, and the connectivity between nodes is determined by the distance between atoms. We consider a cylindrical region with a radius of $r_{c}=16\AA$, centered at the screw dislocation core, and only include atoms within this region as graph nodes. Within this graph, two atoms are connected if their distance is below the cutoff distance of 2.8 $\AA$, which is the distance between nearest neighbors in a perfect crystal of molybdenum. The graph representation consists of 705 nodes and 2610 edges.

The node features contain information about chemicals (solute type) and structural defect (screw dislocation). The structural defect feature is constructed from the $z$-component of the displacement derived from the atomistic simulation in pure Mo. This feature remains identical for all the atomic configurations in the dataset, eleminating the need for atomic relaxation for new configurations during inference. The chemical feature of the node is uniqely determined for each composition represented as the solute types converted to one-hot encoded representations ([1,0,0], [0,1,0], and [0,0,1] for Nb, Mo, Ta, respectively). Moreover, the edge features representing the bond type are encoded as one-hot representations based on the neighboring nodes as shown in Table 1

Table 1: Edge feature representation for bond type connecting neighboring nodes, $i$ and $j$.

node type i	node type j	edge feature
Nb	Nb	[1,0,0,0,0,0]
Mo	Mo	[0,1,0,0,0,0]
Ta	Ta	[0,0,1,0,0,0]
Nb (Mo)	Mo (Ta)	[0,0,0,1,0,0]
Nb (Ta)	Ta (Nb)	[0,0,0,0,1,0]
Mo (Ta)	Ta (Mo)	[0,0,0,0,0,1]

Graph neural networks (GNNs)

The GNN model is developed based on the deep-learning framework PyTorch [64] and its geometric extension library PyTorch Geometric [65]. The GNN architecture is depicted in Figure 1(b). First of all, the input graph is sent to the input block which uses a combination of PNA convolutional (PNAConv) layer, Gated Recurrent Unit cell (GRUCell), and Batch Normalization layer (BatchNorm) to upscale the dimension of node features to 50. Then the graph is passed to the message passing block which contains 10-time repetition of the combined layers. Within the block, nodes communicate with each other by passing the messages given the node features plus connectivity and update their own node features considering the received messages. The last MLP that returns the predicted Peierls barrier or potential energy change has a structure of an input layer of size 30, a hidden layer of size 20, another hidden layer of size 10, and an output layer of one neuron. ReLU is adopted as the activation function in this MLP. We adopt default weight and bias initialization of all of the layers in the model defined by PyTorch Geometric.

Model training and evaluation

All the datasets created for the training/validation compositions (see Figure S1) are split into the train set (90% data) and validation set (10% data). The models were trained with a batch size of 32 using the Adam optimization method [66] for 250 epochs on one NVIDIA Tesla V100s with 32GB memory. Training starts with a learning rate of 0.0005, and a dynamic learning rate scheduler named ReduceLROnPlateau reduces the learning rate by half if no improvement is seen for 10 epochs to minimize the validation MSE. The learning curves of the GNNs are shown in Fig. S4 in the SI which indicate the convergence of training for both models.

Solute-strengthening theory in BCC alloys

We apply Maresca-Curtin screw strengthening theory for non-dilute to high-entropy alloys to compute the yields stress [24]. The theory is established based on the assumption that the initially straight dislocation becomes spontaneously kinked at zero load and zero temperature so as to lower their total energy. Three mechanisms contribute to the screw dislocation strengthening; (I) Peierls-like mechanism, (II) kink glide mechanism, and (III) cross-kink formation and unpinning and the alloy strength at temperature $T$ is expressed by

$$\tau(\dot{\epsilon},T)=\tau_{xk}(\dot{\epsilon},T)+\text{min}[\tau_{k}(\dot{\epsilon},T),\tau_{p}(\dot{\epsilon},T)]$$ (2)

where $\dot{\epsilon}$ is the experimental strain rate and $\tau_{p}$, $\tau_{k}$, and $\tau_{xk}$ are the Peierls strength, kink migration strength, and cross-kink unpinning strength, respectively. [24, 21]

Agent design

We design AI agents using the state-of-the-art all-purpose LLM GPT-4 and dynamic multi-agent collaboration is implemented in AutoGen framework[67], an open-source ecosystem for agent-based AI modeling. Additional agents are introduced as described below.

In our multi-agent system, the human user agent is constructed using UserProxyAgent class from Autogen, and Assistant, Planner, Reviewer, coder agents are created via AssistantAgent class from Autogen, while multi-modal agent is constructed via MultimodalConversableAgent class. Each agent is assigned a role through a profile description as shown in Figures S8-S12 in the Supplementary Material, included as system_message at their creation.

Function and tool design

All the tools implemented in this work are defined as python functions. Each function is characterized by a name, a description, and input properties. The full list of tools and their descriptions can be found in the corresponding codes.

Data and code availability

All data and codes are available on GitHub at https://github.com/lamm-mit/AlloyAgents . Alternatively, they will be provided by the corresponding author based on reasonable request.

Author Contributions: M.J.B and A.G. conceived the overall concept. A.G and M.J.B developed the GNN model and multi-agent system. A.G. curated the training and testing data for the GNN model, performed the tests for various problems, analyzed the results and prepared the first draft of the paper. M.J.B supported the analysis, revised and finalized the paper with A.G.

Supplementary Materials

Additional materials are provided as Supplementary Materials.

Supplementary Materials

Rapid and Automated Alloy Design with Graph Neural Network-Powered LLM-Driven Multi-Agent Systems

Figure S1: Ternary diagram showing compositions for training and validation (red) and new unseen compositions for testing (blue).

Figure S2: Evaluation of the GNN model for potential energy change prediction on the validation set.

Figure S3: Evaluation of the GNN model for Peierls barrier prediction on the validation set.

Figure S4: Loss evolution during training. Learning curves over 250 epochs for the models trained on (left) potential energy changes and (right) Peierls barrier.

Figure S5: Results of the analysis performed by the multi-modal agent powered by GPT-4o for the plot shown in Figure 6(c)

Figure S6: Results of the analysis conducted by the multi-modal agent powered by GPT-4o for the plot shown in Figure 7(c).

Figure S7: Results of the analysis conducted by the multi-modal agent powered by GPT-4o for the plot shown in Figure 7(d).

Figure S8: Profile of the Assistant AI agent.

Figure S9: Profile of the planner AI agent.

Figure S10: Profile of the reviewer AI agent.

Figure S11: Profile of the coder AI agent.

Figure S12: Profile of the multi-modal AI agent.