Zhiyong Xu,Zhongjin He,Xuebo Quan,Delin Sun,Zhaohong Miao,Hai Yu,Shengjiang Yang,Zheng Chen,Jinxiang Zeng,Jian Zhou,*

1 School of Chemistry and Chemical Engineering,Guangdong Provincial Key Lab for Green Chemical Product Technology,South China University of Technology,Guangzhou 510640,PR China

2 National Center for International Research on Deep Earth Drilling and Resource Development,Faculty of Engineering,China University of Geosciences,Wuhan 430074,China

3 Biosciences and Biotechnology Division,Physical and Life Sciences Directorate,Lawrence Livermore National Laboratory,Livermore,CA 94550,USA

Keywords:Molecular simulation Complex fluid Charged system Soft matter Electrostatic interaction

ABSTRACT Molecular simulation plays an increasingly important role in studying the properties of complex fluid systems containing charges,such as ions,piezoelectric materials,ionic liquids,ionic surfactants,polyelectrolytes,zwitterionic materials,nucleic acids,proteins,biomembranes and etc.,where the electrostatic interactions are of special significance.Several methods have been available for treating the electrostatic interactions in explicit and implicit solvent models.Accurate and efficient treatment of such interactions has therefore always been one of the most challenging issues in classical molecular dynamics simulations due to their inhomogeneity and long-range characteristics.Currently,two major challenges remain in the application field of electrostatic interactions in molecular simulations;(i) improving the representation of electrostatic interactions while reducing the computational costs in molecular simulations;(ii)revealing the role of electrostatic interactions in regulating the specific properties of complex fluids.In this review,the calculation methods of electrostatic interactions,including basic principles,applicable conditions,advantages and disadvantages are summarized and compared.Subsequently,the specific role of electrostatic interactions in governing the properties and behaviors of different complex fluids is emphasized and explained.Finally,challenges and perspective on the computational study of charged systems are given.

1.Introduction

The full understanding of interactions and mechanisms at the molecular scale is of critical importance for insights into the thermodynamic and kinetic properties of charged complex fluid systems.Among the various molecular interactions,electrostatic interactions are extremely significant,which appear ubiquitously in complex charged systems.For example,in the application of ion separation and electrodeionization,electrostatic interactions can be used to drive ion transport by external electrical fields [1].Under the external pressure,the internal structure of the piezoelectric materials deforms,the non-centrosymmetric positive and negative ions deviate and polarize with generating a potential difference and vice versa,which realize the mutual conversion of mechanical energy and electrical energy[2,3].In another instance,ionic liquids,composed of a large number of organic cations and inorganic anions,have unique advantages in thermal stability and conductivity compared with traditional conventional solvents;these properties are mainly determined by the strength and directionality of Columbic interactions between the cations and anions[4].Ionic surfactants are composed of a positively or negatively charged head group and a long alkyl chain tail,in which electrostatic interactions also have profound effects on the microscopic aggregation behaviors of ionic surfactants at different phase interfaces[5–9].Zwitterionic polymer is a kind of novel polyelectrolyte materials with both cationic and anionic groups in the same monomer residue,simultaneously maintains the overall charge neutrality [10].Due to their strong hydration capability via electrostatic interactions,zwitterionic polymers are endowed with unique molecular structure and physicochemical properties [11,12].For proteins,the ionizable amino acids and polar groups of proteins play critical roles in regulating protein folding,structural stability as well as enzyme catalysis [13].The surface of nucleic acids always carry negative charges,thus the electrostatic effect caused by cations in the surrounding environment has a strong influence on their folding and packing into a specific conformation to correctly perform their biological functions [14–16].Phospholipids,which consist of hydrophilic headgroups and long hydrophobic alkyl chains,are major components of biomembranes,thus the dynamic process of some functional nanomaterials entering into the interior of cells may be strongly affected by electrostatic interactions[17].Electrostatics,therefore,is of central important in the studies of nature and behavior of charged complex fluid systems,such as ions,ionic liquids,ionic surfactants,zwitterionic materials,nucleic acids,proteins,cell membranes and so on.

However,due to the long-range and multi-body nature of electrostatic interactions,it is rather difficult to develop a robust characterization method to quantify the electrostatic interactions in charged complex fluid systems.Although many experimental studies have made significant progresses on this important topic,it is often uneasy to systematically investigate the detailed interactions owing to the resolution limitations in both time and space[18,19].With its high accuracy,low cost,flexible and adjustable parameters,computer simulation technology has gradually been developed from experimental aids to an independent new branch of science with predictive ability.However,even with the most state-of-art molecular simulation techniques,the treatment of electrostatic interactions is still the computational bottleneck,which occupies the vast majority of running time in the simulation process [20].In particular,the simulation processes are timeconsuming when computing the electrostatic interactions in biomolecular systems;however,the calculation of electrostatic interactions cannot be ignored since their representation has a profound impact on the accuracy of simulation results[13].Therefore,it is necessary to carefully construct accurate algorithms to avoid artifacts and to improve computational efficiency.

Over the past few decades,numerous researches have been performed to find an efficient way to deal with electrostatic interactions.The simplest way to deal with electrostatic interactions is the method that based on the truncation scheme,where electrostatic interactions are completely ignored when beyond the cutoff value[21].Reaction field method,as another computation method is superior to the simple truncation method.It explicitly calculates electrostatic interactions within the cutoff radius,while the electrostatic interactions beyond the cutoff value is described by a mean-field approximation [22].The first standard algorithm to accurately calculate electrostatic interactions is the Ewald summation,proposed by Ewald in 1921 [23].Based on the traditional Ewald summation method,many improved algorithms have been developed in the 1970s,such as the PPPM method based on gridding algorithm proposed by Hockney and Eastwood [24],the PME method proposed by Darden [25],the smooth particle mesh Ewald (SPME) method developed by Pedersen et al.[26]and the Fast Fourier-Poisson method introduced by York and Yang [27].With the enhancement of computer architecture,some new calculation methods also have been proposed [28–30].

In this review,the basic ideal about classical calculation methods for dealing with long-range electrostatic interactions in molecular simulations under explicit solvents is first presented.Then,the effect of electrostatic interactions on the thermodynamic and kinetic properties of charged complex fluid systems,which including ions,piezoelectric materials,ionic liquids,ionic surfactants,zwitterionic materials,nucleic acids,proteins and cell membranes are summarized.Finally,the future directions of the electrostatic interaction representation and their applications in the molecular simulation field of charged complex fluid systems are discussed.

2.Methods for Treating Electrostatic Interactions

Accurate and effective treatment of electrostatic interactions have been a difficult task in molecular simulations due to the inhomogeneity,complicated geometry,multiple scales and the manybody nature of underlying charged complex fluid systems.It accounts for the vast majority of running time in molecular simulations.Considering that most charged systems are embedded in an aqueous environment,the surrounding solvent has a very important influence on the properties of solute molecules.There are two models to describe the electrostatic interactions:the explicit solvent model and the implicit solvent model[31,16,32].For the explicit solvent model,solvent molecules are explicitly included.Therefore,it can realistically describe the specific interactions between solute–solvent,such as hydrogen bonds and salt bonds.However,the random fluctuation of the solvent molecules contributes a lot to the total energy of the system,thus the simulation sampling of long times must be performed to eliminate the effect of fluctuation,which dramatically increases the computational costs.The implicit solvent model does not include solvent molecules in the simulation system,it characterizes solvent effect by its average effect on the solute.Although the implicit solvent model can accelerate the simulation process by orders of magnitude,it is rather difficult to probe the specific interactions between solvent and solute in detail.

This review only presents a brief overview of the electrostatic computational methods with an explicit solvent.The more indepth discussion of these methods based on other solvation models,such as Poisson–Boltzmann and generalized Born methods,can refer to several reviews [13,31,17,33].

2.1.Truncation method

Electrostatic interactions are computed via a sum over N(N -1)/2 pairs for a system with N charged atoms,which makes the calculation complexity up to O(N2).On the other hand,electrostatic interactions are long-ranged which decays with distances as 1/r,thus a large-scale system must be considered in the calculation.To reduce the computational costs,the truncation scheme such as abrupt truncation,switching function and force shifting for calculating electrostatic interactions were proposed.The idea of abrupt truncation is to set the interaction energy between two atoms to 0 when they are computed up to a given truncation radius rcut[21].It is suitable for atom-scale simulations,especially in simulation systems with large sizes,such as polymers or biomolecules.Unfortunately,the abrupt truncation can give rise to artifacts,which may strongly impair the reliability of the dynamics or thermodynamics properties of the system [34–37].There are many methods to alleviate the artifacts arising from abrupt truncation such as the introduction of a switching function or a force shifting function,where the interaction energy between the runcation radius rsand rcutis gradually adjusted to 0 for the former,the interaction energy and its derivative are equal to 0 when r ＞rcutfor the later [38,39].Even then,artifact is still hard to avoid,unless the truncation radius is large.Another group-based cutoff method has been implemented in some software packages (such as GROMACS [40],CHARMM [41]and AMBER [42]).to speed up the simulation for larger systems.In this case,for larger molecules containing several groups,the neighboring atoms are divided into the same group.Only pair interactions between all groups are needed to calculate,thus the number of pair calculations is significantly reduced.

2.2.Reaction field methods

To improve the truncation method,the reaction field method was proposed afterwards.This method treats the pairwise electrostatic interactions within the truncation radius rcut,while the electrostatic interactions outside rcutare regarded as the effect of a specific dielectric constant εRFon charge [43,44,22].According to the reaction field method,the electrostatic potential within r ＜rcutis

Here,the first term is pairwise electrostatic interactions within r ＜rcut,the second term is the interaction formed by reaction field,and the third term makes the potential vanish at r=rcut.It is worth noting that the changes in energy or force at r=rcutare discontinuous,and a switching function(or called weighting factor)needs to be added in the simulation to eliminate this discontinuous problem.

2.3.Ewald-based method

The most popular methods for calculating electrostatic interactions based on the grid summation include Ewald [23,45],PPPM[24,46],PME [25],SPME [26]and FFP [27].

Ewald summation method has become the most commonly used method for handling long-range electrostatic interactions in computer simulations of charged three-dimensionally periodic systems [47,48].The calculation formula of the Ewald summation is expressed as Eq.(2) under the periodic boundary conditions:

In fact,Equation(2)cannot be summed directly because it converges slowly.The ideal of Ewald method is converting the slowly convergent electrostatic energy into two rapidly convergent summations(i.e.,real space and reciprocal space)and a constant term,as shown in Eq.(3):

where rij,m=ri,m-rj,m+nL and erfc are error function complements,which makes electrostatic energy rapidly converge as r increases.The first summation is short-ranged,it corresponds to the calculation of real space,in which the electrostatic energy between all pairs of charged particles within a cutoff distance rcutare calculated.The second summation is long-ranged,it corresponds to the calculation of reciprocal space,in which the summation of the electrostatics energy is obtained via Fourier transform.In this way,the electrostatic energy in Equation (2) can be expressed as

where,Ureal,Urec,and Uselfrepresents real space electrostatic energy,reciprocal space electrostatic energy and self space electrostatic energy,respectively.Uselfis used to eliminate the redundant electrostatic energy calculation of a single charge by Fourier transform when calculating electrostatic energy in reciprocal space.Among them,the calculation for Urealand Uselfdoes not occupy too much time while the calculation for the reciprocal space term Ureccosts the major CPU time,which can be expressed as follows:

The α parameter is used to tune the relative calculation weight of real space and reciprocal space.Large α value makes the erfc in real space decay fast and a lower computational cost is needed,while higher computational cost is needed in reciprocal space.Conversely,small α value makes the opposite effect.Therefore,an appropriate α value should be given to balance the computational cost between real space and reciprocal space to optimize the final computing efficiency,generally people take α=5/L.The S(n) in Equation (6) is structure factor and can be expressed as:

The original Ewald method needs to calculate many summation terms,and it will perform non-uniform Fourier transform (NFFT)multiple times during the summation in the reciprocal space,resulting in the time complexity of the algorithm reaches O (N2).The time complexity still reaches O (N3/2) even after optimizing the Ewald convergence parameters α and the corresponding truncation radius.Thus,the Ewald method is not suitable for huge charged systems.

To further reduce the complexity of electrostatic calculation,the fast Fourier transform (FFT) techniques is used to accelerate the summation of reciprocal space,which can reduce the computational time complexity from O (N2) to scale as O (Nlg N).Hockney and Eastwood proposed the PPPM method in the 1970s [24,46],wherein a general truncation was applied to the real space,while the reciprocal sum part was obtained by solving Poisson’s equation of atomic charge distribution on a grid via FFT under the periodic boundary condition.Darden et al.developed the PME algorithm,which uses Lagrange interpolation method to distribute charges onto the grid [25].Later,Esmmann et al.used cardinal B-splines interpolation to smooth charge distribution on a grid,which was called the SPME algorithm [26].Although the PPPM method is the first method to efficiently handle long-range electrostatic interactions in molecular simulations,this method is largely replaced by the PME/SPME method in the atomic scale simulations.Because the energy of reciprocal space only needs perform single FFT in SPME while it needs 4 FFTs in the PPPM method.In addition,the calculation speed of PME/SPME algorithm is much faster than the original Ewald method for large systems,while the Ewald method may be better for the small systems since it avoids establishing atomic charge distribution on grid.In particular,another different grid-based method is the FFP method proposed by York and Yang[27],which directly samples the compensating Gaussians onto the grid.This method is more suitable for higher accuracy requirements than the above methods.However,the calculation cost is much higher than the other grid-based methods.

It is worth noting that employing the conventional Ewald summation technique becomes more challenging for simulation systems with the slab geometry,such as polymer layers,selfassembled monolayers and charged nanomaterials,since there is no periodicity in the z direction of the simulation box[49].A series of two-dimensional Ewald summations (EW2D) were introduced into slab geometry systems for exploring the electrostatic interactions [50–52].However,the Bessel function used in 2D Fourier transforms cause the efficiency of EW2D not nearly as high as the conventional three-dimensional Ewald summation (EW3D)methods [51].Another widely-spread method is to still use the EW3D algorithm with adding a sufficiently large vacuum layer between periodic replicas in z direction to avoid an artificial influence from the periodic images in the z direction [53].However,Spohr et al.found that the electrostatic potential still cannot converge to a satisfactory result even when the length of the vacuum layer was five times larger than the simulation cell in x or y direction [52].Later,Yeh et al.[54]proposed a modificatory threedimensional Ewald summation with shape-dependent correction term (EW3DC) to calculate the electrostatic interaction of the slab geometry.This EW3DC technique could significantly reduce the computing time and obtain a highly consistent results with that from simulations using the rigorous EW2D.

2.4.Electrostatics in mesoscale modeling

Due to the extreme complexity of biological systems,mesoscopic coarse-grained molecular dynamics (CGMD) simulations are often adopted.In recent years,the MARTINI CG force field developed by Marrink et al.[55]has been widely used in the simulation studies of biomolecular systems,especially biomembranes.In the CG force field,the truncation method with a shift function is used for processing electrostatic interactions,and a dielectric constant of 15 was used to compensate the neglect of explicit polarization of the MARTINI water model,which can implement consistency between simulation results and experimental data.However,in the interactions between biomembranes and nanomaterials whose size exceeds the typical cutoff value,the long-range electrostatic interaction becomes more critical and special treatment is required.In this case,PME method is often adopted to replace the truncation method,which can better realize the accurate calculation of long-range electrostatic interactions.However,it should be noted that the introduction of the PME method will undoubtedly reduce the simulation speed compared with the truncation method.

Dissipative particle dynamics(DPD)method is another popular CG method applied in the studies of biomembranes.In DPD simulations,the particles interact with each other by the soft repulsion potential,which allows the non-bonded particles to overlap.Therefore,it is highly possible for two oppositely point-charge species to form an artificial ion-pair.To avoid the formation of ionpair,one general method is to smear out the charges of the DPD particles by some kinds of distributions.The slater-type charge density distribution is often used [56,57],which can well describe electrostatic interaction between charged species at the mesoscopic scale.

2.5.Other methodological developments

It is worth mentioning the evolutive methodologies which use existing computer architecture tool to accelerate the calculation of electrostatic interactions.One of the recent developments was the extension of LAMMPS’ PPPM algorithm for the secondgeneration Intel Xeon Phi [29].This optimized PPPM method is mainly dedicated to solving the mapping steps including charge distribution and the Coulomb force calculation process by the means of vectorization techniques.It was shown that the PPPM long-range solver could achieve 2–3 times speedup for the across a wide range of cutoff radius.In recent years,the rapid development of graphics processing unit (GPU) has shown stronger computing power,which provides new possibilities to improve the computing power of MD simulations [58,59].At present,mainstream MD software such as GROMACS,NAMD,AMBER,CHARMM etc.,they all use PME algorithm as one of the basic algorithms in processing electrostatic interactions,which provides a superior performance than the traditional CPU-based implementations[60,61].

Another interesting development was the birth of Anton 2,a special-purpose supercomputer designed and built for MD simulations [30].The Anton 2 platform integrates the overlap of computation with communication,in which possesses a wide range of new algorithms,optimized software and hardware to run efficiently.The Anton 2 platform originally used Gaussian split Ewald method [62]to handle long-range electrostatic interactions,a different decomposition called the u-series was implemented on it to gain the additional performance.The most noteworthy is that Anton 2 can achieve simulation rates of multiple microseconds per day for a system containing millions of atoms.

3.Regulation of Electrostatic Interactions on the Behavior of Different Charged Systems

3.1.Electrostatic interactions in ion systems

Electrostatic interactions play a critical role in ion hydration,ion transport and interfacial phenomena in electrolyte solution.Ions are coordinated by a shell of water molecules,driven by the electrostatic interactions between ions and water[63,64],i.e.positively charged cations are directly coordinated by the negatively charged oxygen atoms of water with both hydrogen atoms pointing away,while negatively charged anions are directly coordinated by the positively charged hydrogen atom of water with the oxygen atom pointing away.Thus,the water molecules in ionic hydration shell show specific orientations.In addition to coordination number of ions,molecular dynamics(MD)simulations by Zhou et al.[65]have demonstrated that the orientation of water molecules around ions is another important factor for ionic hydration and proposed a hydration factor to quantitatively characterize this orientation.In series of technological applications,such as ion separation and electro-deionization,electrostatic interactions of ions can be used to drive ion transport by external electrical fields.Recent MD simulations [66–69]indicated that under strong electrical fields,the structure,dynamics and transport properties of ionic solution became anisotropic;ion mobility was significantly enhanced and the first hydration shells of ions were dramatically weakened,as well as ion-water interaction energy was significantly reduced,mainly due to the serverely distrurbed water orientaion in the hydration shells.Anomalous interfacial phenomena of alkali halide solutions at air/water interface were observed[70–72],due to different electrostaic interactions arising from specific ion properties,including ion size,polarizability and hydration strength.The MD simulation work by Feng et al.[73]found that Cl-,Br-,I-and Cs+are expelled to the air/water interface and show higher density than in bulk.Such surface preference of ions is also observed near the interface between solution and the walls of carbon nanotubes(CNTs) [74,75].For fluoride solutions,the preference of cations at the air/water interface was found beacuse F-had a stronger hydration ability than most alkali metal cations [76].

In rencent years,ionic hydration under the nanoconfinement of nanochannels has drawn wide attention for many potential applications,such as water desalination and ion separation.The behaviors of water molecules and ions under nanoconfinement have been summarized in some review papers[77,78].For ion hydration under nanoconfinement,some anomalous phenomena will emerge.The MD simulation by He et al.[79]found that ice-like water structures with well-organized dipole orientation could form in CNT(8,8) at 300 K,once the cations of Na+and K+were hydrated by such ice-like water,the dipole orientions of such ice-like water would concertedly change and induce even more favorable electrostatic interactions with cations than the bulk counterpart.Similarly,the simulation work by Strauss et al.[80]also observed such ultralong polarization chains induced by hydrated ions in confined water monolayer between graphene sheets.In addition,it was reported that in CNTs with diameters around 1 nm,the preferential orientation of water molecules in the hydration shell of Na+and K+presents an anomalous change and causes a diameter-dependent variation in ion-water interaction [81].The MD simulation work by Shao et al.[82]explored the effects of high temperature on ionic hydration in CNTs,they found that the different charges on anions and cations and their hydration strength will affect the ability of ions to resist the disturbance of high temperature.Zhu et al.[83]performed MD simulation to investigate the effect of inner wall modified carbonyl group on the hydration of cations in CNT and revealed that modification groups could enhance or weaken cation hydration depending on CNT diameter.When confined in nanochannels,ions are often dehydrated,namely,the second hydration shell and the outer hydration shell are stripped;when in extremely narrow nanochannels,the first hydration shell would also be partially stripped,consequently,the electrostatic interaction between ions and water was weakened [84,79,85].

The change of ion hydration under nanoconfinement could lead to some potential technological applications.On one hand,in narrow nanochannels,ions were partially dehydrated and the resulting weakened ion-water interaction would impose high energy barriers and then hinder ion transport across nanochannels [86].The energy barriers for anion transport in narrow non-polar pores were quantified by Richards et al.with MD simulations[87].Based on such ion-dehydration mechanism,hydrophobic nanochannels,such as CNTs,were proposed to effectively reject ion and used for efficient water desalination[88].The simulation work by Corry found that modification of charged groups at the entrance of CNTs could help to prevent the passage of ions which is essentially benefitted by the electrostatic interaction between ions and charged groups[89].Inspired by the hydrophobic gating mechanism in biological ion channels,He et al.[90]designed a CNT-based nanogate,where partial dehydration in the continuous water-filled hydrophobic constriction in the tube center could effectively close the channel for ion conduction.Moreover,different ions face different energy barriers when crossing the same narrow nanochannels due to their different dehydration energies,thus,narrow nanochannels,such as CNTs and graphene nanopores have intrinsic ion selectivity and could be used for ion separation [91,92].On the other hand,in functionalized nanochannels,the polar or charged groups could substitute some of water molecules in ion hydration shells to directly coordinate ions and change the transport behavior of ions.Inspired by the selectivity filters of KcsA K+and NavAb Na+channels,He et al.[1]used MD simulations to design bioinspired nanopores in graphene sheets by functionalizing the nanopore rims with carbonyl or carboxylate groups,as shown in Fig.1.Under external electrical fields,these biomimetic nanopores could selectively conduct K+or Na+,where the electrostatic interactions between ions and functional groups in the nanopore play a critical role.In the nanopore modified with carboxylate groups,Na+binds to the negatively charged carboxylate groups and block the nanopore,then another Na+comes to the nanopore and the electrostatic repulsion interactions between them make the Na+leave the nanopore,such ‘‘knock on”ion conduction also occurs in biological channels.Similarly,Lu group designed biomimetic graphene nanopores for Mg2+/Li+separation via mimicking the selectivity filter of Mg2+channel and revealed that the dehydration of the second hydration shell and the electrostatic interactions between ions and functional groups played an essential role[93,94].In addition,the electrostatic interactions from external charges outside of nanochannels have specific effects on the conduction of water and ions across nanochannels.Fang et al.[95]performed MD simulations to devise a electrostatic gating for the single-file water transport in CNT(6,6) by adjusting the distance of external point charge to the nanochannel.In the similar way,Gong et al.[96]designed a molecular switch to tune the selective transport of cations or anions across CNT(9,9)via changing the distance of external point charge.Nevertheless,the outside static charges cannot drive a continuous flow of water molecules through CNTs [97].

Fig.1.(A,D) Selectivity filters of the KcsA K+ and NavAb Na+ channels,which consists of backbone carbonyl groups and carboxylate groups from glutamate residue side chains,respectively.(B,E) The 4CO and 3COO nanopores are functionalized with four carbonyl and four carboxylate groups.(C,F)Ion conduction modes in 4CO and 3COO nanopores [1]

Liquid transport in micro-/nanoscale channels is of fundamental importance.Electroosmotic flow (EOF) is the motion of ionized fluids under externally applied electric fields.The presence of charged surfaces and salt ions in EOF may significantly alter the transport properties such as viscosity and slip length in nanochannels.Through MD simulations,Celebi et al.[98]demonstrated that nanoconfinement induced by walls did not significantly alter the apparent viscosity of the ionic solution,but wall surface charge increased the apparent viscosity compared with the uncharged surfaces;the slip lengths for the charged surfaces are smaller than the neutral surfaces.In addition,they concluded that slip length was normalized with the Debye length which could be adjusted by ionic concentration.Jafari et al.[99]simulated the non-Newtonian EOF through nanochannels,their results shown that the shape and magnitude of such achieved velocity profiles would be remarkably affected by the ratio of one-half of the channel height to Debye length,electric field strength,zeta potential and the flow behaviour index values.

Franz Hofmeister first observed that certain inorganic salts and ions showed different abilities of precipitating proteins,the attendant Hofmeister series of anions were proposed [100].The Hofmeister effect has a great influence on chemical and biological systems,such as colloidal system stability,microemulsion systems,solution viscosity and surface tension,precipitation and crystallization of polymer electrolytes,properties of ionic surfactants,enzyme activity,protein stability,amino acid optical rotation and etc [101].For more details,please refer to the relevant review[102].

3.2.Electrostatic interactions in inorganic non-metallic materials

Piezoelectric materials are a kind of information functional material that can convert mechanical energy into electrical energy and vice versa [2].It plays an important role in medical imaging,communication,non-destructive testing,high-precision positioning etc.[3].The piezoelectric effect means that a potential difference is generated on the surface of the materials under an externally applied mechanical pressure and vice versa,the materials were deformed by an applied external electric field [103].MD simulations technology have been used to predict and calculate the piezoelectricity of semiconductor materials such as ZnO,BN,MoS2,GaN etc.For instance,Momeni et al.[104]used MD simulation to investigate the effect of size scale on the piezoelectricity of ZnO nanobelts.The results indicated that the charge redistribution near the free surface of ZnO nanobelts resulted in the decreased piezoelectric coefficient with the increase of lateral dimensions.

The alternating arrangement of boron and nitrogen atoms in the boron nitride lattice makes it excellent in piezoelectric properties.The MD simulation work by Nan et al.investigated the influence of shape and size of monolayer hexagonal boron nitride nanosheets on its piezoelectric properties [105].It was shown that the piezoelectric constant of monolayer hexagonal boron nitride nanosheets was closely related to its macroscopic shape,but nearly independent of its macroscopic size.The piezoelectric behaviors of the zigzag direction of the boron nitride honeycomb (BNHCs) structures were studied by Xie and his co-workers via MD simulation [103].It was found that the piezoelectric coefficient of the multilayer BNHCs was smaller than that of the monolayer BNHCs as opposite magnetic poles in the multilayer BNHCs material cancel each other more than the monolayer BNHCs.Zhang et al.[106]carried out MD simulations to explore the influence of piezoelectric effect on the resonance frequency of boron nitride nanosheets.The results demonstrated that the resonance frequency of boron nitride nanosheets could be tuned by applying an external electric filed.An interesting phenomenon was observed,the resonance frequency of the boron nitride nanosheets with odd layers moved upward under a negative electric field and downward under a positive electric field;it was reduced under both negative and positive electric field.

Molybdenum disulfide (MoS2) and gallium nitride (GaN) are another two attractive piezoelectric materials that have widely applications in electronics and optoelectronics.MoS2has three different stacking configurations including 1T,2H and 3R-MoS2structure.Tan et al.[107]calculated the piezoelectric coefficient of 3RMoS2multilayer structures.Their results showed that the 5-layer 3R-MoS2structures had the highest piezoelectric coefficient(0.457 C·m-2)of all reported MoS2single and multilayer structures since the 3R-MoS2configuration lacks inversion symmetry and thus leads to a robust large polarization for multilayer structures.The size-and shape-dependent piezoelectric properties of GaN nanowires (cylindrical,square and hexagonal nanowire) and nanotubes (hexagonal nanotube) were studied by Jiang’s research group [108].The results indicated that the piezoelectric material parameter decreased with the increase of the crystal size and shape factor of GaN.For the square and hexagonal GaN nanowires,the influences of size and shape on the piezoelectric properties could be neglected when the size was higher than 41.5 nm and 70.5 nm for square GaN nanowires and hexagonal GaN nanowires,respectively.

The piezoelectric properties of piezoelectric materials depend strongly on the shape,size,temperature,surface properties and electrostatic interaction.The current studies mainly focus on the shape and size of the material.Besides the size,the shape and the temperature,the piezoelectric properties of semiconductor materials can also be tuned by surface modification or alloying.

3.3.Electrostatic interactions in organic acids

Organic acids are widely used in various key areas,such as mineralized solutions in biomineralization[109],electrolytes in batteries [110]or capacitors [111],buffering agent and eluent [112].

In recent years,people pay more attentions to the interaction between citrate and various solid surfaces due to their special application in nanotechnology or biomineralization.Monti et al.[113]studied the adsorption mode of trisodium citrate on gold nanoparticles (AuNPs) by MD simulations,their results showed that citrate ions were able to adsorb on the AuNP surface in a fully deprotonated form with its strong electrostatic interaction and altered the surface charge to confer stability to the metal support.The MD simulation work by Perfilieva et al.[114]investigated the relationship between trisodium citrate and AuNPs with polarization effects.They found that the water molecules closest to the AuNPs caused the redistribution of AuNPs,then the negatively charged citrate ions were attracted to the positively charged edges of AuNPs through electrostatic interactions,formed a network around AuNPs.Citrate preferred to adsorb on the Au (111) face over Au (100) face as the negative field above Au (111) surface was weaker than Au (100).It is worth mentioning that citrate is also a key factor in the regulation of HAP crystal growth.It showed an inhibitory effect on HAP crystal growth due to the adsorption of negatively charged citrate ions.Wang et al.[115]adopted metadynamics simulations to investigate the adsorption mode of citrate on the HAP surface.They found that citrate was able to adsorb Ca2+on the HAP surface with its terminal carboxyl group,the interfacial electrostatic interactions was the key factor,the interaction of citrate with HAP (100) face was also stronger than that with HAP (001) face.

Oxalic acid plays different functions in different living organisms and can adsorb on many minerals.Xue et al.[116]probed the adsorption behavior of organic acid on minerals surface based on MD simulations.The adsorption equilibrium states of oxalic acid on feldspar surface with three different configurations are shown in Fig.2.With the deprotonation of oxalic acid,the driving force for adsorption would transform from weak hydrogen bond(system A)to strong electrostatic force(systems B and C).In addition,Biriukov et al.[117]explored adsorption of oxalic acid ions on the rutile (110) surface.Their results were consistent with the experimental results for all simulated charge densities from neutral (pH=6) to +0.208 C·m-2(pH=3.7);the outer-sphere complexes were the most favorable,while other inner-sphere complexes were unstable.Furthermore,they found that nominal charges exaggerated ion–water interactions in their simulation systems.

Fig.2.Adsorption equilibrium states of oxalic acid with three different configurations on feldspar surface [116].

In most studies,benzoic acid represents one of the model aromatic carboxylates used to investigate the adsorption of emerging contaminants on different environmental surfaces.The recent MD simulation work [118]reported the adsorption of benzoic acid derivatives on charged silica surfaces.The pulling simulations in NaCl solutions disclosed that all benzoic acid derivatives were adsorbed on the silica surface at pH 2–3,while the strength of interactions decreased at pH 7 owing to electrostatic repulsion.More interestingly,with the help of Na+-mediated attractive interactions,the interactions increased again at pH 9–10 in case of 1,2,4-benzenetricarboxylic acid and phthalic acid.The dissolution of organic crystals appeals to the pharmaceutical sector.Toroz et al.[119]investigated the lifetime and stability of representative molecular-clusters of Para-amino benzoic acid (PABA) in aqueous solution,which had two polymorphic forms and constituted an enantiotropic pair with MD simulations.It revealed that the selection of a suitable electrostatic potential was very important for the description of the dissolution process of the nano-sized crystals.

Salicylic acid permeation was studied by employing a nitrocellulose membrane and modelled via DPD simulation regarding drug delivery in the work of Daniel et al.[120]The result showed that charged salicylate always diffused faster than neutral salicylic acid,which demonstrated a consistent behaviour irrespective of the proportion of charged and neutral particles.Moreover,it validated the DPD model as a discriminatory mechanism to study the effect of pH or charge of particles on diffusion rates.To demonstrate the possibility of TiO2nanotubes as drug reservoir in medical applications,Yang et al.[121]performed MD simulation and DPD simulation to explore the release process of aspirin and vitamin C from TiO2nanotubes.They found that aspirin molecules and vitamin C were captured by polylactic acid nanoparticles (the encapsulating material of TiO2nanotubes)due to the electrostatic and hydrophilic/hydrophobic interactions of PLA and aspirin molecules/vitamin C,which retarded the release of drugs.

3.4.Electrostatic interactions in ionic liquids

Ionic liquids(ILs)are recognized as green solvents for industrial and pharmaceutical applications owing to their unique properties such as nonvolatility,high thermal stability and high ionic conductivity [122–124].These properties are mainly determined by the strength and directionality of Columbic interactions between cations and anions.In this section,we briefly review some recent molecular simulation works on ILs,the more detailed simulation reviews have been provided by Zhang group and Voth group[125,126].

Cellulose is abundant on the earth and its conversion to biofuels provides a promising way to satisfy the ever-increasing global energy demand.However,cellulose is not soluble in water or other common organic solvents and this represents one major obstacle to the cellulose-biofuel conversion process.ILs have been demonstrated to be effective solvents for the dissolution of cellulose.A number of experimental and simulation studies have been carried out to investigate the mechanism of cellulose dissolution in ILs.It is generally accepted that the hydrogen bond formation between IL anions and the hydroxyl groups of cellulose is the driving force for cellulose dissolution in ILs[127].However,it remains unclear how the cations and anions of IL affect the dissolution process.Liu et al.developed an atomistic force field for the 1-ethyl-3-methylimidazolium acetate ([C2mim+][OAc-]) IL and they investigated the behavior of a series of (1–4) linked β-D-glucose oligomers in [C2mim+][OAc-].The [OAc-]anions were found to form strong hydrogen bonds with the hydroxyl groups of the cellulose oligomers.Besides,the [C2mim+]cations were found to associate with the cellulose oligomers through hydrophobic interactions.The simulation results support that cations also play significant roles in the dissolution of cellulose in [C2mim+][OAc-]IL [128].Zhao et al.carried out MD simulations to study the dissolution of cellulose oligomer in a series of ILs with varying heterocyclic structures and alkyl chain lengths[129].It was shown that,for the imidazolium based ILs,the shorter the alkyl chain was,the higher the cellulose solubility was.The simulation results also found that the presence of electron-withdrawing allyl groups in the alkyl chain of the cation dramatically enhanced the interaction between the cation and the cellulose,thus promoting the cellulose dissolution.To explore the underlying mechanism of cellulose dissolution only in ILs with unsaturated heterocyclic cations,Li et al.simulated the dissolution of the cellulose bunch in four kinds of ILs composed of unsaturated and saturated cations:1-butyl-3-methylimidazolium,1-butylpyridinium,1-butyl-1-methylpyrrolidinium,and 1-butyl-1-methylpipperidinium and the acetate anion [130].The simulations confirmed that cellulose was only soluble in ILs with unsaturated cations.The π electron delocalization of the unsaturated heterocyclic ring made the cation more active to interact with cellulose and provided more space for acetate anions to form strong hydrogen bonds with cellulose.Besides,the large volume of cations with the saturated heterocyclic ring resulted in a slow transfer of both cations and anions,which inhibited the dissolution of cellulose.Li et al.also performed MD simulations to explore the dissolution of cellulose bunch in 1-ethyl-3-methylimidazolium acetate ([Emim+][Ac-]),1-ethyl-3-methylimidazolium chloride([Emim+][Cl-]),1-butyl-3-methylimidazolum chloride ([Bmim+][Cl-])and water[131].Complete dissolution of the cellulose bunch was observed only in [Emim+][Ac-].They proposed a synergistic role of the cations and anions in dissolving the cellulose.Initially,cations bind to the cellulose bunch and anions insert into the cellulose strands to form hydrogen bonds with the hydroxyl groups.The cations then start to intercalate into cellulose chains due to the strong electrostatic interactions with the inserted anions.The hydrogen bonds formed between Cl-and hydroxyl groups cannot effectively separate the cellulose chain and this explains why[Emim+][Cl-]and [Bmim+][Cl-]ILs dissolve cellulose bunch much more slowly.

CO2capture and storage has been considered as one promising way of mitigating greenhouse emissions.Among all emerging materials for CO2capture,ILs have attracted much attentions.Molecular simulations have been frequently used to study the interactions of CO2with ILs to uncover the molecular determinants for the high solubility of CO2in ILs.For instance,Maggin and coworkers have investigated the factors which determines the high CO2solubility in imidazolium based ILs (i.e.,[Bmim+]and [Bmmim+][132].They found that CO2organized strongly around theanions but was more diffusely distributed around the two different imidazolium rings.Their simulation results showed that the major factor controlling CO2solubility in the imidazolium ILs was the chemical nature of the anions.In the MD simulation study by Yang et al.,the authors studied the solvation of CO2in four ILs:[Emim+][Bmim+][Emim+][Cl-]and [Bmim+][Cl-][133].Their data showed that the van der Waals attractions between CO2and ILs dominated the CO2solubility in ILs,and thus cations with longer alkyl chain and anions with higher hydrophobicity could enhance the CO2solubility.It is predicted from the simulations that cations with long alkyl chain,anions with high hydrophobicity and ILs with smaller surface tension are effective in capturing CO2.To elucidate the molecular determinants for the solubility of CO2in imidazolium based ILs,Klahn et al.simulated different pure and CO2containing imidazolium ILs with varying alkyl chain length and anions of different size and hydrophobicity [134].It was found that the solubility of CO2in ILs was not significantly influenced by the CO2-ion interactions but the amount of free volume in the ILs.ILs with relatively weaker cation–anion electrostatic interactions tended to have larger free volume to accommodate CO2.

Proton exchange membrane fuel cells (PEMFCs) are promising technology for generating clean power for use in electric vehicles and portable devices.The key component of PEMFCs is the proton exchange membrane (PEM) sandwiched between two electrodes.The PEM should meet several requirements to be applied in PEMFCs:(1) high proton conductivity in both dry and wet conditions;(2) good mechanical strength;(3) good chemical and thermal stability;(4) cost effective.The Nafion polymer is the most widely used material in PEMFCs.However,the proton conductivity of Nafion drops at elevated temperatures above 100 °C due to the evaporation of water.ILs are expected to replace water and ILs incorporated Nafion membrane can be used under high temperature conditions.Molecular simulations have been reported to investigate the ILs-Nafion interactions as well as the morphology of the materials.Sunda simulated the Nafion membrane doped with the diethylmethylammonium triflate ([Dema+][TfO-]) IL[135].The hydrogen bonding interactions were found to predominantly exist between the N-H group of the ammonium cation and thegroup of the anion or the Nafion polymer.An increase in the IL concentration and the temperature results in faster diffusion of the ionic liquid and higher proton conductivity.Sun et al.reported all-atom and coarse-grained simulations of Nafion/[Bmim+]composite membranes [136].The CG models can well predict the structural and diffusion properties of the[Bmim+]IL in the membrane.It was found that IL formed nanocluster in the membrane,the continuous IL cluster might act as a pathway for the proton hopping in the membrane.In addition,Mai and Zhou et al.[137]employed DPD simulation to explore the self-assembled meso-structure and related properties of Nafion/[Bmim+][TfO]composite membranes at different ionic liquid concentration and temperatures.The results indicated that the state of ionic liquid transformed from dispersed phase to the continuous IL channels in the Nafion/IL membrane with the increase of the ionic liquid concentration.Moreover,the size of the formed channels became smaller and more concentrated as the temperature increased.The distribution of sulfonic acid groups in Nafion side chain had an important influence on distribution of anion and cation groups of IL at the microphase interface,and thus affected its microscopic phase behavior.

Deep understanding of the interaction mechanism between ILs and material surfaces is of great research significance for improving the application of ILs in electrochemical energy storage and electrocatalysis.Vatamanu et al.investigated potential and temperature dependences of the differential capacitance of ILs at the graphite electrode surface [138].The results showed that both cations and anions near the electrodes surface were primarily perpendicular to the surface as the electrode potential equal to zero.Interestingly,the first layer counterions increasingly parallel to the electrodes if increasing the electrode potential.Their findings indicate that regulating the ion shape and electrostatic interactions between ILs and electrodes surface is of great importance to control the behavior of differential capacitance.Shim et al.studied the charge distribution of the internal layer of carbon nanotubes(CNTs) supercapacitors in the ILs via MD simulations [139].They found that the extended charge distribution of bulky ILs counterions led to the orientational alignment,which formed a charge density of multiple layers with alternating signs inside the electrified CNTs.In addition,Feng et al.also investigated the capacitive behavior of carbon-based supercapacitors with ILs[140].They found that more counterions moved to the charged carbons electrode surface as the temperature increases,whereas the number density of counterions slightly changed for CNTs supercapacitors,and eventually resulted in a temperature sensitive carbon supercapacitor and temperature-resistant CNT supercapacitor.

The high viscosities of ILs may lead to problems in chemical processes,such as the reduction of the transfer rate of mass and heat in the reaction and separation processes.Therefore,viscosity is one of the most important properties to consider when ionic liquids are used in chemical processes.Developing computational approach to predict the IL viscosity has attracted the attention of a few research groups.Yan et al.[141]have developed an electronically polarizable model for the 1-ethyl-3-methylimidazolium([Emim][NO3]) IL.Their MD simulation results showed that,using the classical fixed charge model for [Emim][NO3],the computed viscosity was about 50% higher than the experimental data.However,when using the polarizable model,they derived a viscosity value,which dropped to within about 7% of the experimental value.Micaelo et al.[142]have developed united-atom models for 1-butyl-3-methylimidazolium hexafluorophosphate ([Bmim][PF6]) and 1-butyl-3-methylimidazolium nitrate ([Bmim][NO3])ILs based on the GROMOS96 force field.They run MD simulations of the two ILs at temperatures ranging from 298 to 363 K.A nonequilibrium periodic perturbation method was used to calculate the viscosity.Using the united-atom models,their calculated viscosities agreed very well with experimental data.Maginn and co-workers [143]have utilized non-equilibrium MD simulation approach to predict viscosity of 1-ethyl-3-methylimidazolium bis(trifluoromethanesulfonyl)imide ([Emim][NTF2]) as a function of water content and temperature.Using a fixed-charge model and the simulation method,they obtained very good agreement between predicted viscosity and experimental data.

3.5.Electrostatic interactions in ionic surfactants

Ionic surfactants mainly contain anionic and cationic surfactants,which are generally composed by a positively or negatively charged head group and a long chain alkyl tail.It is widely used in a diverse range of industrial applications such as cosmetics,mining flotation,material manufacture,enhanced oil recovery and pollutant removal owing to their unique amphiphilic molecular structure [9].A large number of experimental studies have reported the changes in the thermodynamics,kinetics,surface tension,zeta potential,critical micelle concentration (CMC) of ionic surfactants adsorbed on different interfaces [5,6,144].The structure of ionic surfactant,mixed surfactants,counterions,salt,ionic concentration and other factors affect the adsorption orientation and adsorption stability of ionic surfactants,but experimental methods cannot explain the behaviors of ionic surfactants on the interface or surface at the atomistic level.Therefore,computer simulation technology is necessary to investigate the structural changes and micro-aggregation behavior of ionic surfactants at different phase interfaces,which is conducive to better design and control industrial production processes.This section reviews the effect of electrostatic interaction on the self-assembly of ionic surfactants in solution or the adsorption of ionic surfactants at different interfaces.

Considerable multi-scale simulation studies have proved that electrostatic interaction plays a critical role in the solution selfassembly or the adsorption of ionic surfactants at different interfaces [7,145,8,146–148].The structures of ionic surfactants have a significant effect on its properties.Zhang et al.[7]studied the influence of the head group structure on the aggregation behavior of four surfactants at the oil/water interface.The simulation results showed that the oxygen atoms,benzene rings and oxyethyl groups in the head groups would further reduce the interfacial tension of the system,which indicated that the three groups were beneficial to enhance the activity of surfactants.Liu et al.[145]explored the aggregation behavior of sodium alkylbenzene sulfonate surfactant molecules with different hydrophobic groups at the oil/water interface to study the effect of the difference in tail chain structure on the properties of surfactants.It indicated that the steric hindrance and electrostatic repulsion between the surfactant molecules increased with the increase of the length of the alkyl chain adjacent to the sulfonic acid group,which in turn made the arrangement at the interface looser and reduced the adsorption at interface.Besides,modifying the structure of ionic surfactants can improve their interfacial activity at the oil–water interface,the mixed surfactants also affect the molecular array behaviors at the interface.Zhao et al.[8]investigated in detail the adsorption behavior of the mixed system containing green surfactant alkyl polyglycoside (APG),anionic surfactant sodium dodecyl sulfonate(SDS) and cationic surfactant dodecyltrimethylammonium bromide (DTAB) at the dodecane-water interface by MD simulations.It was found that the zeta potentials of the emulsified oil droplets changed with the mole ratio of APG:DTAB:SDS.The interface charge was mostly close to neutralization state when the mole ratio was equal to 1:1:1,which reduced the electric repulsion,promoted the dense arrangement of surfactant molecules in the interface layer,thus improved the oil–water interface activity.

In the aqueous solution,ionic surfactants can generally form the larger micelle by self-assemble when its concentration exceeds the CMC.However,the stability of the charged micelle may be significantly reduced because of existing strong electrostatic repulsive interaction between the same type of charged micelles,which further affects severely the stability of the foam layer formed.It was reported that the addition of counterions may change the interaction among ionic surfactants and two immiscible phases by means of electrostatic interactions with the charged headgroup of ionic surfactants and water molecules,which leads to the interfacial adsorption layer’s structure and properties changed,such as interfacial tension,zeta potential and thickness [147].Lima et al.[146]investigated that the effect of six hydrotropic counterions on the structure and properties of cationic micelles composed of DTAB through MD simulations.Their results revealed that the nonpolar group of the counterion,rather than the polar group determined micellar properties.The desolvation of the non-polar group drove the expansion of the adsorption of ions at the micelle interface,hence the structure of the micelle was maintained stable and not destroyed.Moreover,the location of the counterion,and its orientation and the frequency of ion exchange between ‘‘a(chǎn)dsorbed”and ‘‘free”states were also determined by the non-polar group.

To explore the contribution of various chemical species and moieties to the surface tension distribution in these aqueous solutions at various surface coverages.Hantal et al.[147]considered the liquid-vapor interface of five different amphiphilic molecule solutions,representative of anionic,cationic and nonionic (alcoholic) surfactants.The results confirmed that counterions played a key role in contributing to the surface tension in case of two ionic surfactants,namely,the anionic sodium dodecyl sulfate (SDS) and DTAC,whose effect depended on their ‘‘hardness”within the Hofmeister series though there was a large compensation between the ionic surfactants’ charged headgroup and counterions.It is worth noting that counterion ratios and types should be more carefully considered in developing methods to manipulate rheological properties of interfaces and of colloids in general.Riccardi et al.[148]calculate the interfacial tension,Debye length,ζ potential of bis(2-ethylhexyl) succinate sulfonate (AOT) surfactant at dodecane-water interface under constant ionic strength and different Ca2+/Na+molar ratios.They found that the molar ratio of Ca2+/Na+ions was an important parameter for enhancing the interfacial activity of ionic surfactants.Adding Ca2+ions and other counterions could control the viscoelasticity of emulsions containing ionic surfactants.In addition,Sangwai et al.[149]showed accurately micellar assemble process of cetyltrimethylammonium chloride(CTAC) at the existence of sodium salicylate and sodium chloride using CGMD with the MARTINI force field and explicit CG solvent.Their simulation results showed that the phase interface was nearly saturated with adsorbing amphiphilic salicylate ions when the molar ratio of salicylate ions to CTAC was larger than 1,which caused a reduction in micelle/water interfacial tension,and thus induced a sphere to rod transition in the micelle.

The interplay and synergetic effects of many factors involving the structure of ionic surfactants,the surface charge of the substrate,the hydrophobicity,hydrophilicity or the charged solid surface as well as solution pH often render the adsorption of ionic surfactants at the solid/water interface more intricate than at the air/water and oil/water interfaces.It is widely accepted in literatures that the surface-surfactant headgroup electrostatic interaction is a major driving force of adsorption of ionic surfactants on charged substrates.Typically,it seems that ionic surfactants like cetyltrimethylammonium bromide (CTAB) are more inclined to attract the oppositely charged NPs (silica NPs) via electrostatic interaction,thus forming a more stable emulsion.On the contrary,the same charged ionic surfactants and NPs generate electrostatic repulsion,resulting in an unstable emulsion.In addition,it is well known that low salinity aqueous solutions with proper ratio of monovalent to divalent ions can change the charge distribution between two interfaces.By changing the surface charge (zeta potential) of NPs,Wang et al.[150]successfully implemented MD to simulate the detachment of oil from the solid substrate and claimed that nanofluids containing hydrophobic NPs could possess better functionality to alter the wettability of the rock toward more water-wet.Recent research by Aminian and ZareNezhad [151]regarding the molecular simulation of oil detachment from carbonate calcium surface has also demonstrated that a mixture of TiO2NPs,AOT anionic surfactant and aqueous solution of CaCl2,MgSO4,NaCl and KCl salts leads to the maximum synergism such that the lowest contact angle leads to the maximum oil detachment from the calcium carbonate surface.

Amphiphilic surfactant molecules self-assemble and form various micellar structures,the variety of both micellar structures and rheological behavior makes the surfactant solutions have wide applications in industries.MD simulation provides a powerful technique to explore the correlation between rheological properties and micelle structures.Liu et al.[152]studied the rheological behavior in a binary system of cetyltrimethylammonium chloride(CTAC) and counter ion aqueous solution.The simulation results showed that the shear viscosity of the surfactant aqueous solution had a plateau at low shear rates and a shear-thinning regime at high shear rates;the shear viscosity increased with the increase of the concentration ratio of the counter ion salt to CTAC surfactant.Zhao et al.[153]investigated the rheological properties of a pH-sensitive clean fracturing fluid system formed by the anionic surfactant N-hexadecyl iminodiacetic acid (HIA) and salt-free cationic hexadecyl trimethyl ammonium bromide (CTAB).They demonstrated that electrostatic repulsion among the carboxyl groups of HIA was screened by the positively charged heads of CTAB,which could lead to an increase of the solution viscosity.The excellent rheological properties and pH sensitivity can help increase the productivity of wells after fracturing and control the construction cost.

3.6.Electrostatic interactions in polyelectrolytes

Polyelectrolytes are defined as the polymers with ionizable units which dissociate in water,making the polymers charged[154].The polyelectrolytes possess the properties of both electrolytes and polymers,such as electrostatic interaction of electrolytes,viscoelasticity of polymers,electrical conductivity of ions,hydrophilicity of the material etc.,which make polyelectrolytes exhibit unique peculiarities,for instance,polyelectrolyte effect,anti-polyelectrolyte effect,self-assembly etc [155].Owing to the special properties and applications,polyelectrolytes have always been a research hotspot for scientists,not only in experiment for synthesis applications but also in computer simulation for property exploration [156–159].In fact,in the research of polyelectrolytes,the electrostatic interaction which dominates the peculiar behaviors of the polyelectrolytes system is nonnegligible.Thus,in computer simulations,one should pay more attention to the electrostatic interaction in polyelectrolyte systems.Generally,according to the charge of ionized macromolecules,polyelectrolytes can be classified as polyanions electrolytes,polycations electrolytes and polyampholyte.

Polyanion electrolytes ionize in the polar solvent,with numbers of anions on the macromolecular chains,and the corresponding small cations exist in the solution.By contrast,the macromolecular chains of polycation electrolytes possess amounts of cations and the small anions are in the solution.Both of polyanion and polycation electrolytes possess the properties of dispersion,flocculation,drag reduction and thickening etc.,which make them charming in tertiary oil recovery,sewage treatment,coatings,papermaking industries etc.With the prosperity and development of polyanion and polycation electrolytes,the study on them has attracted increasing attention[160,161,159].The inner electrostatic interactions in polyelectrolyte chain itself or electrostatic interaction between different polyelectrolyte chains makes it hold a richer structure.Miao et al.[162]study the electrostatic effect on the surface morphologies of a spherical PE brush in the presence of trivalent counterions under good,theta and poor solvent conditions based on a coarse-grained model by means of MD simulations.They found a nonmonotonic dependence of surface morphology on electrostatic strength,which represented clearly the electrostatic correlation effect mediated by the multivalent counterions.Their simulation results clearly demonstrated that ordered patterns could be induced by the electrostatic correlation effect in the presence of trivalent counterions,which is absent in the system with monovalent ions.Teraob [163]studied the properties of many-body interactions among charged dendrimers with the module of polyamidoamine (PAMAM) by stochastic MD simulations.It was found that the triplet force on charged dendrimers became repulsive,which showed the importance of the nonlinear feature due to the excluded volume effect incorporated with the screened Coulomb interaction.Ni et al.[164]explore the salt effect on polyanionic electrolytes sodium polyacrylate in salt-free and CaCl2solution by MD.The simulation results indicated that polyanionic electrolytes in aqueous solution had an extended chain structure due to the role of strong electrostatic repulsion.While with the introduction of Ca2+to the solution,the polyanionic electrolytes collapsed,and an explanation was offered for the collapse by analysis of the radius of gyration and the radial distribution function with atomistic resolution.

The study of the adsorption behaviors of polyelectrolytes at different interfaces are of great significance for expanding its applications in sewage treatment,biomedicine,and membrane separation.Quesada-Pérez et al.[165]used Monte Carlo simulations with a coarse-grained model to study the electrostatic interactions responsible for the trivalent-counterion-mediated adsorption of polyelectrolytes onto a like-charged planar surface.It was concluded that a delicate balance between charge inversion and screening effects governed the polyelectrolyte adsorption onto like-charged surfaces mediated by trivalent cations by the integrated analysis of the charge density and ionic distributions.It was also demonstrated that small ions favored the polyelectrolyte adsorption.Shiraki et al.[166]investigated the manner in which cationic polyethylene glycol (PEG) and polyelectrolyte chains(PEGylated polyelectrolyte) bound to anionic α-amylase by enzyme kinetic experiments and MD simulations,which was the first MD simulation study to predict the binding behavior between protein and PEGylated polyelectrolyte.In the work,the nonspecific interactions between PEG-b-PAMA and α-amylase were observed by MD simulation and it indicated that the random electrostatic interaction of PEGylated polyelectrolyte caused reversible inactivation of the α-amylase.Liu et al.[167]explored the interaction of anionic polyelectrolyte with cationic gemini surfactant by CGMD simulation.In fact,polyelectrolyte facilitates the oppositely charged ionic surfactants to aggregate by suppressing the electrostatic repulsion between ionic head groups and leading to the formation of micellar complex.They found that the amount of surfactant affected the conformation of polyion chain and more free micelles were formed by surfactants in mixtures at the same time.Increasing the length of spacer or tail chain in gemini surfactant will weaken its interaction with polyelectrolyte and simultaneously strengthen its tendency to self-assemble.It was revealed that the electrostatic interaction played an important role in the interaction of polyelectrolyte with gemini surfactant.Shi et al.[168]employed Monte Carlo simulations to investigate the interaction between an adsorbing linear flexible cationic polyelectrolyte and a binary fluid membrane.The result indicated that the cooperativity effect and the electrostatic interaction of the polyelectrolyte beads could significantly affect the segregation extent and the concentration gradients of the phosphatidylinositol molecules,and further cooperate to induce the complicated hierarchical mobility behaviors of PIP2 molecules.Yuan and Xu et al.[169]studied the adsorption of polyelectrolyte–surfactant mixture of sodium poly(acrylic acid) (NaPAA) and dodecyl trimethyl ammonium bromide(C12TAB)at the air/water interface by MD simulations,which indicated that the electrostatic interaction was the main driving force for the binding of surfactants to the polyelectrolyte.

With the continuous development of industrial production and biological nanomaterials etc.,the application of polyanion electrolytes and polycation electrolytes will be more extensive.Although there has been a lot of works done by scientists,the research on the polyelectrolytes must be constantly carried out,not only on experiments but also on simulations.As the electrostatic interaction affects the behaviors of polyelectrolytes,it should be explored carefully and thoroughly.

3.7.Electrostatic interactions in zwitterionic materials

Zwitterionic material is a kind of novel polyelectrolyte materials with both cationic and anionic groups in the same monomer residue,and has the characteristics of high dipole moment and high charge density [10].Carboxybetaine methacrylate (CBMA),sulfobetaine methacrylate (SBMA)and 2-methacryloxyethyl phosphorylcholine(MPC)are three commonly-used zwitterionic materials.They have been widely used in antifogging [11],antifouling[170],drug delivery [171],fuel cells [172]and underwater selfcleaning [173]due to its strong hydration capability via electrostatic interactions [12].

Previous studies have shown that the excellence performance of zwitterionic materials depends on the structure of zwitterionic moieties via MD simulation,including the charge density and the carbon spacer lengths (CSL) between cationic and anionic groups[174].Shao et al.[175]adopted MD simulations to investigate the hydration,self-association and protein interactions of a series of zwitterionic moieties with different charge densities,including three negatively groups(carboxylic(CO2),sulfonate(SO3),and sulfate (OSO3)) and four positively groups (quaternary ammonium(NC4),tertiary ammonium (NC3),secondary ammonium (NC2),and primary ammonium (NC1)).They found that all these zwitterionic moieties had strong hydration ability,but the structural and dynamic properties of water molecules near the charged groups were related to the types of positively and negatively groups.The coordination numbers of water molecules near the charged groups followed the order COO ＜OSO3＜SO3and NC1 ＜NC2 ＜NC3 ＜NC4;however,the order of the residence time of these water molecules is OSO3＜SO3＜COO when the NC4 group is the cationic group.Moreover,the self-association and protein interactions of these zwitterionic moieties also showed certain correlations with the charge density of charged groups.These results indicated that the charge density of charged groups played an important role in excellent performance of zwitterionic materials.It is known that the carbon spacer lengths between oppositely charged groups of zwitterionic molecules could influence the partial charges of the charged groups and further affect the hydration of zwitterionic moieties.Shao et al.[176]combined MD simulations with quantum chemical calculations to investigate the effect of carbon spacer lengths(CSL=0,1,2,3 and 4)on the molecular properties of zwitterionic carboxybetaine molecules.Simulation results indicated that the charged groups could form strong interplay between oppositely charged groups with shorter CSLs,leading to the weakness of the hydration of carboxybetaine molecules.For CSL ≥3,the hydration of carboxybetaine molecule strengthened with the increase of CSL and reached a steady state.Recently,Chen et al.[177]investigated the effect of CSL on the hydration and underwater oleophobicity of sulfobetaine-terminated self-assembled monolayers (SB-SAMs) by MD.They found that the hydration of SB-SAMs is positively dependent on CSL,the underwater oleophobicity is strengthened and then weakened with the increase of CSL,reaching optimal performance when CSL=3.These results provided a rich pool for understanding the roles of charged groups in determining the properties of zwitterionic materials.

Zwitterionic materials are highly resistant to nonspecific protein adsorption due to their strong hydration properties[178,179].He et al.[180]studied the interactions between lysozyme and phosphorylcholine-terminated self-assembled monolayers (PC-SAMs) to reveal the mechanism of resist nonspecific protein adsorption from the molecular level.The zwitterionic PCSAM surface could form a strong hydration layer to produce strong repulsive force when the lysozyme was close to it.Compared with OEG-SAM surface,PC-SAM surface was more tightly bound to water molecules.Liao et al.[181]found similar evidence that PC-SAM has lower underwater oil contact angle and oil adhesion than OEG-SAM.Zwitterionic monomers were also usually polymerized into polymer brushes and copolymers to modify different substrate materials.Liu et al.[182]investigated the surface hydration and surface resistance to a lysozyme of three zwitterionic polymer brushes (pCBMA,pSBMA and pMPC) by using MD and SMD simulations via CHARMM force field.For the three zwitterionic brushes and PEG brushes,they found that all brushes had strong hydration and surface resistance to lysozyme adsorption;the order of the hydration and repulsive force was pCBMA ＞pMPC ＞pSBMA ＞PEG.Moreover,the grafting density of zwitterionic brushes was a key factor affecting their structures and antifouling properties.At a higher grafting density,zwitterionic brushes exhibited a more organized structure,which could form a closely combined hydrated layer at the interface to resist pollutant adsorption;at a lower grafting density,the brushes exhibited a randomly oriented structure and the antifouling of brush was achieved via the deformation of the zwitterionic branches [183].A zwitterionic copolymer poly(ether sulfone)-blo ck-polycarboxybetaine methacrylate (PES-b-PCBMA) has been used as the blending additive to modify PES membrane by the nonsolvent induced phase separation(NIPS)method[184].During the exchange process of solvent/nonsolvent in NIPS process,the hydrophilic PCBMA segments of the copolymer migrate and enrich on the membrane surface via surface segregation,resulting in the improvement of hydrophilicity and antifouling properties.

Recently,some pH-responsive polymers containing zwitterionic polymer segments have been applied in drug delivery [185,186].The DPD simulation work by Liao et al.studied the selfassembled morphologies of poly(lactic acid)/poly(ethylene glycol)(PLA-PEG) and poly(lactic acid)/poly(carboxybetaine) (PLA-PCB)in aqueous solutions[186].It was found that the PLA-PCB system could form homogeneous shell layers while the PLA-PEG system had uneven shell layers,as shown in Fig.3(a).This result is due to the fact that the PEG shell layer is amphiphilic,while the PCB shell layer is strongly hydrophilic.Hao et al.[187]developed a novel pH-responsive triblock copolymer docosahexaenoic acid-bpoly(γ-benzyl-L-glutamate)-b-poly (carboxybetaine methacrylate)(DHA-PBLG-PCB)to control the release of doxorubicin by DPD simulations.They found that doxorubicin could be released effectively by adjusting the composition,concentration of DHA-PBLG-PCB copolymer and doxorubicin content.Under the physiological pH condition,the DHA-PBLG-PCB copolymer could self-assemble into a stable spherical micelle to protect doxorubicin core,while the micelle could be decomposed under acidic conditions and then doxorubicin further released.Zwitterionic monomers are also used as surfactants in a variety of consumer and industrial products,which have the advantages of hard water resistance,temperature stability,low irritation and good biodegradability.A series of methods to produce of zwitterionic surfactants have been developed,many of which contain positively charged quaternary ammonium groups and negatively charged groups,such as carboxylic,sulfonate and phosphate.Mafi et al.[188]developed an all-atom type force field by combining the TEAM force field with the GAFF force field to investigate the interfacial behavior of zwitterionic surfactant N-dodecyl-N,N-dimethyl-3-ammonio-1-propanesulfonate(DDAPS) on water surface.It was revealed that the head groups of DDPAS were nearly parallel to the water surface and the increase of DDPAS surface coverage had little influence on the orientation of the head groups.

In addition,positively and negatively charged moieties are evenly distributed through electrostatic interaction to form mixed-charged materials,which are equivalent to zwitterionic materials and have strong resistance to the adsorption of nonspecific proteins and oil droplet[189,190].Cheng et al.[190]investigated the underwater oleophobicity of nonionic hydrophilic selfassembled monolayers (SAMs) and mixed-charged SAMs by MD simulations.Simulation results indicated that mixed-charged SAMs had stronger hydration capacity and underwater oleophobicity than nonionic SAMs,the schematic diagram of underwater oleophobicity is shown in Fig.3(b).Underwater oleophobicity of mixed-charged SAMs could be modulated from superoleophilic to superoleophobic by adjusting the alkyl chain length and ionic strength [191].

Fig.3.The self-assembled structures of the PLA-PEG and PLA-PCB copolymers(a)[186];the underwater-oil contact angles of Nonionic SAM and zwitterionic SAM(b)[190].

3.8.Electrostatic interactions in nucleic acids

The study of biomolecules at the molecular scale is mainly accomplished by simulating the interactions between biomolecules.Understanding the role of electrostatic interactions is the central key point to gain deep knowledge of the folding and structural stability of protein,gene regulation function of nucleic,catalytic activity of enzymes and the transport of substances across biomembranes.

Nucleic acids (NAs),including deoxyribonucleic acid (DNA)[192]and ribonucleic acid (RNA) [193],play fundamental roles in biology.NAs are related to different gene regulation functions(i.e.,gene silencing and protein translation) as well as carriers of genetic information.Furthermore,NAs are biopolymers and also belong to polyelectrolytes due to their negatively charged backbones (each phosphate group on the backbone of a NA carries a negative charge) [194].To exert its biological function properly,NAs need to fold and pack into specific conformation,in which DNA folds into B-form double helix[192]and RNA exists as a single chain [193].It has been demonstrated that the neutralizing counterions around NAs,which is referred to as an ion atmosphere,have a close association with kinds of properties of NAs [14,15]and the electrostatic effect is very important for the stability structure,folding and function of highly charged NAs [195].

In the lack of counterions,the folding of a～400 nucleotide RNA will encounter about 600 kcal·mol-1(1 cal=4.1868 J)electrostatic repulsion[14,196]and the electrostatic repulsion could be reduced with the presence of counterions due to the screening effect.The folding process is hard to accomplish for most functional RNA molecules under low salt conditions and the Mg2+is usually used to induce folding of RNA.Inversely,electrostatic interactions are reduced at higher salt concentrations,which is unfavourable for the formation of DNA–protein and RNA–protein complexes [197].For experiments,it is difficult to provide a microscopic description for the ion-NA interactions because it is highly dynamic.As an alternative,computer MD simulations were effective way to reveal the interaction nature of ions,water and surface with NAs[198,14,199,200],in which not only the electrostatic interaction(i.e.,the ‘‘passive”screening effects) could be well described,but also the ‘‘a(chǎn)ctive”effects of ions on NAs dynamics at the atomic details could be obtained [201].Fischer et al.[200]investigated the influence of Na+and Mg2+ions on RNA structures by MD simulations,they found that Mg2+conserved the experimental structure of RNA better than Na+and more Mg2+ions were closer to RNA than Na+.Li et al.[201]studied the effect of Mg2+and Na+ions on structure and conformational dynamics of DNA,they found that the binding of Mg2+ions made the DNA duplex more rigid and the overall conformational entropy was reduced significantly for Mg2+-bound DNA compared with that of Na+-bound DNA;in addition,the local bending amplitudes were decreased with the addition of Mg2+ions.

Furthermore,the adsorption of NAs on different surfaces has been widely researched in recent years.Muraru et al.[202]investigated the adsorption behavior of DNA on graphene oxide (GO)and reduced graphene oxide functionalized with aminated polyethylene glycol (rGO-PEG-NH2) surface in the presence of Mg2+ions.They demonstrated that the adsorption process of DNA on the two surfaces were promoted by the electrostatic force and hydrogen bonds.It is worth pointing out that DNA adsorbed on GO was mediated by Mg2+ions,where Mg2+ions acted as bridges between the GO surface and DNA through electrostatic interactions.Similarly,Martín-Molina et al.[203]researched the adsorption mechanism of DNA on anionic lipid surfaces by MD simulations,they revealed that the adsorption of DNA onto anionic lipid monolayers was mediated by Ca2+and the optimal conditions for cation-mediated adsorption of DNA onto negatively charged surfaces were established.

3.9.Electrostatic interactions in proteins

Protein is the expression vector of life information.The folding and stability mechanisms of proteins are key to comprehend the structure and function of proteins.Proteins are biological macromolecules formed by the polymerization of amino acids,in which charged and hydrophobic residues impart important properties to proteins.It is well acknowledged that hydrophobic interactions between nonpolar residues are the driving force for protein folding,while the electrostatic interactions between charged residues can significantly modulate the folding stability [204,205].Duan and Kollman successfully simulated the folding process about 1 μs of a small protein for the first time,sparking an upsurge in protein folding research [206].

The charged amino acids contained in proteins have different charges under different pH environments.The protonation or deprotonation of ionizable groups leads to large repulsion between similar charges,which leads to changes in their structure.Azia et al.[205]employed CGMD simulations to understand the effect of charge–charge interactions on protein folding involved in three protein systems.They found that the charge mutation mainly affected the unfolded-state,while the folded state remained similar to the wild type for both N-terminal L9 protein domain and RNase Sa protein.In addition,the structurally homologues of βlactalbumin and hen egg-white lysozyme exhibited different folding kinetics due to the different distribution of charged residues,which demonstrates that electrostatic interactions can significantly modulate protein folding kinetics.Intrinsically disordered proteins (IDPs) are functional proteins without a specific tertiary structure,which fold into a specific conformation that could impose a kinetic bottleneck,but they can achieve efficient folding upon binding of specific targets [207].Ganguly and co-workers[208]combined experiments and CGMD simulations to investigate how IDPs bind with target and fold into a specific conformation.The results demonstrated that electrostatic interactions played an extremely important role in this process,it not only accelerated the encounter rate but also promoted folding-competent topologies in the encounter complexes by modulating the early stages of coupled binding and folding.The effects of charge mutation,pH and salts on the protein folding,unfolding and assembly structures have been extensively studied.Srivastav et al.[209]gave an insight into the structural stability of bovine serum albumin at five different ionic strengths and pH via MD simulations.Their simulation results indicated that the protein structure was highly charged at both low or high pH and high ionic strength conditions;the strong intramolecular electrostatic repulsion would induce denaturation of the protein structure.

The fully understanding of protein adsorption behavior on charged surfaces is of fundamental importance for the development of applications including biocatalysis,biosensors,biofuel cells and biomaterials [210–214],since the adsorption orientation and conformation of adsorbed proteins greatly impacts their biological functions.Zhou and his co-workers have done systematic researches to reveal how interface charge density (SCD),interface topology,solution ionic strength (IS),pH value and external electric field affect the adsorption behavior of different proteins[210–229].They found that van der Waals and electrostatic interactions co-determined protein adsorption on surfaces,while the overall electric dipole of a protein(the electric dipole definition of a protein is miu=where qiis the local charge of each atom,and rCOMis the center of mass position of the protein)played a dominant role in controlling the adsorption orientation of the protein on charged surfaces [215–218].Consequently,the adsorption of proteins on charged surfaces with a highly ordered orientation can be regulated based on the nature of the protein’s electric dipole and the charged characteristics of the surfaces.

Fig.4.The final snapshots of α-ChT adsorbed on pristine CNTs (a),low charged CNTs (b),and high charged CNTs (c),respectively [221].

Liu et al.[219,220]investigate the adsorption orientation of prototypical and mutated protein G B1 on positively and negatively charged SAMs via the combined parallel tempering Monte Carlo(PTMC) and all-atom molecular dynamics (AAMD) simulation approach.The simulation results showed that electrostatic interactions dominated the protein adsorption orientation of both kinds of protein,the prototype protein absorbed on a negatively charged SAMs surface was beneficial to improve the efficiency of antibody binding.They also found that the total net charge of protein was not the main factor for determining its adsorption behavior;the distribution of charged residues on the protein surface had a significant impact on the adsorption behavior.Zhao et al.[221]explored the interaction mechanism between carboxylated CNTs and a-chymotrypsin (α-ChT) at the molecular level.It was found that α-ChT interacted with pristine CNTs through vdW interaction,and the active site was towards the solution,which was conducive to its catalytic function;while the main driving force was electrostatic interaction when it bound onto carboxylated CNTs with the active site facing toward the CNTs surface,which caused enzymatic activity inhibition.The final simulation snapshots and the key adsorbed are shown in Fig.4.

For electrostatically dominated protein adsorption processes,it strongly depends on solution conditions,such as solution pH and IS[230–232].In general,higher IS will shield the electrostatic interaction between the protein and the surface,making the adsorption strength of the system with the opposite charge weakened,while adsorption can occur in systems that were mutually exclusive.When both SCD and IS are high,a counterion layer may be formed near the surface,thereby mediating the adsorption of proteins with the same charge [218].Liu et al.[222]investigated the effects of SCD and IS on the adsorption behavior of feruloyl esterase from Aspergillus niger (AnFaeA) on the charged SAMs surface.It was found that the electrostatic interactions between AnFaeA and the SAMs surface became weakened as the IS increased.The counterion layer formed near the surface played a key role in the adsorption of AnFaeA on the positively NH2-SAM surface.The adsorption loading of proteins or the final adsorption state will also be affected by the different types of ions in solution[233].Peng et al.[223]investigated the influence of anion type (chlorine and phosphate ions)and IS on the adsorption behavior of Cyt-c on positively charged SAM surfaces.It was shown that positively-charged Cyt-c could not be stably adsorbed on NH2-SAM surface regardless of chloride ion under high or low IS as chloride ions did not have obvious tendency to aggregate near the surface of NH2-SAM,which could not shield the electrostatic repulsion between Cyt-c and the NH2-SAM surface.In contrast,phosphate ions could form a distinct counterion layer near the surface at both low and high IS,thus promoted the adsorption of Cyt-c.

The environmental pH value determines the surface charged state of the protein,the protein will carry net negative charge when the pH value is higher than its isoelectric point;while it will carry net positive charge when the pH value is below the isoelectric point[224].In addition,changes in the pH of the solution cause protonation or deprotonation of some charged surfaces,and therefore affect the interaction between the protein and the surface.Johonson et al.[234]pointed out that the adsorption orientation of protein G can be adjusted by changing the pH value of the solution.By CGMD simulations,Yu et al.showed that lysozyme adsorbed onto the hydrophobic charged induction chromatography (HCIC) surfaces was driven by the hydrophobic interaction when solution pH=7 [224].In contrast,the lysozyme would desorb from the HCIC surfaces due to the electrostatic repulsion interaction between lysozyme and the surface when solution pH=4.Liu et al.[189]conducted a simulation work on the interaction mechanism between human gamma fibrinogen and mixed charged SAM(which could be regulated by solution pH environment)before and after hydrolysis.The results showed that negative charged human gamma fibrinogen proteins could stably adsorb on the positively charged surface through electrostatic interactions before hydrolysis.After hydrolysis,however,the surface charge changed from a positively charged state to a pseudo-zwitterionic state,which led to the protein detached from the surface due to the disappear of electrostatic interaction between the protein and pseudozwitterionic SAM surface.

In addition to controlling protein behavior through factors such as protein mutation,designing new surface materials,or adjusting the microenvironment of solutions,the regulation of protein adsorption could also be controlled by the external electric field.Xie et al.[225]studied the effect of electric field on the adsorption of lysozyme on the COOH-SAM surface.The simulation results showed that a positive electric field could promote the adsorption of lysozyme,while a negative electric field would weaken or even resist its adsorption.However,the promotion or inhibition of protein adsorption by the electric field did not change monotonically with the strength of the electric field.They believed that the competitive adsorption of counterion ions and lysozyme on negatively charged surfaces had significant effect on the process of protein adsorption.Later,they also investigated the effect of external electric field on the adsorption of Cyt-c on PC-SAM surface [226].Because the PC-SAM surface contained both positively charged choline groups and negatively charged phosphate groups,changing the size and direction of the electric field can make the surface appear three different charged states,thereby achieving the purpose of regulating protein adsorption.The adsorption behavior of lysozyme on electrically responsive mixed COOH/OH-SAM surface was also explored by MD[235].It was found that the positive and negative electric fields affect the adsorption behavior of lysozyme by regulating the hydrophilicity and hydrophobicity of the mixed COOH/OH-SAM surface.

Metalloproteins represent a large portion of proteins,it has been estimated that about one-third of all proteins are associated with a metal,i.e.,Mg,Zn,Ca,Mn,Fe,Pb,Co,K,Na and Ni [236].The metal ions play a significant role in biological activity of the proteins.Recently,MD simulations were widely applied to study metalloproteins,despite the challenge introduced by the presence of metal ions.Zhao et al.[237]investigated the adsorption behavior of Cyt-c on graphene and GO surfaces,they found that Cyt-c adsorbed on graphene surface through hydrophobic interaction,while the adsorption of Cyt-c on the GO surface was dominated by electrostatic forces.Furthermore,it had also been revealed that the adsorption orientation of Cyt-c on GO surface facilitated the electron transfer (ET) with the heme plane slightly deviated from the normal direction to the surface and the axial ligand Met80 close to the surface,it was not conducive to ET when Cyt-c adsorbed on graphene surface with heme plane horizontally oriented and far away from the graphene surface.Yang et al.[238]employed Monte Carlo and all-atom MD simulation methods to studied the adsorption mechanism of the heme-contain myoglobin(Mb) on rutile (1 1 0) and (0 0 1) surfaces.The simulation results indicated that Mb adsorbed on the rutile (0 0 1) with the heme and the electron transfer pathway of Mb close to the surface,which was favorable for faster ET.However,the heme was far away from surface when Mb was adsorbed on the rutile (1 1 0) surface because of electrostatic interactions between the heme and the outmost atoms of rutile.They concluded that the rutile(0 0 1)surface was more suitable for Mb immobilization.The adsorption behavior of the standard Ni-Fe hydrogenase on aminoterminated alkanethiol SAMs with different levels of protonation was revealed through MD simulation combined with surfaceenhanced infrared absorption (SEIRA) measurements by Utesch et al.[239].The results showed that hydrogenase presented a strong adsorption on SAMs with a high degree of protonation;however,on lowly ionized SAM surface,a weak immobilization was observed.It indicated that SCD was a crucial parameter for the electrostatic binding of the target enzyme on charged surfaces.

Enzymes,as an efficient biocatalyst,have been widely applied in food,medicine,energy and environmental protection.The immobilization of enzymes can not only significantly improve their catalytic activity and stability,but also can facilitate recovery and recycling,overwhelmingly reduce the cost of biocatalytic processes [240,241].In enzyme fuel cell,biosensors or biochips,the orderly immobilization of an enzyme on an electrode can effectively promote the direct electron transfer of the electrode,which is beneficial to the direct conversion of chemical energy and electrical energy [242–245].Liu et al.[246]studied the direct electron transfer(DET)efficiency of laccase(TvL)on charged SAMs surfaces.The simulation results showed that compared with the negatively charged surface,the positively charged surface could regulate the adsorption orientation of T1 copper site of TvL closer to the surface,which made the path of the electron transfer shorter,thereby achieving higher direct electron transfer efficiency.A similar study was conducted by Yang et al.[247]to explore the electron transfer behavior between bilirubin oxidase(MvBOx)and SAM surface with different SCDs.It was indicated that the COOH-SAM with lower SCD was more suitable for the quick DET process of MvBOx immobilized on the electrode surface.Recently,they performed a research work on a novel TcAChE-based amperometric biosensor for the detection of organophosphorus pesticides.They investigated the adsorption details of Torpedo californica (TcAChE) on the charged SAMs with different SCDs [248].It was found that TcAChE could stably adsorb on both COOH-SAM and NH2-SAM surface,the driving force for adsorption was electrostatic interactions.Compared with adsorption on negatively charged COOH-SAM surface,TcAChE adsorbed on positively charged NH2-SAM surface could gain an efficient bio-catalytic reaction and quick DET for the reason that the active-site gorge and active center was close to the surface,which shortened the path for electron and substrate transfer.For enzyme catalysis,enzymes adsorbed on surfaces in a specific orientation and with the active part of the enzyme toward the reaction system could effectively improve the enzyme utilization and catalytic efficiency[249].Zhao et al.[250]investigated the adsorption behavior of lipase at four different surfaces (i.e.,graphite,TiO2,NH2-SAM,COOH-SAM).It was shown that lipase was stably adsorbed on the graphite surface through π-π stacking and hydrophobic interaction,with the active center towards the solvent,which was beneficial to the binding of the enzyme and substrate.Due to the blocking effect of the two hydration layers adsorbed at TiO2surface,the interaction between lipase and the TiO2surface was weak.Additionally,the interaction between lipase and the NH2-SAM surface was weak and unstable,it underwent adsorption and desorption processes.Although lipase carries a negative charge,it could be strongly anchored to the COOH-SAM surface via a single Arg249 residue through electrostatic interaction.The catalytic activity center of the lipase eventually faced to the solution,which was beneficial to improve the catalytic activity of lipase.Liu et al.[251]investigated the interaction mechanism of CalB at the GO surface with different oxygen contents.It was indicated that the electrostatic force dominated the adsorption,and it was enhanced with the increase of oxygen content of GO surface.

3.10.Electrostatic interactions in biomembranes

Phospholipids,which consist of charged hydrophilic headgroups and long hydrophobic alkyl chains,are main components of biomembranes.It has been demonstrated that the charge properties of phospholipids can significantly influence the structure of biomembranes,as well as the dynamic entrance process of some functional nanomaterials into the interior of cells.Therefore,understanding the electrostatic interaction between phospholipids and nanomaterials is particularly crucial in the studies of nanomaterial-biomembrane interactions.

In recent years,increasing number of studies have been devoted to investigate the interactions between nanomaterials and biomembranes through computer simulations.Heikkila et al.[252]simulated the adsorption behaviors of cationic and AuNPs on lipid bilayers by AAMD simulations,they indirectly inferred the membrane permeation mechanism of AuNPs.However,due to the time limitation of AAMD,the specific transmembrane process is difficult to be observed by AAMD.Therefore,large-scale mesoscopic CGMD simulations have been performed to study nanomaterial-biomembrane interactions,which mainly include DPD simulations and the MARTINI force filed-based CG simulations.

Through the introduction of charge smearing on DPD particles,the interactions between charged NPs and lipid bilayers have been widely studied by DPD simulations[28,253–256].Ding et al.[253]systematically investigated the mechanism of cellular uptake of cationic polymeric nano-vectors containing DNA molecules through DPD simulations.The results showed that the property of polyelectrolyte chains grafted to nano-vector and DNA molecules had important impacts on the endocytosis.Li et al.[255]explored the mechanism of the cooperative endocytosis between like-charged NPs;they demonstrated that it was definitely caused by the interplay of particle size,the charged ligand density on particle surface and local concentration of NPs.Diffusion dynamics of charged NPs on the lipid membrane is of essential importance to cellular functioning.Chen et al.[256]presented a computational investigation to uncover the pivotal role of electrostatics in the diffusion dynamics of charged NPs on the lipid membrane,the results demonstrated that the diffusive behavior and directional transport of a charged NP significantly depended on the sign and spatial distribution of charges on its surface.Wang et al.proposed an effective ENUF–DPD method to compute the electrostatic interaction based on Ewald summation method,which was successfully used in simulating the interactions between PAMAM dendrimers and DMPC lipid bilayers [28].

Based on the MARTINI force field,Lee et al.explored the effect of electrostatic interaction treatments on the properties of pure DPPC lipid bilayers and the interplay between PAMAM dendrimers and lipid bilayers [257,258].It has been demonstrated that the area per lipid in the bilayer obtained with PME summation compared well with that obtained with the truncation method.However,for an experimental phenomenon in which dendrimers was inserted into lipid bilayers and induced pore formation,only the simulation with the PME method could reappear.Tian et al.[259]explored the interactions between PAMAM dendrimers and lipid bilayers under different pH conditions,they revealed the importance of electrostatic interaction on the transmembrane delivery of dendrimers.Hu et al.[258]investigated the effect of hydrophobicity and surface charge of NPs on the membrane translocation process,they demonstrated that charged NPs generally translocated quickly across the pulmonary surfactant film,while hydrophobic NPs were trapped by the surfactant film and encapsulated in lipid protrusions upon film compression.

There are also several CGMD simulation works about the interactions between AuNPs and biomembranes.Based on the structure and dynamical properties of AuNPs in experiments,Lin et al.[260]established a CG model of monolayer-protected AuNPs under the framework of the MARTINI force field.Furthermore,they investigated the interactions between lipid membranes and monolayerprotected AuNPs with different surface charge properties.Through the PME method,they successfully observed the penetration process of AuNPs into membranes.The CGMD simulation work conducted by Quan et al.[261]systematically explored the interactions between AuNPs with different surface chemistries and lipid membranes with different composition distribution(Fig.5).Simulation results showed that AuNPs with different signs of charge could spontaneously adhere to the membrane surface or penetrate into the membrane core.Besides,the asymmetric distribution of charged lipids in membranes could facilitate the penetration of cationic AuNPs.Increasing the surface charge density of AuNPs could not only improve the penetration efficiency,but also lead to more severe disruption of membrane structure.This work revealed the driving force for AuNPs penetrating into lipid membranes and built a connection between transmembrane efficiency and cytotoxicity of AuNPs.Besides,they also studied the transmembrane behaviors of pH-responsive zwitterionic AuNPs under different dissociation degrees by CGMD simulations [262].Simonelli et al.[263]investigated the membrane penetration of anionic AuNPs,they found that the interactions between AuNPs and phospholipid membranes could be divided into three stages:the adsorption of AuNP onto the membrane surface;the hydrophobic contact formed between hydrophobic ligands of AuNPs and protruding solvent-exposed lipid tails;the embedding of AuNP into the membrane core with a snorkeling configuration.Gupta et al.[264]systematically investigated effects of particle size and surface charge on the transdermal delivery of AuNPs.Shimizu et al.[265]simulated the membrane penetration behaviors of cationic AuNPs under an applied electric field,they found that AuNPs could fully cross the lipid bilayer through a membrane pore induced by the electric field.

Fig.5.Representative snapshots of equilibrated states for AuNPs with different signs of charge interacting with the three kinds of lipid membranes [261].

In the overwhelming majority of cases,most attention has been focused on the interactions between NPs and single-phase,homogeneous lipid bilayers.However,cell membranes,especially plasma membranes,are multicomponent and structurally heterogeneous membranes which are separated into many nanosized domains such as lipid rafts.NPs may show different interactions with different domains and influence the lateral organization of biomembranes.Recent simulation works have gradually been focused on this point.Chen et al.[266]investigated the interactions between NPs coated with neutral and/or charged ligands and phase-separated lipid bilayers using CGMD simulations.They demonstrated that both penetration and adsorption processes as well as the final distribution of the NPs could be readily modulated by varying the ligand density and the surface charge of nanoparticles.Noh et al.[267]explored the structural properties of mixed phospholipid bilayers in the presence of ligand-functionalized NPs.Simulation results showed that charged NPs aggregated near unsaturated phospholipid regions,the phase-separation of phospholipid bilayer could be promoted in the presence of NPs,the heterogeneous components of a phospholipid bilayer played a significant role in the lateral organization of NPs.

4.Conclusions and Perspectives

In the fields of chemistry,chemical engineering,material science,biochemistry and biomolecular science and engineering etc.,the long-range electrostatic interaction,plays a dominant role in determining the properties of charged complex fluid systems including inorganic electrolytes,organic electrolytes,polyelectrolytes and bio-electrolytes.It is interesting and challenging to understand the relationships between the macroscopic properties of complex charged systems with the microscopic structures induced by electrostatic interactions.

(A) By modulating the electrostatic interactions between ions and the charged groups of nanomaterials,the structure of the ion hydration layer can be changed,and the behavior of water molecules and ions can be indirectly controlled to achieve water desalination and ion separation.The piezoelectric properties of piezoelectric materials depend strongly on the shape,size,temperature,surface property and electrostatic interaction.Moreover,free charge in the piezoelectric materials can be tuned by an externally applied electric field and thus achieve electromechanical coupling for nanoscale actuation,sensing,detection and energy harvesting.

(B) The electrostatic interactions mediated by ions in the solution environment make the organic electrolyte exhibit different behaviors at different interfaces,thus exerting multiple functions in biomineralization,battery electrolyte,buffering agent and eluent.In addition,electrostatic interactions regulate the adsorption and micro-aggregation behaviors of ionic surfactants at the gas–liquid,oil–water,and solid–liquid interfaces,making it widely used in cosmetics,mining flotation,material manufacture,enhanced oil recovery and pollutant removal.

(C) The change in the interplay between ions and the polyelectrolyte mediated by the counterion and the ionic strength of the solution makes the hydration and micro-aggregation behavior of the polyelectrolyte an incredible change.Understanding the structure–property relationships of zwitterionic materials,such as the effects of ion species,hydrophilicity and hydrophobicity,as well as main chain structure on material properties,is useful for expanding its applications in anti-fouling coatings,protein modification,drug delivery,membrane separation etc.

(D)Electrostatic interactions also make a pivotal effect in studying the folding pathway and structural stability of biomolecules.By electrostatic interactions,the specific adsorption orientation of enzymes on surfaces could be regulated,which is useful for the application of enzymes in biodiesel,supercapacitors,fuel cells,electrode materials and etc.Understanding the electrostatic interaction between phospholipids and nanomaterials is particularly crucial in the studies of nanomaterial-biomembrane complex system.

It needs to be emphasized that the regulation and prediction of the above-mentioned physical and biochemical characteristics of complex fluids will only be realized once we are able to establish reasonable,accurate and effective electrostatic interaction treatment methods.However,there are always irreconcilable contradiction between the calculation cost and accuracy of electrostatic interactions.After decades of development,despite significant progress have been achieved in the qualitative understanding and quantitative evaluation of electrostatic interactions,several aspects need further attention.

(i) To make full use of the existing and emerging computer architectures such as graphic processing units(GPUs);to construct new,faster as well as more efficient and accurate algorithms for the computation of electrostatics.

(ii) To master the qualitative understanding and quantitative evaluation of electrostatic interactions;to further gain a deep understanding the role of electrostatics in the realization of the quantitative predictive physical/biochemical characteristics of complex fluids,and thus to better serve related chemical process and product development.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (21776093,21376089,41976203,21506178,21908066).