ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/chrisDissertation/Electrostatics.tex
(Generate patch)

Comparing trunk/chrisDissertation/Electrostatics.tex (file contents):
Revision 3016 by chrisfen, Thu Sep 21 23:21:37 2006 UTC vs.
Revision 3025 by chrisfen, Tue Sep 26 16:02:25 2006 UTC

# Line 1 | Line 1
1   \chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES}
2  
3 < In molecular simulations, proper accumulation of the electrostatic
4 < interactions is essential and is one of the most
3 > In molecular simulations, proper accumulation of electrostatic
4 > interactions is essential and one of the most
5   computationally-demanding tasks.  The common molecular mechanics force
6   fields represent atomic sites with full or partial charges protected
7 < by repulsive Lennard-Jones interactions.  This means that nearly
8 < every pair interaction involves a calculation of charge-charge forces.
7 > by repulsive Lennard-Jones interactions.  This means that nearly every
8 > pair interaction involves a calculation of charge-charge forces.
9   Coupled with relatively long-ranged $r^{-1}$ decay, the monopole
10   interactions quickly become the most expensive part of molecular
11   simulations.  Historically, the electrostatic pair interaction would
# Line 32 | Line 32 | Wolf {\it et al.},\cite{Wolf99} which we further exten
32  
33   In this chapter, we focus on a new set of pairwise methods devised by
34   Wolf {\it et al.},\cite{Wolf99} which we further extend.  These
35 < methods along with a few other mixed methods (i.e. reaction field) are
36 < compared with the smooth particle mesh Ewald
35 > methods, along with a few other mixed methods (i.e. reaction field),
36 > are compared with the smooth particle mesh Ewald
37   sum,\cite{Onsager36,Essmann99} which is our reference method for
38   handling long-range electrostatic interactions. The new methods for
39   handling electrostatics have the potential to scale linearly with
40 < increasing system size since they involve only a simple modification
40 > increasing system size, since they involve only a simple modification
41   to the direct pairwise sum.  They also lack the added periodicity of
42   the Ewald sum, so they can be used for systems which are non-periodic
43   or which have one- or two-dimensional periodicity.  Below, these
44 < methods are evaluated using a variety of model systems to
45 < establish their usability in molecular simulations.
44 > methods are evaluated using a variety of model systems to establish
45 > their usability in molecular simulations.
46  
47   \section{The Ewald Sum}
48  
# Line 92 | Line 92 | $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the
92   where $\alpha$ is the damping or convergence parameter with units of
93   \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to
94   $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric
95 < constant of the surrounding medium. The final two terms of
96 < equation (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term
97 < for interacting with a surrounding dielectric.\cite{Allen87} This
98 < dipolar term was neglected in early applications in molecular
99 < simulations,\cite{Brush66,Woodcock71} until it was introduced by de
100 < Leeuw {\it et al.} to address situations where the unit cell has a
101 < dipole moment which is magnified through replication of the periodic
102 < images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the
103 < system is said to be using conducting (or ``tin-foil'') boundary
104 < conditions, $\epsilon_{\rm S} = \infty$. Figure
105 < \ref{fig:ewaldTime} shows how the Ewald sum has been applied over
106 < time.  Initially, due to the small system sizes that could be
107 < simulated feasibly, the entire simulation box was replicated to
108 < convergence.  In more modern simulations, the systems have grown large
109 < enough that a real-space cutoff could potentially give convergent
110 < behavior.  Indeed, it has been observed that with the choice of a
111 < small $\alpha$, the reciprocal-space portion of the Ewald sum can be
112 < rapidly convergent and small relative to the real-space
113 < portion.\cite{Karasawa89,Kolafa92}
95 > constant of the surrounding medium. The final two terms of equation
96 > (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term for
97 > interacting with a surrounding dielectric.\cite{Allen87} This dipolar
98 > term was neglected in early applications of this technique in
99 > molecular simulations,\cite{Brush66,Woodcock71} until it was
100 > introduced by de Leeuw {\it et al.} to address situations where the
101 > unit cell has a dipole moment which is magnified through replication
102 > of the periodic images.\cite{deLeeuw80,Smith81} If this term is taken
103 > to be zero, the system is said to be using conducting (or
104 > ``tin-foil'') boundary conditions, $\epsilon_{\rm S} = \infty$.
105  
106   \begin{figure}
107   \includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf}
# Line 121 | Line 112 | convergence for the larger systems of charges that are
112   convergence for the larger systems of charges that are common today.}
113   \label{fig:ewaldTime}
114   \end{figure}
115 + Figure \ref{fig:ewaldTime} shows how the Ewald sum has been applied
116 + over time.  Initially, due to the small system sizes that could be
117 + simulated feasibly, the entire simulation box was replicated to
118 + convergence.  In more modern simulations, the systems have grown large
119 + enough that a real-space cutoff could potentially give convergent
120 + behavior.  Indeed, it has been observed that with the choice of a
121 + small $\alpha$, the reciprocal-space portion of the Ewald sum can be
122 + rapidly convergent and small relative to the real-space
123 + portion.\cite{Karasawa89,Kolafa92}
124  
125   The original Ewald summation is an $\mathcal{O}(N^2)$ algorithm.  The
126   convergence parameter $(\alpha)$ plays an important role in balancing
# Line 155 | Line 155 | Several studies have recognized that the inherent peri
155   bringing them more in line with the cost of the full 3-D summation.
156  
157   Several studies have recognized that the inherent periodicity in the
158 < Ewald sum can also have an effect on three-dimensional
159 < systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
160 < Solvated proteins are essentially kept at high concentration due to
161 < the periodicity of the electrostatic summation method.  In these
162 < systems, the more compact folded states of a protein can be
163 < artificially stabilized by the periodic replicas introduced by the
164 < Ewald summation.\cite{Weber00} Thus, care must be taken when
158 > Ewald sum can have an effect not just on reduced dimensionality
159 > system, but on three-dimensional systems as
160 > well.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
161 > As an example, solvated proteins are essentially kept at high
162 > concentration due to the periodicity of the electrostatic summation
163 > method.  In these systems, the more compact folded states of a protein
164 > can be artificially stabilized by the periodic replicas introduced by
165 > the Ewald summation.\cite{Weber00} Thus, care must be taken when
166   considering the use of the Ewald summation where the assumed
167   periodicity would introduce spurious effects.
168  
# Line 174 | Line 175 | short-ranged in condensed phase systems and that neutr
175   periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.}
176   observed that the electrostatic interaction is effectively
177   short-ranged in condensed phase systems and that neutralization of the
178 < charge contained within the cutoff radius is crucial for potential
178 > charges contained within the cutoff radius is crucial for potential
179   stability. They devised a pairwise summation method that ensures
180   charge neutrality and gives results similar to those obtained with the
181   Ewald summation.  The resulting shifted Coulomb potential includes
182   image-charges subtracted out through placement on the cutoff sphere
183   and a distance-dependent damping function (identical to that seen in
184 < the real-space portion of the Ewald sum) to aid convergence
184 > the real-space portion of the Ewald sum) to aid convergence:
185   \begin{equation}
186   V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}
187   - \lim_{r_{ij}\rightarrow R_\textrm{c}}
# Line 240 | Line 241 | and showed that this potential does fairly well at cap
241   \label{eq:ZahnPot}
242   \end{equation}
243   and showed that this potential does fairly well at capturing the
244 < structural and dynamic properties of water compared the same
244 > structural and dynamic properties of water compared with the same
245   properties obtained using the Ewald sum.
246  
247   \section{Simple Forms for Pairwise Electrostatics}\label{sec:PairwiseDerivation}
# Line 311 | Line 312 | v(r) = \frac{q_i q_j}{r},
312   v(r) = \frac{q_i q_j}{r},
313   \label{eq:Coulomb}
314   \end{equation}
315 < then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et
316 < al.}'s undamped prescription:
315 > then the {\sc sp} form will give Wolf {\it et al.}'s undamped
316 > prescription:
317   \begin{equation}
318   V_\textrm{SP}(r) =
319   q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad
# Line 334 | Line 335 | simulations.
335   forces at the cutoff radius which results in energy drift during MD
336   simulations.
337  
338 < The shifted force ({\sc sf}) form using the normal Coulomb potential
338 < will give,
338 > The {\sc sf} form using the normal Coulomb potential will give,
339   \begin{equation}
340   V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}}
341   + \left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right]
# Line 349 | Line 349 | This formulation has the benefits that there are no di
349   \label{eq:SFForces}
350   \end{equation}
351   This formulation has the benefits that there are no discontinuities at
352 < the cutoff radius, while the neutralizing image charges are present in
352 > the cutoff radius and the neutralizing image charges are present in
353   both the energy and force expressions.  It would be simple to add the
354   self-neutralizing term back when computing the total energy of the
355   system, thereby maintaining the agreement with the Madelung energies.
# Line 361 | Line 361 | insufficient for accurate determination of the energy
361   Wolf \textit{et al.} originally discussed the energetics of the
362   shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was
363   insufficient for accurate determination of the energy with reasonable
364 < cutoff distances.  The calculated Madelung energies fluctuated around
365 < the expected value as the cutoff radius was increased, but the
364 > cutoff distances.  The calculated Madelung energies fluctuated wildly
365 > around the expected value, but as the cutoff radius was increased, the
366   oscillations converged toward the correct value.\cite{Wolf99} A
367 < damping function was incorporated to accelerate the convergence; and
367 > damping function was incorporated to accelerate this convergence; and
368   though alternative forms for the damping function could be
369   used,\cite{Jones56,Heyes81} the complimentary error function was
370   chosen to mirror the effective screening used in the Ewald summation.
# Line 374 | Line 374 | v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
374   v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
375   \label{eq:dampCoulomb}
376   \end{equation}
377 < the shifted potential (Eq. (\ref{eq:SPPot})) becomes
377 > the {\sc sp} potential function (Eq. (\ref{eq:SPPot})) becomes
378   \begin{equation}
379   V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}
380   - \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right)
# Line 389 | Line 389 | + \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r
389   \quad r\leqslant R_\textrm{c}.
390   \label{eq:DSPForces}
391   \end{equation}
392 < Again, this damped shifted potential suffers from a
393 < force-discontinuity at the cutoff radius, and the image charges play
394 < no role in the forces.  To remedy these concerns, one may derive a
395 < {\sc sf} variant by including the derivative term in
396 < equation (\ref{eq:shiftingForm}),
392 > Again, this damped shifted potential suffers from a discontinuity in
393 > the forces at the cutoff radius, and the image charges play no role in
394 > the forces.  To remedy these concerns, one may derive a {\sc sf}
395 > variant by including the derivative term present in
396 > equation~(\ref{eq:shiftingForm}),
397   \begin{equation}
398   \begin{split}
399   V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&
# Line 423 | Line 423 | If the damping parameter $(\alpha)$ is set to zero, th
423   \end{split}
424   \end{equation}
425   If the damping parameter $(\alpha)$ is set to zero, the undamped case,
426 < equations (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly
427 < recovered from equations (\ref{eq:DSPPot} through \ref{eq:DSFForces}).
426 > equations (\ref{eq:SPPot}) through (\ref{eq:SFForces}) are correctly
427 > recovered from equations (\ref{eq:DSPPot}) through (\ref{eq:DSFForces}).
428  
429 < This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot}
429 > This new {\sc sf} potential is similar to equation (\ref{eq:ZahnPot})
430   derived by Zahn \textit{et al.}; however, there are two important
431   differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from equation
432   (\ref{eq:shiftingForm}) is equal to equation (\ref{eq:dampCoulomb})
# Line 436 | Line 436 | would be expected to have sudden jumps as particle dis
436   portion is different.  The missing $v_\textrm{c}$ term would not
437   affect molecular dynamics simulations (although the computed energy
438   would be expected to have sudden jumps as particle distances crossed
439 < $R_c$).  The sign problem is a potential source of errors, however.
440 < In fact, it introduces a discontinuity in the forces at the cutoff,
441 < because the force function is shifted in the wrong direction and
442 < doesn't cross zero at $R_\textrm{c}$.
439 > $R_c$); however, the sign problem is a potential source of errors.  In
440 > fact, equation~(\ref{eq:ZahnPot}) introduces a discontinuity in the
441 > forces at the cutoff, because the force function is shifted in the
442 > wrong direction and does not cross zero at $R_\textrm{c}$.
443  
444   Equations (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an
445   electrostatic summation method in which the potential and forces are
446   continuous at the cutoff radius and which incorporates the damping
447   function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of
448 < this paper, we will evaluate exactly how good these methods ({\sc sp},
449 < {\sc sf}, damping) are at reproducing the correct electrostatic
448 > this chapter, we will evaluate exactly how good these methods ({\sc
449 > sp}, {\sc sf}, damping) are at reproducing the correct electrostatic
450   summation performed by the Ewald sum.
451  
452  
# Line 455 | Line 455 | classical molecular mechanics simulations: Monte Carlo
455   As mentioned in the introduction, there are two primary techniques
456   utilized to obtain information about the system of interest in
457   classical molecular mechanics simulations: Monte Carlo (MC) and
458 < Molecular Dynamics (MD).  Both of these techniques utilize pairwise
458 > molecular dynamics (MD).  Both of these techniques utilize pairwise
459   summations of interactions between particle sites, but they use these
460   summations in different ways.
461  
# Line 476 | Line 476 | cumulative, one should expect greater deviation at lon
476   electrostatic summation techniques, the dynamics in the short term
477   will be indistinguishable.  Because error in MD calculations is
478   cumulative, one should expect greater deviation at longer times,
479 < although methods which have large differences in the force and torque
479 > and methods which have large differences in the force and torque
480   vectors will diverge from each other more rapidly.
481  
482   \subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
483  
484 + \begin{figure}
485 + \centering
486 + \includegraphics[width = 3.5in]{./figures/dualLinear.pdf}
487 + \caption{Example least squares regressions of the configuration energy
488 + differences for SPC/E water systems. The upper plot shows a data set
489 + with a poor correlation coefficient ($R^2$), while the lower plot
490 + shows a data set with a good correlation coefficient.}
491 + \label{fig:linearFit}
492 + \end{figure}
493   The pairwise summation techniques (outlined in section
494   \ref{sec:ESMethods}) were evaluated for use in MC simulations by
495   studying the energy differences between conformations.  We took the
# Line 494 | Line 503 | Sample correlation plots for two alternate methods are
503   correlation (slope) and correlation coefficient for these regressions
504   indicate perfect agreement between the alternative method and {\sc spme}.
505   Sample correlation plots for two alternate methods are shown in
506 < Fig. \ref{fig:linearFit}.
498 <
499 < \begin{figure}
500 < \centering
501 < \includegraphics[width = 3.5in]{./figures/dualLinear.pdf}
502 < \caption{Example least squares regressions of the configuration energy
503 < differences for SPC/E water systems. The upper plot shows a data set
504 < with a poor correlation coefficient ($R^2$), while the lower plot
505 < shows a data set with a good correlation coefficient.}
506 < \label{fig:linearFit}
507 < \end{figure}
506 > figure \ref{fig:linearFit}.
507  
508   Each of the seven system types (detailed in section \ref{sec:RepSims})
509   were represented using 500 independent configurations.  Thus, each of
510   the alternative (non-Ewald) electrostatic summation methods was
511   evaluated using an accumulated 873,250 configurational energy
512 < differences.
512 > differences. Results for and discussions regarding the individual
513 > analysis of each of the system types appear in appendix
514 > \ref{app:IndividualResults}, while the cumulative results over all the
515 > investigated systems appear below in section~\ref{sec:EnergyResults}.
516  
515 – Results and discussion for the individual analysis of each of the
516 – system types appear in appendix \ref{app:IndividualResults}, while the
517 – cumulative results over all the investigated systems appear below in
518 – sections \ref{sec:EnergyResults}.
519 –
517   \subsection{Molecular Dynamics and the Force and Torque
518   Vectors}\label{sec:MDMethods} We evaluated the pairwise methods
519   (outlined in section \ref{sec:ESMethods}) for use in MD simulations by
# Line 528 | Line 525 | forces (and torques) on each molecule in each configur
525   comparing $\Delta E$ values.  Instead of a single energy difference
526   between two system configurations, we compared the magnitudes of the
527   forces (and torques) on each molecule in each configuration.  For a
528 < system of 1000 water molecules and 40 ions, there are 1040 force
529 < vectors and 1000 torque vectors.  With 500 configurations, this
530 < results in 520,000 force and 500,000 torque vector comparisons.
531 < Additionally, data from seven different system types was aggregated
532 < before the comparison was made.
528 > system of 1000 water molecules and 40 ions, there are 1040 force and
529 > 1000 torque vectors.  With 500 configurations, this results in 520,000
530 > force and 500,000 torque vector comparisons.  Additionally, data from
531 > seven different system types was aggregated before comparisons were
532 > made.
533  
534   The {\it directionality} of the force and torque vectors was
535   investigated through measurement of the angle ($\theta$) formed
# Line 547 | Line 544 | between two different electrostatic summation methods,
544   unit sphere.  Since this distribution is a measure of angular error
545   between two different electrostatic summation methods, there is no
546   {\it a priori} reason for the profile to adhere to any specific
547 < shape. Thus, gaussian fits were used to measure the width of the
547 > shape. Thus, Gaussian fits were used to measure the width of the
548   resulting distributions. The variance ($\sigma^2$) was extracted from
549   each of these fits and was used to compare distribution widths.
550   Values of $\sigma^2$ near zero indicate vector directions
# Line 610 | Line 607 | crystals), so the systems studied were:
607   \item a high ionic strength solution of NaCl in water (1.1 M), and
608   \item a 6~\AA\  radius sphere of Argon in water.
609   \end{enumerate}
613 –
610   By utilizing the pairwise techniques (outlined in section
611   \ref{sec:ESMethods}) in systems composed entirely of neutral groups,
612   charged particles, and mixtures of the two, we hope to discern under
# Line 699 | Line 695 | inverted triangles).}
695   inverted triangles).}
696   \label{fig:delE}
697   \end{figure}
702 –
698   The most striking feature of this plot is how well the Shifted Force
699   ({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy
700   differences.  For the undamped {\sc sf} method, and the
# Line 721 | Line 716 | shows an excellent correlation and quality of fit with
716   For the {\sc sp} method, inclusion of electrostatic damping improves
717   the agreement with Ewald, and using an $\alpha$ of 0.2~\AA $^{-1}$
718   shows an excellent correlation and quality of fit with the {\sc spme}
719 < results, particularly with a cutoff radius greater than 12~\AA\.  Use
719 > results, particularly with a cutoff radius greater than 12~\AA . Use
720   of a larger damping parameter is more helpful for the shortest cutoff
721   shown, but it has a detrimental effect on simulations with larger
722   cutoffs.
# Line 735 | Line 730 | limitations, primarily that it was developed for use i
730  
731   The reaction field results illustrates some of that method's
732   limitations, primarily that it was developed for use in homogeneous
733 < systems; although it does provide results that are an improvement over
734 < those from an unmodified cutoff.
733 > systems. It does, however, provide results that are an improvement
734 > over those from an unmodified cutoff.
735  
736   \section{Magnitude of the Force and Torque Vector Results}\label{sec:FTMagResults}
737  
# Line 759 | Line 754 | inverted triangles).}
754   inverted triangles).}
755   \label{fig:frcMag}
756   \end{figure}
757 <
758 < Again, it is striking how well the Shifted Potential and Shifted Force
759 < methods are doing at reproducing the {\sc spme} forces.  The undamped and
760 < weakly-damped {\sc sf} method gives the best agreement with Ewald.
761 < This is perhaps expected because this method explicitly incorporates a
767 < smooth transition in the forces at the cutoff radius as well as the
757 > Again, it is striking how well the {\sc sp} and {\sc sf} methods
758 > reproduce the {\sc spme} forces.  The undamped and weakly-damped {\sc
759 > sf} method gives the best agreement with Ewald.  This is perhaps
760 > expected because this method explicitly incorporates a smooth
761 > transition in the forces at the cutoff radius as well as the
762   neutralizing image charges.
763  
764   Figure \ref{fig:frcMag}, for the most part, parallels the results seen
# Line 775 | Line 769 | method is generating usable forces.  Further increases
769  
770   With moderate damping and a large enough cutoff radius, the {\sc sp}
771   method is generating usable forces.  Further increases in damping,
772 < while beneficial for simulations with a cutoff radius of 9~\AA\ , is
772 > while beneficial for simulations with a cutoff radius of 9~\AA\ , are
773   detrimental to simulations with larger cutoff radii.
774  
775   The reaction field results are surprisingly good, considering the poor
776   quality of the fits for the $\Delta E$ results.  There is still a
777 < considerable degree of scatter in the data, but the forces correlate
778 < well with the Ewald forces in general.  We note that the reaction
779 < field calculations do not include the pure NaCl systems, so these
780 < results are partly biased towards conditions in which the method
781 < performs more favorably.
777 > considerable degree of scatter in the data, but in general, the forces
778 > correlate well with the Ewald forces.  We note that the pure NaCl
779 > systems were not included in the system set used in the reaction field
780 > calculations, so these results are partly biased towards conditions in
781 > which the method performs more favorably.
782  
783   \begin{figure}
784   \centering
# Line 798 | Line 792 | inverted triangles).}
792   inverted triangles).}
793   \label{fig:trqMag}
794   \end{figure}
801 –
795   Molecular torques were only available from the systems which contained
796   rigid molecules (i.e. the systems containing water).  The data in
797 < fig. \ref{fig:trqMag} is taken from this smaller sampling pool.
797 > figure \ref{fig:trqMag} is taken from this smaller sampling pool.
798  
799 < Torques appear to be much more sensitive to charges at a longer
800 < distance.   The striking feature in comparing the new electrostatic
801 < methods with {\sc spme} is how much the agreement improves with increasing
802 < cutoff radius.  Again, the weakly damped and undamped {\sc sf} method
803 < appears to be reproducing the {\sc spme} torques most accurately.  
799 > Torques appear to be much more sensitive to charge interactions at
800 > longer distances.  The most noticeable feature in comparing the new
801 > electrostatic methods with {\sc spme} is how much the agreement
802 > improves with increasing cutoff radius.  Again, the weakly damped and
803 > undamped {\sc sf} method appears to reproduce the {\sc spme} torques
804 > most accurately.
805  
806   Water molecules are dipolar, and the reaction field method reproduces
807   the effect of the surrounding polarized medium on each of the
# Line 816 | Line 810 | performs best of all of the methods on molecular torqu
810  
811   \section{Directionality of the Force and Torque Vector Results}\label{sec:FTDirResults}
812  
813 < It is clearly important that a new electrostatic method can reproduce
814 < the magnitudes of the force and torque vectors obtained via the Ewald
815 < sum. However, the {\it directionality} of these vectors will also be
816 < vital in calculating dynamical quantities accurately.  Force and
817 < torque directionalities were investigated by measuring the angles
818 < formed between these vectors and the same vectors calculated using
819 < {\sc spme}.  The results (Fig. \ref{fig:frcTrqAng}) are compared through the
820 < variance ($\sigma^2$) of the Gaussian fits of the angle error
821 < distributions of the combined set over all system types.
813 > It is clearly important that a new electrostatic method should be able
814 > to reproduce the magnitudes of the force and torque vectors obtained
815 > via the Ewald sum. However, the {\it directionality} of these vectors
816 > will also be vital in calculating dynamical quantities accurately.
817 > Force and torque directionalities were investigated by measuring the
818 > angles formed between these vectors and the same vectors calculated
819 > using {\sc spme}.  The results (figure \ref{fig:frcTrqAng}) are compared
820 > through the variance ($\sigma^2$) of the Gaussian fits of the angle
821 > error distributions of the combined set over all system types.
822  
823   \begin{figure}
824   \centering
# Line 839 | Line 833 | and 15~\AA\ = inverted triangles).}
833   and 15~\AA\ = inverted triangles).}
834   \label{fig:frcTrqAng}
835   \end{figure}
842 –
836   Both the force and torque $\sigma^2$ results from the analysis of the
837   total accumulated system data are tabulated in figure
838 < \ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc
839 < sp}) method would be essentially unusable for molecular dynamics
840 < unless the damping function is added.  The Shifted Force ({\sc sf})
841 < method, however, is generating force and torque vectors which are
842 < within a few degrees of the Ewald results even with weak (or no)
850 < damping.
838 > \ref{fig:frcTrqAng}. Here it is clear that the {\sc sp} method would
839 > be essentially unusable for molecular dynamics unless the damping
840 > function is added.  The {\sc sf} method, however, is generating force
841 > and torque vectors which are within a few degrees of the Ewald results
842 > even with weak (or no) damping.
843  
844   All of the sets (aside from the over-damped case) show the improvement
845   afforded by choosing a larger cutoff radius.  Increasing the cutoff
# Line 929 | Line 921 | However, at larger values of $\alpha$, it is possible
921   The complimentary error function inserted into the potential weakens
922   the electrostatic interaction as the value of $\alpha$ is increased.
923   However, at larger values of $\alpha$, it is possible to over-damp the
924 < electrostatic interaction and to remove it completely.  Kast
924 > electrostatic interaction and remove it completely.  Kast
925   \textit{et al.}  developed a method for choosing appropriate $\alpha$
926   values for these types of electrostatic summation methods by fitting
927   to $g(r)$ data, and their methods indicate optimal values of 0.34,
928   0.25, and 0.16~\AA$^{-1}$ for cutoff values of 9, 12, and 15~\AA\
929   respectively.\cite{Kast03} These appear to be reasonable choices to
930 < obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on
930 > obtain proper MC behavior (figure \ref{fig:delE}); however, based on
931   these findings, choices this high would introduce error in the
932 < molecular torques, particularly for the shorter cutoffs.  Based on our
933 < observations, empirical damping up to 0.2~\AA$^{-1}$ is beneficial,
934 < but damping may be unnecessary when using the {\sc sf} method.
932 > molecular torques, particularly for the shorter cutoffs.  Based on the
933 > above observations, empirical damping up to 0.2~\AA$^{-1}$ is
934 > beneficial, but damping may be unnecessary when using the {\sc sf}
935 > method.
936  
937  
938   \section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}\label{sec:ShortTimeDynamics}
# Line 948 | Line 941 | that a method similar (but not identical with) the dam
941   using equations (\ref{eq:ZahnPot}) and
942   (\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated
943   that a method similar (but not identical with) the damped {\sc sf}
944 < method resulted in properties very similar to those obtained when
944 > method resulted in properties very close to those obtained when
945   using the Ewald summation.  The properties they studied (pair
946   distribution functions, diffusion constants, and velocity and
947   orientational correlation functions) may not be particularly sensitive
948   to the long-range and collective behavior that governs the
949   low-frequency behavior in crystalline systems.  Additionally, the
950 < ionic crystals are the worst case scenario for the pairwise methods
950 > ionic crystals are a worst case scenario for the pairwise methods
951   because they lack the reciprocal space contribution contained in the
952   Ewald summation.
953  
954 < We are using two separate measures to probe the effects of these
954 > We used two separate measures to probe the effects of these
955   alternative electrostatic methods on the dynamics in crystalline
956 < materials.  For short- and intermediate-time dynamics, we are
957 < computing the velocity autocorrelation function, and for long-time
958 < and large length-scale collective motions, we are looking at the
959 < low-frequency portion of the power spectrum.
956 > materials.  For short- and intermediate-time dynamics, we computed the
957 > velocity autocorrelation function, and for long-time and large
958 > length-scale collective motions, we looked at the low-frequency
959 > portion of the power spectrum.
960  
961   \begin{figure}
962   \centering
# Line 977 | Line 970 | are stiffer than the moderately damped and {\sc spme}
970   are stiffer than the moderately damped and {\sc spme} methods.}
971   \label{fig:vCorrPlot}
972   \end{figure}
973 <
981 < The short-time decay of the velocity autocorrelation function through
973 > The short-time decay of the velocity autocorrelation functions through
974   the first collision are nearly identical in figure
975   \ref{fig:vCorrPlot}, but the peaks and troughs of the functions show
976   how the methods differ.  The undamped {\sc sf} method has deeper
977 < troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than
977 > troughs (see inset in figure \ref{fig:vCorrPlot}) and higher peaks than
978   any of the other methods.  As the damping parameter ($\alpha$) is
979   increased, these peaks are smoothed out, and the {\sc sf} method
980   approaches the {\sc spme} results.  With $\alpha$ values of 0.2~\AA$^{-1}$,
# Line 995 | Line 987 | important.
987   important.
988  
989   \section{Collective Motion: Power Spectra of NaCl Crystals}\label{sec:LongTimeDynamics}
998 –
999 – To evaluate how the differences between the methods affect the
1000 – collective long-time motion, we computed power spectra from long-time
1001 – traces of the velocity autocorrelation function. The power spectra for
1002 – the best-performing alternative methods are shown in
1003 – fig. \ref{fig:methodPS}.  Apodization of the correlation functions via
1004 – a cubic switching function between 40 and 50~ps was used to reduce the
1005 – ringing resulting from data truncation.  This procedure had no
1006 – noticeable effect on peak location or magnitude.
990  
991   \begin{figure}
992   \centering
# Line 1015 | Line 998 | functions of NaCl crystals at 1000~K while using {\sc
998   100~cm$^{-1}$ to highlight where the spectra differ.}
999   \label{fig:methodPS}
1000   \end{figure}
1001 + To evaluate how the differences between the methods affect the
1002 + collective long-time motion, we computed power spectra from long-time
1003 + traces of the velocity autocorrelation function. The power spectra for
1004 + the best-performing alternative methods are shown in
1005 + figure \ref{fig:methodPS}.  Apodization of the correlation functions via
1006 + a cubic switching function between 40 and 50~ps was used to reduce the
1007 + ringing resulting from data truncation.  This procedure had no
1008 + noticeable effect on peak location or magnitude.
1009  
1010   While the high frequency regions of the power spectra for the
1011   alternative methods are quantitatively identical with Ewald spectrum,
1012   the low frequency region shows how the summation methods differ.
1013   Considering the low-frequency inset (expanded in the upper frame of
1014 < figure \ref{fig:dampInc}), at frequencies below 100~cm$^{-1}$, the
1014 > figure \ref{fig:methodPS}), at frequencies below 100~cm$^{-1}$, the
1015   correlated motions are blue-shifted when using undamped or weakly
1016   damped {\sc sf}.  When using moderate damping ($\alpha =
1017 < 0.2$~\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly
1018 < identical correlated motion to the Ewald method (which has a
1017 > 0.2$~\AA$^{-1}$), both the {\sc sf} and {\sc sp} methods produce
1018 > correlated motions nearly identical to the Ewald method (which has a
1019   convergence parameter of 0.3119~\AA$^{-1}$).  This weakening of the
1020   electrostatic interaction with increased damping explains why the
1021   long-ranged correlated motions are at lower frequencies for the
1022   moderately damped methods than for undamped or weakly damped methods.
1023  
1033 – To isolate the role of the damping constant, we have computed the
1034 – spectra for a single method ({\sc sf}) with a range of damping
1035 – constants and compared this with the {\sc spme} spectrum.
1036 – Fig. \ref{fig:dampInc} shows more clearly that increasing the
1037 – electrostatic damping red-shifts the lowest frequency phonon modes.
1038 – However, even without any electrostatic damping, the {\sc sf} method
1039 – has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode.
1040 – Without the {\sc sf} modifications, an undamped (pure cutoff) method
1041 – would predict the lowest frequency peak near 325~cm$^{-1}$.  {\it
1042 – Most} of the collective behavior in the crystal is accurately captured
1043 – using the {\sc sf} method.  Quantitative agreement with Ewald can be
1044 – obtained using moderate damping in addition to the shifting at the
1045 – cutoff distance.
1046 –
1024   \begin{figure}
1025   \centering
1026   \includegraphics[width = \linewidth]{./figures/increasedDamping.pdf}
# Line 1056 | Line 1033 | motions.}
1033   motions.}
1034   \label{fig:dampInc}
1035   \end{figure}
1036 + To isolate the role of the damping constant, we have computed the
1037 + spectra for a single method ({\sc sf}) with a range of damping
1038 + constants and compared this with the {\sc spme} spectrum.  Figure
1039 + \ref{fig:dampInc} shows more clearly that increasing the electrostatic
1040 + damping red-shifts the lowest frequency phonon modes.  However, even
1041 + without any electrostatic damping, the {\sc sf} method has at most a
1042 + 10 cm$^{-1}$ error in the lowest frequency phonon mode.  Without the
1043 + {\sc sf} modifications, an undamped (pure cutoff) method would predict
1044 + the lowest frequency peak near 325~cm$^{-1}$, an error significantly
1045 + larger than that of the undamped {\sc sf} technique.  This indicates
1046 + that {\it most} of the collective behavior in the crystal is
1047 + accurately captured using the {\sc sf} method.  Quantitative agreement
1048 + with Ewald can be obtained using moderate damping in addition to the
1049 + shifting at the cutoff distance.
1050  
1051   \section{An Application: TIP5P-E Water}\label{sec:t5peApplied}
1052  
# Line 1114 | Line 1105 | p(t) =         \frac{1}{\textrm{d}V}\sum_\mu
1105          + \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right],
1106   \label{eq:MolecularPressure}
1107   \end{equation}
1108 < where $V$ is the volume, $\mathbf{P}_{\mu}$ is the momentum of
1109 < molecule $\mu$, $\mathbf{R}_\mu$ is the position of the center of mass
1110 < ($M_\mu$) of molecule $\mu$, and $\mathbf{F}_{\mu i}$ is the force on
1111 < atom $i$ of molecule $\mu$.\cite{Melchionna93} The virial term (the
1112 < right term in the brackets of equation \ref{eq:MolecularPressure}) is
1113 < directly dependent on the interatomic forces.  Since the {\sc sp}
1114 < method does not modify the forces (see
1115 < section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} will
1116 < be identical to that obtained without an electrostatic correction.
1117 < The {\sc sf} method does alter the virial component and, by way of the
1118 < modified pressures, should provide densities more in line with those
1119 < obtained using the Ewald summation.
1108 > where d is the dimensionality of the system, $V$ is the volume,
1109 > $\mathbf{P}_{\mu}$ is the momentum of molecule $\mu$, $\mathbf{R}_\mu$
1110 > is the position of the center of mass ($M_\mu$) of molecule $\mu$, and
1111 > $\mathbf{F}_{\mu i}$ is the force on atom $i$ of molecule
1112 > $\mu$.\cite{Melchionna93} The virial term (the right term in the
1113 > brackets of equation
1114 > \ref{eq:MolecularPressure}) is directly dependent on the interatomic
1115 > forces.  Since the {\sc sp} method does not modify the forces (see
1116 > section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp}
1117 > will be identical to that obtained without an electrostatic
1118 > correction.  The {\sc sf} method does alter the virial component and,
1119 > by way of the modified pressures, should provide densities more in
1120 > line with those obtained using the Ewald summation.
1121  
1122   To compare densities, $NPT$ simulations were performed with the same
1123   temperatures as those selected by Rick in his Ewald summation
# Line 1147 | Line 1139 | greater force on the central particle. The error bars
1139   Ewald summation, leading to slightly lower densities. This effect is
1140   more visible with the 9~\AA\ cutoff, where the image charges exert a
1141   greater force on the central particle. The error bars for the {\sc sf}
1142 < methods show plus or minus the standard deviation of the density
1143 < measurement at each temperature.}
1142 > methods show the average one-sigma uncertainty of the density
1143 > measurement, and this uncertainty is the same for all the {\sc sf}
1144 > curves.}
1145   \label{fig:t5peDensities}
1146   \end{figure}
1154 –
1147   Figure \ref{fig:t5peDensities} shows the densities calculated for
1148   TIP5P-E using differing electrostatic corrections overlaid on the
1149   experimental values.\cite{CRC80} The densities when using the {\sc sf}
# Line 1176 | Line 1168 | As a final note, all of the above density calculations
1168   important role in the resulting densities.
1169  
1170   As a final note, all of the above density calculations were performed
1171 < with systems of 512 water molecules. Rick observed a system sized
1171 > with systems of 512 water molecules. Rick observed a system size
1172   dependence of the computed densities when using the Ewald summation,
1173   most likely due to his tying of the convergence parameter to the box
1174   dimensions.\cite{Rick04} For systems of 256 water molecules, the
# Line 1231 | Line 1223 | identical.}
1223   identical.}
1224   \label{fig:t5peGofRs}
1225   \end{figure}
1234 –
1226   The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc
1227   sf} technique with a various parameters are overlaid on the
1228 < $g_\textrm{OO}(r)$ while using the Ewald summation. The differences in
1229 < density do not appear to have any effect on the liquid structure as
1230 < the $g_\textrm{OO}(r)$s are indistinguishable. These results indicate
1231 < that the $g_\textrm{OO}(r)$ is insensitive to the choice of
1232 < electrostatic correction.
1228 > $g_\textrm{OO}(r)$ while using the Ewald summation in figure
1229 > \ref{fig:t5peGofRs}. The differences in density do not appear to have
1230 > any effect on the liquid structure as the $g_\textrm{OO}(r)$s are
1231 > indistinguishable. These results indicate that the $g_\textrm{OO}(r)$
1232 > is insensitive to the choice of electrostatic correction.
1233  
1234   \subsection{Thermodynamic Properties}\label{sec:t5peThermo}
1235  
# Line 1368 | Line 1359 | property trends with temperature seen when using the E
1359  
1360   As observed for the density in section \ref{sec:t5peDensity}, the
1361   property trends with temperature seen when using the Ewald summation
1362 < are reproduced with the {\sc sf} technique. Differences include the
1363 < calculated values of $\Delta H_\textrm{vap}$ under-predicting the Ewald
1364 < values. This is to be expected due to the direct weakening of the
1365 < electrostatic interaction through forced neutralization in {\sc
1366 < sf}. This results in an increase of the intermolecular potential
1367 < producing lower values from equation (\ref{eq:DeltaHVap}). The slopes of
1368 < these values with temperature are similar to that seen using the Ewald
1369 < summation; however, they are both steeper than the experimental trend,
1370 < indirectly resulting in the inflated $C_p$ values at all temperatures.
1362 > are reproduced with the {\sc sf} technique. One noticable difference
1363 > between the properties calculated using the two methods are the lower
1364 > $\Delta H_\textrm{vap}$ values when using {\sc sf}. This is to be
1365 > expected due to the direct weakening of the electrostatic interaction
1366 > through forced neutralization. This results in an increase of the
1367 > intermolecular potential producing lower values from equation
1368 > (\ref{eq:DeltaHVap}). The slopes of these values with temperature are
1369 > similar to that seen using the Ewald summation; however, they are both
1370 > steeper than the experimental trend, indirectly resulting in the
1371 > inflated $C_p$ values at all temperatures.
1372  
1373   Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$
1374   values all overlap within error. As indicated for the $\Delta
# Line 1398 | Line 1390 | the onset of a more frustrated or glassy behavior for
1390   indicate a more pronounced transition in the supercooled regime,
1391   particularly in the case of {\sc sf} without damping. This points to
1392   the onset of a more frustrated or glassy behavior for TIP5P-E at
1393 < temperatures below 250~K in these simulations. Because the systems are
1394 < locked in different regions of phase-space, comparisons between
1395 < properties at these temperatures are not exactly fair. This
1396 < observation is explored in more detail in section
1397 < \ref{sec:t5peDynamics}.
1393 > temperatures below 250~K in the {\sc sf} simulations, indicating that
1394 > disorder in the reciprical-space term of the Ewald summation might act
1395 > to loosen up the local structure more than the image-charges in {\sc
1396 > sf}. Because the systems are locked in different regions of
1397 > phase-space, comparisons between properties at these temperatures are
1398 > not exactly fair. This observation is explored in more detail in
1399 > section \ref{sec:t5peDynamics}.
1400  
1401   The final thermodynamic property displayed in figure
1402   \ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy
# Line 1435 | Line 1429 | relation using the mean square displacement (MSD),
1429   self-diffusion constants ($D$) were calculated with the Einstein
1430   relation using the mean square displacement (MSD),
1431   \begin{equation}
1432 < D = \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
1432 > D = \lim_{t\rightarrow\infty}
1433 >    \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
1434   \label{eq:MSD}
1435   \end{equation}
1436   where $t$ is time, and $\mathbf{r}_i$ is the position of particle
# Line 1446 | Line 1441 | regions:
1441   \begin{enumerate}[itemsep=0pt]
1442   \item parabolic short-time ballistic motion,
1443   \item linear diffusive regime, and
1444 < \item poor statistic region at long-time.
1444 > \item a region with poor statistics.
1445   \end{enumerate}
1446   The slope from the linear region (region 2) is used to calculate $D$.
1447   \begin{figure}
# Line 1515 | Line 1510 | easier comparisons in the more relevant temperature re
1510   easier comparisons in the more relevant temperature regime.}
1511   \label{fig:t5peDynamics}
1512   \end{figure}
1513 < Results for the diffusion constants and reorientational time constants
1513 > Results for the diffusion constants and orientational relaxation times
1514   are shown in figure \ref{fig:t5peDynamics}. From this figure, it is
1515   apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using
1516   the Ewald sum are reproduced with the {\sc sf} technique. The enhanced
# Line 1524 | Line 1519 | diffuse a little faster than with the Ewald sum; howev
1519   insight into differences between the electrostatic summation
1520   techniques. With the undamped {\sc sf} technique, TIP5P-E tends to
1521   diffuse a little faster than with the Ewald sum; however, use of light
1522 < to moderate damping results in indistinguishable $D$ values. Though not
1523 < apparent in this figure, {\sc sf} values at the lowest temperature are
1524 < approximately an order of magnitude lower than with Ewald. These
1522 > to moderate damping results in indistinguishable $D$ values. Though
1523 > not apparent in this figure, {\sc sf} values at the lowest temperature
1524 > are approximately an order of magnitude lower than with Ewald. These
1525   values support the observation from section \ref{sec:t5peThermo} that
1526   there appeared to be a change to a more glassy-like phase with the
1527   {\sc sf} technique at these lower temperatures.
# Line 1540 | Line 1535 | calculation from a 10~ps trajectory with {\sc sf} with
1535   for this deviation between techniques. The Ewald results were taken
1536   shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick
1537   calculation from a 10~ps trajectory with {\sc sf} with an $\alpha$ of
1538 < 0.2~\AA$^-1$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$, placing
1539 < the result more in line with that obtained using the Ewald sum. These
1540 < results support this explanation; however, recomputing the results to
1541 < meet a poorer statistical standard is counter-productive. Assuming the
1542 < Ewald results are not the product of poor statistics, differences in
1543 < techniques to integrate the orientational motion could also play a
1544 < role. {\sc shake} is the most commonly used technique for
1545 < approximating rigid-body orientational motion,\cite{Ryckaert77} where
1546 < as in {\sc oopse}, we maintain and integrate the entire rotation
1547 < matrix using the {\sc dlm} method.\cite{Meineke05} Since {\sc shake}
1548 < is an iterative constraint technique, if the convergence tolerances
1549 < are raised for increased performance, error will accumulate in the
1550 < orientational motion. Finally, the Ewald results were calculated using
1551 < the $NVT$ ensemble, while the $NVE$ ensemble was used for {\sc sf}
1538 > 0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$,
1539 > placing the result more in line with that obtained using the Ewald
1540 > sum. These results support this explanation; however, recomputing the
1541 > results to meet a poorer statistical standard is
1542 > counter-productive. Assuming the Ewald results are not the product of
1543 > poor statistics, differences in techniques to integrate the
1544 > orientational motion could also play a role. {\sc shake} is the most
1545 > commonly used technique for approximating rigid-body orientational
1546 > motion,\cite{Ryckaert77} where as in {\sc oopse}, we maintain and
1547 > integrate the entire rotation matrix using the {\sc dlm}
1548 > method.\cite{Meineke05} Since {\sc shake} is an iterative constraint
1549 > technique, if the convergence tolerances are raised for increased
1550 > performance, error will accumulate in the orientational
1551 > motion. Finally, the Ewald results were calculated using the $NVT$
1552 > ensemble, while the $NVE$ ensemble was used for {\sc sf}
1553   calculations. The additional mode of motion due to the thermostat will
1554   alter the dynamics, resulting in differences between $NVT$ and $NVE$
1555   results. These differences are increasingly noticeable as the
# Line 1565 | Line 1561 | consider an extension of these techniques to include p
1561   neutralizing the cutoff sphere with charge-charge interaction shifting
1562   and by damping the electrostatic interactions. Now we would like to
1563   consider an extension of these techniques to include point multipole
1564 < interactions. How will the shifting and damping need to develop in
1564 > interactions. How will the shifting and damping need to be modified in
1565   order to accommodate point multipoles?
1566  
1567 < Of the two techniques, the least to vary is shifting. Shifting is
1567 > Of the two techniques, the easiest to adapt is shifting. Shifting is
1568   employed to neutralize the cutoff sphere; however, in a system
1569   composed purely of point multipoles, the cutoff sphere is already
1570   neutralized. This means that shifting is not necessary between point
# Line 1586 | Line 1582 | than considering only the interactions between single
1582   replacing $r^{-1}$ with erfc$(\alpha r)\cdot r^{-1}$ in the multipole
1583   expansion.\cite{Hirschfelder67} In the multipole expansion, rather
1584   than considering only the interactions between single point charges,
1585 < the electrostatic interactions is reformulated such that it describes
1585 > the electrostatic interaction is reformulated such that it describes
1586   the interaction between charge distributions about central sites of
1587   the respective sets of charges. This procedure is what leads to the
1588   familiar charge-dipole,

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines