CAVIAR: A simulation package for charged particles in environments surrounded by conductive boundaries

We introduce CAVIAR, new scientific software that is developed for molecular simulation of ionic-liquids or charged colloids inside conductive boundaries. CAVIAR imports computer-aided-design geometries and uses them to simulate boundary walls. Then based on this geometry, a finite element mesh is generated and utilized for solving the Poisson equation. To avoid the complexity due to the singularity of point charges, we propose a new method, using the advantages of the superposition theorem of the linear partial differential equations. Within this paper, the CAVIAR structure, its features, and basic scientific algorithms are discussed.


I. INTRODUCTION
Soft-matter physics is of high technological importance. This area of physics mainly contains the study of liquids, colloids, polymers, and biological and granular materials. 1 Systems with charged components are an important subset of soft-matter, which present uncommon behaviors. Charged particles, both in the dielectric as well as the conductive environments, are of scientific interest, due to their important rules in the physical behavior of such systems. In particular, charged particles surrounded by conductive boundaries that are observed in many physical systems are both of technological and theoretical interest. As a special field of interest, one can refer to conductive porous electrodes in porous super-capacitors, which have gained much technological attention as a candidate for future energy storage units. [2][3][4][5] Furthermore, capacitive deionization by porous electrodes as a method of water desalination is another example of the above-mentioned systems. [6][7][8] Any progress in these two fields of research can have high impact on human's quality of life. Due to the complexity of the charged systems, simulation is an effective and efficient way to investigate these systems; it enables us to figure out the details of physical events in different scales of time and space up to the available computational resources.
The Molecular Dynamics (MD) simulation is a well-known and popular method for investigating many-body systems. [9][10][11][12] However, for systems with long-range interactions such as electrostatic and gravitational interactions, the complexity of the standard MD algorithms is on the order of O(N 2 ), which is computationally inefficient. Furthermore, in charged systems under partial or full periodic boundary conditions, infinite summations appear that may not converge easily. To deal with this problem, Ewald sum 13 based algorithms have been created. [14][15][16][17][18][19][20][21] In these algorithms, the potential is divided into two short-range and long-range parts in a way that the long-ranged part always converges.
However, the situation can even be worse when charged particles are surrounded by conductive boundaries; here, the complexity arises from the existence of induced charges on the boundaries. In order to address this problem, some methods and algorithms such as Induced Charge MMM2D (ICMMM2D) 22 for 2D, ELCIC 23 for 2D + h, and induced charge computation (ICC) * 24-26 for 3D periodicity have been developed recently. The ICMMM2D (Induced Charge MMM2D) method 22 is an extension of MMM2D, which is able to treat planar dielectric interfaces using image charge particles (ICPs). The ICC * , which was first proposed by Tyagi et al., 24 is a generic method that can be used for dielectric interfaces of arbitrary shapes and is based on the induced charge computation 27 (ICC) method. The ICC * needs a Poisson solver suitable for the same geometry in the absence of the dielectric interface, and the star denotes a place holder for any Poisson solver. For example, the ICC * is used for simulating carbon electrodes with the shape of slit-pore containing ionic-liquids. 28,29 There is another method introduced by Reed et al., 30 which is based on the work of Siepmann and Sprik, 31 which uses a variational procedure for variable charges to set the constant potential condition. This method has also been used in some other super-capacitor simulations. [32][33][34] In this paper, we introduce a new method for simulating charged particles in conductive media. A simple version of this algorithm was introduced and used for simulation of a partially periodic system confined between two metallic plates. 35 This algorithm utilized the P 3 D method, which is a Laplace solver for slab geometries that utilizes both plane waves and finite element method. 17 Based on this algorithm, we develop CAVIAR, which is a new, multi-purpose, and cross-platform scientific software package. The name CAVIAR is from "Charged pArticles in enVironments surrounded by conduc-tIve boundaRies." CAVIAR was originally developed as a molecular dynamics simulation code for charged particles near conductive boundaries with non-trivial geometries. It can also be used to simulate problems with less complexity such as MD systems without conductive or periodic boundaries. CAVIAR has been developed in C++ with an extensible design. The rest of the paper is as follows: In Sec. II, we first present all the details about induced charge calculations and talk about the interactions of particles with geometries. In Sec. III, we introduce the package structure, interpreter, implementation of codes, parallel simulation methods, external libraries, and other open-source software that should be used as an additional tool. In Sec. IV, we discuss the practical uses of the CAVIAR with a simulation example. Section V has the results of some test simulations that compares the performance of the implemented algorithm with another software. At last, a conclusion is provided in Sec. VI.

A. Poisson to Laplace transformation (PLT) algorithm
When dealing with conductive boundaries, induced charges appear on the boundaries, which in turn make the computation of the electrical potential generally a difficult task. It is well known that the relation between charge distribution ρ(⃗ x) and electric potential φ(⃗ x) obeys Poisson's equation ∇ 2 φ = −ρ/ϵ 0 . Even for regular geometries such as a point charge in a cubic container with conductive walls, the analytical methods of solving the Poisson equation, like image charges method, results in an infinite series with dramatically slow convergence. 36 Besides, typical methods of numerical solution of the Poisson equation fail in the case of a system of point charges because of the singularity of the electrical potential at the position of point charges. In order to avoid these singularities, we use the superposition theorem and divide the electric potential into two, singular (φsi) and smooth (φsm) parts, φ(⃗ x) = φsi(⃗ x) + φsm(⃗ x), where the smooth part is the solution of the Laplace equation under modified boundary conditions, whereas φsi in each point is simply The modified boundary condition for the φsm is given by where φ| ⃗ x bc is the original electrical potential on the boundaries, as schematically illustrated in Fig. 1.
To check the accuracy of the algorithm and our package, we compare the results of the simulation of a point charge q in a conductive spherical shell of radius a at a distance of r from its center. The analytical solution to this problem is This case is simulated in CAVIAR-1.0 for a range of point charge distance from the center of the sphere, and the numeric result is compared with its analytical solution. Figure 2 is a plot of the relative error of computed electrical potential on the point charge as well as total induced charge on the sphere as a function of relative distance of the charge from the center of the sphere, for three different levels of mesh refinement. By approaching the conductive boundaries, the gradient of the electrical potential increases and as a result, the relative error of the computed electrical potential also increases. As shown in Fig. 2(f), the computational error is reduced with the mesh refinement. Figure 3 is the flowchart of the CAVIAR package. It contains two major parts: the main loop of the Molecular Dynamics (MD) algorithm and the force field computation of long-range electrostatic forces using the Poisson to Laplace transformation (PLT) algorithm. The MD loop consists of computing the force fields of all short-range interactions such as particle-particle as well as particle-wall interactions, the PLT procedure for evaluating the electrostatic long-range interactions, and the integration of the resulting equation of motions. The PLT procedure, flowcharts  The inset plot of (f) is the total induced electrical charge on the shell. As it can be expected, with a higher level of refinement, a more accurate answer can be achieved.
of which are sketched on the right side of Fig. 3, has the following steps: 1. PLT initialization: This step, which is done just once, includes generating finite element mesh, refining it, and importing it into the code, then assembling it with the weak (integral) form of the Laplace equation.

Apply boundary values: Calculate and apply boundary values
φsm| ⃗ x bc as constraints on the system. 3. Solve Finite Element (FE) equations: In this step, the resulting discrete Laplace equation is solved, then based on it, the electrical field resulted from the smooth part of the Laplace equation is computed.

Compute singular force:
In this step, the singular electrostatic force is computed by summing over all Coulomb interactions between electrical charges. 5. Compute total electrostatic force: By summing the smooth and singular forces applied to each electrical charge, the total electrostatic force is computed in this step.
Taking periodic systems into account, one can utilize the Ewald based methods 13 to compute the singular part of the potential [Eq. (2)] (see the work of Lee and Cai 37 for more information).
The value of induced charge, q, on each 2D mesh element of the electrodes can be calculated by means of the Gauss law, 36 i.e., ∮ S ⃗ E ⋅ ⃗ nda = q/ϵ, for a closed surface S that includes the mesh The short-range interactions include all short-range particle to particle and particle to wall interactions, which are computed using standard methods of MD. The electrostatic force field calculation, on the other hand, utilizes the PLT algorithm for evaluation, the sketched flowchart of which is on the right. element. Here, ⃗ n is the normal vector of the mesh element, ⃗ E is the electric field, and da is its area. The area of each mesh element and its normal is calculated once at the initialization step of the simulation. The electric field is computed by numerical differentiation of the total electric potential computed at each PLT step.
Order of the PLT algorithm depends on the number of 2D boundary cells, NB, the number of charged particles, N, and the order of the finite element solver of the Laplace equation. As discussed before, PLT is divided into different steps with different orders. Setting up boundary values requires the calculation of electric potential on all of the 2D boundary cells, which is in the order of O(NBN). Solving the Laplace equation using the finite element methods depends on the number of mesh points (Nm), which is a function of the geometry. Furthermore, it depends on the solver employed by FE; the best solvers usually are in the order of O(Nm log Nm). Calculating a portion of electrical forces on particles with the smooth potential is done by interpolation in cells that surround them and have the order of O(N). Calculating those with singular potential is done by the Ewald method, O (N log N), as well.

B. Particle-boundary interaction algorithm
A solid boundary or wall in an MD simulation can be created using fixed particles. This is the most common way in granular or atomic simulations. 38,39 However, this is a timeconsuming method since it imposes extra calculation needed to compute forces of these numbers of fixed atoms. Furthermore, for simulating meso-scale and macro-scale systems such as charged colloidal particles and granular materials, the surface of boundaries will be rough if it is created by particles, against those in the real world. As an alternative, geometric shapes can be used to construct the boundary walls. In this case, the nearest distance D from a Cartesian point P to a geometric object can be calculated efficiently. D has a simple form, for convenience, geometry shapes such as a sphere, cylinder, or an unlimited plane. For more complicated cases, however, it is possible to approximate the boundary's geometry with many connected flat polygons, i.e., a polyhedron. In this case, D is the minimum of distances between particles and each polygon. To compute the collision force efficiently, as schematically illustrated in Fig. 4(a), we first compute the compression vector ⃗ ξ from all boundary polygons using the following procedure: 1. Surface: First, we calculate the orthogonal distance D of a particle from the surface of the polygon. If the projected point is inside the polygon and D is less than the size of particle R, the compression vector is equal to ⃗ ξ = (D − R)n. 2. Edge: If the above step does not happen, we check the orthogonal distance D of the particle from the polygon edges. If this In each step, if the distance is less than the criteria distance R, which might be the particle radius or the force field cutoff, the algorithm returns the compression vector. When there is more than one contact point, as schematically represented in (b), vector summation of all compression vectors is used to compute the net force imposed by the wall.
is less than R and the projected point is on the edge, we set the compression vector to ⃗ ξ = (D − R)ê. 3. Vertex: Finally, if the previous step does not happen, we calculate the distance from all of the polygon vertices. The nearest vertex, which has a distance D less than R, is used to return In addition, if none of the above conditions are satisfied, there is not any collision between the particle and polygon.
As schematically illustrated in Fig. 4(b), this method works properly even when there is more than one single polygon for collision at the corner of the boundaries. The above procedure is an optimized version of the Gilbert-Johnson-Keerthi algorithm 40,41 for contact detection in MD. This type of contact algorithm for nonspherical elements is categorized in the branch of discrete-elementmethod algorithms, which is used alongside MD or FE simulations in physics and engineering. [42][43][44] This method is applied to many areas of charged-particle simulation, from the atomic scale to the macro-scale simulations, as well as for a range of force fields. As an example, it can be applied to the Lennard-Jones (LJ) as an atomic level force field, 39 or the Hookean or Hertzian interaction in the simulation of the granular materials. 38 In the case of Hookean or Hertzian forces fields, the compression ξ of the particle into the wall is needed to be calculated. For the LJ force field, ξ is the overlapping distance between the wall and particle potentials. The compression vector is in the direction of the imposed force from the wall to the particle. This procedure is simple and efficient; if the force field parameters and time step are set correctly, it will not let the particles pass through the boundaries.
The order of this algorithm depends on the multiplication of the number of particles N and the number of triangles Nt, resulting in a performance downgrade specially in the case of complicated geometries. In order to increase the performance, we define a boundary cell-list. To do that, we first decompose the simulation box into cubic cells and then dedicate the index of all boundary triangles, which intersect with each cell, to the list of that cell. Furthermore, in order to take all possible cutoffs into account, we add the index of those triangles in a distance of less than or equal to R from cell boundaries to the list of that cell. Now for calculating the particle-boundary interactions, the algorithm primarily finds the cell containing the particle and then computes the interaction of the particle with all triangles stored in the cell-list. This algorithm is also implemented in CAVIAR.
To set the initial position of particles and atoms, the famous ray-casting algorithm 45,46 is used to determine whether a point is inside or outside the simulation volume. In the mathematical view, a polyhedron separates the space into two volumes, meaning inside and outside. 47 In the ray-casting algorithm, one needs to emit a straight ray from the point to infinity and sums the number of surface contacts. An odd number means that the point belongs to the inside of volume, and even numbers indicate that the point belongs to the outside one. Because particles should not have any contact with boundaries at the initial stage of the simulation, we mixed this algorithm with the above-mentioned procedure when calculating the initial positions.

A. Implementation
CAVIAR is a software with its own scripting interpreter. It has been developed with object-oriented programming paradigm 48 and modern C++11 standard. 49 The package uses CMake 50 for its crossplatform configuration tool. The only external dependency is deal. II library. 51,52 The core of the package is a combination of two major components: interpreter and objects. All the physics and mathematics tools of CAVIAR are called "objects." Adding a new object such as a force field to the package is straightforward, the procedure is the duplication of a similar object file, then renaming the file name and the class within it. The "interpreter" handles the creation, calls, and the deletion of the objects. It provides meaningful run-time warnings and errors as well as interprets the CAVIAR scripting language (CASL). An example of CASL is illustrated in Fig. 5.
The CAVIAR main directory is divided into subdirectories, which are listed as follows with a short description: The red words are the base class. The words next to them are the child object type and the object name. The script obeys the object-oriented structure of the code. One can create unlimited numbers of any objects and use them. For example, here we have created two "Neighborlist" objects and only the first one is used in the "Md_simulator" object.
• include: C++ header files for the interpreter, objects, and utility. • src: C++ source files including the package main function "main_CAVIAR.cpp" and the other source files. • doc: Package documentations.
• example: Working scripts for the tutorials and references.
• test: Working scripts that has been used to test the code features. • scripts: Bash-shell scripts for users and developers.
• CMake: Additional CMake script for the package configuration. • tools: Some preprocessing and post-processing C++ codes alongside other tools.
CAVIAR objects are divided into 10 base classes. They are listed as follows: • Atom_data: The container of all the atoms and molecules data that are active in a simulation. • Constraint: It has different constraint algorithms such as NVE, 53 Berendsen, 54 and Nosé-Hoover 55,56 thermostats. In addition, bond fixes with RATTLE, 57 SHAKE, 58 and M-SHAKE 59 algorithms are also implemented.
• Domain: It has the specification of the simulation box. It handles the domain decomposition in the case of parallel runs. • Force_field: The base class of all the force fields. The most important implemented force fields are Lennard-Jones (LJ), geometry with LJ, geometry spring, PLT, dissipative particle dynamics (DPD), electrostatic, periodic electrostatic with Ewald sums, spring bond, and spring angle. • Integrator: Different integration schemes such as Euler, leap-frog, and velocity-verlet. • Md_simulator: This object connects all the objects related to an MD simulation. • Neighborlist: It has a cell-list and a verlet-list as its child objects. Short-range force fields need a pointer to this type of object. • Shape: It has different objects that handle the boundary shapes. The geometric force fields need a pointer to this type of object. • Unique: The objects that do not belong to any base classes are put here. For example "Distribution" generates initial positions of atoms.

ARTICLE scitation.org/journal/adv
• Writer: This object is responsible for writing simulation results to the output files.
PLT force field classes are implemented in a way that it does not have any potential formula. It uses other force field potential such as simple electrostatic or periodic electrostatic to calculate the potential values. In this way, the user does not need to create a new PLT class for a new type of potential. In addition, the PLT algorithm needs a finite element routine. CAVIAR currently uses deal. II, 51,52 which is a C++ library for FE simulations on hexahedral meshes. CAVIAR can be configured without deal. II library too. In this case, it can do MD simulations without the PLT algorithm. There are some examples provided in the package and described in the user guide that works without deal. II, including a simple point charge extended (SPC/E) water model simulation and a polymer with fluctuations.

B. Parallelization
CAVIAR runs in both sequential and parallel mode using message passing interface (MPI). While using PLT, every time step has MD and FE calculations and each of them needs a different method of parallelization. To make the MD procedure parallel, the simulation box decomposes into smaller boxes according to the number of processors, Np. This decomposition could be 1, 2, or 3-dimensional in a way that domains have the minimum common interface between each other. This method is called the spatial-decomposition algorithm and is discussed in detail by Plimpton. 9 For FE parallelization, the mesh is divided into Np parts. This division depends on the deal. II library installation. The current MPI implementation of the PLT algorithm in CAVIAR uses PETSC [60][61][62][63] or TRILINOS 64-68 library, depending on deal. II configurations. Figure 6 illustrates three different ways of parallelization for a PLT simulation with two processors. Here, for the sake of illustration, we consider just four particles in a simulation box with a rough spherical mesh. Figure 6(a) is an illustration of sequential simulation, which is used as a reference for other cases. Here, the MD domain is shown by a black solid box, whereas the FE mesh is shown by a color-filled mesh. In Figs. 6(b) and 6(c) the MD procedure is done in parallel, while FE is sequential. In this case, all the processors have to solve the same FE mesh, but the particles are divided between processors depending on their position. Setting boundary values for FE requires having information on all the charged particles. To do that, all processes calculate the singular potential due to their dedicated charged particles on the boundary points and, then, proceed with an MPI collective summation over all processes. In Figs. 6(d) and 6(e), the FE procedure is parallel while MD is not. Therefore, the FE mesh is divided between processes while data of particles is the same for the processes. We have implemented this in a way that only the root process does the MD calculations, then it broadcasts the data to all processes. Using these complete data of the particles, each process evaluates the potential on the boundary points and then solves the partial mesh according to its implicit algorithm in the library. Subsequently, the electric field is interpolated and sent to the root process for MD integration. The parallel cases of Figs. 6(f) and 6(g) are much more complex; in the previous cases, processes have all the information of MD or FE, while here no process has all of them. Besides, MD domain decomposition does not necessarily coincide with the FE mesh decomposition. This complexity raises a parallel efficiency problem. Both cases of Figs. 6(b) and 6(c) and Figs. 6(d) and 6(e) have been implemented in CAVIAR, while the last case will appear in the next version.

IV. PACKAGE ABILITIES
The CAVIAR software package is designed to be a general MD simulator. The physical systems that can be simulated with the first version of CAVIAR are atomic or molecular MD, such as the dynamics and rheology of uncharged or charged polymer, dynamics of massive particles with gravity, electric double-layer capacitors or supercapacitors, porous media, water droplet simulations on flat or curved surfaces, and water desalination with polarized electrodes. Besides, it is possible to simulate granular materials with CAVIAR. In the granular material simulations, one can use polyhedrons to simulate realistic flat boundaries instead of rough walls that are generated by grains. In addition, the grains can be permanently charged and affected by electromagnetic fields, or they can electrostatically interact with metallic surfaces using the PLT algorithm.
As an example of the CAVIAR package, we have simulated a simple model of ionic-liquid inside a pore with an asymmetric shape. The pore has 1D periodicity in the z-axis direction. Figures 7(a) and 7(b) are snapshots of the system at the initial stage and at the t = 12 ns after running the simulation, respectively. Particles are arranged on a simple cubic lattice at the initial stage of the simulation, and the yellow/red particle color indicates the positive/negative electrical charges. The yellow/red electrode demonstrates an anode/cathode. In this example, the number of particles is 1276 and the time step is equal to Δt = 6 fs. There is a potential difference between the yellow and red boundaries in which its numeric value is equal to ΔV = 4 V, whereas the temperature is set to be equal to 400 K. To make the governing equations of motion dimensionless, we proceeded as follow: first, we set the following numeric values for fundamental parameters; the length scale is the ion diameter σ = 5 Å, the mass unit is the mass of the hexafluorophosphatem = 144 g/mol, and the energy is scaled withε = 1 kJ/mol. With these parameters, temperature scales asT = ϵ/kB and the time is dimensionless with τ = √m σ 2 /ε. The Langevin thermostat with a damping constant equal to ξ = 10 ps −1 is used to control the overall temperature of the simulation. The particles are interacting through the LJ potential, 28 which is cut at r cutoff = 2 (1/6) σ and shifted by u 0 to fix the discontinuity due to cutting. This assures a fully repulsive LJ potential.

ARTICLE scitation.org/journal/adv
In addition to the LJ, electrostatic interactions are also taken into account. There is another pure repulsive LJ interaction with ϵ wall = 1 kJ/mol and σ wall = 3.37 Å between the particles and the electrode walls (Sec. II B). Figure 7(d) is a plot of the generated FE mesh to compute the smooth electrostatic potential. The color-bar at the top of this mesh indicates the voltage difference across the mesh. It is worth noticing that the numeric values of this color-bar belongs to the dimensionless potential and differs from the SI values. In this simulation, 1 V equals 96.4923 in dimensionless units. Figure 8(a) is a plot of the dynamical behavior of induced charge density on different boundaries. These boundaries are schematically plotted in Fig. 8(b) and distinguished using different colors. A unique ID has been set to each boundary. All boundaries on the left (ID < 6) are equipotential surfaces and have a potential difference with that on the right (ID = 6). Figure 8(a) shows that nearly 15 ns is required to charge the capacitor at 4 V. At the equilibrium stage, different areas of the positive electrode have the same value of induced charge density. One can see in Fig. 7(b) that the positive (yellow) particles are attracted to the negative (red) electrode and vice-versa. This simulation is provided in the package's user guide as an example.

V. PACKAGE PERFORMANCE
In this section, some performance facts of CAVIAR simulation with the PLT algorithm are presented. As for the only available alternative to the PLT algorithm, the ICC * algorithm is used to check the PLT performance. Currently, ICC * is only implemented in the ESPResSo software package. 69 At the time of writing this work, the latest release of ESPResSo is its version 4.1.2. In this version, ICC * can only be used alongside a P3M Laplace solver. However, the CAVIAR-1.0 does not have P3M yet; instead, it has a basic Ewald algorithm. The ICC * algorithm only sets the zero potential on the electrodes. To apply a potential difference in ICC * simulations, one needs to create a mesh for solving the Laplace equation with the desired potential as its boundary condition. In order to make the comparison straightforward, we conducted all simulations in zero potential difference.
The following simulations are done in dimensionless units. The simulation domain is a cubic box of length L = 20, in which the planes Z = 20 and Z = 0 represent a metallic surface. In the case of the PLT algorithm, these surfaces are a part of the FE mesh that has been indicated with zero potential. The initial FE mesh is a simple cube of length 20 that encounters the whole simulation box, and four levels of its global refinements are used for the simulations (Fig. 9). For each level of the refinement, an amount of 200 time steps is averaged and recorded as the central processing unit (CPU)-time for an MD step. In the case of the ICC * algorithm, the electrode surfaces are made of static charged particles that are called ICC charges. The simulations are done for different numbers of ICC charges, and the CPU-time for one MD time step is calculated in the same way as before.
The simulation box contains some freely moving charges. We have set up two configurations for the simulations as follows:  However, in the PLT algorithm, we do not need to set particles to generate walls, but it is needed in the ICC * method. It should be noted that the optimal number of the mesh refinement needs to be investigated for each problems before long production runs.
The CAVIAR software package is developed to assist the scientists and researchers in their research. It is a new package but it has some advantages that can make it interesting for the users.
• It is easy to install on any computer, it has the least number of external dependency. • It is a general purpose MD software, from MD atomic simulations up to granular and gravitational particles. • It has a particle-geometry interaction algorithm. Comparing with other algorithms, there is no need of particle distribution on the boundaries. This makes simulations of complicated geometries possible by utilizing computer-aideddesign (CAD) software. • CAVIAR has the capability to simulate any kind of boundary conditions, including non-periodic or periodic in one or more dimensions. • It has an automatic initial particle distribution generator, which is able to generate any desired distribution inside the simulation domain or in a given polyhedron. It also can mix multiple polyhedrons by using boolean operators to generate complex initial particle distributions. • Some of the MD constraint algorithms are designed only for one type of the MD integrator, i.e., depending on the implemented integrator on the MD software, only some of the constraint can be added. As an example, one can refer to the SHAKE and RATTLE algorithm, both of which are used to simulate water molecules. The former is compatible with verlet and the latter uses velocity-verlet for integration. 57 The CAVIAR is designed completely modular, so the choice of the integrator is user dependent and according to it, the suitable constraint can be selected by the user. • It is easy to modify the CAVIAR code by the researchers and developers to add new features. It is completely objectoriented, and it has a developer-friendly software design. • The CAVIAR structure has been completely re-designed many times to ensure a flexible and long-life code. • It can be used as a C++ library in other codes.

VI. CONCLUSION
In this paper, the PLT algorithm for the MD simulations of charged particles under metallic boundary conditions is introduced. This algorithm was implemented with an FE Poisson solver; then, a test simulation for a point charge inside a conductive spherical shell is conducted. It was shown that PLT is highly accurate and the errors can easily be reduced with mesh refinement. In addition, another algorithm for the calculation of the interactions of MD particles with polyhedrons is described. With this algorithm, the polyhedrons can be used as MD walls. The advantage of these algorithms is that there is no need to put static MD particles as walls, which reduces the computation cost. It was shown that calculating the initial position of particles inside complex geometries can be done with a ray-casting algorithm.
These algorithms, as well as other standard MD methods, are implemented in the CAVIAR software package. This package can be used for a wide range of MD simulations. An example of its usage is given in this paper. CAVIAR is an object-oriented and highly extensible package. The package's basic structure and the implemented methods are described in this paper. Different parallelization schemes for mixing FE and MD are also introduced.
This package has the capability to utilize other Poisson solving methods, such as the boundary element (BE) method. The BE method has higher efficiency of solving the Poisson equation in some special geometries. This feature will be added to the next version of the CAVIAR package. In addition, novel methods of the Poisson equation solving such as Isogeometric Formulation (IF) can be utilized in the PLT algorithm. 70,71 In the IF method, instead of discretizing the boundaries into a polyhedron mesh, a continuous description of the boundary surfaces is used, 72 which in turn helps to improve the performance of the PLT algorithm.