Learning physics-consistent particle interactions

Abstract Interacting particle systems play a key role in science and engineering. Access to the governing particle interaction law is fundamental for a complete understanding of such systems. However, the inherent system complexity keeps the particle interaction hidden in many cases. Machine learning methods have the potential to learn the behavior of interacting particle systems by combining experiments with data analysis methods. However, most existing algorithms focus on learning the kinetics at the particle level. Learning pairwise interaction, e.g., pairwise force or pairwise potential energy, remains an open challenge. Here, we propose an algorithm that adapts the Graph Networks framework, which contains an edge part to learn the pairwise interaction and a node part to model the dynamics at particle level. Different from existing approaches that use neural networks in both parts, we design a deterministic operator in the node part that allows to precisely infer the pairwise interactions that are consistent with underlying physical laws by only being trained to predict the particle acceleration. We test the proposed methodology on multiple datasets and demonstrate that it achieves superior performance in inferring correctly the pairwise interactions while also being consistent with the underlying physics on all the datasets. While the previously proposed approaches are able to be applied as simulators, they fail to infer physically consistent particle interactions that satisfy Newton’s laws. Moreover, the proposed physics-induced graph network for particle interaction also outperforms the other baseline models in terms of generalization ability to larger systems and robustness to significant levels of noise. The developed methodology can support a better understanding and discovery of the underlying particle interaction laws, and hence, guide the design of materials with targeted properties.


Introduction
Interacting particle systems play a key role in nature and engineering as they govern planetary motion (1), mass movement processes (2) such as landslides and debris flow, bulk material packaging (3), magnetic particle transport for biomaterials (4), and many more. Since the macroscopic behavior of such particle systems is the result of interactions between individual particles, knowing the governing interaction law is required to better understand, model, and predict the kinetic behavior of these systems. Particle interactions are determined by a combination of various factors including contact, friction, electrostatic charge, gravity, and chemical interaction, which affect the particles at various scales. The inherent complexity of particle systems inhibits the study of the underlying interaction law. Hence, they remain largely unknown and particle systems are mostly studied in a stochastic framework or with simulations based on simplistic laws.
Recent efforts on developing machine learning (ML) methods for the discovery of particle interaction laws have shown great potential in overcoming these challenges (5)(6)(7)(8)(9)(10)(11). These ML methods, such as the Graph Network-based Simulators (GNS) (12) for simulating physical processes, Dynamics Extraction From cryo-em Map (DEFMap) (13) for learning atomic fluctuations in proteins, the SchNet (14,15) which can learn the molecular energy, and the neural relational inference model (NRI) (16) developed for inferring heterogeneous interactions, can be applied on various types of interacting particle systems such as water particles, sand, and plastically deformable particles. They allow implicit and explicit learning of the mechanical behavior of particle systems without prior assumptions and simplifications of the underlying mechanisms. A commonly applied approach is to predict directly the kinetics of the particles without explicitly modeling the interactions (12,14,(17)(18)(19)(20)(21). The neural networks, then, map directly the input states to the particle acceleration, occasionally by virtue of macroscopic potential energy (12,14,(17)(18)(19). While these approaches give an accurate prediction of the particle system as it evolves, they do neither provide any knowledge about the fundamental laws governing the particle interactions nor are they able to extract the particle interactions precisely.
Recent work (22) proposed an explicit model for the topology of particle systems, which imposes a strong inductive bias and, hence, provides access to the individual pairwise particle interactions. Reference (22) demonstrated that their Graph Network (GN) framework predicts well the kinetics of the particles system. However, as we will show, the inferred particle interactions violate Newtonian laws of motion, such as the action-reaction property, which states that two particles exert the same but opposed force onto each other. Therefore, the extracted pairwise particle interactions do not correspond to the real underlying particle interaction force or potential, which are the fundamental properties of a physical system. The origin of these discrepancies lies in the design of the GN approach, which does not sufficiently constrain the output space, and clearly demonstrates the need for a physicsconsistent Graph Neural Network (GNN) framework for particle interactions.
Besides interacting particle systems, recent works applied GN models to learn microlevel dynamics in other physical systems. For example, the work of Moon et al. (23) integrated domain knowledge into the GNN to learn drug-target binding affinities in a supervised way. Ha and Jeong (24) trained the GN model as a simulator for interacting discrete physical systems other than interacting particles, including cellular automata, the Vicsek model, the active Ornstein-Uhlenbeck particle model, and the movement of a flock of birds. Another follow-up research line aimed to improve GNNs to learn particle dynamics in more complicated scenarios. For example, the work of Huang et al. (25) learns particle dynamics in interacting particle systems with geometrical constraints where some particles are connected by bounds. However, to the best of our knowledge, no other previous work aimed to reveal physics-consistent particle interactions from observed particles trajectories without any direct supervision.
Here, we propose a GNN framework that incorporates universal physical laws, specifically Newton's second law of motion, to learn the interaction potential and force of any physical particle system. The proposed algorithm, termed physics-induced graph network for particle interaction (PIG'N'PI), combines the GNN methodology with deterministic physics operators to guarantee physics consistency when indirectly inferring the particle interaction forces and potential energy (Fig. 1). We will show that PIG'N'PI learns the pairwise particle potential and force by only being trained to predict the particle acceleration (without providing any supervision on the particle interactions). Moreover, we will show that the inferred interactions by PIG'N'PI are physically consistent (contrary to those inferred by purely data-driven approaches). We will further demonstrate that predictions provided by PIG'N'PI are more accurate, generalize better to larger systems, and are more robust to noise than those provided by purely data-driven graph network approaches. Moreover, we will demonstrate on a case study that is close to real applications that the proposed algorithm is scalable to large systems and is applicable to any type of particle interaction laws.

PIG'N'PI: physics-induced graph network for particle interaction
We propose a framework that is able to infer pairwise particle forces or potential energy by simply observing particle motion in time and space. In order to provide physics-consistent results, a key requirement is that the learnt particle interactions need to satisfy Newtonian dynamics. One of the main challenges in developing such a learning algorithm is that only information on the particle position in time and space along with particle properties (e.g., charge and mass) can be used for training the algorithm and no ground truth information on the interactions is available since it is very difficult to measure it in real systems.
The proposed framework comprises the following distinctive elements ( Fig. 1): (1) a graph network with a strong inductive bias representing the particles, their properties, and their pairwise interactions; and (2) physics-consistent kinetics imposed by a combination of a neural network for learning the edge function and a deterministic physics operator for computing the node function within the graph network architecture. In addition, the proposed framework consists of the two steps: (1) training the network to predict the particle motion in time and space; and (2) extraction of the pairwise forces or the pairwise potential energy from the edge functions of the trained network.
Particle systems: We consider particle systems that are moving in space and time and are subject to Newtonian dynamics without any external forces. A particle system in this research is represented by the directed graph G = (V, E), where nodes V = {v 1 , v 2 , …, v |V| } correspond to the particles and the directed edges E = {e ij : v i , v j ∈ V, i = j} correspond to their interactions. The graph is fully connected if all particles interact with each other. Each particle i, represented by a node v i , is characterized by its time-invariant properties, such as charge q i and mass m i , and time-dependent properties such as its position r t i and its velocityṙ t i . We use η t i to denote the features of particle i at time t, η t i = [r t i ,ṙ t i , q i , m i ]. We limit our evaluations to particle systems comprising homogeneous particle types. This results in particles exhibiting only one type of interaction with all its neighboring particles, leading to |E| = |V|(|V| − 1). We further assume that the position r t i of each particle i is observed at each time step t and that this information is available during training. Based on the position information r t i , velocityṙ t i , and accelerationr t i are computed. Proposed framework: The proposed PIG'N'PI framework extends the general GN framework proposed by (26), which is a generalized form of message-passing GNNs. The architecture of the proposed framework is illustrated in Fig. 2. We use a directed graph to represent the interacting particle system where nodes correspond to the particles and edges correspond to their interactions. The framework imposes a strong inductive bias and enables to learn the position-invariant interactions across the entire particle system network. Given the particle graph structure, the input is then defined by the node features η i (representing particle's characteristics). The target output is the accelerationr t i of each node at time step t. The standard GNs block (26), typically, comprises two neural networks: an edge neural networkĜ E (·; θ E ) and a node neural net-workĜ V (·; θ V ), where θ E and θ V are the trainable parameters. Here, we propose to substitute the node neural networkĜ V (·; θ V ) by a deterministic node operator G V ( · ) to ensure that the learned particle interactions are consistent with the underlying physical laws. The main novelty compared to the standard GN framework is that we impose known basic physical laws to ensure that the inferred pairwise force or potential energy corresponds to the real force or potential energy while only being trained on predicting the Figure 1 Framework of the proposed model to learn pairwise force or pairwise potential energy. (A) The interacting particle system contains three particles that evolve over time. At every time step, each particle is described by multiple features, which include position, velocity, charge, and mass (represented by the bar). Position and velocity evolve with time whereas other properties remain constant. (B to C) The proposed method learns physics-consistent pairwise force or pairwise potential at every time step t. The model has two components: the edge part (B) and the node part (C). In the edge part (B), two nodes' vectors are concatenated as edge feature (process 1 ). An edge neural networkĜ E (·; θ E ) (θ E represents the trainable parameters) takes the edge feature as input (process 2 ) and outputs a learnt vector on that edge representing the pairwise force or potential energy. In the node part (C), the output vectors by the edge neural network and the raw node feature are aggregated on each node (process 3 ). We design the deterministic node operator G V ( · ) by incorporating physics knowledge to derive the net acceleration on nodes (process 4 ). By minimizing the loss on node-level accelerations, the edge neural networkĜ E (·; θ E ) will output pairwise force or potential energy exactly. The workflow where the edge neural networkĜ E (·; θ E ) takes edge features as input. The corresponding output message M i j is the predicted pairwise force or potential energy, depending on the physics operator (B) or (C) in the node part. Parameters θ E inĜ E (·; θ E ) are trained by minimizing the loss on particle acceleration. acceleration of the particles. Whether to use the force-based PIG'N'PI (blue color in Fig. 2) or the potential-based PIG'N'PI (purple color in Fig. 2) depends on different demands or downstream applications in practice, i.e., we can use force-based PIG'N'PI if we are only interested in the pairwise force without caring about the pairwise potential energy; or we can use potential-based PIG'N'PI if we want to know the pairwise potential. Furthermore, note that the force-based PIG'N'PI method is not limited to energy-conserving systems, but can also be extended to dissipative systems (e.g., friction).
It is important to emphasize that only information on particle positions is used for training the algorithm and the ground-truth information on the forces and the potential energy is not available during training. For each edge, the property vectors η i of two nodes connected by an edge are concatenated as the edge feature vector. The edge neural networkĜ E (·; θ E ) outputs a message on every edge that corresponds to the pairwise force or potential energy. The output dimension is set to be the same as the spatial dimension d (two or three) ifĜ E (·; θ E ) is targeted to learn the pairwise force or one to learn the pairwise potential energy. Edge messages are aggregated on nodes and the node part operator computes the output corresponding to the acceleration of nodes, imposing physics-consistency on edge messages. Trained to predict the node-level acceleration, once applied to a new particle system, the GN predicts the particle motion at consecutive time steps. The pairwise forces or the pairwise potential energy can then be extracted from the edge function for each time step.
Contributions of the present work compared to precious research: Here, we propose a methodology to learn the pairwise force or pairwise potential energy from the observed particle trajectories. This focus distinguishes our work from many previous works such as (14,(27)(28)(29)(30)(31) that learn the energy of the system and then derive the per-particle dynamics from the global energy. Moreover, as outlined above, our proposed approach does not have access to any ground truth information during training but rather learns to infer the force and potential energy indirectly. This is contrary to the previously proposed approaches that rely on such information (17,21).
While our proposed framework has several similarities with two previously proposed frameworks that are also aiming to infer pairwise force and pairwise potential energy using also only the particle accelerations for training (22), none of the previously proposed methods is able to infer the underlying particle interactions that are consistent with the underlying physical laws. We demonstrate in our experiments that the learnt particle interactions of the previously proposed approaches are not consistent with the underlying physical laws and do not correspond to the real forces or potential energy.
In fact, the proposed algorithm has also similarities to the Physics-informed neural networks (PINNs) (32), which aim to solve partial differential equations. Both PINNs and PIG'N'PI integrate known physical laws. While PIG'N'PI integrates Newton's second law, PINNs enforce the structure imposed by partial differential equation at a finite set of collocation points.

Performance evaluation metrics
We evaluate the performance of the proposed PIG'N'PI framework on synthetic data generated from two-(d = 2) and three-(d = 3) dimensional numerical simulations. The key distinctive property of the generated datasets is the definition of the inter-particle potential energy P, which defines the inter-particle pairwise force by F = −∂P/∂r. The selected cases, which have also been used in prior work (22) and can be considered as a benchmark case study, cover a wide range of particle interaction features, including dependence on particle properties, e.g., mass and charge, dependence on interaction properties, e.g., stiffness, and varying degrees of smoothness (see Table 1 and Fig. S2 in SI Appendix J for visualization).
The method developed by ref. (22), which applies multilayer perceptrons [MLPs (33)] in the edge and node part, serves as the baseline for comparison. We do not change the architecture of the baseline except for changing the output dimension of its edge part MLPs when learning the pairwise force or potential energy. The output dimension is a d-dimensional vector for learning the pairwise force and a one-dimensional scalar for the potential energy. Besides the baseline, we compare the performance of PIG'N'PI to Table 1. The force and potential energy equations for different datasets, where F i j is the force from particle j to particle i, P ij is the potential incurred by particle j on particle i, r ij is the Euclidean distance between particle i and particle j, n i j is the unit vector pointing from particle i to particle j, and q i and m i are the electric charge and mass of particle i. k, L, c, and are constants.

Dataset
Pairwise force an alternative method proposed by ref. (34) that is also based on GN and was specifically designed to infer pairwise forces. We denote this method as GN+ and the details are introduced in the "Details of the method to learn pairwise force introduced by (34)" section.
We split each dataset into training, validation, and testing datasets and use T train , T valid , and T test to indicate the corresponding simulation time steps for these different splits. Details regarding the dataset generation are provided in the "Details about simulation and experiments" section. The baseline algorithm and PIG'N'PI are trained and evaluated on the same training and testing datasets from simulations with an 8-particle system. Further, the evaluation of the generalization ability uses a 12-particle system.
It should be noted that ref. (22) measures the quality of the learnt forces by quantifying the linear correlation between each dimension of the learnt edge message and all dimensions of the ground-truth pairwise force. This is a necessary but not sufficient condition to claim the correspondence of the learnt edge message with the pairwise interactions and to evaluate the performance of the indirect inference of the pairwise interactions. Instead, we evaluate the proposed methodology with a focus on two key aspects: (1) supervised learning performance, and (2) consistency with underlying physics. For all the evaluations, the mean absolute error on the testing dataset of various particle and interaction properties is used and is defined as follows and respectively, where the superscript hat indicates the predicted values. Here,φ i j and φ ij are the predicted and corresponding groundtruth, respectively, of a physical quantity between particle j and particle i (e.g., pairwise force), andφ i = jφi j and φ i = j φ ij are the aggregated prediction and the corresponding ground-truth, respectively, on particle i (e.g., net force). l 1 (x, y) computes the sum of absolute differences between each element in x and y, l 1 (x, y) = i |x i − y i |, if x and y are vectors or the absolute difference, l 1 (x, y) = |x − y|, if x and y are scalars. Hence, MAE part measures the averaged error of the physical quantity on particles over T test , and MAE inter is the averaged error of the inter-particle physical quantity over T test . The supervised learning performance is evaluated on the prediction of the acceleration MAE acc = MAE part (r,r). The true acceleration values serve as target values during training. The physical consistency is evaluated on two criteria. First, we evaluate the ability of the proposed framework to infer the underlying physical quantities that were not used as target during training (e.g., pairwise force), and second, we evaluate physical consistency by verifying whether Newton's action-reaction property is satisfied.
The following metrics are used to evaluate the consistency with the true pairwise interaction. For pairwise force, we use MAE ef = MAE inter (F, F); and for potential energy case, we evaluate the increment in potential energy MAE ep = MAE inter (P −P 0 , P − P 0 ), where superscript 0 refers to the initial configuration.
For the second part of the evaluation of the physical consistency, we verify whether Newton's action-reaction property is satisfied. For that, we evaluate the symmetry in either inter-particle forces with MAE F symm = 1

Performance evaluation between PIG'N'PI, GN+, and baseline for pairwise force
First, we analyze PIG'N'PI for application on particle systems with interactions given by pairwise forces. We start with evaluating the supervised learning performance by evaluating the prediction of the acceleration using MAE acc . The results show that PIG'N'PI provides slightly better predictions than the GN+ and the To verify the physical consistency, we first evaluate if the implicitly inferred pairwise forces are consistent with the true physical quantity. PIG'N'PI is able to infer the force field around a particle correctly, while the baseline model fails to predict the force field (see Fig. 3A for the spring dataset). A force field needs to be precise in both amplitude and direction. The error of the magnitude (see Fig. 3B) and angle (see Fig. 3C) demonstrate unambiguously the superior performance of PIG'N'PI compared to the baseline model. We quantitatively summarize the performance of the pairwise force inference with MAE ef , which shows PIG'N'PI is able to infer the pairwise force correctly, while both GN+ and the baseline fail to infer the pairwise force [two to three orders of magnitude worse inference performance for the spring dataset (see Fig. 3D) and all other datasets (see Fig. 4A and SI Appendix B)].
Secondly, we verify the consistency of the implicitly inferred pairwise forces with Newton's action-reaction law by evaluating the symmetry of the inter-particle forces with MAE F symm . Our results demonstrate that PIG'N'PI satisfies the symmetry property. However, GN+ and the baseline model are not able to satisfy the underlying Newton's laws for the spring dataset (see Fig. 3D) and the other datasets (see Fig. 4B and SI Appendix B).
Furthermore, we test the robustness of PIG'N'PI and the baseline model to learn from noisy data. We impose noise to the measured positions and then compute the noisy velocities (firstorder derivative of position) and noisy accelerations (second-order derivative of position). The noisy accelerations serve then as the target values for the learning tasks of all the models. The performance of PIG'N'PI decreases with increasing noise level (see SI Appendix I; Tables S11, S12, and S13). This is to be expected given that adding noise makes the training target (particle accelerations) less similar to the uncorrupted target that is associated with particle interactions. However, PIG'N'PI can still learn reasonably well the particle interactions despite the corrupted data. The performance of the baseline model fluctuates; however, with different noise levels significantly. This is due to the fact that the baseline model does not learn the particle interactions but rather the particle kinematics and is, therefore, more sensitive to noise.
Finally, we note that the proposed algorithm is also able to generalize well when trained on an 8-particle system and applied to a 12-particle system for all datasets (see Figs. 3D and 4, and Table S4).
Overall, the results demonstrate that the proposed algorithm learns correctly the pairwise force (that is consistent with the underlying physics) without any direct supervision, i.e., without access to the pairwise force in the first place, and that the inferred forces are consistent with the imposed underlying physical laws.

Performance evaluation between PIG'N'PI and baseline for pairwise potential energy
Besides learning the pairwise force, the proposed methodology is extended to learn the pairwise potential energy (see node part design for learning potential in the "PIG'N'PI: physics-induced graph network for particle interaction" section). In this case, the physics operator computes the pairwise force via partial derivative. Since GN+ was solely designed for learning the pairwise force, it is not possible to apply it to infer the pairwise potential. Therefore, for the task of pairwise potential energy inference, we only compare PIG'N'PI with the baseline model. Our results show that PIG'N'PI performs well in the supervised learning of the acceleration (MAE acc ). Here again, its performance is considerably better compared to the baseline model (see Fig. 5B). Moreover, the performance is similar to that in the force-based version of the algorithm. However, when comparing the performance of the baseline model on the supervised learning task between the potentialbased version and the force-based version of the model, the performance reduces significantly in the potential-based implementation (compare to Fig. 3D). This drop of performance is potentially explained by the adjustment of the output dimension of the edge neural network in the baseline model to enable the extraction of the potential energy.
Further, our results demonstrate a superior performance of the PIG'N'PI algorithm on consistency with underlying physics. Firstly, it infers well the increment of the potential energy (see Fig. 5A). It clearly infers the inrement of the potential field correctly. On the contrary, the baseline model fails to infer the potential energy. This can be quantitatively assessed with MAE ep . The results (see Figs. 5B and 6) show that PIG'N'PI is able to infer the potential energy, while the performance of the baseline model indicates that it is not able to learn the potential energy only from observing particles movements.
It is important to note that our algorithm cannot predict the absolute value of the potential energy; only the increment (see SI  Fig. S1). The reason is that the model is trained on the acceleration, which is computed from the derivative of the potential energy (i.e., the force). Hence, the model only constrains the derivative, and the constant of integration remains unknown. This limitation can be overcome by one of the two following options: either the potential energy is constrained by a spatial boundary condition or by an initial condition. In the former, we can impose a known value for a given value of r ij , e.g., we could use the assumption that the potential energy for a charge interaction approaches zero with increasing particle distance. The alternative (but less likely) approach consists in knowing the potential energy at a given time, e.g., at the beginning of the observation and add the inferred increment to the known initial value. Nevertheless, knowing the absolute value of the potential energy is, in fact, not crucial as only its derivative determines the dynamics of a particle system. This is also confirmed in our experiments by the accurate prediction of the acceleration (see MAE acc ).
Secondly, similar to the pairwise force prediction, PIG'N'PI also provides a superior performance on satisfying Newton's actionreaction property, while again, the baseline model fails to satisfy the underlying Newton's law. The performance is quantified by the symmetry of the inter-particle potential energies (MAE P symm ). (Figs. 5B and 6, and Table S3).
Finally, we test the generalization ability of the learning algorithms in a similar way as in the pairwise force case study. We apply the models trained on an 8-particles system to a new particle system comprising 12 particles. The results (see SI Appendix C; Table S5) show that PIG'N'PI predicts well the pairwise force and potential energy, and outperforms considerably the baseline model (see Figs. 5B and 6). This demonstrates that the PIG'N'PI model provides a general model for learning particle interactions.

Case study under more realistic conditions: learning pairwise interactions for an Lennard-Jones (LJ)-argon system
To evaluate the performance of the proposed framework on a more realistic system with a larger particle interaction system (to evaluate the scalability), we apply PIG'N'PI on a large LJ system.
We adopt the dataset introduced in ref. (35). This dataset simulates the movements of liguid argon atoms governed by the LJ potential. The LJ potential, which is given by V(r) = 4 {(σ /r) 1 2 − (σ /r) 6 }, is an extensively used governing law for two nonbonding atoms (36). The simulation contains 258 particles in a cubic box whose length is 27.27Å. The simulation is run at 100 K with periodic boundary conditions. The potential well depth is set to 0.238 kilocalorie/mole, the van der Waals radius σ is set to 3.4Å, and the interaction cutoff radius for the argon atoms is set to 3σ . The mass of argon atom is 39.9 dalton. The dataset is run for 10 independent simulations. Each simulation contains 1000 time steps with randomly initialized positions and velocities. The position, velocity, and acceleration of all particles are recorded at each time step. Figure 7 summarizes the learning pipeline. Contrary to the previous case study where for a small number of particles, a fully connected graph is considered, in this case study, we construct the graph of neighboring particles at every time step (Fig. 7). We connect the particles within the defined interaction cutoff radius while taking the periodic boundary conditions into consideration. Particles in the LJ-argon system are characterized by their position, velocity, and mass. The charge is not part of the particle properties, which is different from the particle systems considered in the previous case study (section "Performance evaluation between PIG'N'PI, GN+, and base-line for pairwise force" and section "Performance evaluation between PIG'N'PI and baseline for pairise potential energy"). Moreover, we compute the position difference under the periodic boundary condition and use it as an edge feature. This edge feature is required because the distance between two particles in this simulation does not correspond to the Euclidean distance in the real world due to the periodic boundary conditions. The node features and the edge features are then concatenated and are used as the input to the edge part of PIG'N'PI ( Fig. 7C) or the baseline model. Similarly to the previous case study, the learning target is the accelerations of the particles. The pairwise force and pairwise potential energy are then inferred from the intermediate output of edge part.
We evaluate the performance of PIG'N'PI on inferring pairwise interactions of the LJ-argon particle system with the same performance metrics as in the previous case study. The results are reported in Fig. 8 and Table S6. Because the particles in this dataset have the same mass, we also test a variant of GN+ such that we assign all nodes with a unique learnable scalar. We denote this variant as GN+ uni . The results confirm the very good performance of PIG'N'PI as observed in the previous case study. Generally, GN+ uni outperforms the GN+, but PIG'N'PI still surpasses GN+ uni and the baseline. On the one hand, PIG'N'PI performs better than the baseline, GN+, and GN+ uni on the supervised prediction task of predicting the acceleration [achieving less than half of the MAE acc compared to the baseline and the GN+ uni (Fig. 8A)]. On the other hand, PIG'N'PI is also able to infer the learn pairwise force correctly. Again, the baseline model is not able to  infer the pairwise force (PIG'N'PI outperforms the baseline by more than two orders of the magnitude on the MAE ef ). Moreover, particle interactions inferred by PIG'N'PI are consistent with Newton's action-reaction law (MAE F symm ). The bad performance of the baseline model indicates that the learnt interactions do not satisfy Newton's law.
To summarize, similar to the cases discussed in the "Performance evaluation between PIG'N'PI, GN+, and base-line for pairwise force" section, PIG'N'PI learns the pairwise force well without any direct supervision for this complex and large system. Besides, we test PIG'N'PI to learn the pairwise potential energy for this LJ system. Results are reported in Fig. 9 and Table S7. We first examine the MAE acc that is the learning target. The MAE acc of PIG'N'PI is similar to that in the force-based version of the algorithm. PIG'N'PI performs significantly better than the baseline model with more than two orders of magnitude (Fig. 9A). And, similar to the cases in the "Performance evaluation between PIG'N'PI and baseline for pairwise potential energy" section, we again observe the performance drop of the baseline model in this potential-based version with the force-based version. Then, we evaluate the MAE ep and MAE ef that are the two metrics for measuring the quality of the learnt pairwise potential energy. Results show that PIG'N'PI provides again superior performance on inferring the increment of the potential energy MAE ep (Fig. 9B). The bad performance of the baseline model clearly shows that it is not able to infer the potential energy correctly. Finally, we check Newton's action-reaction property in the potential energy by MAE P symm . Here again, the potential energy inferred by PIG'N'PI follows Newton's laws while the baseline model fails to infer the underlying interaction laws correctly (Fig. 9D). All evaluations demonstrate that the predicted pairwise potential energy by PIG'N'PI is consistent with the LJ potential used in the simulation, even though PIG'N'PI does not access the ground truth information.
The results on this case study demonstrate the scalability of PIG'N'PI to larger systems and the applicability to more realistic case studies. Moreover, the results confirm the results obtained in the previous case study that PIG'N'PI is able to infer the underlying interaction laws correctly. While the other baseline methods are able to precisely predict the particle dynamics, they fail to infer the underlying pairwise interaction laws correctly.

Comparison of PIG'N'PI to alternative hyperparameter choices and an alternative regularized architecture
We compare the performance of the proposed approach first to alternative hyperparameter choices, in particular to different activation functions, and second, to an alternative way of imposing physical consistency in the network architecture.
First, we evaluate different choices of activation functions following the observations made in previous studies (18,21,37) that confirmed their significant influence on the performance of MLPs in approximating physical quantities. The performance of PIG'N'PI with different activation functions is reported in Table S8 in SI Appendix F and Table S9 in SI Appendix G. The results demonstrate that PIG'N'PI with SiLU activation function (which was in fact used in all case studies) performs consistently best on most test datasets compared to PIG'N'PI with other commonly used activation functions, such as ReLU or LeakyReLU. Based on this observation, the performance of the baseline with the SiLU activation function was evaluated (SI Appendix B). The results show that the SiLU activation improves the learning performance of the baseline model to some degree (when only evaluating the prediction performance MAE acc ). However, it still performs consistently worse than PIG'N'PI and, more importantly, the consistency with underlying physics in terms of the inferred force (or potential) and interaction symmetry worsens even considerably.
Second, we compare the performance of PIG'N'PI to an alternative way of imposing physical consistency: we add a regularization into the baseline model to enforce the symmetry property onto the output messages of the edge function. The goal of imposing the symmetry regularization term is to ensure that the model satisfies the action-reaction physical consistency requirement. It is expected that by satisfying this symmetry constraint, the model performance on learning physics-consistent pairwise forces and potential energy can be improved. We add a symmetry regularization term on the learnt pairwise corresponding messages to enforce the action-reaction property. The details on this regularization term can be found in the "Details on imposing a symmetry regularization on the base-line model" section.
While the performance is improved compared to the baseline model without any regularization, the results demonstrate that PIG'N'PI still performs considerably better on inferring physical meaningful quantities for the pairwise force and potential energy than the symmetry-regularized baseline model (see SI Appendix H; Table S10 for detailed results and SI Appendix E for the evaluation on the LJ-argon system).

Conclusions
In this paper, we propose the PIG'N'PI algorithm to learn particle interactions that are consistent with the underlying physical laws.
The main novelty of the proposed algorithm is in the design of the physics operator in the node part. The designed physics operator on nodes guides the edge neural network to learn the pairwise force or pairwise potential energy exactly. This design also reduces the model complexity for this machine learning algorithm by reducing the number of tunable parameters.
While our method shows a similar performance on the supervised learning task of predicting accelerations compared to the other baseline models (purely data-driven GNs), it is also able to infer the particle interactions correctly that follow the underlying physical laws. However, while the baseline models are able to be applied as simulators, they fail to infer physically consistent particle interactions that satisfy Newton's laws. Moreover, the proposed PIG'N'PI also outperforms the other baseline models in terms of generalization ability to larger systems and robustness to significant levels of noise.
The proposed methodology can generalize well to larger particle systems. However, we have to point out that the trained model cannot extrapolate the data arbitrarily far from the training distribution. In our experiments, we found that the edge neural network converges to linear functions outside the training input space. This observation matches the discussion in ref. (38), which is an inherent limitation of MLPs.
The developed methodology will help to make a step forward in developing a flexible and robust tool for the discovery of physical laws in material mechanics. Such tools will be able to support, for example, additive manufacturing with heterogeneous materials that are particularly subject to highly varying material properties, e.g., sustainable or recycled materials (39).

Notations and formal task description
We use a fully connected directed graph G = (V, E) to represent the interacting particle system, where nodes V = {v 1 , v 2 , …, v |V| } correspond to the particles and the directed edges E = {e ij : v i , v j ∈ V, i = j} correspond to particle interactions. Under this notation, v i refers to the i-th particle, and e ij is the directed edge from v j to v i . We use {η t i } i,t to denote the observation of particle states at different time steps, where η t i is a vector describing the state of particle v i at time t. We note that η t i ∈ R 2d+2 (d is the space dimension) includes position r t i ∈ R d , velocityṙ t i ∈ R d , electric charge q i ∈ R, and mass m i ∈ R. The velocityṙ t i and accelerationr t i at time t are computed from the position series of particle v i . We use M i j to denote the message from v j to v i learnt by the neural networkĜ E (·; θ E ) with parameters θ E . Our goal is to infer the pairwise force F t i j and the potential energy P t i j on every edge e ij at each time t given the observation of particle trajectories. All notations are summarized in Table S1.

PIG'N'PI details
PIG'N'PI contains an edge part to learn the pairwise interaction and a node part to aggregate the interactions to derive node accelerations (see Fig. 1). In the edge part, we use MLPs as universal approximators (40,41) to learn the pairwise force or pairwise potential energy. We denote this edge neural network asĜ E (·; θ E ). G E (·; θ E ) takes the vectors η i and η j of two nodes as input. The output M i j ofĜ E (·; θ E ) is the inferred pairwise forceF i j or potential energyP i j on edge e ij , depending on the operator in the node part. We design the physics operator G N ( · ) to aggregate the edge messages in the node part and derive the accelerationr t i for every particle v i at time t. We optimize parameters θ E by minimizing the mean absolute error between the predicted acceleration and the true acceleration. The objective function is given by arg min In the following, we explain the design of the edge neural net-workĜ E (·; θ E ) and the node part operator G N ( · ) in two cases: inferring the pairwise force and inferring the pairwise potential energy.

Learning pairwise force
We use an MLP as the edge neural networkĜ E (·; θ E ) to learn the pairwise force from v j to v i on each edge e ij . The output dimension ofĜ E (·; θ E ) is the same as the spatial dimension d. We first concatenate η t i and η t j which is the input ofĜ E (·; θ E ). We denote the corresponding output as M t i j ∈ R d , e.g., According to Newton's Second law, the net acceleration of every particle is equal to the net force divided by its mass. Hence, in the node part of PIG'N'PI, we first sum up all incoming messages M i = j =i j M i j of every particle v i , and then divide it by the mass of the particle m i . The output of G N ( · ) is the predicted accelerationr We optimize the parameters θ E inĜ E (·; θ E ) by minimizing the objective function Eq. (3). Through this process, the node part operator G N ( · ) guides the edge neural networkĜ E (·; θ E ) to predict the pairwise force exactly, e.g., This is illustrated in Block (B) of Fig. 2.

Learning pairwise potential energy
For the pairwise potential energy case, the edge neural network G E (·; θ E ) is designed to output the pairwise potential energy. Here, the output dimension ofĜ E (·; θ E ) is one because the potential energy is a scalar. We still first concatenate η t i and η t j as the input ofĜ E (·; θ E ) and use MLPs asĜ E (·; θ E ). The corresponding output M t i j ∈ R is denoted as We know that the net force of every particle equals to the negative partial derivative of the potential energy with respect to its position. Hence, in the node part, we first sum up all incoming messages M i = j =i j M i j for every particle i, then compute the negative derivative with respect to the input position and finally divide it by the mass. The final output corresponds then to the predicted acceleration. The node part operator G N ( · ) for the potential energy case is given bŷr Analogously to the force-based case, we optimize for the parameters θ E inĜ E (·; θ E ) by minimizing the acceleration loss (Eq. 3). The node part operator G N ( · ) here guides the edge neural network G E (·; θ E ) to learn the pairwise potential energy exactly. The learnt message on each edge corresponds to the predicted pairwise potential energy, and the negative partial derivative is the predicted pairwise force, e.g.,P This is illustrated in Block (C) of Fig. 2. We note that the commonly used ReLU activation function is not suitable as activation function inĜ E (·; θ E ) for learning the potential energy. The reason is that we compute the partial derivative of M i j =Ĝ E (concat(η i , η j ); θ E ) to derive the predicted accelerations for every particle. The derivative should be continuous and even smooth considering physical forces. However, ReLU approximates the underlying function by piece-wise linear hyper-planes with sharp boundaries. The first-order derivative is, thus, piecewise constant that does not change with input (21). Details on selecting the activation function inĜ E (·; θ E ) are explained in the "Details about simulation and experiments" section.

Details on imposing a symmetry regularization on the baseline model
As mentioned in the "Comparison of PIG'N'PI to alternative hyperparameter choices and an alternative regularized architecture" section, to ensure that the model satisfies the action-reaction physical consistency requirement, we also test an extension of the baseline model by imposing a symmetry regularization on the corresponding pairwise messages in the baseline model. This can be considered as an alternative way of imposing physical consistency. In details, let M i j be the message from v j to v i which is the output of the edge neural network of the baseline model. In our experimental setup, the message M i j corresponds to the force from v j to v i . We impose the symmetry regularization by adding a regularization term on the learnt messages in the objective function [Eq. (3)]. This results in the following objective function: Acceleration loss on nodes Symmetry regularization loss on edges ), where α is a weight parameter. The original baseline model can be considered as the special case with α = 0 in Eq. (10). In our experiments, we evaluate the impact of the regularization term with different weights (α = 0.1, 1.0, 10, 100). The results are reported in Table S10.

Details of the method to learn pairwise force introduced by ref. (34)
Reference (34) proposed a method that has a similar goal to the proposed PIG'N'PI applied to pairwise force prediction. The authors also impose Newton's second law in the standard GN block by dividing the aggregated messages by the node property. We denote this method as GN+. The main difference between GN+ and PIG'N'PI is that GN+ treats the node property as a learnable parameter. It assigns an individual learnable scalar w i for each particle v i and predicts the acceleration of v i by dividing the aggregated incoming messages by 10 w i . The learnable scalars on all nodes representing the pairwise force are learnt together with all other parameters. It is important to point out that GN+ was designed solely for learning the pairwise force while PIG'N'PI can be applied both: to infer the pairwise forces and also the pairwise potential energy. The detailed results of GN+ are reported in Tables S2, S4, and S6, and Figs. 3D, 4, and 8.

Details about simulation and experiments
Here, we summarize the different force functions used in our simulation. Please note that in this work, we used the same case studies as in previous work (22). However, we adapted the parameters of the particle systems slightly to make the learning more challenging.
r Spring force: We denote the spring constant as k and balance length as L. The pairwise force from v i to v j is k(r i j − L)n i j and its potential energy is 0.5k(r ij − L) 2 , where r i j = r j − r i is the Euclidean distance and n i j = r j −r i r j −r i is the unit vector pointing from v i to v j . We set k = 2.0 and L = 1.0 in our simulations.
r Charge force: The electric charge force from v i to v j is −cq i q j n i j /r 2 i j and the potential energy is cq i q j /r ij , where c is the charge constant, and q i , q j are the electric charges. We set c = 1.0 in the simulation. Furthermore, to prevent any zeros in the denominator, we add a small number δ (δ = 0.01) when computing distances.
r Orbital force: The orbital force from v i to v j equals to m i m j n i j /r i j and the potential energy is m i m j ln(r ij ), where m i , m j are the masses of v i and v j . We again add a small number δ (δ = 0.01) when computing distances to prevent zeros in the denominator and logarithm. r Discontinuous force: We set threshold constant = 2.0 such that the pairwise force is 0 if the Euclidean distance r ij is strictly smaller than this threshold and (r i j − 1)n i j otherwise. The corresponding potential is 0 if r ij is strictly smaller than this threshold and 0.5(r ij − 1) 2 otherwise.
We intentionally omit units for variables because the simulation data can be at arbitrary scale. Moreover, the presented cases serve as proof of concept to learn the input-output relation. Further, we note that m i is sampled from the log-uniform distribution within the range [−1] (ln(m i ) ∼ U (−1, 1)). q i is uniformly sampled from the range [−1, 1]. Initial location and velocity of particles are both sampled from the normal Gaussian distribution N (0, 1). Each simulation contains eight particles. Each particle is associated with the corresponding features including position, velocity, mass, and charge. The target for prediction is node accelerations. Every simulation contains 10,000 time steps with step size 0.01. We randomly split the simulation steps into training dataset, validation dataset, and testing dataset with the ratio 7:1.5:1.5. We use T train , T valid , and T test to indicate the simulation time steps corresponding to training split, validation split, and testing split. We train the model on the training dataset [optimizing the parameters θ E inĜ E (·; θ E )] by optimizing Eq. (3), fine-tune hyperparameters, and select the best trained model on the validation dataset and report the performance of the selected trained model on the testing dataset. For generalization tests, we re-run each simulation on 12 particles with 1500 time steps (same size as original testing dataset). The previously selected trained model with eight particles is tested on the new testing dataset.
We only fine-tune hyperparameters on the spring validation dataset and use the same hyperparameters in all experiments. We set the learning rate to 0.001, the number of hidden layers in the edge neural network to 4, the units of hidden layers to 300, max training epochs to 200. The dimension of the output layer in the edge neural network is d to learn the force or one to learn the potential energy. We use the Adam optimizer with the mini-batch size of 32 for the force case study and eight for the potential case study to train the model. The SiLU activation function is used in all PIG'N'PI evaluations.