Molecular Simulation with Discrete Fast Fourier Transform

Molecular simulation is widely used as a computer experiment to study molecular systems based on the first principle and accumulated knowledge. It extends our reach to the atomic level and to extreme conditions by providing information that is not available from experiment. As a complement to experimental studies, molecular simulation has played important roles in understanding the physics and chemistry of molecular systems. The rapid increase in computing power has facilitated extensive applications of molecular simulation in the study of complicated systems, such as ion channels, and difficult problems, such as protein folding.


Introduction
Molecular simulation is widely used as a computer experiment to study molecular systems based on the first principle and accumulated knowledge.It extends our reach to the atomic level and to extreme conditions by providing information that is not available from experiment.As a complement to experimental studies, molecular simulation has played important roles in understanding the physics and chemistry of molecular systems.The rapid increase in computing power has facilitated extensive applications of molecular simulation in the study of complicated systems, such as ion channels, and difficult problems, such as protein folding.
Molecular simulations are often expensive, because the sizes of simulation systems and the time scales of simulations are many orders of magnitude smaller than studied in experiments.For example, many biological molecular assemblies have millions of atoms and take milliseconds or longer to function.These are beyond the reach of current simulations without size reduction, or the use of additional assumptions to simplify and reduce the scope of the problem.To be relevant to real experiments, simulators tend to maximize the sizes of their simulated systems and/or the time scales of their simulations to the limit they can afford with their available computing resources.Therefore, improving calculation efficiency is always a focus of the development in molecular simulation methods.
Molecular interactions are the basis of all macroscopic properties.Accurate and efficient calculation of molecular interactions is the key for a successful simulation study.Long-range interactions, such as the electrostatic interaction and van der Waals interaction, play very important roles in the properties of molecular systems.However, they reach far beyond the size of a typical simulation system and are the most expensive part in molecular simulations.
Improving the calculation efficiency for long-range interactions has long been the goal of method development.The Ewald sum (Ewald 1921) is a well known method for calculating electrostatic interactions without the need to deal with a vacuum boundary interface by approximating large systems as small systems with periodicity.Recently, a method called the isotropic periodic sum (IPS) was developed as a general approach to the calculation of long-range interactions of all types of potentials (Wu and Brooks 2005;Takahashi, Yasuoka et al. 2007;Wu and Brooks 2008;Wu and Brooks 2009;Takahashi, Narumi et al. 2010;Takahashi, Narumi et al. 2011;Takahashi, Narumi et al. 2011) and has been applied in many Here, k is the Boltzmann constant and T is temperature.A transition from conformation m Ω to conformation n Ω has a probability of The commonly used Metropolis Monte Carlo method (Metropolis and Ulam 1949) propagates conformations according to eq. ( 2) to sample the canonical ensemble distribution.As can be seen from eq. ( 2), Monte Carlo simulations need calculate only energies.Further details of the Metropolis Monte Carlo method can be found elsewhere (Allen and Tildesley 1987).

Molecular dynamics
Molecular dynamics generates trajectories of particles by numerically solving the equation of motion (Allen and Tildesley 1987;Brooks, Brooks et al. 2009).The motion of simulation systems produces trajectories in the conformational space so that the ensemble properties as well as dynamics properties can be studied.Particles of a molecular system move according to the Newtonian equation of motion: Here, i f is the net force acting on particle i.The vectors, i p and i p $ , are the momentum and its time derivative of particle i.As can be seen from eq. (3), molecular dynamics simulations rely on forces to generate trajectories.Molecular dynamics simulations can be performed in many types of ensembles (Allen and Tildesley 1987).

Other sampling methods
Many other methods have been developed to sample the conformational space, such as the genetic algorithm (Dandekar and Argos 1992;Le Grand and Merz 1994;Ogata, Akiyama et al. 1995;Beckers, Buydens et al. 1997;Jones, Willet et al. 1997).In addition, there are many improved methods based on above mentioned methods or their combinations.For example, the self-guided molecular dynamics (SGMD) and self-guided Langevin dynamics (SGLD) were developed based on molecular dynamics or Langevin dynamics.they accelerate conformational search through enhancing low frequency motions of simulation systems (Wu and Wang 1998;Wu and Wang 1999;Wu and Brooks 2003;Wu and Brooks 2011a;Wu and Brooks 2011b) and has been applied in many simulation studies such as protein folding (Wu and Sung 1999;Wu and Wang 2000;Wu and Wang 2001;Wu, Wang et al. 2002;Wen, Hsieh et al. 2004;Wen and Luo 2004;Wu and Brooks 2004;Lee and Chang 2010;Lee and Olson 2010), ligand binding (Varady, Wu et al. 2002), docking (Chandrasekaran, Lee et al. 2009), conformational transition (Damjanovic, Miller et al. 2008;Damjanovic, Wu et al. 2008;Damjanovic, Garcia-Moreno E et al. 2009;Pendse, Brooks et al. 2010), and surface absorption (Abe and Jitsukawa 2009;Sheng, Wang et al. 2010;Sheng, Wang et al. 2010;Tsuru, Yosuke et al. 2010).

Molecular interactions
Molecular interactions are essential for all microscopic and macroscopic properties.The energy and forces of molecular interactions are the basis for conformational searching and sampling.In molecular simulations, molecular interactions are represented by force fields, which describe quantitatively relations of the interaction energy with system conformations.For computational feasibility, molecular interactions are often described as pairwise terms and total interactions are the sum of all pair interactions.For example, the CHARMM force field takes the following form (Brooks, Brooks et al. 2009): The potential energy is a sum of individual terms representing the internal and nonbonded contributions.The internal terms include contributions from bond (b), valence angle (θ), Urey-Bradley (UB, S), dihedral angle (ϕ), improper dihedral angle (ω), and backbone torsional correction (CMAP, φ, ψ).The parameters b k , θ k , UB k , ϕ k , and ω k are the respective force constants and the variables with the subscript "0" are the respective equilibrium values.The nonbonded terms include Coulombic interactions between the point charges ( i q and j q ) and the Lennard-Jones (LJ) 6-12 term, which is used for the treatment of the core- core repulsion and the attractive van der Waals dispersion interaction.
To further improve the accuracy, additional terms can be added.Polarizable force fields have conformational dependent atomic properties such as charges or dipole moments.In large molecule simulations, implicit solvation terms can be used to replace solvent molecules.There are also many attempts to directly calculate molecular interactions using quantum mechamics.

Short-range interactions
Interactions with limited ranges are categorized into short-range interactions.All internal terms shown in eq. ( 4) are limited within a molecule and are typical short-range interactions.Some non-bonded interactions, such as the repulsion between atoms, are also considered to be in this category.Short-range interactions occur within a short distance range and the number of interactions is very small as compared to long-range interactions.Therefore, the computing cost for short-range interactions is relatively low and no special treatment is usually needed.

Long-range interactions
Long-range interactions are interactions that can reach far beyond molecular sizes.Typically, long-range interactions have infinite interaction ranges.In molecular simulation, electrostatic and van der Waals interactions are two major long-range interactions.Because long-range interactions play a crucial role in system properties and are expensive to calculate, their accurate and efficient calculation is the focus of many method developments.

Methods to calculate long-range interactions
Many methods have been developed for the calculation of long-range interactions.In the early stage of molecular simulations, due to the limit in computing resources, molecular interactions are limited to a certain range by methods like the minimum image method or the cutoff based methods.When accurate long-range interactions are desired, typically for electrostatic interactions, more sophisticated methods are developed.Below we briefly describe four typical methods for the calculation of long-range interactions.Fig. 1 illustrates the concepts of these methods.

Cutoff methods
Cutoff based methods consider only interactions within a certain distance range.Fig. 1 (a) illustrates this concept.The cutoff based methods assume that a particle interacts with only the particles or their images within a certain distance (Steinbach and Brooks 1994).This distance is called the cutoff distance.To avoid discontinuities in energy and/or forces in the cutoff region, interaction potentials are shifted or switched as below: Here, r c is the cutoff distance, and ) , ( c cut r r S represents a shift or switch function.Under periodic boundary conditions (PBC), the minimum image convention must be used when applying the cutoff based methods.A pair of particles can have an infinite number of interaction pairs between the two particles or their images.The minimum image convention considers only the pair with the minimum distances for each pair of particles.The cutoff based methods save computing cost by neglecting interacting pairs with distances larger than r c .However, this approximation often causes large errors in simulation results, especially for electrostatic interactions.

Reaction field method
The reaction field method describes a simulation system with a local cavity and a long-range reaction field, as illustrated in Fig. 1 (b).For each atom, there is a local spherical cavity centred on it.The central particle interacts directly with all particles within the cavity.Interactions with particles beyond the cavity are replaced by a reaction field.The reaction field acting on the particle is proportional to the moment of the cavity surrounding it: Here, s ε is the dielectric constant of the surrounding environment and j μ is the dipole moment of molecule j.This method acts like a cutoff method but retains a certain physical basis.Its use requires an artificial switching function to overcome the discontinuity crossing the cavity boundary.

The Ewald sum
The Ewald sum (Ewald 1921) is an efficient method for the calculation of interactions with the lattice images created by the periodic boundary conditions (PBC), as shown in Fig. 1 (c).
For electrostatic interactions, the total energy of a system of N particles is: Here, q i i s t h e c h a r g e o f p a r t i c l e i.The first summation runs over all cell vectors, . Here, www.intechopen.comMolecular Simulation with Discrete Fast Fourier Transform 143 Here, V is the volume of the PBC box.The parameter β controls the relative contributions from the real space, the first term, and the reciprocal space, the second term.The complementary error function, erfc (x), takes the following form: The Ewald sum is a method of order N 2 and is expensive for large systems.

The Isotropic Periodic Sum
The IPS method was developed to overcome some artefacts of the lattice sum methods.The lattice sum methods, like the Ewald sum, assume that the remote regions can be represented by images created by periodic boundary conditions (PBC).The PBC images distribute discretely in lattice points throughout the space and are anisotropic in nature, as shown in Fig. 1 (c).When the PBC size is comparable to macromolecules, which is often the case due to the consideration of computing cost, macromolecules have detectable orientation bias due to their interaction with their images (Wu and Brooks 2005).The calculation of lattice sums is expensive and special techniques like the Ewald method and the particle meshed Ewald method are needed to save computing time.As shown in Fig. 1 (d), the IPS long-range interactions are isotropic so that the orientation bias caused by image interactions can be overcome.The IPS long-range interactions are represented by equivalent short-range functions, which can be calculated as efficiently as the cut-off based methods and the calculation is especially efficient for parallel computing.
The concept of the IPS method is using isotropic distributed images of a local region to represent remote environment to calculate long-range interactions.The difference between the IPS method and the lattice sum methods such as the Ewald sum lies in the shape and distribution of remote images.The images for the lattice sum are generated from the periodic boundary conditions and are discretely and anisotropically positioned at the lattice points in space.The images for the IPS calculation are imaginary, which means they do not explicitly exist in a simulation system, and are distributed in an isotropic and periodic way around each particle.The IPS images are distributed equally in all homogeneous dimensions.Analytic solutions for IPS are available for many potential types and the calculation of the IPS potentials is straightforward and very efficient.
Fig. 2 illustrates the definition of the local region and its isotropic periodic images in a square periodic boundary system.The local region of particle 1 is enclosed by a dashed circle of radius R c .The isotropic periodic images of the local region and their particles are shown as dotted circles and dotted particles labeled correspondingly.There are an infinite number of image shells around the local region.The image regions of the first layer are bounded with the local region and occupy the area with a radius from R c to 3R c .In this layer, the isotropic periodic images of particle 1 are distributed on image shells with a radius of 2R c .The image regions on an image shell are statistical representation of conformations around this image shell and can overlap with each other.Because the image regions are translation of the local region, the images of each particle will distribute on its own image shells centered at this particle.As shown in Fig. 2, the isotropic periodic images of particle 2 distribute on image shells centered at particle 2. Particle 1 only interacts with particles within its local region, i.e., particles 2, 3, and 4, and their isotropic periodic images, including itself.Similarly, all particles in the local region will interact with the isotropic periodic images of particle 1.All other particles, such as, particles 5, 6, 7, and all images generated by the periodic boundary condition that are outside the local region are not seen by particle 1 and are represented by the isotropic periodic images of the local particles in the calculation of long-range energies.Particle 4 is at the boundary of the local region of particle 1 and has the same distance, R c , to particle 1 and its nearest isotropic periodic image on the first image shell.Due to the periodicity, the total force on particle 4 from particle 1 and its images is zero.Please note that the total interaction between particle 1 and all images of particle 2 will be the same as that between particle 2 and all images of particle 1.If we assume that the structure beyond the local region can be represented by the isotropic periodic images of the local region, the summation over all particles beyond the local region becomes a function of the local region structure: The IPS potentials have analytic forms for many commonly used potentials (Wu and Brooks 2005).For electrostatic interaction, where 1 q and 2 q are the charges of the two interacting particles, its IPS image interaction is where γ is the Euler's constant, , and The IPS analytic solutions are often very complicated and are time consuming to compute directly.Instead, we use numerical functions that fit these analytic solutions for efficient calculations in simulations.To avoid numerical mistakes in implementation, we use the following simplified polynomial with rational coefficients: Eq. ( 14) has an average deviation of c j i R q q 4 10 7.7 − × from eq. ( 13).We call eq. ( 13) or eq.( 14) the non-polar IPS electrostatic potential.
For polar systems, the opposite charges play a screening effect in charge-charge interactions.After considering the opposite charge screening effect, we obtain the following expression for the polar IPS electrostatic potential: The Lernnard-Jones potential can be separated into repulsion and dispersion terms: where, 0 ε and * r are the energy minimum and the minimum distance.are the constants for repulsion and dispersion interactions, respectively.We obtained the following polynomials with rational coefficients by fitting into the analytic solutions of the dispersion and repulsion IPS potentials (Wu and Brooks 2005;Wu and Brooks 2009): The average deviations from analytic solutions are With eqs.( 14), ( 15), ( 17), and ( 18), the IPS method calculates long-range interactions as the short-range functions with c R r < .Obviously, the computating cost is comparable to the cutoff method with c c R r = .Therefore, the IPS method is of order N, and can scale very well for parallel computing.
The IPS method described above assumes that a system is fully homogeneous in all dimensions and is called the three-dimensional (3D) IPS method.For partial homogeneous systems like membrane or fibril systems, two-dimensional (2D) and one dimensional (1D) IPS methods have been developed.Please refer to the original paper (Wu and Brooks 2005) for details about the 2D-IPS and 1D-IPS methods.

DFFT for long-range interactions
For pairwise molecular interactions, the interaction energy can be expressed as a convolution of a momentum distribution with potential functions.For example, electrostatic interaction energy can be written as the charge distribution convolved with the electrostatic potential: Here, Q is the charge distribution.
) (r ψ is the potential function at position r from a unit charge at origin and its all images.self ψ is the self-potential to correct the self-interaction included in the convolution.Fourier transform provides a convenient way to deal with such convolution operations: After separating direct interactions, the electrostatic interaction energy can be written as: Eq. ( 21) separates potentials to a short-range part calculated by a direct pairwise summation and a remaining long-range part calculated as a convolution.Other long-range interactions such as the Lennard-Jones energy can also be treated like eq. ( 21).

Particle Mesh Ewald (PME)
The PME method assigns the charge density to a finely spaced mesh in a simulation box so that the reciprocal part of the Ewald sum can be calculated as convolutions that can be accelerated with DFFT.Here we describe a smooth particle mesh Ewald method based on the b-spline interpolation (Essmann, Perera et al. 1995).If represents the b-spline function of the n-th order, we have: Here u is the coordinate to be interpolated.In a simulation box, atomic charges can be distributed over a set of predefined grid points, (k 1 , k 2 , k 3 ), where k 1 =1,2,…, K 1 , k 2 =1,2,.., K 2 , k 3 =1,2,…,K 3 .Assume the charges of a simulation system of N particles are { } N q q q ,..., , 2 1 = q .After spreading q on the grid points we have a distribution, Q, (Essmann, Perera et al. 1995): We use to represent the inverse discrete Fourier transform: The structure factor can be expressed as: www.intechopen.com Fourier Transform -Materials Analysis 148 Where ( ) The electrostatic long-range sum can be calculated in the following way: Here, and The force can be calculated by The PME algorithm is of order Nln (N).As compared with the Ewald sum of order N 2 , PME is more suitable for large systems.

IPS/DFFT method
The underlying assumption of the IPS method is that the simulation system is homogenous and isotropic over the size of the local region (defined by the local region radius, R c ), and for convenience the local region in a homogeneous system is defined to be the same as the region within a cutoff distance, r c ., i.e., R c= r c .However, in many cases a simulation system is not homogenous in such a length scale (r c ~ 10 Å).To accurately describe the long-range interaction of heterogeneous systems, IPS need use a local region large enough to cover the heterogeneous range, up to the size of the simulation system or the periodic boundary box.Obviously, it is highly time consuming to do direct pair-wise calculation with such a large cutoff distance.To efficiently calculate interactions within such a large local region, the IPS/DFFT method split long-range interactions into two parts, a cutoff part and a longrange part.The cutoff part is calculated by summing atom pairs within a cutoff range (about 10 Å), and the long-range part is treated as convolutions and is calculated efficiently using DFFT.
A simple case of a heterogeneous system is a two-ion periodic system as shown in Fig. 3.With a small cutoff distance of r c , one ion often falls out of the sight of the other, therefore, www.intechopen.comsuch a region is hard to be representative of the system.In other words, this system is hetereougeneous in a size scale of r c .If we look at the system with the images created by the PBC, we can see that the PBC forces the system to be repeated throughout the space.With a local region radius, R c , larger than the PBC size, a local region has a reasonable number of particles and can be a good representative of the system.That is, the system is homogeneous in a size scale of R c .With such a large local region radius, IPS energies can well approximate the long range interactions.Therefore, to use IPS to describe this system, it is desirable to use a local region (defined by R c ) larger than the cutoff regions (defined by r c ).When R c is larger than the PBC size, the lattice symmetry from the PBC will be imposed into long-range interactions.Obviously, the IPS energies will approach that of the lattice sum when R c becomes infinitely large., and the long-range (LR) part, The cutoff part goes to zero at the cutoff distance, r c .
The long-range part is the difference between the IPS energy and the cutoff part, eq. ( 34).Because IPS potentials have a non-zero boundary energy, to avoid energy discontinuities when particles move across the boundary, boundary energies are subtracted from the IPS potentials and calculated separately.Therefore, the actual long-range part used in the calculation is: And the boundary energy is defined as: The boundary energies are summed separately in a simulation: Here N is the number of particles and M is the number of PBC image particles.V is the volume of the simulation system.The total IPS energy is a sum over all particle pairs within R c : The cutoff sum, C E , can be calculated in a pair wise way like the cutoff methods, and the long range sum, LR E Δ , can be calculated using the DFFT technique.
For electrostatic potential, eq. ( 9), we use the following smoothing function: Alternatively, the polar IPS potential, eq. ( 15), can be used as the smooththing function: The advantage to use eq.( 41) is that the long-range part, eq. ( 36), becomes a measure of the electrostatic hetereogeneity of a simulation system in a length scale of r c .If a simulation system is homogeneous in the length scale of r c , For the Lernnard-Jones potential, eq. ( 16), the repulsion part is short-ranged and its IPS is calculated only within the cutoff distance, c r , using eq. ( 18 with c c r R = .For the dispersion part, the smoothing function is: For the long-range sum, because the dispersion term is atom type dependent, to simplify the calculation, we transform the atom pair dependent dispersion parameters, C ij , to a transferable quantity, d i , which is defined as the dispersion momentum to measure the contribution of each atom to long range dispersion interactions.

∑∑
In a simulation box, like atomic charges, atomic dispersion momentums are distributed over a set of predefined grid points, (k . After spreading q and d on the grid points we have grid distributions Q and D. Like the charge spreading in eq. ( 24), we again use Cardinal b-spline to do the d spreading (Essmann, Perera et al. 1995).
The electrostatic long-range sum can be calculated in the following way: ( ) The energy type subscript is dropped here to indicate that eq. ( 45) applies to all energy types.The Fourier transform of q and d can be approximated by Where is defined by eq. ( 28).We have is calculated with eq. ( 30).
The long-range sum of electrostatic interaction, eq. ( 32), can be rewritten to: Similarly, the dispersion long-range sum is: Pressure tensors are important quantities in molecular simulation.The contributions from these long-range sums are calculated directly as a summation in the Fourier space: Here, α, β stand for either x, y, or z, and The forces acting on each particle can be derived from the long-range energies, eqs.( 51) and ( 52).
In summary, the IPS/DFFT method uses eq. ( 34) to calculate the cutoff part directly for all atom pairs within the cutoff distance, c r , and uses eq. ( 51) and ( 52) to calculate the long- range part through DFFT.The total interaction is a sum of these two parts as shown in eq.(33).

IPS/DFFT for finite systems
Finite systems like proteins in vacuum are often the objects of simulation studies.They are hetereogeneous in nature.With some modifications, the IPS/DFFT method can be extended to finite systems.
The first modification is to create a virtual periodic boundary for DFFT.The size of the boundary box is twice the maximum dimensions of the actual system plus b-spline widths.
Here n is the order of the b-spline and x Δ , y Δ , and z Δ are the grid sizes in x, y, and z directions, respectively.The charges and dispersion momentums are spread over the grid points in the PBC box according to eqs. ( 24) and ( 45).The second modification is that the local region radius will be infinity, . The 3D IPS interaction with an infinity local region radius becomes purely the original pair potential: . To avoid interaction with images created by the virtual PBC, the long-range sum must be limited to a half of the box size in each dimension: Here, x ij , y ij , and z ij are the differences of coordinates in x, y, and z directions, respectively.Fig. 4 illustrates the virtual PBC box defined for a protein in vacuum.With these modifications, we can get rid of the artificial image interactions due to the virtual periodic boundary condition while taking advantage of DFFT, and a finite system can be treated in the same way as a periodic boundary system described above.
x max L x =2x max +nΔx y max L y =2y max +nΔy Fig. 4. The virtual periodic boundary for a finite system as defined by eq. ( 54).

Applications
Due to efficient calculation of convolutions through DFFT, PME and IPS/DFFT can efficiently and accurately calculate long-range electrostatic interactions for large systems.In addition, IPS/DFFT can be applied to all types of potentials, including the Lennard-Jones potential.Here we present several examples to demonstrate the application of these methods.

Water interfaces
A water interface system is created by enlarging a cubic PBC box length along the z-axis, so that a gas phase is produced above and below the water liquid.There are 2180 TIP3P water molecules in the 40×40×80 Å 3 orthorhombic periodic boundary box.Fig. 5 shows a snapshot of this system.This is a typical heterogeneous system involving phase equilibrium.In this system, 2 ns MD simulations at constant temperature (300K), constant volume (40×40×80 Å 3 ) are performed with PME, 3D IPS, 2D IPS, and the IPS/DFFT method.Here, 2D IPS is a method designed specifically for two dimensional partial homogenous systems (Wu and Brooks 2005;Klauda, Wu et al. 2007).
An important property of interface systems is the surface tension.Because the surface tension is very sensitive to long range interactions, an accurate calculation of surface tension is often time consuming.The surface tension is evaluated from, where L z is the size of the simulation box normal to the interface, P zz is the normal component of the internal pressure tensor and P xx and P yy are the tangential components.Since the MD simulations here contain two interfaces (see Fig. 5), the prefactor 0.5 is required to obtain γ on a per interface basis.Pressure tensors are sensitive to long-range structures.In the case that the box size is small as compared to the homogeneity scale, i.e., the thickness of the water layer in this case, it is recommended to set the R c to twice of the longest box side or larger to equally consider all images in every direction.Fig. 6 shows the surface tension results from different methods.As can be seen, the results from PME and 3D IPS strongly depend on the cutoff distance, while 2D IPS and the IPS/DFFT method produce results showing little dependence on the cutoff distance.Obviously, as the cutoff distance increases, the results from PME and 3D IPS approach that from 2D IPS and the IPS/DFFT methods.We can see that to calculate surface tension, the IPS/DFFT method with a normal cutoff distance, r c =10 Å, is more efficient in the interface system simulation.
Another property that is sensitive to the long-range interaction is the electrostatic potential profile across the layer, which is calculated by a double integration of the Poisson's equation,  These results indicate that both the IPS/DFFT method and the 2D IPS method can accurately describe the long-range interactions of this interface system.The IPS/DFFT method is much faster than the 2D IPS because its short-range cutoff contains much fewer atom pairs than the cylinder cutoff in the 2D IPS method (Wu and Brooks 2005).

A sodium aqueous system
Solvation of an ion involves interactions far beyond the simulation box.The small number of ions in a simulation box makes the system highly heterogeneous.The solvation energy of an ionic solution can reflect the enthalpy effect as well as the entropy effect of charged ions onto the solvent, which provides a good case to examine the IPS/DFFT method.
We performed MD simulations of a sodium aqueous solution in both charged and neutral states to examine the energy difference, defined as the electrostatic solvation energy: ( ) Here aq E and vac E represent the average potential energies of the system with and without water in the same periodic boundary box, and 1 = q and 0 = q represent the charged and neutral states of the solvated ion, respectively.The system contains one sodium ion and 265 TIP3P water molecules.Simulations of 20 ns are performed at 300 K and in a 20×20×20 Å 3 cubic periodic boundary box for both the charged and neutral states with PME, 3D IPS, and the IPS/DFFT method.A cutoff distance of 10 Å is used for all simulations.Clearly, the results are almost independent of grid sizes.The cpu times of the IPS/DFFT method are comparable to those of the PME method.

Proteins in vacuum
Systems without periodic boundary conditions are clearly heterogeneous by nature.By imposing a virtual periodic boundary and avoiding interactions with images, the IPS/DFFT method can be applied to such finite systems.
We chose the x-ray structure of acetylcholine binding protein (ACHBP) (Brejc, van Dijk et al. 2001) (PDB code: i9b) to examine the energy calculation for non-periodic systems.We use its monomer and pentamer as examples of systems of small and large sizes.Fig. 8 shows the image of this protein in its monomer and pentamer form.than the cutoff method at small cutoff distances, they both show significant force deviations with the cutoff distances up to 50 Å.By contrast, the IPS/DFFT result is better than these two methods by an order of magnitude.
For the monomer, the IPS/DFFT costs more cpu time than the cutoff method at a given cutoff distance.However, at 10 Å, the IPS/DFFT method has a force deviation, f δ , of 0.019 www.intechopen.com Fourier Transform -Materials Analysis 158 kcal/mol·Å, while for the cutoff (force switch) method, = f δ 1.92 kcal/mol·Å.Even with a cutoff distance of 50 Å, the cutoff method has = f δ 0.022 kcal/mol·Å.Therefore, the IPS/DFFT method can reach a better accuracy with a 10 Å cutoff than the force switch method with a cutoff distance of 50 Å.In other words, to reach the same accuracy, the IPS/DFFT method needs less cpu time.This is obvious for large systems.For the pentamer, at a 10 Å cutoff, the IPS/DFFT method can reach = f δ 0.022 kcal/mol·Å with 0.158 seconds of cpu time.The force switch method with a cutoff distance of 55 Å can only reach an accuracy of = f δ 0.289 kcal/mol·Å, but uses 2.03 seconds of cpu time.Clearly, for small systems, accurate forces can be calculated by summing directly over all atom pairs, while for large systems, the IPS/DFFT method is a superior way to efficiently get accurate forces.

Lipid bilayer system
Venable et al. applied IPS/DFFT and PME methods for simulations of lipid bilayers and monolayers (Venable, Chen et al. 2009).The method is demonstrated to be highly accurate for simple bulk fluids, liquid/liquid and liquid/vapor interfaces, and lipid bilayers and monolayers (Klauda, Wu et al. 2007).Values for r C (the cutoff distance for direct evaluation of pairs) and R C (the local region radius) equal to 10 Å and twice the longest edge of the periodic cell, respectively, provide excellent efficiency and accuracy.Dimyristoylphosphatidylcholine (DMPC) monolayers and bilayers are simulated with the CHARMM (Chemistry at HARvard Molecular Mechanics) C27r lipid parameter set using IPS/DFFT and PME.This picture shows that even highly heterogeneous systems such as monolayers can be accurately treated as homogenous when particles from many periodic replicates are included.With the IPS/DFFT method, interactions within a 10-12 Å cutoff distance r C are calculated directly, like the "real-space" part in PME.The remaining IPS interactions within R C are then evaluated by DFFT on a grid, and the direct interaction is subtracted to correct for overcounting.This is analogous to the splitting of real and k-space terms in PME that allows the Ewald equations to be solved in N ln N rather than N 2 time, where N is the number of particles.
Fig. 10.The long (RC, equal to twice the longest edge length) and short (r C =10 Å) cutoffs of the 3D-IPS/DFFT method illustrated for a DMPC monolayer.Coloring is as follows: water, blue; hydrocarbon chains, grey; carbonyl oxygens, red; phosphate groups, green; quaternary amines, purple.The primary cell in the center is darker than the images, and the vapor phase between the chains is white.Top panel shows the substantial number of image atoms within R C (leading to homogeneity of the region) while bottom panel highlights the highly anisotropic distribution of atoms within r C (and why the longer cutoff is necessary) (Venable, Chen et al. 2009) One important property for interfacial systems is the surface tension, which is very sensitive to long-range interactions, both electrostatic and Lennard-Jones interactions.Proper consideration of long-range interactions is essential to obtain accurate surface tension results.Table 2 compares the surface tensions for DMPC obtained with several methods.The most important and clear cut result is that γ for the monolayer for IPS/DFFT is www.intechopen.comsubstantially, and statistically significantly higher than for PME.The difference is 10 dyn/cm for both cutoffs, and the trend holds for other surface areas of the DMPC isotherm (Figure 11), where the average difference is 8 dyn/cm.This parallels the results obtained for octane/vapor interfaces reported in the study, where very long cutoffs on the LJ interactions are required.In other words, the inclusion of long-range LJ interactions raises the surface tension.Consequently, the very good agreement of experimental monolayer surface tensions and those obtained from C27r simulated with PME arise from cancellation of errors.Table 2. Surface tensions of bilayers (dyn/cm/leaflet) and monolayers (dyn/cm) of DMPC (T=310 K, area/per lipid = 59.6 Å 2 ) simulated with the force field C27r, and of DPPC (T=323 K, area/per lipid = 64 Å 2 ) with a modified version (Sonne, Jensen et al. 2007) of C27 (Venable, Chen et al. 2009).a (Lide 2000), b (Small 1986)

Conclusion
Long-range molecular interactions play a crucial role in molecular systems.Because longrange interactions require the most of the computing time in molecular simulations, improving their calculation efficiency is a major focus in molecular simulation method development.One approach to an efficient calculation of long-range interactions is treating long-range interactions as convolutions between related property distributions and potential functions so that the calculation can be efficiently handled with DFFT.Particle-mesh-Ewald (PME) is a widely used method that utilizes DFFT to perform Ewald summation for electrostatic interaction.IPS/DFFT is a recently developed method for the calculation of potentials of all kinds.For applications where additional long-range interactions, such as Lennard-Jones interactions, are crucial, IPS/DFFT is a better choice for calculating longrange interactions.

Fig. 1 .
Fig. 1.Calculation methods for long-range interactions.(a) the cutoff methods: interactions with particle within a cutoff distance; (b) the reaction field method: in addition to the cutoff interaction, the environment has an induced field on the cutoff region; (c) the lattice sum (the Ewald sum): interactions with all lattice images; (d) the isotropic periodic sum: interactions with all isotropic periodic images

Fig. 2 .
Fig. 2. The local region and the isotropic periodic images in a square periodic boundary system.The local region of particle 1 is enclosed by the dashed circle.The image shells of particles 1 and 2 are shown as dotted-dashed circles around particles 1 and 2, respectively.The image shells of other particles are not shown for clarity.The isotropic periodic images of the local region shown as dotted circles distribute around the local region and can overlap with each other.Particle 1 interacts with particles 2, 3, and 4 in its local region and the isotropic periodic image particles, shown as dotted particles.Particle 4 is at the boundary of the local region and has the same distance to particle 1 and the nearest image of particle 1

Molecular
Simulation with Discrete Fast Fourier Transform 145 Here, the summation, ∑ ≤ c ij R r, runs over all particles, including any PBC image particles, that are within the range of R c from particle i. long-range contribution as a function of ij r and c R .
IPS potential, which is the sum of the pair interactions within the local region and the image interactions:

Fig. 3 .
Fig.3.A two-ion periodic system with a square periodic boundary condition.A region within a small radius of r c is highly heterogeneous.However, a region within a large radius of R c would be much closer to a homogeneous system due to the PBC To work with a large R c , we define a smoothing function, ) , ( c r r ϑ to split the IPS potential into two smooth parts, the cutoff (C) part, ) , ( c C r r ε, and the long-range (LR) part, smoothing function can be the IPS potential, eq.(17), with a local region radius of r c : Lennard-Jones hetereogeneity of a simulation system in a length scale of r c .

Fourier
based on eqs.(24) and (43) from the property of the b-spline functions

Fig
Fig. 5.A water interface system with 2108 TIP3P water molecules in a 40×40×80 Å3 orthorhombic boundary box density and electrostatic potential along the z direction, respectively.

Fig. 7
Fig. 7 shows the electrostatic potential profiles calculated from simulations using the PME, 3 D I P S , 2 D I P S , a n d t h e I P S / D F F T m e t h o d s .C l e a r l y , 3 D I P S c a n n o t p r o d u c e c o r r e c t electrostatic potential profiles with small cutoff distances.However, as the cutoff distance increases, the result from 3D IPS becomes closer to the PME result.The 2D IPS, the IPS/DFFT, as well as the PME method, produce almost identical results.

Fig. 7 .
Fig. 7.The electrostatic potential profile cross the water interface system from simulations with different methods , and K 3 are the grid numbers along the three sides of the 20×20×20 Å 3 cubic box

Fig. 8 .
Fig. 8. Acetylcholine binding protein (ACHBP) (PDB code: i9b) in its monomer and pentamer forms.The backbones are shown as ribbons.For clarity, atoms are not shown in the top views Figs. 9 (a) and (b) show the force root-mean-square deviations and cpu times at different cutoff distances for the monomer and pentamer.Even though 3D IPS shows better results than the cutoff method at small cutoff distances, they both show significant force deviations with the cutoff distances up to 50 Å.By contrast, the IPS/DFFT result is better than these two methods by an order of magnitude.

Fig. 9 .
Fig.9.The root-mean-square deviations (rmsd) of the forces and cpu times from the cutoff method, 3D IPS, and the IPS/DFFT method with different cutoff distances.The rmsd is calculated against the forces calculated with no cutoff.(a) the ACHBP monomer, i9b; (b) the ACHBP pentamer, (i9b)

Fig. 10
Fig.10illustrates a lipid monolayer and periodic images with R C equal to twice the longest edge length.This picture shows that even highly heterogeneous systems such as monolayers can be accurately treated as homogenous when particles from many periodic replicates are included.With the IPS/DFFT method, interactions within a 10-12 Å cutoff distance r C are calculated directly, like the "real-space" part in PME.The remaining IPS interactions within R C are then evaluated by DFFT on a grid, and the direct interaction is subtracted to correct for overcounting.This is analogous to the splitting of real and k-space terms in PME that allows the Ewald equations to be solved in N ln N rather than N 2 time, where N is the number of particles.

Table 1 .
Electrostatic solvation energies of a sodium ion calculated with different methods.K 1 , K 2 Table1lists the electrostatic solvation energies from different simulations.As can be seen, 3D IPS underestimates the solvation energy by ca. 3 kcal/mol as compared to the PME result.This difference indicates that the homogenous approximation with a cutoff of 10 Å causes significant error for the solvation energy calculation.By contrast, the IPS/DFFT method produces results very close to the PME results, supporting the idea that with a cutoff larger than the homogeneity scale, 3D IPS can well approximate a heterogeneous system.Table 1 also lists the result of the IPS/DFFT method with different grid sizes.