Reflecting boundary conditions for classical molecular dynamics simulations of nonideal plasmas

The influence of boundary conditions on results of the classical molecular dynamics simulations of nonideal electron-ion plasma is analyzed. A comprehensive study is performed for the convergence of per-particle potential energy and pressure with the number of particles using both the nearest image method (periodic boundaries) and harmonic reflective boundaries. As a result an error caused by finiteness of the simulation box is estimated. Moreover the electron oscillations given by the spectra of the current autocorrelation function are analyzed for both types of the boundary conditions. Nonideal plasmas with the nonideality parameter in range 0.26-2.6 is considered. To speed up the classical molecular dynamics simulations the graphics accelerators code is used.


Introduction
Theoretical studies of strongly coupled Coulomb systems such as electron-ion nonideal plasma [1][2][3][4][5], ion-ion plasma [6], warm dense matter [7][8][9], ion liquids [10], dusty plasmas [11], etc. often rely on atomistic simulations by the methods of classical molecular dynamic (MD) and Monte-Carlo (MC). In this paper we concentrate on the MD simulations of fully ionized nonideal electron-ion plasmas, although most of our results are still valid for other types of the above mentioned systems of charged particles.
The general feature of MD simulation of plasma is long-range Coulomb interaction. It determines the time and distance scales for MD plasma simulations that are several orders of magnitude lower than those for atomic liquids and condensed matter. In contrast to simulations of atomic systems where times and distances measured in the scales of nanoseconds and micrometers, the plasma simulations are limited by picoseconds and tens of nanometers. Nevertheless, there are plasma objects that fit within these scales, for example, a nanoplasma created by rapid ionization of nanosized atomic clusters by femtosecond laser pulses [12,13]. In this case no boundaries are needed except for those at a very large distance to avoid run away particles. Simulation of Coulomb system require much more computational effort due to several reasons. First, a relatively small timestep (10 −19 -10 −17 s) is needed for an explicit electron dynamics. In the classical MD electrons are represented as point-like charges and the electron-ion interaction is given by a pseudopotential which accounts for quantum effects to a certain extent. In the case of the wave packet molecular dynamics (WPMD) [14] or electron force field (eFF) [8] simulations, the electron wave function is represented by a single Gaussian wave packet with changing width. This provide better description of electron-ion bound states, ionization and recombination processes. The use of multiple Gaussians per electron (SplitWP method) [15][16][17] increases the accuracy of the low energy bound states and allows one to simulate quantum branching (e.g., tunneling).
The second reason for the slowdown of plasma simulations compared to atomic systems is related to the need for a more accurate account for the Coulomb interaction that is responsible for electron correlations, plasma waves, transport properties, etc. The approximate electrostatic solvers such as the particle-particle-particle-mesh (PPPM) are used to avoid this problem in atomic simulation. Unfortunately, this method is hardly applicable to plasma studies. Although there exist techniques such as MicPIC [18,19] and TreeMD [20] which reduce computational time for large systems (over 10 5 particles), many simulations still rely on the full "all-to-all" interaction model.
Direct MD simulations of larger objects are not feasible and require a multiscale approach where MD is used for obtaining properties of a small system fragment which can be used then in a higher level hydrodynamic codes or semi-empirical models. The computational time in classical MD grows as t simul ∼ N 2 ∼ L 6 , where N is the number of particles and L is the system size parameter. In this case the MD system size is to be taken as small as possible and governed by a characteristic scale of the phenomena under study. One should note, that there is no reason to increase the system size for better averaging, as it can be replaced by averaging over time for an equilibrium MD trajectory or over ensemble of nonequilibrium runs from a macroscopically identical initial states [21].
A typical choice for MD simulations of an extended matter is periodic boundary conditions (PBC) [22,23] where each particle interacts either with the nearest images of other particles (nearest image method) or with all periodic images within the Ewald summation approach [22,24]. The latter method is more appropriate for ordered systems like crystals, although sometimes it is used as well for a disordered nonideal plasma [2]. In our view, the applicability of the Ewald summation for a disorderd system is doubtful and needs additional argumentation. In some works [25,26] the authors insist on using the reflecting walls (RW) for MD simulations of nonideal plasmas. In this paper we compare the identical systems under the PBC provided by nearest image method (NI) and RW conditions and reveal possible artefacts caused by each type of boundary conditions. All types of boundary conditions introduce an implicit cut off of the interparticle interaction potential and in system with long-range interaction this phenomena can be critical. Fortunately in the case of two-component (electron-ion) disordered plasmas the potential is effectively screened and seen by particles as a Debye-type interaction. It should be noted that the Debye interaction is not introduced by the MD model but arises naturally in simulations due to the particle correlation build up.
In determining a suitable model system size from the screening length r scr one has, however, to take into account that the screening length is underestimated [3] by the standard Debye formula r D = (k B T ) 1/2 (4πn e e 2 ) −1/2 in the case of the nonideal plasma with the nonideality parameter Γ = e 2 (4πn e /3) 1/3 (k B T e ) −1 0.1, where e and n e are the electron charge and the electron number density (the case of singly charged ions), T is the temperature. For Γ > 1 it is close to the mean interparticle distance r S ≈ (3Γ) 1/2 r D . In WPMD simulations there is another issue related to the boundary conditions. An application of this method to rather dilute nondegenerate extended systems leads to the known problem of unrestricted wave packet spreading [27]. Although this spreading is a natural effect for the Gaussian wave function in an unrestricted space, the model of an extended system should also provide correct limitations on the density of allowed quantum states, which remains arbitrary when the wave packets are allowed to spread. The reflective boundary conditions for simulation box, which we investigate in this work, may be used as a direct method to provide the finiteness of the partition function. The paper is organized as follows. In Section 2 we describe our molecular dynamics simulation method for the confined systems both in classical and quantum cases. In Section 3 we present and discuss our simulation results. In particular, for the classical system the convergence of the internal energy and pressure is studied as a function of the number of particles and excitations of the system are analyzed. The final section contains conclusions and suggestions for using confined systems in extended medium simulations.

Simulation method 2.1. Classical molecular dynamics
We consider a classical system of electrons and singly charged ions. The Hamiltonian is given by where K i and K e are the kinetic energies of ions and electrons and U = U ii + U ei + U ee is the potential energy N i and N e are the numbers of ions and electrons (N = N i +N e ), R i and q i are the ion coordinates and charges, r k are the electron coordinates, σ ei = σ ee = 0.318 nm are the parameters responsible for the potential correction at short distances, U ext is an external field potential which may include wall boundaries as described below. This form of the potential was used is a number of classical nonideal plasma simulations [13,28,29], in particular in the simulations of ionized sodium nanoclusters [13,29]. The same form of the potential is obtained in WPMD simulations with the frozen wave packet width. The ion to electron mass ratio is taken to be 100. The real ion to electron mass ratio is lower but as shown earlier [3], setting such a low mass ratio does not influence the electron dynamics and thermodynamical properties. This allow us to integrate motion equation for ion and electron with same time step 2 × 10 −18 s. The equations of motion are integrated using the second order Leap-frog scheme.
Simulation start from a random particle distribution and apply the Langevin thermostat to equilibrate the system with a given temperature T 0 . A special procedure is used to check if the equilibrium is established. Then the thermostat is switched off and the main MD run is performed for an adiabatic system. In this paper we report on the results for T 0 = 3 × 10 4 K. We considered systems with electron densities n e = 2.42 × 10 19 cm −3 , 2.89 × 10 21 cm −3 , and 2.42 × 10 22 cm −3 which correspond to the nonideality parameters of Γ = 0.26, 1.28 and 2.6, θ = 2 (3π 2 n e ) 2/3 /(2m e k B T ) = 0.011, 0.29 and 1.19 respectively. The equilibrium trajectory length ranges from 5×10 6 steps (10 ps) for N i = 16 to 1.32×10 5 steps (0.27 ps) for N i = 32000. To increase accuracy of the measuments the results are averaged over 3− 6 statistically independent trajectories.

Equilibration and obtaining the required temperature
Running the main part of simulations without using a thermostat avoids an unphysical distortion of the velocity distribution function. The drawback is that the the actual mean temperature T actual obtained by averaging over the trajectory may slightly differ from the target temperature T 0 . In our simulations this difference was up to 5% (drift of the total energy was less then 0.2%). As it could influence precise calculations of the thermodynamic quantities we determined the derivative of a required quantity A over the temperature using a series of simulations and then used the correction:  In this way, we applied the corrections to the mean per-particle potential energy and pressure. The temperature dependencies in the vicinity of T 0 can be approximated by linear equation ( figure 1) with a good accuracy as given below for Γ = 1.

Boundary conditions
The simulations are performed for two types of the boundary conditions: the periodic boundaries within the nearest image approach and the reflecting walls. The pure reflecting walls means immediate reflection of the particle when it crosses the wall and it implementation without essential sophistication of the integration scheme causes violation of the total energy conservation. Instead of it we use a softer harmonic-like wall-particle interaction potential for each direction: where L = (N e /n e ) 1/3 is the simulation box length, a is a steepness parameter, r kα is the α-component of k-th particle coordinate vector. The dependence of the results for thermodynamical properties on the steepness parameter a is shown in figure 2. It is seen that the results converge well with the increase of the wall steepness. With respect to this dependence we took a = 267 J/m 2 for all simulations discussed below, that provide thermodynamical value close to their limits and at the same time provides acceptable total energy conservation for the Leap-Frog scheme. As seen from figure 2a the internal region is slightly affected by the boundary parameters.  For precise measurements one should take into account an effective expansion of the system and lowering of the mean particle density compared to PBC where the system volume is fixed. An extra volume caused by the transient presence of particles behind the wall can be eliminated by shifting the walls inwards the cell: This shift is shown in figure 3b as the difference between solid and dashed vertical lines. As example, for the plasma parameters given above, this volume correction changes the mean potential energy by 1.2% for N i = 32000. Due to finiteness of the simulated system, a layer of particles with different energy is formed near boundaries. This effect is observed on the energy distribution (figure 3). The higher is the plasma density the greater is the mean energy difference between the boundary layer and the rest plasma body (see figure 3b). The width of the boundary layer is about the screening lenght r scr for all Γ. One should note, that the boundary layer width slightly changes with the system size.
To achieve the required accuracy of thermodynamical values, long MD trajectories are needed. As the long range Coulomb interactions embarrasses standard parallel techniques such as the domain decomposition, we run simulations on graphics accelerators (GPU) with particle (force) decomposition for speed up calculation. To increase the GPU efficiency we use several threads per particle and a partial force reduction in the GPU shared memory.

Boundary effects to the mean potential energy and pressure
In this paper we study the convergence of the mean potential energy and pressure with the system size growth in the classical MD with different types of the boundary conditions ( figure 4). The results for the plasma with three different values of the nonideality parameter are considered. Note, that the results for Γ = 1.28 and Γ = 2.6 are obtained with better accuracy due to better averaging and more pronounced boundary effects.
The pressure is obtained using the common virial equation [5,23]: where V = L 3 is the system volume and f ij are interparticle forces. As seen from the figure in the case of PBC (black circles) the energy and pressure converge to their limiting values U ∞ and P ∞ faster than in the case of reflecting walls (triangles). It is also seen that due to screening one can use rather small systems for nonideal plasma simulations with acceptable accuracy of the obtained thermodynamical properties. For instance, the error in pressure determination at Γ = 1.28 is only 1.0% for N i = 16 and 0.2% for N i = 128, the error of the mean potential energy is 2.7% for N i = 16 and 0.5% for N i = 128. It opens a way to use more complicated and time consuming simulation techniques like AWPMD where the number of particles becomes a crucial parameter. As stated above, the confined system is preferred for WPMD simulations because it avoids the wave packet spreading problem. At the same time the boundaries produce a layer with modified per-particle potential energy as seen in figure 3 where the energy distribution along a single axis is visualized. The width of this layer is governed by the screening length r scr ≈ 0.5 nm that is almost independent of the system size. Thus the boundary effects should scale as N −1/3 which is demonstrated by inverse cubic root fits to the deviations of energy and pressure from their limiting values (see triangles and solid curves in figure 4): here ζ and η are fitting parameters. Moreover if the potential energy is calculated only for the internal region of the cell, it gives even better results for RW than for PBC (see rhombus in figures 4a, 4c, 4e).

Electron oscillation modes for periodic and constrained systems
In addition to the static properties we compared the spectrum of electron oscillations for both types of boundary conditions using the current autocorrelation functions (ACF) [3,4]: where v i are particle velocities and J is the total current density. According to the linear response theory [30] the current ACF is related to the dynamical plasma conductivity σ(ω) = (V k B T ) −1 J; J ω . As seen from figure 5a, using PBC leads to a monotonic spectrum (dashed curve) with no resonances in accordance with earlier simulations [3,4,31]. When the wall boundaries are applied the spectra is changed substantially demonstrating well pronounced resonances for both Γ. These resonances are related to the Mie-like oscillations that are known, for example, for ionized clusters [12,13]. They can not be observed for a periodic system as a shift of all electrons does not produce charge separation and a corresponding returning force. In the constrained system the charge separation appears at the boundaries and this fundamental mode of electron oscillations do exist. As pointed out in [31] in this case the value of J; J ω is related to the so called external conductivity in contrast to the internal conductivity for PBC results. Spectrum for Γ = 0.26 are not presented due to the lack of statistics.  The standard equation for the Mie frequency reads ω Mie = (4πe 2 n i /(3m e )) 1/2 = ω pl / √ 3 where ω pl is the plasma frequency and n i = n e is the ion number density. This equation has been obtained for a spherical system while we a using a cubic box cell. The Mie-like oscillation for the cubic geometry are studied in [32]. There are more than one resonance in this case and the frequencies of the strongest resonances are shown in figure 5a by dashed vertical lines (0.43 ω pl , 0.50 ω pl and 0.79 ω pl ). They are in agreement with the MD results.
Another type of electron oscillations in plasmas is the Langmuir plasma waves [3,33]. In a finite system the cell size should contain an integer number of the wavelengths λ and therefore such oscillations are not seen in the total current ACF spectra (positive and negative contributions of electron current cancel each other). However, if the system is large enough the current ACF can be calculated for a subcell (figure 5b). In this case the plasma waves are revealed near ω pl for both PBC and RW cases. As the influence of the boundaries is smaller for the subcell the ACF spectra look closer to each other. For the different nonideality parameters size of inner part where Mie-like oscillation dumped are different.

Conclusions
Classical molecular dynamics simulations were performed for the nonideal singly charged electron-ion plasma with different types of boundary conditions and different nonideality parameters. It was shown that the periodic boundary conditions (the nearest image method) provide the fastest convergence of the calculated mean potential energy and pressure with the system size growth in studied range of systems. The results for an extremely small system of 16 ions and 16 electrons deviate only within 0.5 − 1% from the limiting values. It enables to use small systems for advanced quantum simulation methods like AWPMD where large numbers of particles are hard to reach.
For simulations with the reflecting boundaries the boundary layer was found to be of the order of the mean interparticle distance. It leads to an additional error proportional to N −1/3 for both thermodynamical quantities. Moreover an effective expansion of the system and the mean particle density decrease due to soft boundaries should be taken into account. The use of a non-periodic system also leads to the appearance of Mie-like oscillations seen at the current ACF spectrum. It affects computing of the dynamical conductivity and similar properties. The ACF spectra for an internal subcell are close to each other and reveal the Langmuir plasma oscillations. With all these corrections we can conclude that the results for the reflecting walls and periodic boundaries are in a good agreement.