Physics informed neural network for charged particles surrounded by conductive boundaries

Molecular dynamics of charged particles in porous conductive media have received considerable attention in recent years due to their application in cutting-edge technologies such as batteries and supercapacitors. Due to the presence of long-range electrical interactions, induced charges present at the boundary, and the influence of boundary conditions, the simulation of these systems is more challenging than the simulation of typical molecular dynamic systems. Simulating these kinds of systems typically involves using a numerical solver to solve the Poisson equation, which is a very time-consuming procedure. Recently, Physics-Informed Neural Networks (PINNs) have been introduced as an alternative to numerical solutions of PDEs. In this paper, we present a new PINN-based model for predicting the potential of point-charged particles surrounded by conductive walls. As a result of the proposed PINN model, the mean square error is less than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$7\%$$\end{document}7% and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{2}$$\end{document}R2 score is more than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90\%$$\end{document}90% for the corresponding example simulation. Results have been compared with typical neural networks and random forest as standard machine learning algorithms. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{2}$$\end{document}R2 score of the random forest model was \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$70\%$$\end{document}70%, and a standard neural network could not be trained well.


Introduction
Computational Electromagnetic Simulation (CES) plays a significant role in many areas of science and engineering, such as soft matter, electrical engineering, biomedical engineering and chemistry.In addition, it has numerous applications in industry.For example, it is one of the main tools in investigating and designing the process of supercapacitors, which are porous energy storage devices with many applications in industry, especially when high power consumption or transfer is neededMiller and Simon [2008].Here, studying the physical mechanisms arising from charge storage in supercapacitors is essential for further technological developmentSalanne et al. [2016], Simon and Gogotsi [2008].Solving Maxwell's equation, especially the Poisson equation in this study, is an essential part of computational electromagnetic algorithmsJackson [1962].Solving the Poisson equation can help scientists to calculate the potential of electrical sources in any system.However, many difficulties arise in practice due to the long-range nature of electrical interactions.In particular, estimating the potential of point-charged components in an environment with conductive walls is challenging because of the induced charges presented on the boundaries.
Generally, there are two approaches to solving the Poisson equation: analytical solutionJackson [1962] and numerical methods.There are limited techniques for solving analytically, like image charges methods applicable for cases with regular geometries; however, there is no guarantee to achieve practical results.If, for example, a particle is placed in a cubic conductive container, the image charges method will produce an infinite series.On the other hand, numerical methods lead to approximate solutions based on discretizing space and/or time domains.One of the typical numerical methods is the Finite Element Method (FEM)Jin [2014] which discretizes the continuous partial differential equations (PDEs) and forms a linear set of algebraic equationsS.et al. [1991].Nevertheless, even FEM fails in calculating the potential in a charged particle's position since the electrical potential is singular at the place of charges.There are a number of methods and algorithms that have been developed to address this problem, including Induced Charge MMM2D (ICMMM2D) Tyagi et al. [2007] for 2D, ELCICTyagi et al. [2008] for 2D + h, Induced Charge Computation (ICC * ) Tyagi et al. [2010], Kesselheim et al. [2010], Arnold et al. [2013] for 3D periodicity, and a method introduced by Reed et al.Reed et al. [2007] have been developed.In addition, recently, there has been another algorithm named PLT.It was first demonstrated for a partially periodic system constrained between two metallic plates in Rostami et al. [2016], and then it was applied to CAVIAR Biagooi et al. [2020], a molecular dynamics simulation package for charged particles surrounded by non-trivial conductive boundaries.Numerical solving of these problems with the CAVIAR package is accurate; moreover, it took less time than ICC * Biagooi et al. [2020] but is still time and memory-consuming.
Recently another data-driven approach to solving the PDEs based on deep machine learning is also of great current interest.For instance, Shan et al.Shan et al. [2020] present a CNN to predict the electric potential with different excitations and permittivity distribution in 2D and 3D models.It is fast and efficient compared with FEMJin [2014].However, a couple of problems prevent it from utilizing as a Poisson solver in the MD simulation process; first, it could not work in the case of discrete density functions such as those of point charges, and second, it is a physics-free approach which makes it hard to consider boundary conditions.To overcome the first problem, one can use the PLT algorithm.Additionally, Raissi et al. introduced the physics-informed neural network (PINN) that the loss function defined by Raissi et al. [2019] is an excellent alternative to the conventional deep learning method because of the governing equations, boundary conditions, and initial conditions used in its definition.
In this paper, we applied a new PINN-based model to predict the potential of point-charged particles surrounded by conductive walls.We then compared the results with typical neural networks and random forests as a standard machine learning algorithm.For instance, we tried to implement these models for a charged particle in a spherical container.The reason for utilizing this simple example was that there is an exact solution to this problem through the analytical method, the image charges method.As a starting point, we used the PLT algorithm to transfer the Poisson equation into the Laplace equation with modified boundary conditions.Then we trained the model to solve the Laplacian equation with new boundary conditions.The input data is included the position in which we want to evaluate the potential on it and the modified boundary conditions; the output data is the corresponding electrical potential of that position.

Methods
This article aims to build a machine-learning model (ML-Model) to predict the potential of point-charged particles surrounded by conductive walls.The potential of charged particles is calculated by solving the Poisson equation, which can be written as Jackson [1962]: where φ is the potential and ρ is a charge distribution.The first and straightforward ML-Model that jumps to mind is a model that includes x q and x as an input, and φ(x q , x) as an output.Here x q is the position of point-charged particle, x is the position in which we want to calculate the potential on it, and φ(x q , x) is the corresponding potential.So the number of input features depends on the number of charged particles; for instance, in 3 dimensions, if there is N charged particles, the input features have to be 3 + 3 × N .Therefore, this kind of model could only predict the potential of fix number of charged particles.In many applications of this method, such as the molecular dynamic simulation, this number is not fixed and could even increase or decrease during the simulation.We use the PLT algorithm to transpose the Poisson equation into the Laplace equation with new boundary conditions to overcome this problem.This algorithm will be discussed in more detail in 2.1.So we can train a model which includes x and modified boundary conditions as input features and φ(φ b , x) as an output.We define the boundary conditions only on N b fixed points on the boundary {φ 1 , φ 2 , . . ., φ N b }.In this case, with the PLT algorithm, we can build a model with a fixed number of input features that can predict any charged particles' potential.

Poisson to Laplace Transformation (PLT)
According to the PLT algorithm, the electrical potential is divided into two parts: singular potential (φ si ) and smooth potential (φ sm ); φ( x) = φ si ( x) + φ sm ( x).It is important to note that the smooth part here is the solution of the Laplace equation with modified boundary conditions, while φ si obeys the famous Columb laws It can be seen that the modified boundary condition for φ sm is represented by where φ x bc corresponds to the initial electrical potential on the boundaries.Finally, with the PLT algorithm, we could transfer the Poisson to the Laplace equation with new modified boundary conditions, then train an ML-Model with these modified boundary conditions as an input parameter and the smooth potential as an output.Afterward with the summation of singular and predicted smooth potential, we can reach the total potential.Advantage of utilizing PLT algorithm is that it leads to having a fixed number of input data since the number of input data would be independent of the number of point-charged particles.

Data Engineering
For training a highly accurate model, having a nice train set is crucial.In this work, the reference set is , where the input is concluded {x, y, z, ϕ bc = {φ 1 , φ 2 , ..., φ N b }} and φ is the target.
x, y, z is the coordinate of a point in the container on which we want to calculate their potential, φ is the numeric value of the potential at this point, and {φ 1 , φ 2 , ..., φ N b } is the boundary condition on N b points on the boundary.First, in the container, N p positions are chosen to predict the potential in their situations.
In fact, for each boundary condition, {φ 1 , φ 2 , ..., φ N b }, there are N p points that we want to calculate the potential at their positions.Then, the reference set could be created for N q different boundary conditions.So the reference set consists of N = N q × N p samples which could split to train and test set.In this case, our container is a sphere, we also set N p = 78, N b = 26, and N q = 100.So, our reference consists of 100 different boundary conditions and for each boundary condition {φ 1 , φ 2 , ..., φ 26 } there are 78 points in the sphere on which we want to calculate the potential on it.We use the solution of the image charges method (Eq.5) to calculate targets of the reference set: where a is the conductive spherical shell radius and r is the distance of a point charge q from its center.The numeric value of potential is minimal (∼ 10 −9 ), which conducts to significant rounding error during computation; therefore, the potential of an electron in a 1m distance of it, 1.44 × 10 −9 [V ], is used as a unit to make equation 5 dimensionless.We randomly chose 5000 and 1000 samples from the reference set to create a train and test set.The train and test set have no samples in common.In addition, the best model could adequately predict the potential of test samples and samples with different boundary conditions from the train and test set.So to evaluate the model better, we prepare an extrapolation set that includes 1000 samples with 55 distinct boundary conditions.

ML Algorithms
In this work, three different supervised learning methods have been used, and their regression accuracy, based on the metrics presented in 2.4, has been evaluated.Mainly, we stick to Physics-Informed Neural Networks (PINN, Raissi et al. [2019]), but to compare our results with other ML algorithms, we use Random Forest (RF, Breiman [2001]) and Artificial Neural Networks (ANN).All the models are briefly introduced, the hyperparameters are fine-tuned, and their performance is reported.Scikit-learn Pedregosa et al. [2011], Tensorflow Abadi [2016], Keras Chollet [2015], and NumPy Walt et al. [2011] are all the python libraries that have been used in this project.

Random Forest(RF)
RF is one of the most popular machine learning algorithms in regression problems for many reasons, but in this project, this model has been chosen since a) It is speedy to learn.b) It is robust against over-fitting.over-fitting is detected when the performance of train samples is perfect while the performance of test samples is poor.RF is an ensemble model in which an average of many uncorrelated trees determines the predicted potential for the target data set.Although each tree is a weak learner, they make a strong learner when many trees are grouped.The RF randomizes the trees by choosing a subset of training data and features for each tree.Here we use scikit-learn Pedregosa et al. [2011] RF implementation.

ANN
Typical neural network architecture consists of the input layer, multiple hidden layers, and the output layer with several neurons in each layer.Totally: • Input layer: The neurons in the input layer are the input features.
• Hidden layers: The value of every neuron in the hidden layers is a linear combination of the neurons in the previous layer followed by the implementation of an activation function(Eq.6); in most cases, the activation function is non-linear.
n is the layer number, w and b are the model parameters, weights and bias respectively, and σ l is the activation function based on Goodfellow et al. [2016].
• Output layer: The neurons in the output layer are the model targets and they are calculated with Eq.6 with linear activation function.
• Loss function: There is a function in all neural networks that must be minimized over the model parameters during the training stage via back-propagation, typically the loss function is the mean square error between the true and the predicted values.
U and T are the predicted output and true target values, respectively, X is the input data, and w is the parameter of neural networks, weights, and biases.The first sentence in Eq.7 is a mean square error, and The second sentence exists to prevent over-fitting, namely L2 regularizationa.K. Connect et al. [1992], that is used in order to reduce the effects of the large weights.

PINN
PINNRaissi et al. [2019] enforces the Laplace equation, a physical law of the electromagnetic system, as a constraint on the neural network.This study proposes a PINN-based approach to solve the Laplace equation with changeable boundary conditions.Fig. 3 shows a schematic of the neural network layout for this approach.PINN-based models are neural networks with modified loss functions: The first and the last term of Eq.9 are the same as typical neural networks in Eq.7.The second term corresponds to the governing physical equation, i.e., the is Laplace, and the third term corresponds to the boundary conditions; and Here we define f (x, u; w) with Dirichlet boundary conditions λ 1 , λ 2 , λ 3 in Eq.9 correspond to the weight coefficients for the data contributions, Laplace equation, and boundary losses.We use the weight coefficient by motivating from the study of Kag et al.Kag et al. [2022].
The last sentence is the L2 regularizationa.K. Connect et al. [1992].Notice that the model with λ 2 = λ 3 = 0.0 is exactly a typical neural network described in the previous subsection.

Evaluating Metrics
The performance evaluation of different algorithms for potential estimation depends on different metrics, ∆φ = φ T rue − φ P red , (15) Where φ T rue is the true potential, φ P red is the predicted potential, and φT rue is the mean true potential of a given test sample.In this study we used scatter σ, R 2 score and M SE as our evaluating metrics.

Result
In this paper, we predict the smooth potential of a point-charged particle in a spherical conductive container.First, we set the train and test set with 5000 and 1000 samples (2.2), then we train our models to predict smooth potential.We can calculate total potential by summating smooth and singular potential (more detailed in 2.1).However, in this work, to compare our results with CAVIARBiagooi et al. [2020], we investigate the smooth potential.

hyperparameter for RF
We optimize over the only hyperparameter, the number of trees in the forest that influences the fitting of the random forest model.In Fig. 4, we plot M SE (the left panel) and R 2 score (the right one) as a function of the number of trees for the test set to determine the optimal hyperparameter, which we find 100 trees since progress after 100 trees are negligible.Afterward, we trained the RF model using 100 trees.The RF method is relatively fast; however, it works when the predicted potential is smooth and relatively small, it is not suitable in the case of point-charged particles (where a gradient of potential as well as its numeric value is high at the position of the charge).Furthermore, it fails to predict the potential near the boundaries since the gradient of the potential is considerable.

PINN based model and NN
By setting both λ 2 and λ 3 to zero in PINN-based model, one can get exactly NN model.Therefore we investigate both models together and report their reulsts at the same time in the following section.

hyperparameter for PINN and NN
Unlike the RF model, we define several hyperparameters: a number of neurons, a number of layers, λ 2 , and λ 3 .To tune all hyperparameters, we train the model up to 100000 epochs, using the L-BFGS-B optimizerLiu and Nocedal [1989], until the model's tolerance reaches the level of machine epsilon.For all layers except the last one, we use a tanh activation function.Table 1 is reported the M SE between the predicted and the same potential for a different value of hyperparameters; λ 2 = [0, 0.1, 0.2, 0.3], λ 3 = [0, 0.1, 0.2, 0.3, 0.4], number of hidden layers= [1,3,5,7] and number of neurons per hidden layer= [10,30,50] for 1000 samples of the test set.As shown in Table 1, we observe that a model with one hidden layer could not predict the potential well.Also, a model with ten neurons per layer could not work well.So, Table 2 and Table 3 reported the results for just [3,5,7] layers as well as for [30,50] neurons per layer.We chose λ 4 = 0.0001 to prevent over-fitting.

PINN and NN Prediction
In Contrast with Table 1, Table 2 is reported not only the M SE but also the R 2 score of the test set.As can be seen in Table 2, the model with seven layers and 50 neurons per layer resulted better when λ 2 and λ 3 are 0.3, 0.3, or 0.2, 0.4, respectively.When λ 2 and λ 3 are zero, a standard neural network, the model has not worked well; it is observed from Tabel 2 and Fig. 4.    1: M SE between the predicted and the exact potential φ(x) for a different value of λ 2 , and λ 3 , and the different number of hidden layers, and neurons per hidden layer in PINN for 1000 different sample of the Test set.Here, λ 4 = 0.0001 is fixed and λ 1 = 1 − (λ 2 + λ 3 + λ 4 ).In this table, the bold number means M SE < 0.1.

Comparison
We evaluate RF, NN, and PINN to estimate the potential of point-charged particles surrounded by conductive walls.According to Fig. 6, NN was not trained well, while RF and PINN-based models could predict potential precisely.However, RF did not work well to estimate φ T rue > 0.3.Apart from this, the best model could estimate not only the potential of the train and the test sets but also the potential of point-charged particles that are not in the train or test set.So we evaluate the best tuned-PINN model and RF on the extrapolation samples; the results are reported in Table 3.As can be seen in Fig. 8 PINN-based model could predict the potential of newly charged particles better than the RF model, where PINN could predict φ T rue > 0.3 by far better than RF.

Generalization (Multi charged particles)
For generalization, we test the PINN-based model with λ 2 = 0.3 and λ 3 = 0.3 for the case of more than one charged particle surrounded with conductive boundaries.Since the Laplace equation is a linear function, we predict the potential of each charged particle and then calculate the total potential with a superposition of the corresponding predicted potential.After that, we report the M SE between the predicted smooth potential and the exact smooth solution, which is calculated by the image charges method.Fig. 9 shows the relation between M SE and N , the number of charged particles.As expected, the M SE is independent of the number of charged particles.Therefore, it leads to the fact that we can also use this method for problems with any desired particles.

Conclusion
In this study, we have trained a machine to predict the smooth potential of charged components surrounded by conductive boundaries.In this scene, the total potential could be easily calculated by the summation of predicted smooth potential with singular potential due to the PLT algorithm.The reference set consists of analytic solutions, the solution of the image charge method, which is split into a train set with 5000 samples and a test set with 1000 samples.To check the accuracy of our model, we set another data set called extrapolation set consisting of 1000 samples with different boundary conditions which were not in the train or even test set.Our main conclusion can be summarized as follows: • We find that the PINN-based model trained better than RF and NN models.RF could not predict high potential; on the other hand, the NN could not be trained well at all.
• our PINN-based model could predict the potential of the test set with M SE = 0.069, R 2 score = 0.902, and scatter σ = 0.02.It also could predict the potential of the extrapolation set with M SE = 0.089, R 2 score = 0.851, and scatter σ = 0.02.• Since the Laplace equation is a linear equation, the trained model could predict the potential of more than one charged particle by summating every particle's predicted potential.Besides, we show that the M SE of more than one particle is independent of a number of particles.

Figure 2 :
Figure 2: Schematic of the PLT method: a) the main system which had point charges inside of it b) the new system without any point charges and the boundaries were modified c) N b points on the boundary are shown to be used as our model input.
Figure 3: Physics-informed neural network scheme for solving Laplace equation with variable boundaries

Figure 4 :
Figure 4: M SE (left), and R 2 score (right) for 1000 different sample of Test set

Figure 5 :
Figure 5: Potential estimation of RF model: a) show the train data set with 5000 samples with scatter of σ = 0.02, b) show the test data set with 1000 samples with scatter of σ = 0.07.The dashed red line shows where the predicted potential equals the true potential.The pink-shaded region marks 1σ scatter of potential errors.

Figure 6 :
Figure 6: Potential estimation of best-tuned (a) NN with scatter of σ = 0.07 and (b) PINN model on the 5000 samples of the train set with scatter of σ = 0.01.The dashed red line shows where the predicted potential equals the true potential.The pink-shaded region marks 1σ scatter of potential errors.Both plots in Fig.6compare the true and predicted potentials for the best-tuned NN and the best-tuned PINN model with 7 layers and 50 neurons per layer, λ 2 = 0.3 and λ 3 = 0.3-on the train set.As can be seen, the NN model is not trained well, while the PINN-based model could predict the potential precisely with a scatter of 0.01.Although the PINN-based model predicts the train set well, aiming to clarify that over-fitting has not been accrued, we also evaluate the model on the test set, Fig.7.

Figure 7 :
Figure 7: Potential estimation of best-tuned PINN model on the 1000 samples of the test set with scatter of σ = 0.02, M SE = 0.069 and R 2 score = 0.851.The dashed red line shows where the predicted potential equals the true potential.The pink-shaded region marks 1σ scatter of potential errors.

Figure 8 :
Figure 8: Potential estimation of best-tuned (a) RF with scatter of σ = 0.07, and (b) PINN-based model with scatter of σ = 0.02 on the 1000 samples of the Extrapolation set.The dashed red line shows where the predicted potential equals the true potential.The pink-shaded region marks 1σ scatter of potential errors.

Figure 9 :
Figure 9: M SE between true and predicted potential as a function of charged particles number; in 100 different problems

Table 2 :
M SE, and R 2 score between the predicted and the exact potential φ(x) for a different value of λ 2 , and λ 3 , and the different number of hidden layers, and neurons per hidden layer in PINN for 1000 different samples of the Test set.Here, λ 4 = 0.0001 is fixed and λ 1 = 1 − (λ 2 + λ 3 + λ 4 ).In this table bold numbers show cases with M SE < 0.1 and R 2 score > 0.9.

Table 3 :
M SE and R 2 score between the predicted and the exact potential φ(x) for a different value of λ 2 , and λ 3 , and the different number of hidden layers, and neurons per hidden layer in PINN for 1000 different samples of the Extrapolate set.Here, λ 4 = 0.0001 is fixed and λ 1 = 1 − (λ 2 + λ 3 + λ 4 ).In this table bold number shows the best hyperparameters for our PINN-based model.