Particle Mesh Ewald Method Research Articles

PME Switching in Confinement Simulations of Charged Solutes.

The confinement method is a reaction coordinate-free enhanced sampling method for the calculation of conformational free energy differences. We show that in explicit solvent, artifacts occur when treating charged solutes. These artifacts are resolved by switching off the particle mesh Ewald (PME) method at the start of confinement. Calculations of the free energy cost of this switching converge rapidly, with small statistical error. The effectiveness and accuracy of confinement with PME switching is demonstrated by its application to a series of solutes of different charge; its ability to treat complex systems is illustrated by evaluating the free energy difference between B and Z-DNA.

Molecular dynamics simulation of the effect of temperature on the conformation of ubiquitin protein.

Ubiquitin, a ubiquitous small protein found in all living organisms, is crucial for tagging proteins earmarked for degradation and holds pivotal importance in biomedicine. Protein functionality is intricately linked to its structure. To comprehend the impact of diverse temperatures on ubiquitin protein structure, our study delved into the energy landscape, hydrogen bonding, and overall structural stability of ubiquitin protein at varying temperatures. Through meticulous analysis of root mean square deviation and root mean square fluctuation, we validated the robustness of the simulation conditions employed. Within our simulated system, the bonding energy and electrostatic potential energy exhibited linear augmentation, while the van der Waals energy demonstrated a linear decline. Additionally, our findings highlighted that the α-Helix secondary structure of the ubiquitin protein gradually transitions toward helix destabilization under high-temperature conditions. The secondary structure of ubiquitin protein experiences distinct changes under varying temperatures. The outcomes of our molecular simulations offer a theoretical framework that enhances our comprehension of how temperature impacts the structural stability of ubiquitin protein. These insights contribute not only to a deeper understanding of iniquity's behavior but also hold broader implications in the realm of biomedicine and beyond. All the MD simulations were performed using the GROMACS software with GROMOS96 force field and SPC for water. The ubiquitin protein was put in the center of a cubic box with a length of 8nm, a setting that allowed > 0.8nm in the minimal distance between the protein surface and the box wall. To remove the possible coordinate collision of the configurations, in the beginning, the steepest descent method was used until the maximum force between atoms was under 100kJ/mol·nm with a 0.01nm step size. Minimization was followed by 30ps of position-restrained MD simulation. The protein was restrained to its initial position, and the solvent was freely equilibrated. The product phase was obtained with the whole system simulated for 10ns without any restraint using an integral time step of 1fs with different temperatures. The cutoff for short-range electronic interaction was set to 1.5nm. The long-range interactions were treated with a particle-mesh Ewald (PME) method with a grid width of 1.2nm.

Biopolymers as a Potential Alternative for the Retention of Pollutants from Vinasse: An In Silico Approach.

Vinasse, a waste from the bioethanol industry, presents a crucial environmental challenge due to its high organic matter content, which is difficult to biodegrade. Currently, no sustainable alternatives are available for treating the amount of vinasse generated. Conversely, biopolymers such as cellulose, carboxymethylcellulose, and chitosan are emerging as an interesting alternative for vinasse control due to their flocculating capacity against several organic compounds. This study seeks to determine the thermodynamic behavior of in silico interactions among three biopolymers (cellulose, carboxymethylcellulose, and chitosan) regarding 15 organic compounds found in vinasse. For this, the Particle Mesh Ewald (PME) method was used in association with the Verlet cutoff scheme, wherein the Gibbs free energy (ΔG) was calculated over a 50 ns simulation period. The findings revealed that cellulose showed a strong affinity for flavonoids like cyanidin, with a maximum free energy of -84 kJ/mol and a minimum of -55 kJ/mol observed with phenolic acids and other flavonoids. In contrast, chitosan displayed the highest interactions with phenolic acids, such as gallic acid, reaching -590 kJ/mol. However, with 3-methoxy-4-hydroxyphenyl glycol (MHPG), it reached an energy of -70 kJ/mol. The interaction energy for flavonoid ranged from -105 to -96 kJ/mol. Finally, carboxymethylcellulose (CMC) demonstrated an interaction energy with isoquercetin of -238 kJ/mol, while interactions with other flavonoids were almost negligible. Alternatively, CMC exhibited an interaction energy of -124 kJ/mol with MHPG, while it was less favorable with other phenolic acids with minimal interactions. These results suggest that there are favorable interactions for the interfacial sorption of vinasse contaminants onto biopolymers, indicating their potential for use in the removal of contaminants from the effluents of the bioethanol industry.

Enhanced sampling in explicit solvent by deep learning module in FSATOOL.

FSATOOL is an integrated molecular simulation and data analysis program. Its old molecular dynamics engine only supports simulations in vacuum or implicit solvent. In this work, we implement the well-known smooth particle mesh Ewald method for simulations in explicit solvent. The new developed engine is runnable on both CPU and GPU. All the existed analysis modules in the program are compatible with the new engine. Moreover, we also build a complete deep learning module in FSATOOL. Based on the module, we further implement two useful trajectory analysis methods: state-free reversible VAMPnets and time-lagged autoencoder. They are good at searching the collective variables related to the conformational transitions of biomolecules. In FSATOOL, these collective variables can be further used to construct a bias potential for the enhanced sampling purpose. We introduce the implementation details of the methods and present their actual performances in FSATOOL by a few enhanced sampling simulations.

Drude Polarizable Lipid Force Field with Explicit Treatment of Long-Range Dispersion: Parametrization and Validation for Saturated and Monounsaturated Zwitterionic Lipids.

Accurate empirical force fields of lipid molecules are a critical component of molecular dynamics simulation studies aimed at investigating properties of monolayers, bilayers, micelles, vesicles, and liposomes, as well as heterogeneous systems, such as protein-membrane complexes, bacterial cell walls, and more. While the majority of lipid force field-based simulations have been performed using pairwise-additive nonpolarizable models, advances have been made in the development of the polarizable force field based on the classical Drude oscillator model. In the present study, we undertake further optimization of the Drude lipid force field, termed Drude2023, including improved treatment of the phosphate and glycerol linker region of PC and PE headgroups, additional optimization of the alkene group in monounsaturated lipids, and inclusion of long-range Lennard-Jones interactions using the particle-mesh Ewald method. Initial optimization targeted quantum mechanical (QM) data on small model compounds representative of the linker region. Subsequent optimization targeted QM data on larger model compounds, experimental data, and dihedral potentials of mean force from the CHARMM36 additive lipid force field using a parameter reweighting protocol. The use of both experimental and QM target data during the reweighting protocol is shown to produce physically reasonable parameters that reproduce a collection of experimental observables. Target data for optimization included surface area/lipid for DPPC, DSPC, DMPC, and DLPC bilayers and nuclear magnetic resonance (NMR) order parameters for DPPC bilayers. Validation data include prediction of membrane thickness, scattering form factors, electrostatic potential profiles, compressibility moduli, surface area per lipid, water permeability, NMR T1 relaxation times, diffusion constants, and monolayer surface tensions for a variety of saturated and unsaturated lipid mono- and bilayers. Overall, the agreement with experimental data is quite good, though the results are less satisfactory for the NMR T1 relaxation times for carbons near the ester groups. Notable improvements compared to the additive C36 force field were obtained for membrane dipole potentials, lipid diffusion coefficients, and water permeability with the exception of monounsaturated lipid bilayers. It is anticipated that the optimized polarizable Drude2023 force field will help generate more accurate molecular simulations of pure bilayers and heterogeneous systems containing membranes, advancing our understanding of the role of electronic polarization in these systems.

Quantifying electrostatics from (bio)molecular dynamics trajectories using g-elpot

Electrostatics plays a pivotal role in many physiological processes, defining ion selectivity of ion channels, driving formation of supramolecular complexes, and determining energetics of transmembrane transporters, to name a few. Atomistic molecular dynamics (MD) simulations can serve as a rich source of structural and dynamic information on such processes. However, tools for robust high-resolution quantification of biomolecular electrostatics from MD trajectories are still limited. Here we present g_elpot, a GROMACS-based tool for calculating electrostatic potentials from MD simulations using the Smooth Particle Mesh Ewald (SPME) method.

A hierarchical matrix approach for computing hydrodynamic interactions

For simulations of large numbers of small, spherical particles in a Stokes flow, the long-range hydrodynamic interactions approximated by the Rotne–Prager–Yamakawa (RPY) kernel can be summed rapidly using, for example, the fast multipole method (FMM) or the particle-mesh Ewald (PME) method. In this paper, we develop new fast methods for computing these sums using the H2 hierarchical matrix representation, for open and for periodic boundary conditions. To the best of our knowledge, the method for infinite periodic sums using the H2 hierarchical matrix representation is the first such method developed. We also consider a more general RPY kernel that handles polydisperse particle radii, and show analytically and experimentally that the proxy surface method for efficiently constructing the H2 hierarchical matrix representation remains effective in this case. Numerical tests demonstrate the well-controlled accuracy of the H2 summation methods and their linear-scaling computation and storage cost. We find that the H2 matrix approach has lower cost for computing the summations compared to FMM and PME, but higher precomputation cost (required for each particle configuration). This precomputation cost can be amortized over several summations when computing Brownian displacements or forces in Brownian and Stokesian dynamics simulations with very large numbers of particles.

Validating the CHARMM36m protein force field with LJ-PME reveals altered hydrogen bonding dynamics under elevated pressures

The pressure-temperature phase diagram is important to our understanding of the physics of biomolecules. Compared to studies on temperature effects, studies of the pressure dependence of protein dynamic are rather limited. Molecular dynamics (MD) simulations with fine-tuned force fields (FFs) offer a powerful tool to explore the influence of thermodynamic conditions on proteins. Here we evaluate the transferability of the CHARMM36m (C36m) protein force field at varied pressures compared with NMR data using ubiquitin as a model protein. The pressure dependences of J couplings for hydrogen bonds and order parameters for internal motion are in good agreement with experiment. We demonstrate that the C36m FF combined with the Lennard-Jones particle-mesh Ewald (LJ-PME) method is suitable for simulations in a wide range of temperature and pressure. As the ubiquitin remains stable up to 2500 bar, we identify the mobility and stability of different hydrogen bonds in response to pressure. Based on those results, C36m is expected to be applied to more proteins in the future to further investigate protein dynamics under elevated pressures.

G_elpot: A Tool for Quantifying Biomolecular Electrostatics from Molecular Dynamics Trajectories.

Electrostatic forces drive a wide variety of biomolecular processes by defining the energetics of the interaction between biomolecules and charged substances. Molecular dynamics (MD) simulations provide trajectories that contain ensembles of structural configurations sampled by biomolecules and their environment. Although this information can be used for high-resolution characterization of biomolecular electrostatics, it has not yet been possible to calculate electrostatic potentials from MD trajectories in a way allowing for quantitative connection to energetics. Here, we present g_elpot, a GROMACS-based tool that utilizes the smooth particle mesh Ewald method to quantify the electrostatics of biomolecules by calculating potential within water molecules that are explicitly present in biomolecular MD simulations. g_elpot can extract the global distribution of the electrostatic potential from MD trajectories and measure its time course in functionally important regions of a biomolecule. To demonstrate that g_elpot can be used to gain biophysical insights into various biomolecular processes, we applied the tool to MD trajectories of the P2X3 receptor, TMEM16 lipid scramblases, the secondary-active transporter GltPh, and DNA complexed with cationic polymers. Our results indicate that g_elpot is well suited for quantifying electrostatics in biomolecular systems to provide a deeper understanding of its role in biomolecular processes.

Multilevel summation for periodic electrostatics using B-splines.

Fast methods for calculating two-body interactions have many applications, and for molecular science and cosmology, it is common to employ periodic boundary conditions. However, for the 1/r potential, the energy and forces are ill-defined. Adopted here is the model given by the classic Ewald sum. For the fast calculation of two-body forces, the most celebrated method is the fast multipole method and its tree-code predecessor. However, molecular simulations typically employ mesh-based approximations and the fast Fourier transform. Both types of methods have significant drawbacks, which, in most respects, are overcome by the less well-known multilevel summation method (MSM). Presented here is a realization of the MSM, which can be regarded as a multilevel extension of the (smoothed) particle mesh Ewald (PME) method, but with the Ewald softening replaced by one having a finite range. The two-level (single-grid) version of MSM requires fewer tuning parameters than PME and is marginally faster. Additionally, higher-level versions of MSM scale well to large numbers of processors, whereas PME and other two-level methods do not. Although higher-level versions of MSM are less efficient on a single processor than the two-level version, evidence suggests that they are more efficient than other methods that scale well, such as the fast multipole method and tree codes.

Analytical Hessians for Ewald and particle mesh Ewald electrostatics.

The particle mesh Ewald (PME) method has become ubiquitous in the molecular simulation community due to its ability to deliver long range electrostatics accurately with ON ⁡log(N) complexity. Despite this widespread use, spanning more than two decades, second derivatives (Hessians) have not been available. In this work, we describe the theory and implementation of PME Hessians, which have applications in normal mode analysis, characterization of stationary points, phonon dispersion curve calculation, crystal structure prediction, and efficient geometry optimization. We outline an exact strategy that requires O(1) effort for each Hessian element; after discussing the excessive memory requirements of such an approach, we develop an accurate, efficient approximation that is far more tractable on commodity hardware.

CHARMM36 Lipid Force Field with Explicit Treatment of Long-Range Dispersion: Parametrization and Validation for Phosphatidylethanolamine, Phosphatidylglycerol, and Ether Lipids.

Long-range Lennard-Jones (LJ) interactions have been incorporated into the CHARMM36 (C36) lipid force field (FF) using the LJ particle-mesh Ewald (LJ-PME) method in order to remove the inconsistency of bilayer and monolayer properties arising from the exclusion of long-range dispersion [Yu, Y.; Semi-automated Optimization of the CHARMM36 Lipid Force Field to Include Explicit Treatment of Long-Range Dispersion. J. Chem. Theory Comput. 2021, 10.1021/acs.jctc.0c01326. (preceding article in this issue)]. The new FF is denoted C36/LJ-PME. While the first optimization was based on three phosphatidylcholines (PCs), this work extends the validation and parametrization to more lipids including PC, phosphatidylethanolamine (PE), phosphatidylglycerol (PG), and ether lipids. The agreement with experimental structure data is excellent for PC, PE, and ether lipids. C36/LJ-PME also compares favorably with scattering data of PG bilayers but less so with NMR deuterium order parameters of 1,2-dimyristoyl-sn-glycero-3-phospho-(1'-rac-glycerol) (DMPG) at 303.15 K, indicating a need for future optimization regarding PG-specific parameters. Frequency dependence of NMR T1 spin-lattice relaxation times is well-described by C36/LJ-PME, and the overall agreement with experiment is comparable to C36. Lipid diffusion is slower than C36 due to the added long-range dispersion causing a higher viscosity, although it is still too fast compared to experiment after correction for periodic boundary conditions. When using a 10 Å real-space cutoff, the simulation speed of C36/LJ-PME is roughly equal to C36. While more lipids will be incorporated into the FF in the future, C36/LJ-PME can be readily used for common lipids and extends the capability of the CHARMM FF by supporting monolayers and eliminating the cutoff dependence.

Recent Updates to the Charmm Lipid Force Fields

Since its release in the early 1990s, the CHARMM lipid force field (FF) has been updated frequently to agree with experimental data and to work with state-of-art simulation protocols. The most recent all-atom version, C36, has been heavily used due to the diversity of lipids it covers and the well-designed parameters, but it is not suitable for monolayer simulations due to the inconsistency of surface tension with experiment. To fix this problem, we introduced the long-range Lennard-Jones (LJ) interactions using the LJ Particle Mesh Ewald (LJ-PME) method and optimized the C36 parameter set to work with LJ-PME.

A GPU-Accelerated Fast Multipole Method for Gromacs. Performance and Accuracy

The evaluation of Coulomb forces is the computationally most expensive task in molecular dynamics (MD) simulations. Most simulation codes rely on the Particle Mesh Ewald (PME) method to evaluate these forces. This is for a good reason - PME is extremely fast. However, due to the underlying Fourier transformations, PME suffers from a communication bottleneck which prohibits efficient scaling to large numbers of compute nodes. The fast multipole method (FMM) is an alternative to PME that by construction needs significantly less communication volume in parallel.

(PEFC&E 2nd Place Best Poster Winner) Molecular Analysis of Cerium Ion Transport Properties in Polymer Electrolyte Membrane

In order to improve the durability of Polymer electrolyte fuel cell (PEFC), to clarify the cause of deterioration and take measures is necessary. Chemical degradation of polymer electrolyte membrane (PEM) is considered as one of the causes of degradation of PEFC. Hydroxyl radical (OH ∙) is generated by decomposition of hydrogen peroxide generated by oxygen reduction reaction (ORR) during PEFC operation. This OH ∙ has high reactivity, and OH ∙ would degrade the polymer membrane. Based on this degradation mechanism, a method of adding a substance that inactivates OH ∙ before reacting with the polymer (radical scavenger) has been put into practical use. One of the most useful radical scavengers is cerium ions. Cerium ions in the PEM have been reported to migrate during operation and the distribution becomes non-uniform, and degradation occurs at points where there are only a few cerium ions in the membrane. Understanding the cerium ion transport mechanism in PEM is important for controlling cerium ion migration. Because investigating the cerium ion transport in operating PEFC is limited in experiments, nanoscale analysis by numerical simulation is useful to understand the properties of cerium ion transport. In this study, we analyze the distribution and transport properties of cerium ions in the membrane by modeling the PEM containing cerium ions using molecular dynamics (MD) simulations. Diffusion due to the concentration gradient, migration due to the ionomer potential gradient, and convective ion transport due to water crossover are considered as the transport factors of the cerium ions in the membrane. In the present study, the effect of water contents and water flux on the convection ion transport is analyzed. Simulation system is composed of polymer chains, cerium ions, and solvent molecules (water molecules and hydronium ions). Nafion® model was used as polymer and the equivalent weight (EW) was set at 1146. The potential based on the DREIDING force field was adopted for the molecular interaction of Nafion®, and particle mesh Ewald method was used to calculate Coulomb force.The cutoff distance of the LJ potential and Coulomb potential was set at r c=12 Å. The positive charge of the hydronium ion and the cerium ion was balanced with the negative charge of the sulfonic acid group in the system. The number of water molecules was determined by the water content λ which represents the number of water molecules per sulfonate group. The water content was set at λ = 4, 8, 12 and 16, which correspond to the water content in the bulk Nafion® membrane at approximately RH = 30%, 70%, 85% and 100%, respectively. Cerium ion concentration is equivalent to 1.2 wt% cerium ion relative to the polymer mass in the membrane, which is the typical concentration for Ce-Nafion composite membranes. Three-dimensional periodic boundary condition was applied to simulation box. The annealing process was applied in order to equilibrate the system. After annealing, the production run was performed for 15 ns using the NVT ensemble. The temperature in the system was 350 K. The time step was set at 1 fs, and the sampling interval was 10,000 steps. The simulation system is shown in Fig. 1. To evaluate the convection ion transport, an external force was applied to water molecules along the uniaxial direction to obtain flux of water. The carbon atom at the base of the side chain in Nafion® was fixed to prevent the flowing of Nafion®. The convection ion transport coefficient, K convection, is defined as the ratio of the mass flux of cerium ions and the mass flux of water molecules. The mass flux can be obtained by multiplying the concentration of X and the mean velocity of X. In addition, structural analysis will be performed to identify the factors that are dominant in cerium ion transportation. Figure 1

Velocity jump processes: An alternative to multi-timestep methods for faster and accurate molecular dynamics simulations

We propose a new route to accelerate molecular dynamics through the use of velocity jump processes allowing for an adaptive time step specific to each atom-atom pair (two-body) interactions. We start by introducing the formalism of the new velocity jump molecular dynamics, ergodic with respect to the canonical measure. We then introduce the new BOUNCE integrator that allows for long-range forces to be evaluated at random and optimal time steps, leading to strong savings in direct space. The accuracy and computational performances of a first BOUNCE implementation dedicated to classical (non-polarizable) force fields are tested in the cases of pure direct-space droplet-like simulations and of periodic boundary conditions (PBC) simulations using Smooth Particle Mesh Ewald method. An analysis of the capability of BOUNCE to reproduce several condensed-phase properties is provided. Since electrostatics and van der Waals two-body contributions are evaluated much less often than with standard integrators using a 1 fs time step, up to a 400% direct-space acceleration is observed. Applying the reversible reference system propagator algorithms [RESPA(1)] to reciprocal-space (many-body) interactions allows BOUNCE-RESPA(1) to maintain large speedups in PBC while maintaining precision. Overall, we show that replacing the BAOAB standard Langevin integrator by the BOUNCE adaptive framework preserves a similar accuracy and leads to significant computational savings.

Exact long-range Coulombic energy calculation for net charged systems neutralized by uniformly distributed background charge using fast multipole method and its application to efficient free energy calculation.

In molecular dynamics (MD) calculations of the free energies of ions and ionic molecules, we often encounter net charged molecular systems where the electrical neutrality condition is broken. This charge causes a problem in the evaluation of long-range Coulombic interactions under periodic boundary conditions. A standard remedy for this problem is to consider a hypothetical homogeneous background charge density to neutralize the total system. Here, we present a new expression for the evaluation of Coulombic interactions for such systems including background charge using the fast multipole method (FMM). Furthermore, an efficient scheme is developed to evaluate solute-solvent interaction energies using the FMM, reducing the computational burden for the far-field part. We calculate the hydration free energies of Mg2+, Na+, and Cl- ions dissolved in a neutral solvent using the new expression. The calculated free energies show good agreement with the results obtained using the well-established particle mesh Ewald method. This demonstrates the validity of the proposed expression. This work should make a contribution to highly parallelized MD calculations for large-scale charged systems (particularly, those with over million particles).

A New Lipid Force Field (FUJI).

To explore inhomogeneous and anisotropic systems such as lipid bilayers, the Lennard-Jones particle mesh Ewald (LJ-PME) method has been applied without a conventional isotropic dispersion correction. As the popular AMBER and CHARMM lipid force fields were developed using a cutoff scheme, their lipid bilayers unacceptably shrink when using the LJ-PME method. In this study, a new all-atom lipid force field (FUJI) was developed on the basis of the AMBER force-field scheme including the Lipid14 van der Waals parameters. Point charges were calculated using the restrained electrostatic potentials of many lipid conformers. Further, torsion energy profiles were calculated using high-level ab initio molecular orbitals (LCCSD(T)/Aug-cc-pVTZ//LMP2/Aug-cc-pVTZ), following which the molecular mechanical dihedral parameters were derived through a fast Fourier transform. By incorporation of these parameters into a new lipid force field without fitting experimental data, the desired lipid characteristics such as the area per lipid and lateral diffusion coefficients were obtained through GROMACS molecular dynamics simulations using the LJ-PME method and virtual hydrogen sites. The calculated area per lipid and lateral diffusion coefficients showed satisfactory agreement with experimental data. Furthermore, the electron-density profiles along the membrane normal were calculated for pure lipid bilayers, and the resulting membrane thicknesses agreed well with the experimental values. As the new lipid force field is compatible with FUJI for protein and small molecules, the new FUJI force field will offer accurate modeling for complex systems consisting of various membrane proteins and lipids.

Modifying the CHARMM36 Lipid Force Field for LJ-PME Simulations

Long-range Lennard-Jones (LJ) interactions were historically excluded from energy calculations in Molecular Dynamics (MD) simulations due to high computational cost. However, these interactions have a significant impact on the properties of nonpolar systems. The LJ particle-mesh Ewald (LJ-PME) method offers a solution to this problem and retains good efficiency. Here we introduce a modification to the additive all-atom CHARMM lipid force field (FF) which makes it compatible with the LJ-PME method. The last generation of the CHARMM36 (C36) lipid FF is suitable for NPT simulations of lipid bilayer systems using the PME method for the nonbonded interactions but systematically yields lower values of surface area per lipid compared to experimental estimates when the long-range LJ interactions are included. Here we develop an automated workflow to optimize the parameters of the nonbonded interactions (partial atomic charges and Lennard-Jones parameters). 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) is used as the prototype and various structural properties (surface area per lipid, NMR deuterium order parameters, radial distribution functions of water for selected atoms) are used as targets. In each optimization step, free energy perturbation theory is used to estimate the properties from trajectories generated by the optimal parameters in the previous step. QM calculations of partial atomic charges are also used to guide the optimization. The parametrization improves the above-mentioned properties significantly with small changes to the C36 lipid FF and allows MD simulations of lipid systems to use the LJ-PME method.

Fast multipole method for three-dimensional systems with periodic boundary condition in two directions.

We derived a new expression for the electrostatic interaction of three-dimensional charge-neutral systems with two-dimensional periodic boundary conditions (slab geometry) using a fast multipole method (FMM). Contributions from all the image cells are expressed as a sum of real and reciprocal space terms, and a self-interaction term. The reciprocal space contribution consists of two parts: zero and nonzero terms of the absolute value of the reciprocal lattice vector. To test the new expressions, electrostatic interactions were calculated for a randomly placed charge distribution in a cubic box and liquid water produced by molecular dynamics calculation. The accuracy could be controlled by the degree of expansion of the FMM. In the present expression, the computational complexity of the electrostatic interaction of N-particle systems is order N, which is superior to that of the conventional two-dimensional periodic Ewald method for a slab geometry and the particle mesh Ewald method with a large empty space at an interface of the unit cell. © 2020 Wiley Periodicals, Inc.