| 86 |
|
\centering |
| 87 |
|
\includegraphics[width=0.50 \textwidth]{chargesystem.pdf} |
| 88 |
|
\caption{NaCl crystal showing (a) breaking of the charge ordering in the direct spherical truncation, and (b) complete $(NaCl)_{4}$ molecule interacting with the central ion. } |
| 89 |
< |
\label{figure1} |
| 89 |
> |
\label{fig:NaCl} |
| 90 |
|
\end{figure} |
| 91 |
|
|
| 92 |
< |
Any charge in a NaCl crystal is surrounded by opposite charges. Similarly for each pair of charges, there is an opposite pair of charge to its adjacent as shown in Fig 1a. Furthermore for each group of four charges, there should be an oppositely aligned group of four charges as shown in Fig 1b. If we consider any group of charges, suppose $(NaCl)_4$, far from the central charge, they have little electrostatic interaction with the central charge (acts like point octopole when it is far from the center ). But if the cutoff sphere passes through the $(NaCl)_4$ molecule leaving behind net positive or negative charge, the electrostatic contribution due to these broken charges going to be very large (for point charge radial function $1/r_c$ and for point octupole $1/r_c$). Because of this reason, although the nature of electrostatic interaction short ranged, the hard cutoff sphere creates very large fluctuation in the electrostatic energy for the perfect crystal. In addition, the charge neutralized potential proposed by Wolf et al. converged to correct Madelung constant but still holds oscillation in the energy about correct Madelung energy.\cite{Wolf99}. This oscillation in the energy around its fully converged value should be due to the non-neutralized value of the dipole and higher order moments within the cutoff sphere. Recently, \textit{Ikuo Fukuda} used neutralization of the higher order moments for the calculation of the electrostatic interaction of the point charges system.\cite{Fukuda13} |
| 92 |
> |
Any charge in a NaCl crystal is surrounded by opposite charges. Similarly for each pair of charges, there is an opposite pair of charge to its adjacent as shown in Figure ~\ref{fig:NaCl}. Furthermore for each group of four charges, there should be an oppositely aligned group of four charges as shown in Fig 1b. If we consider any group of charges, suppose $(NaCl)_4$, far from the central charge, they have little electrostatic interaction with the central charge (acts like point octopole when it is far from the center ). But if the cutoff sphere passes through the $(NaCl)_4$ molecule leaving behind net positive or negative charge, the electrostatic contribution due to these broken charges going to be very large (for point charge radial function $1/r_c$ and for point octupole $1/r_c$). Because of this reason, although the nature of electrostatic interaction short ranged, the hard cutoff sphere creates very large fluctuation in the electrostatic energy for the perfect crystal. In addition, the charge neutralized potential proposed by Wolf et al. converged to correct Madelung constant but still holds oscillation in the energy about correct Madelung energy.\cite{Wolf99}. This oscillation in the energy around its fully converged value should be due to the non-neutralized value of the dipole and higher order moments within the cutoff sphere. Recently, \textit{Ikuo Fukuda} used neutralization of the higher order moments for the calculation of the electrostatic interaction of the point charges system.\cite{Fukuda13} |
| 93 |
|
|
| 94 |
|
\subsection{Disordered system} |
| 95 |
< |
The $r ^{-5}$ convergence behaviors is not only limited to the perfect crystals but also applied in the highly disordered crystal.\cite{Wolf99} At high temperature there should be local ordering of the charge and higher multipole moments in the liquids (To form the structure which is electrostaticaly neutral) but this ordering disappears at the long range. As in ionic crystal, even for liquid positive ion tends to be surrounded by the negative ion and vice versa, so the spherical truncation breaks the short range charge ordering present in the liquid system which results in oscillation (smaller amplitude in electrostatic energy of liquid as compared to crystal).\cite{Wolf99} This idea can also be generalized in molecule with multipole moments assuming local ordering is even true for multipoles, which is supported by the presence of the oscillation of the electrostatic energy as it plotted against the cutoff radius for dipolar liquid 2a. For quadrupolar liquid oscillation damped pretty quickly as seen in Fig 2c because of short range nature of the quadrupole-quadrupole interaction. |
| 95 |
> |
The $r ^{-5}$ convergence behaviors is not only limited to the perfect crystals but also applied in the highly disordered crystal.\cite{Wolf99} At high temperature there should be local ordering of the charge and higher multipole moments in the liquids (To form the structure which is electrostaticaly neutral) but this ordering disappears at the long range. As in ionic crystal, even for liquid positive ion tends to be surrounded by the negative ion and vice versa, so the spherical truncation breaks the short range charge ordering present in the liquid system which results in oscillation (smaller amplitude in electrostatic energy of liquid as compared to crystal).\cite{Wolf99} This idea can also be generalized in molecule with multipole moments assuming local ordering is even true for multipoles, which is supported by the presence of the oscillation of the electrostatic energy as it plotted against the cutoff radius for dipolar liquid in Figure ~\ref{fig:rcutConvergence}. For quadrupolar liquid oscillation damped pretty quickly as seen in Figure ~\ref{fig:rcutConvergence_hardQuadrupole} because of short range nature of the quadrupole-quadrupole interaction. |
| 96 |
|
\begin{figure} |
| 97 |
|
\centering |
| 98 |
|
\includegraphics[width=0.45 \textwidth]{rcutConvergence_hard_dipolar.pdf} |
| 100 |
|
\begin{figure} |
| 101 |
|
\centering |
| 102 |
|
\includegraphics[width=0.45 \textwidth]{rcutConvergence_hard_quadrupole.pdf} |
| 103 |
< |
\caption{Total energy per molecule against cutoff radius, rc for (a) dipolar liquid (b) dipolar crystal, (c) quadrupolar liquid and (b) quadrupolar crystal to compare the oscillation in the electrostatic energy. The crystalline system shows larger oscillation as compared to liquid. Also the fluctuation in the dipolar system is very large as compared to quadupolar system. } |
| 104 |
< |
\label{figure1} |
| 103 |
> |
|
| 104 |
> |
\caption{Total energy per molecule against cutoff radius, $r_c$ for (a) dipolar liquid (b) dipolar crystal, (c) quadrupolar liquid and (b) quadrupolar crystal to compare the oscillation in the electrostatic energy. The crystalline system shows larger oscillation as compared to liquid. Also the fluctuation in the dipolar system is very large as compared to quadupolar system. } |
| 105 |
> |
\label{fig:rcutConvergence} |
| 106 |
|
\end{figure} |
| 107 |
|
\subsection{Oscillation in the electrostatic energy} |
| 108 |
|
The oscillation of the electrostatic potential energy per molecule for a direct truncation method is associated with the charge neutrality.\cite{Wolf99} The electrostatic energy of a central molecule due to all other molecules within cutoff sphere is plotted against the cutoff radius for (i)dipolar liquid (i) perfect dipolar crystal, (iii) quadrupolar liquid, and (iv)quadrupolar cyrstal in Figure 2a, 2b, 2c and 2d. The larger amplitude in the oscillation of electrostatic energy for the perfect crystal is because of the long range of multipolar ordering.\cite{Wolf99} Moreover, the oscillation damped much faster for the system of higher order multipoles (Compare range of oscillation between the dipolar and quadrupolar system in figure 2a, 2b, 2c and 2d). As in the case of the charge system, this oscillating nature of the electrostatic energy of the central molecule should be due to the net charge-multipole within the cut off sphere. If the amplitude of the oscillation is very large such as in ionic\cite{Wolf99} and dipolar crystal (Figure 3) then it will introduce huge error in the calculation of the absolute energy. On the other hand, if the oscillation is small, it can contribute error in the evaluation of the energy differences between configurations. |
| 114 |
|
Any force field associated with MD simulation should address two major issues in the electrostatic interaction. First, it should deal with the breaking of the charge or multipole ordering due to direct spherical truncation. Second, the electrostatic energy, force and the torque between central molecule and any other molecule should smoothly approaches to zero as $r$ tends to $r_c$. The first issue is associated with the oscillation of the total electrostatic potential energy of the central molecule due to all other molecules within cutoff sphere and second issue is related with the continuous nature of the electrostatic interaction at the cutoff radius, which eventually related with the conservation of total energy in the MD simulation. The mathematical detail for the SP, GSF and TSF has already been discussed in detail in previous paper I.\cite{PaperI} |
| 115 |
|
|
| 116 |
|
\subsection{Shifted potential (SP) } |
| 117 |
< |
As it shown in Fig 2, a discontinuous truncation of the electrostatic potential at the cutoff sphere introduces severe artifact(Oscillation in the electrostatic energy) even for molecules with the higher-order multipoles. This artifact is due to the existence of higher order moments within the cutoff spheres.The net multipole moment within cutoff sphere is contributed by the breaking of the multipole ordering in direct truncation of the cutoff sphere. The multipole moments of the cutoff sphere can be neutralized by placing image multipole on the surface, for every each multipole within it. The electrostatic potential between multipoles for the SP method is given by, |
| 117 |
> |
As it shown in Figure ~\ref{fig:rcutConvergence}, a discontinuous truncation of the electrostatic potential at the cutoff sphere introduces severe artifact(Oscillation in the electrostatic energy) even for molecules with the higher-order multipoles. This artifact is due to the existence of higher order moments within the cutoff spheres.The net multipole moment within cutoff sphere is contributed by the breaking of the multipole ordering in direct truncation of the cutoff sphere. The multipole moments of the cutoff sphere can be neutralized by placing image multipole on the surface, for every each multipole within it. The electrostatic potential between multipoles for the SP method is given by, |
| 118 |
|
\begin{equation} |
| 119 |
< |
U_{SF}(\vec r)=U(\vec r) - U(\vec r_c) |
| 119 |
> |
U_{SP}(\vec r)=U(\vec r) - U(\vec r_c) |
| 120 |
> |
\label{eq:SP} |
| 121 |
|
\end{equation} |
| 122 |
|
The SP method compensates the artifact created by truncation of the multipole ordering by placing image on the cutoff surface. Also, the potential energy between central multipole and other multipole within sphere approaches smoothly to zero as $r$ tends to $r_c$. But the force and torque obtained from the shifted potential are discontinuous at $r_c$. Therefore, the MD simulation will still have the total energy drift for a longer simulation. If we derive the force and torque from the direct shifting about $r_c$ like in shifted potential then inconsistency between the force, torque, and potential fails the energy conservation in the dynamic simulation. |
| 123 |
|
\subsection{Taylor-shifted force(TSF)} |
| 129 |
|
(C_b - D_{b \alpha }\frac{\partial}{\partial r_{b \alpha}}+Q_{b \alpha \beta }\frac{\partial}{\partial r_{b \alpha}\partial r_{b \beta}})\\ |
| 130 |
|
[(\frac{1}{r}-[\frac{1}{r_c}-(r-r_c)\frac{1}{r_c^2}+(r-r_c)^2\frac{1}{r_c^3}+...)] |
| 131 |
|
\end{split} |
| 132 |
+ |
\label{eq:TSF} |
| 133 |
|
\end{equation} |
| 134 |
|
|
| 135 |
< |
where $C$ , $D_{\alpha,\beta}$, and $Q_{\alpha,\beta}$ stands for charge, component of the dipole and quadrupole moment respectively (detail in paperI\cite{PaperI}). The electrostatic force and torque acting on the central molecule due to a molecule within cutoff sphere are derived from the equation (3) with the account of appropriate number of terms. This method is developed on the basis of using kernel potential due to the point charge ($1/r$) and their image charge potential ($1/r_c$) with its Taylor series expansion and considering that the expression for multipole-multipole interaction can be obtained operating the modified kernel by their corresponding operators. |
| 135 |
> |
where $C$ , $D_{\alpha,\beta}$, and $Q_{\alpha,\beta}$ stands for charge, component of the dipole and quadrupole moment respectively (detail in paperI\cite{PaperI}). The electrostatic force and torque acting on the central molecule due to a molecule within cutoff sphere are derived from the equation ~\ref{eq:TSF} with the account of appropriate number of terms. This method is developed on the basis of using kernel potential due to the point charge ($1/r$) and their image charge potential ($1/r_c$) with its Taylor series expansion and considering that the expression for multipole-multipole interaction can be obtained operating the modified kernel by their corresponding operators. |
| 136 |
|
\subsection{Gradient-shifted force (GSF)} |
| 137 |
|
As we mentioned earlier, in the MD simulation the electrostatic energy, force and torque should approach to zero as r tends to $r_c$. Also, the energy, force and torque should be consistent with each other for the total energy conservation. The GSF method is developed to address both the issues of consistency and convergence of the energy, force and the torque. Furthermore, the compensating of charge or multipole ordering breakage in the SP method due to direct spherical truncation will remain intact for large $r_c$. The electrostatic potential energy between central molecule and any molecule inside cutoff radius is given by, |
| 138 |
|
\begin{equation} |
| 139 |
|
U_{SF}(\vec r)=U(\vec r) - U(\vec r_c)-(\vec r-\vec r_c)\cdot\vec \nabla U(\vec r)|_{r=r_c} |
| 140 |
+ |
\label{eq:GSF} |
| 141 |
|
\end{equation} |
| 142 |
|
where the third term converges more rapidly as compared to first two terms hence the contribution of the third term is very small for large $r_c$ value. Hence the GSF method similar to SP method for large $r_c$. Moreover, the force and torque derived from equation 3 are consistent with the energy and approaches to zero as $r$ tends to $r_c$. |
| 143 |
|
Both GSF and TSF methods are the generalization of the original DSF method to higher order multipole-multipole interactions. These two methods are same up to charge-dipole interaction level but generate different expressions in the energy, force and torque for the higher order multipole-multipole interactions. |
| 144 |
|
\section{Test} |
| 145 |
|
\subsection{Test with Ewald} |
| 146 |
|
We have compared the electrostatic force and torque of each molecule from SP, TSF and GSF method with the multipolar-Ewald method. Furthermore, total electrostatic energies of a molecular system from the different methods have also been compared with total energy from the Ewald. In Mote Carlo (MC) simulation, the energy difference between different configurations of the molecular system is important, even though absolute energies are not accurate. The Metropolis Monte Carlo algorithm states that the new configuration of the molecular system is accepted if the energy difference between the new and previous configuration $(\triangle E) < 0$ or if any random number $R< exp(-\triangle E/kT)$, where R is between 0 to 1, for the case $\triangle E >0$. We have compared the change in electrostatic potential energy $\triangle E$ of 250 different configurations of the various multipolar molecular systems (Section IV B) calculated from the Hard, SP, GSF, and TSF methods with the well-known Ewald method. In MD simulations, the force and torque acting on the molecules drives the whole dynamics of the molecules in a system. The magnitudes of the electrostatic force, torque and their direction for each molecule of the all 250 configurations have also been compared against the Ewald’s method. |
| 147 |
< |
We have used least square regression analyses for six different molecular systems to compare $\triangle E$ from Hard, SP, GSF, and TSF with the reference method. Molecular systems were run longer enough to explore various configurations and 250 independent configurations were recorded for comparison. The total numbers of 31,125 energy differences from the proposed methods have been compared with the Ewald. Similarly, the magnitudes of the forces and torques have also been compared by using least square regression analyses. In the forces and torques comparison, the magnitudes of the forces acting in each molecule for each configuration were evaluated. For example, our dipolar liquid simulation contains 2048 molecules and there are 250 different configurations for each system thus there are 512,000 force and torque comparisons. The correlation coefficient and correlation slope varies from 0 to 1, where 1 is the best agreement between the two methods. The directions of the forces and torques were compared by using formula, , where and are unit vectors along direction of the force (or torque) evaluated by Ewald and new methods respectively and is angle between them. The standard deviation about the zero for the angle between the forces and torques was computed for every molecule in each configuration. |
| 147 |
> |
We have used least square regression analyses for six different molecular systems to compare $\triangle E$ from Hard, SP, GSF, and TSF with the reference method. Molecular systems were run longer enough to explore various configurations and 250 independent configurations were recorded for comparison. The total numbers of 31,125 energy differences from the proposed methods have been compared with the Ewald. Similarly, the magnitudes of the forces and torques have also been compared by using least square regression analyses. In the forces and torques comparison, the magnitudes of the forces acting in each molecule for each configuration were evaluated. For example, our dipolar liquid simulation contains 2048 molecules and there are 250 different configurations for each system thus there are 512,000 force and torque comparisons. The correlation coefficient and correlation slope varies from 0 to 1, where 1 is the best agreement between the two methods. |
| 148 |
> |
|
| 149 |
> |
R.A. Fisher has developed a probablity density function to analyse directional data sets is expressed as below, |
| 150 |
> |
\begin{equation} |
| 151 |
> |
p_f(\theta) = \frac{\kappa}{2 \sinh\kappa}\sin\theta \exp(\kappa \cos\theta) |
| 152 |
> |
\label{eq:pdf} |
| 153 |
> |
\end{equation} |
| 154 |
> |
where $\kappa$ measures directional dispersion of the data about mean direction can be estimated as a reciprocal of the circular variance for large number of directional data sets.\cite{vector_statistics} In our calculation, the unit vector from the Ewald method was considered as mean direction and the angle between the vectors from Ewald and the purposed method were evaluated.The total displacement of the unit vectors from the purposed method was calculated as, |
| 155 |
> |
\begin{equation} |
| 156 |
> |
R = \sqrt{(\sum\limits_{i=1}^N \sin\theta_i)^2 + (\sum\limits_{i=1}^N \sin\theta_i)^2} |
| 157 |
> |
\label{eq:displacement} |
| 158 |
> |
\end{equation} |
| 159 |
> |
where N is number of directional data sets and $theta_i$ are the angles between unit vectors evaluated from the Ewald and the purposed methods. The circular variance is defined as $ Var(\theta) = 1 -R/N$. The value of circular variance varies from 0 to 1. The lower the value of $Var{\theta}$ is higher the value of $\kappa$, which expresses tighter clustering of the direction sets around Ewald direction. |
| 160 |
> |
|
| 161 |
|
\subsection{Modeled systems} |
| 162 |
|
We studied the comparison of the energy differences, forces and torques for six different systems; i) dipolar liquid, ii) quadrupolar liquid, iii) dipolar crystal, iv) quadrupolar crystal v) dipolar-quadrupolar liquid(SSDQ), and vi) ions in dipolar-qudrupolar liquid(SSDQC). To simulate different configurations of the crystals, the body centered cubic (BCC) minimum energy crystal with 3,456 molecules was taken and translationally locked in their respective crystal sites. The thermal energy was supplied to the rotational motion so that dipoles or quadrupoles can freely explore all possible orientation. The crystals were simulated for 10,000 fs in NVE ensemble at 50 K and 250 different configurations was taken in equal time interval for the comparative study. The crystals were not simulated at high temperature and for a long run time to avoid possible translational deformation of the crystal sites. |
| 163 |
|
For dipolar, quadrupolar, and dipolar-quadrupolar liquids simulation, each molecular system with 2,048 molecules simulated for 1ns in NVE ensemble at 300 K temperature after equilibration. We collected 250 different configurations in equal interval of time. For the ions mixed liquid system, we converted 48 different molecules into 24 $Na^+$ and $24 Cl^-$ ions and equilibrated. After equilibration, the system was run at the same environment for 1ns and 250 configurations were collected. While comparing energies, forces, and torques with Ewald method, Lennad Jone’s potentials were turned off and purely electrostatic interaction had been compared. |
| 164 |
|
\subsection{Summation methods} |
| 165 |
< |
The Ewald summation for charge, dipole, and quadurpole was performed by using multipolar Ewald’s code in OpendMD/2.1. For different types of multipolar systems, damping alpha (α) for Ewald’s method was derived by plotting total electrostatic energy versus cutoff radius for the different values of $\alpha$ and compared the converged potential energy with the converged electrostatic energy from the pure cutoff method for very large cutoff radius (20 $A^o$). We found 0.3 $A^o^{-1}$ damping alpha for 12 $A^o$ cutoff radius is suitable for all kind of multipolar systems. |
| 165 |
> |
The Ewald summation for charge, dipole, and quadurpole was performed by using multipolar Ewald’s code in OpendMD/2.1. For different types of multipolar systems, damping alpha (α) for Ewald’s method was derived by plotting total electrostatic energy versus cutoff radius for the different values of $\alpha$ and compared the converged potential energy with the converged electrostatic energy from the pure cutoff method for very large cutoff radius (20 $A^o$). We found 0.3 $(A^o)^{-1}$ damping alpha for 12 $A^o$ cutoff radius is suitable for all kind of multipolar systems. |
| 166 |
|
The energies, forces and torques for all methods i) Hard, ii) SP, iii) GSF, and iv) TSF are evaluated for different cutoff radii i) 9, ii) 12, and iii) 15 $A^o$ with damping parameter ($\alpha$) 0.0, 0.1, 0.2, and 0.3 then compared with Ewald’s method. All the simulations for the various systems were conducted in the OpenMD/2.1. |
| 167 |
|
\subsection{Energy conservation and efficiency} |
| 168 |
|
To test conservation of the energy, the mixed molecular system of 2000 dipolar-quadrupolar molecules with 24 $Na^+$, and 24 $Cl^-$ was run for 1ns in the microcanonical ensemble at 300 K temperature for different cutoff methods (Ewald, Hard, SP, GSF, and TSF). The molecular system was run in 12 parallel computers and started with same initial positions and velocities for all cutoff methods. The slope and Standard Deviation of the energy about the slope (SD) were evaluated in the total energy versus time plot, where the slope evaluates the total energy drift and SD calculates the energy fluctuation in MD simulations. Also, the time duration for the simulation was recorded to compare efficiency of the purposed methods with the Ewald. |
| 170 |
|
\section{RESULTS} |
| 171 |
|
\subsection{Electrostatic energy and configurational energy differences} |
| 172 |
|
The magnitude of the fluctuation in the total electrostatic energy per molecule for a dipolar crystal is very high as shown in (Fig … paper I).\cite{PaperI}As soon as, the net dipole moment within a cutoff radius is neutralized in the SP method, the magnitude of the fluctuation in the total electrostatic energy per molecule reduced significantly and rapidly converged to the correct energy constant (Refer figure … Paper I).\cite{PaperI} The GSF potential energy also converged to the correct energy constant for the cutoff radius rc = 6a for the undamped case. The potential energy from the TSF method converges towards the correct value for a very large cutoff radius. The speed of convergence for the all the cutoff methods can be increased by using damping function as shown in figure … Paper I\cite{PaperI}. For the quadrupolar crystals, the fluctuation in the total electrostatic energy for the hard cutoff method is small and short ranged as compared to the dipolar crystals (figure … in the paper I).\cite{PaperI} Similar to the dipolar crystals, the net quadrupolar neutralization of the cutoff sphere in the SP method reduces oscillation rapidly and converge electrostatic energy to the correct energy constant. |
| 173 |
< |
The oscillation in the the electrostatic energy for the hard cutoff method is even true for the dipolar liquids as shown in Figure 3. As we placed image on the surface of the cutoff sphere in SP method, the oscillation in the energy is reduced. The fluctuation in the energy in liquid is much smaller as compared to the crystal (This result is similar to the results observed by \textit{Wolf et al.} in the case of ionic crystal and Mgo melt). The large magnitude in the fluctuation of the electrostatic energy in the crystal is because of large range of multipole ordering in the crystal. When the energy is evaluated by the direct truncation, it breaks up large number of multipolar ordering leaving behind net multipole moment. But in the case of liquid, there is only local ordering of the multipoles and their ordering disappears in the long range. Therefore, the direct truncation results a small oscillation in the electrostatic energy (which is smaller than deviation SP energy from the Ewald -refer figure 3) in the case of liquid. Although, the oscillation in the energy is very small for the case of liquid, this affects the change in potential energy ($\nabla E$), which is observed when $\nabla E$ evaluated from the SP method compared with Ewald as shown in figure 4a and 4b. |
| 173 |
> |
The oscillation in the the electrostatic energy for the hard cutoff method is even true for the dipolar liquids as shown in Figure ~\ref{fig:rcutConvergence_dipolarLiquid}. As we placed image on the surface of the cutoff sphere in SP method, the oscillation in the energy is reduced. The fluctuation in the energy in liquid is much smaller as compared to the crystal (This result is similar to the results observed by \textit{Wolf et al.} in the case of ionic crystal and Mgo melt). The large magnitude in the fluctuation of the electrostatic energy in the crystal is because of large range of multipole ordering in the crystal. When the energy is evaluated by the direct truncation, it breaks up large number of multipolar ordering leaving behind net multipole moment. But in the case of liquid, there is only local ordering of the multipoles and their ordering disappears in the long range. Therefore, the direct truncation results a small oscillation in the electrostatic energy (which is smaller than deviation SP energy from the Ewald Figure ~\ref{fig:rcutConvergence_dipolarLiquid}) in the case of liquid. Although, the oscillation in the energy is very small for the case of liquid, this affects the change in potential energy ($\triangle E$), which is observed when $\triangle E$ evaluated from the SP method compared with Ewald as shown in figure 4a and 4b. |
| 174 |
|
\begin{figure}[h!] |
| 175 |
|
\centering |
| 176 |
< |
\includegraphics[width=0.50 \textwidth]{rcutConvergence_dipolarLiquid.pdf} |
| 176 |
> |
\includegraphics[width=0.50 \textwidth]{rcutConvergence_dipolarLiquid-crop.pdf} |
| 177 |
|
\caption{The energy per molecule plotted against cutoff radius, rc for i) Hard ii) SP iii) GSF, and iv) TSF method. The hard cutoff method shows fluctuation in the electorstatic energy and it disappers in all other methods. } |
| 178 |
< |
\label{figure1} |
| 178 |
> |
\label{fig:rcutConvergence_dipolarLiquid} |
| 179 |
|
\end{figure} |
| 180 |
< |
In MC simulations, the electrostatic differences between the molecules are important parameter for sampling. We have compared $\nabla E$ from the different methods (Hard, SP, GSF, and TSF) with the Ewald using linear regression analysis. The correlation coefficient ($R^2$) of the regression line measures the deviation of the evaluated quantities from the mean slope. We know that Ewald method evaluates accurate value of the electrostatic energy. Hence, if the proposed methods can quantify the electrostatic energy as good as Ewald then the correlation coefficient is 1.The correlation coefficient is 0 for the completely random result for any physical quantity measured by the proposed method. The slope is a measure of the accuracy of the average of a physical quantity obtained from the proposed methods. If the slope is 1 then we can conclude that the average of the physical quantity measured by the method is as good as Ewald. The deviation of the slope from 1 state that the method used in quantifying physical quantities is statistically biased as compared to Ewald. |
| 180 |
> |
In MC simulations, the electrostatic differences between the molecules are important parameter for sampling. We have compared $\triangle E$ from the different methods (Hard, SP, GSF, and TSF) with the Ewald using linear regression analysis. The correlation coefficient ($R^2$) of the regression line measures the deviation of the evaluated quantities from the mean slope. We know that Ewald method evaluates accurate value of the electrostatic energy. Hence, if the proposed methods can quantify the electrostatic energy as good as Ewald then the correlation coefficient is 1.The correlation coefficient is 0 for the completely random result for any physical quantity measured by the proposed method. The slope is a measure of the accuracy of the average of a physical quantity obtained from the proposed methods. If the slope is 1 then we can conclude that the average of the physical quantity measured by the method is as good as Ewald. The deviation of the slope from 1 state that the method used in quantifying physical quantities is statistically biased as compared to Ewald. |
| 181 |
|
\begin{figure} |
| 182 |
|
\centering |
| 183 |
|
\includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf} |
| 184 |
+ |
\label{fig:barGraph1} |
| 185 |
|
\end{figure} |
| 186 |
|
\begin{figure} |
| 187 |
|
\centering |
| 188 |
|
\includegraphics[width=0.45 \textwidth]{slopeComparision_undamped.pdf} |
| 189 |
< |
\caption{The energy per molecule plotted against cutoff radius, rc for i) Hard ii) SP iii) GSF, and iv) TSF method. The hard cutoff method shows fluctuation in the electorstatic energy and it disappers in all other methods. } |
| 190 |
< |
\label{figure1} |
| 189 |
> |
\caption{} |
| 190 |
> |
|
| 191 |
> |
\label{fig:barGraph2} |
| 192 |
|
\end{figure} |
| 193 |
|
The correlation coefficient ($R^2$) and slope of the linear regression plots for the energy differences for all six different molecular systems is shown in figure 4a and 4b.The plot shows that the correlation coefficient improves for the SP cutoff method as compared to the undamped hard cutoff method in the case of SSDQC, SSDQ, dipolar crystal, and dipolar liquid. For the quadrupolar crystal and liquid, the correlation coefficient is almost unchanged and close to 1. The correlation coefficient is smallest (0.696276 for $r_c$ = 9 $A^o$) for the SSDQC liquid because of the presence of charge-charge and charge-multipole interactions. Since the charge-charge and charge-multipole interaction is long ranged, there is huge deviation of correlation coefficient from 1. Similarly, the quarupole–quadrupole interaction is short ranged ($\sim 1/r^6$) with compared to interactions in the other multipolar systems, thus the correlation coefficient very close to 1 even for hard cutoff method. The idea of placing image multipole on the surface of the cutoff sphere improves the correlation coefficient and makes it close to 1 for all types of multipolar systems. Similarly the slope is hugely deviated from the correct value for the lower order multipole-multipole interaction and slightly deviated for higher order multipole – multipole interaction. The SP method improves both correlation coefficient ($R^2$) and slope significantly in SSDQC and dipolar systems. The Slope is found to be deviated more in dipolar crystal as compared to liquid which is associated with the large fluctuation in the electrostatic energy in crystal. The GSF also produced better values of correlation coefficient and slope with the proper selection of the damping alpha (Interested reader can consult accompanying supporting material). The TSF method gives good value of correlation coefficient for the dipolar crystal, dipolar liquid, SSDQ, and SSDQC (not for the quadrupolar crystal and liquid) but the regression slopes are significantly deviated. |
| 194 |
|
\begin{figure} |
| 195 |
|
\centering |
| 196 |
< |
\includegraphics[width=0.50 \textwidth]{energy_combined_slope_correlation.pdf} |
| 196 |
> |
\includegraphics[width=0.50 \textwidth]{energyPlot_slopeCorrelation_combined-crop.pdf} |
| 197 |
|
\caption{The correlation coefficient and regression slope of configurational energy differences for a given method with compared with the reference Ewald method. The value of result equal to 1(dashed line) indicates energy difference is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle)} |
| 198 |
< |
\label{figure1} |
| 198 |
> |
\label{fig:slopeCorr_energy} |
| 199 |
|
\end{figure} |
| 200 |
< |
The combined correlation coefficient and slope for all six systems is shown in Figure 5. The correlation coefficient for the undamped hard cutoff method is does not have good agreement with the Ewald because of the fluctuation of the electrostatic energy in the direct truncation method. This deviation in correlation coefficient is improved by using SP, GSF, and TSF method. But the TSF method worsens the regression slope stating that this method produces statistically more biased result as compared to Ewald. Also the GSF method slightly deviate slope but it can be alleviated by using proper value of damping alpha and cutoff radius. The SP method shows good agreement with Ewald method for all values of damping alpha and radii. |
| 200 |
> |
The combined correlation coefficient and slope for all six systems is shown in Figure ~\ref{fig:slopeCorr_energy}. The correlation coefficient for the undamped hard cutoff method is does not have good agreement with the Ewald because of the fluctuation of the electrostatic energy in the direct truncation method. This deviation in correlation coefficient is improved by using SP, GSF, and TSF method. But the TSF method worsens the regression slope stating that this method produces statistically more biased result as compared to Ewald. Also the GSF method slightly deviate slope but it can be alleviated by using proper value of damping alpha and cutoff radius. The SP method shows good agreement with Ewald method for all values of damping alpha and radii. |
| 201 |
|
\subsection{Magnitude of the force and torque vectors} |
| 202 |
< |
The comparison of the magnitude of the combined forces and torques for the data accumulated from all system types are shown in Figure 6. The correlation and slope for the forces agree with the Ewald even for the hard cutoff method. For the system of molecules with higher order multipoles, the interaction is short ranged. Moreover, the force decays more rapidly than the electrostatic energy hence the hard cutoff method also produces good results. Although the pure cutoff gives the good match of the electrostatic force, the discontinuity in the force at the cutoff radius causes problem in the total energy conservation in MD simulations, which will be discussed in detail in section D. The correlation coefficient for GSF method also perfectly matches with Ewald but the slope is slightly deviated (due to extra term obtained from the angular differentiation). This deviation in the slope can be alleviated with proper selection of the damping alpha and radii ($\alpha = 0.2$ and $r_c = 12 A^o$ are good choice). The TSF method shows good agreement in the correlation coefficient but the slope is not good as compared to the Ewald. |
| 202 |
> |
The comparison of the magnitude of the combined forces and torques for the data accumulated from all system types are shown in Figure ~\ref{fig:slopeCorr_force}. The correlation and slope for the forces agree with the Ewald even for the hard cutoff method. For the system of molecules with higher order multipoles, the interaction is short ranged. Moreover, the force decays more rapidly than the electrostatic energy hence the hard cutoff method also produces good results. Although the pure cutoff gives the good match of the electrostatic force, the discontinuity in the force at the cutoff radius causes problem in the total energy conservation in MD simulations, which will be discussed in detail in section D. The correlation coefficient for GSF method also perfectly matches with Ewald but the slope is slightly deviated (due to extra term obtained from the angular differentiation). This deviation in the slope can be alleviated with proper selection of the damping alpha and radii ($\alpha = 0.2$ and $r_c = 12 A^o$ are good choice). The TSF method shows good agreement in the correlation coefficient but the slope is not good as compared to the Ewald. |
| 203 |
|
\begin{figure} |
| 204 |
|
\centering |
| 205 |
< |
\includegraphics[width=0.50 \textwidth]{force_combined_slope_correlation.pdf} |
| 205 |
> |
\includegraphics[width=0.50 \textwidth]{forcePlot_slopeCorrelation_combined-crop.pdf} |
| 206 |
|
\caption{The correlation coefficient and regression slope of the magnitude of the force for a given method with compared to the reference Ewald method. The value of result equal to 1(dashed line) indicates, the magnitude of the force from a method is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle). } |
| 207 |
< |
\label{figure1} |
| 207 |
> |
\label{fig:slopeCorr_force} |
| 208 |
|
\end{figure} |
| 209 |
< |
The torques appears to be very influenced because of extra term generated when the potential energy is modified to get consistent force and torque. The result shows that the torque from the hard cutoff method has good agreement with Ewald. As the potential is modified to make it consistent with the force and torque, the correlation and slope is deviated as shown in figure 7 for SP, GSF and TSF cutoff methods. But the proper value of the damping alpha and radius can improve the agreement of the GSF with the Ewald method. The TSF method shows worst agreement in the slope as compared to Ewald even for larger cutoff radii. |
| 209 |
> |
The torques appears to be very influenced because of extra term generated when the potential energy is modified to get consistent force and torque. The result shows that the torque from the hard cutoff method has good agreement with Ewald. As the potential is modified to make it consistent with the force and torque, the correlation and slope is deviated as shown in Figure~\ref{fig:slopeCorr_torque} for SP, GSF and TSF cutoff methods. But the proper value of the damping alpha and radius can improve the agreement of the GSF with the Ewald method. The TSF method shows worst agreement in the slope as compared to Ewald even for larger cutoff radii. |
| 210 |
|
\begin{figure} |
| 211 |
|
\centering |
| 212 |
< |
\includegraphics[width=0.50 \textwidth]{torque_combined_slope_correlation.pdf} |
| 212 |
> |
\includegraphics[width=0.5 \textwidth]{torquePlot_slopeCorrelation_combined-crop.pdf} |
| 213 |
|
\caption{The correlation coefficient and regression slope of the magnitude of the force for a given method with compared to the reference Ewald method. The value of result equal to 1(dashed line) indicates, the magnitude of the force from a method is indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle).} |
| 214 |
< |
\label{figure1} |
| 214 |
> |
\label{fig:slopeCorr_torque} |
| 215 |
|
\end{figure} |
| 216 |
|
\subsection{Directionality of the force and torque vectors} |
| 217 |
< |
The accurate evaluation of the direction of the force and torques are also important for the dynamic simulation. In our research, the direction computed from the new methods was compared with the Ewald by measuring angle between them. And, the standard deviation of the angles between the force vectors from the proposed methods and Ewald evaluates the accuracy of the new methods with compared to Ewald. Figure 8 shows that the standard deviation of the combined force and torque vectors for all molecules of the 250 configurations in all system types. |
| 218 |
< |
The result shows that hard cutoff and SP methods show good directional agreement with Ewald in the case of force. The GSF method generate force within few degrees of Ewald method, for example alpha = 0.2 and $r = 12A^o$ the standard deviation is equal to 4.51$^o$. The TSF method is the poorest in evaluating accurate direction as compared to Hard, SP and GSF methods. The direction of the torques has larger deviation of the angle as compared to the direction of the forces. For the same damping $\alpha = 0.2$ and $rc = 12 A^o$, the standard deviation of the torque is found to be 12.01$^o$ for the GSF method. The standard deviation of the both force and torque can also be minimized by varying damping alpha. |
| 217 |
> |
The accurate evaluation of the direction of the force and torques are also important for the dynamic simulation.In our research, the direction data sets were computed from the purposed method and compared with Ewald using Fisher statistics and results are expressed in terms of circular variance ($Var(\theta$).The force and torque vectors from the purposed method followed Fisher probability distribution function expressed in equation~\ref{eq:pdf}. The circular variance for the force and torque vectors of each molecule in the 250 configurations for all system types is shown in Figure~\ref{fig:slopeCorr_circularVariance}. The direction of the force and torque vectors from hard and SP cutoff methods showed best directional agreement with the Ewald. The force and torque vectors from GSF method also showed good agreement with the Ewald method, which can also be improved by varying damping alpha and cutoff radius.For $\alpha = 0.2$ and $r_c = 12 A^o$, $ Var(\theta) $ for direction of the force was found to be 0.002061 and corresponding value of $\kappa $ was 485.20. Integration of equation \ref{eq:pdf} for that corresponding value of $\kappa$ showed that 95\% of force vectors are with in $6.37^o$. The TSF method is the poorest in evaluating accurate direction with compared to Hard, SP, and GSF methods. The circular variance for the direction of the torques is larger as compared to force. For same $\alpha = 0.2, r_c = 12 A^o$ and GSF method, the circular variance was 0.01415, which showed 95\% of torque vectors are within $16.75^o$.The direction of the force and torque vectors can be improved by varying $\alpha$ and $r_c$. |
| 218 |
> |
|
| 219 |
|
\begin{figure} |
| 220 |
|
\centering |
| 221 |
< |
\includegraphics[width=0.50 \textwidth]{angle_forceNtorque_modified.pdf} |
| 222 |
< |
\caption{The standard deviation of the angle between force and torque vectors from a given method and the reference Ewald method. The value of result equal to 0 (dashed line) indicates direction of the vectors are indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle)} |
| 223 |
< |
\label{figure1} |
| 221 |
> |
\includegraphics[width=0.5 \textwidth]{Variance_forceNtorque_modified-crop.pdf} |
| 222 |
> |
\caption{The circular variance of the data sets of the direction of the force and torque vectors obtained from a given method about reference Ewald method. The result equal to 0 (dashed line) indicates direction of the vectors are indistinguishable from the Ewald method. Here different symbols represent different value of the cutoff radius (9 $A^o$ = circle, 12 $A^o$ = square 15 $A^o$ = inverted triangle)} |
| 223 |
> |
\label{fig:slopeCorr_circularVariance} |
| 224 |
|
\end{figure} |
| 225 |
|
\subsection{Total energy conservation and efficiency } |
| 226 |
< |
We have tested the conservation of energy in the SSDQC liquid system by running system for 1ns in the Hard, SP, GSF and TSF method. The Hard cutoff method shows very high energy drifts 433.53 KCal/Mol/ns/particle and energy fluctuation 32574.53 Kcal/Mol (measured by the SD from the slope) for the undamped case, which makes it completely unusable in MD simulations. The SP method also shows large value of energy drift 1.289 Kcal/Mol/ns/particle and energy fluctuation 0.01953 Kcal/Mol. The energy fluctuation in the SP method is due to the non-vanishing nature of the torque and force at the cutoff radius. We can improve the energy conservation in some extent by the proper selection of the damping alpha but the improvement is not good enough, which can be observed in Figure 9a and 9b .The GSF and TSF shows very low value of energy drift 0.09016, 0.07371 KCal/Mol/ns/particle and fluctuation 3.016E-6, 1.217E-6 Kcal/Mol respectively for the undamped case. Since the absolute value of the evaluated electrostatic energy, force and torque from TSF method are deviated from the Ewald, it does not mimic MD simulations appropriately. The electrostatic energy, force and torque from the GSF method have very good agreement with the Ewald. In addition, the energy drift and energy fluctuation from the GSF method is much better than Ewald’s method for reciprocal space vector value ($k_f$) equal to 7 as shown in Figure 9a and 9b. We can improve the total energy fluctuation and drift for the Ewald’s method by increasing size of the reciprocal space, which extremely increseses the simulation time. In our current simulation, the simulation time for the Hard, SP, and GSF methods are about 5.5 times faster than the Ewald method. |
| 226 |
> |
We have tested the conservation of energy in the SSDQC liquid system by running system for 1ns in the Hard, SP, GSF and TSF method. The Hard cutoff method shows very high energy drifts 433.53 KCal/Mol/ns/particle and energy fluctuation 32574.53 Kcal/Mol (measured by the SD from the slope) for the undamped case, which makes it completely unusable in MD simulations. The SP method also shows large value of energy drift 1.289 Kcal/Mol/ns/particle and energy fluctuation 0.01953 Kcal/Mol. The energy fluctuation in the SP method is due to the non-vanishing nature of the torque and force at the cutoff radius. We can improve the energy conservation in some extent by the proper selection of the damping alpha but the improvement is not good enough, which can be observed in Figure 9a and 9b .The GSF and TSF shows very low value of energy drift 0.09016, 0.07371 KCal/Mol/ns/particle and fluctuation 3.016E-6, 1.217E-6 Kcal/Mol respectively for the undamped case. Since the absolute value of the evaluated electrostatic energy, force and torque from TSF method are deviated from the Ewald, it does not mimic MD simulations appropriately. The electrostatic energy, force and torque from the GSF method have very good agreement with the Ewald. In addition, the energy drift and energy fluctuation from the GSF method is much better than Ewald’s method for reciprocal space vector value ($k_f$) equal to 7 as shown in Figure~\ref{fig:energyDrift} and ~\ref{fig:fluctuation}. We can improve the total energy fluctuation and drift for the Ewald’s method by increasing size of the reciprocal space, which extremely increseses the simulation time. In our current simulation, the simulation time for the Hard, SP, and GSF methods are about 5.5 times faster than the Ewald method. |
| 227 |
|
\begin{figure} |
| 228 |
|
\centering |
| 229 |
|
\includegraphics[width=0.45 \textwidth]{log(energyDrift).pdf} |
| 230 |
+ |
\label{fig:energyDrift} |
| 231 |
|
\end{figure} |
| 232 |
|
\begin{figure} |
| 233 |
|
\centering |
| 234 |
|
\includegraphics[width=0.45 \textwidth]{logSD.pdf} |
| 235 |
< |
\caption{The plot showing a) standard deviation, and b) total energy drift in the total energy conservation plot for different values of the damping alpha for different cut off methods. } |
| 236 |
< |
\label{figure1} |
| 235 |
> |
\caption{The plot showing (a) standard deviation, and (b) total energy drift in the total energy conservation plot for different values of the damping alpha for different cut off methods. } |
| 236 |
> |
\label{fig:fluctuation} |
| 237 |
|
\end{figure} |
| 238 |
|
\section{CONCLUSION} |
| 239 |
|
We have generalized the charged neutralized potential energy originally developed by the Wolf et al.\cite{Wolf99} for the charge-charge interaction to the charge-multipole and multipole-multipole interaction in the SP method for higher order multipoles. Also, we have developed GSF and TSF methods by implementing the modification purposed by Fennel and Gezelter\cite{Gezelter06} for the charge-charge interaction to the higher order multipoles to ensure consistency and smooth truncation of the electrostatic energy, force, and torque for the spherical truncation. The SP methods for multipoles proved its suitability in MC simulations. On the other hand, the results from the GSF method produced good agreement with the Ewald's energy, force, and torque. Also, it shows very good energy conservation in MD simulations. |
| 240 |
< |
The direct truncation of any molecular system without multipole neutralization creates the fluctuation in the electrostatic energy. This fluctuation in the energy is very large for the case of crystal because of long range of multipole ordering (Refer paper I).\cite{PaperI} This is also significant in the case of the liquid because of the local multipole ordering in the molecules. If the net multipole within cutoff radius neutralized within cutoff sphere by placing image multiples on the surface of the sphere, this fluctuation in the energy reduced significantly. Also, the multipole neutralization in the generalized SP method showed very good agreement with the Ewald as compared to direct truncation for the evaluation of the $\nabla E$ between the configurations. |
| 240 |
> |
The direct truncation of any molecular system without multipole neutralization creates the fluctuation in the electrostatic energy. This fluctuation in the energy is very large for the case of crystal because of long range of multipole ordering (Refer paper I).\cite{PaperI} This is also significant in the case of the liquid because of the local multipole ordering in the molecules. If the net multipole within cutoff radius neutralized within cutoff sphere by placing image multiples on the surface of the sphere, this fluctuation in the energy reduced significantly. Also, the multipole neutralization in the generalized SP method showed very good agreement with the Ewald as compared to direct truncation for the evaluation of the $\triangle E$ between the configurations. |
| 241 |
|
In MD simulations, the energy conservation is very important. The conservation of the total energy can be ensured by i) enforcing the smooth truncation of the energy, force and torque in the cutoff radius and ii) making the energy, force and torque consistent with each other. The GSF and TSF methods ensure the consistency and smooth truncation of the energy, force and torque at the cutoff radius, as a result show very good total energy conservation. But the TSF method does not show good agreement in the absolute value of the electrostatic energy, force and torque with the Ewald. The GSF method has mimicked Ewald’s force, energy and torque accurately and also conserved energy. Therefore, the GSF method is the suitable method for evaluating required force field in MD simulations. In addition, the energy drift and fluctuation from the GSF method is much better than Ewald’s method for finite-sized reciprocal space. |
| 242 |
+ |
\bibliographystyle{rev4-1} |
| 243 |
|
\bibliography{references} |
| 244 |
|
\end{document} |
| 245 |
|
|
