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

Comparing trunk/multipole/multipole_2/multipole2.tex (file contents):
Revision 4203 by gezelter, Wed Aug 6 19:10:04 2014 UTC vs.
Revision 4207 by gezelter, Thu Aug 7 21:13:00 2014 UTC

# Line 220 | Line 220 | negaitve ions and vice versa.  The reversed-charge ima
220   limited to perfect crystals, but can also apply in disordered fluids.
221   Even at elevated temperatures, there is local charge balance in an
222   ionic liquid, where each positive ion has surroundings dominated by
223 < negaitve ions and vice versa.  The reversed-charge images on the
223 > negative ions and vice versa.  The reversed-charge images on the
224   cutoff sphere that are integral to the Wolf and DSF approaches retain
225   the effective multipolar interactions as the charges traverse the
226   cutoff boundary.
# Line 325 | Line 325 | U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1
325   expressed as the product of two multipole operators and a Coulombic
326   kernel,
327   \begin{equation}
328 < U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r}  \label{kernel}.
328 > U_{ab}(r)= M_{a} M_{b} \frac{1}{r}  \label{kernel}.
329   \end{equation}
330 < where the multipole operator for site $\bf a$, $\hat{M}_{\bf a}$, is
331 < expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf
332 <    a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object
333 < $\bf a$, etc.
330 > where the multipole operator for site $a$, $M_{a}$, is
331 > expressed in terms of the point charge, $C_{a}$, dipole, ${\bf D}_{a}$, and quadrupole, $\mathsf{Q}_{a}$, for object
332 > $a$, etc.
333  
334   % Interactions between multipoles can be expressed as higher derivatives
335   % of the bare Coulomb potential, so one way of ensuring that the forces
# Line 385 | Line 384 | radius.  For example, the direct dipole dot product
384   contributions to the potential, and ensures that the forces and
385   torques from each of these contributions will vanish at the cutoff
386   radius.  For example, the direct dipole dot product
387 < ($\mathbf{D}_{\bf a}
388 < \cdot \mathbf{D}_{\bf b}$) is treated differently than the dipole-distance
387 > ($\mathbf{D}_{a}
388 > \cdot \mathbf{D}_{b}$) is treated differently than the dipole-distance
389   dot products:
390   \begin{equation}
391 < U_{D_{\bf a}D_{\bf b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left(
392 <  \mathbf{D}_{\bf a} \cdot
393 < \mathbf{D}_{\bf b} \right) v_{21}(r) +
394 < \left( \mathbf{D}_{\bf a} \cdot \hat{r} \right)
395 < \left( \mathbf{D}_{\bf b} \cdot \hat{r} \right) v_{22}(r) \right]
391 > U_{D_{a}D_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left(
392 >  \mathbf{D}_{a} \cdot
393 > \mathbf{D}_{b} \right) v_{21}(r) +
394 > \left( \mathbf{D}_{a} \cdot \hat{\mathbf{r}} \right)
395 > \left( \mathbf{D}_{b} \cdot \hat{\mathbf{r}} \right) v_{22}(r) \right]
396   \end{equation}
397  
398   For the Taylor shifted (TSF) method with the undamped kernel,
# Line 418 | Line 417 | U_{D_{\bf a}D_{\bf b}}(r)  = U_{D_{\bf a}D_{\bf b}}(r)
417   which have been projected onto the surface of the cutoff sphere
418   without changing their relative orientation,
419   \begin{equation}
420 < U_{D_{\bf a}D_{\bf b}}(r)  = U_{D_{\bf a}D_{\bf b}}(r) -
421 < U_{D_{\bf a} D_{\bf b}}(r_c)
422 <   - (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot
423 <  \nabla U_{D_{\bf a}D_{\bf b}}(r_c).
420 > U_{D_{a}D_{b}}(r)  = U_{D_{a}D_{b}}(r) -
421 > U_{D_{a}D_{b}}(r_c)
422 >   - (r_{ab}-r_c) ~~~\hat{\mathbf{r}}_{ab} \cdot
423 >  \nabla U_{D_{a}D_{b}}(r_c).
424   \end{equation}
425 < Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{\bf
427 <  a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance
425 > Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{a}$ and $\mathbf{D}_{b}$, are retained at the cutoff distance
426   (although the signs are reversed for the dipole that has been
427   projected onto the cutoff sphere).  In many ways, this simpler
428   approach is closer in spirit to the original shifted force method, in
# Line 445 | Line 443 | and any multipolar site inside the cutoff radius is gi
443   In general, the gradient shifted potential between a central multipole
444   and any multipolar site inside the cutoff radius is given by,
445   \begin{equation}
446 <  U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
447 <    U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{\mathbf{r}}
448 <    \cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
446 > U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) -
447 > U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) - (r-r_c)
448 > \hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right]
449   \label{generic2}
450   \end{equation}
451   where the sum describes a separate force-shifting that is applied to
452   each orientational contribution to the energy.  In this expression,
453   $\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles
454 < ($a$ and $b$) in space, and $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$
454 > ($a$ and $b$) in space, and $\mathsf{A}$ and $\mathsf{B}$
455   represent the orientations the multipoles.
456  
457   The third term converges more rapidly than the first two terms as a
# Line 480 | Line 478 | U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf
478   interactions with the central multipole and the image. This
479   effectively shifts the total potential to zero at the cutoff radius,
480   \begin{equation}
481 < U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
482 < U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
481 > U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) -
482 > U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right]
483   \label{eq:SP}
484   \end{equation}          
485   where the sum describes separate potential shifting that is done for
# Line 526 | Line 524 | in the test cases are given in table~\ref{tab:pars}.
524   used the multipolar Ewald sum as a reference method for comparing
525   energies, forces, and torques for molecular models that mimic
526   disordered and ordered condensed-phase systems.  The parameters used
527 < in the test cases are given in table~\ref{tab:pars}.
527 > in the test cases are given in table~\ref{tab:pars}.
528  
529   \begin{table}
530   \label{tab:pars}
# Line 596 | Line 594 | approximations.\cite{Smith82,Smith98} In all cases, th
594   with a reciprocal space cutoff, $k_\mathrm{max} = 7$.  Our version of
595   the Ewald sum is a re-implementation of the algorithm originally
596   proposed by Smith that does not use the particle mesh or smoothing
597 < approximations.\cite{Smith82,Smith98} In all cases, the quantities
598 < being compared are the electrostatic contributions to energies, force,
599 < and torques.  All other contributions to these quantities (i.e. from
600 < Lennard-Jones interactions) are removed prior to the comparisons.
597 > approximations.\cite{Smith82,Smith98} This implementation was tested
598 > extensively against the analytic energy constants for the multipolar
599 > lattices that are discussed in reference \onlinecite{PaperI}.  In all
600 > cases discussed below, the quantities being compared are the
601 > electrostatic contributions to energies, force, and torques.  All
602 > other contributions to these quantities (i.e. from Lennard-Jones
603 > interactions) are removed prior to the comparisons.
604  
605   The convergence parameter ($\alpha$) also plays a role in the balance
606   of the real-space and reciprocal-space portions of the Ewald
# Line 958 | Line 959 | temperature of 300K.  After equilibration, this liquid
959   in this series and provides the most comprehensive test of the new
960   methods.  A liquid-phase system was created with 2000 water molecules
961   and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a
962 < temperature of 300K.  After equilibration, this liquid-phase system
963 < was run for 1 ns under the Ewald, Hard, SP, GSF, and TSF methods with
964 < a cutoff radius of 12\AA.  The value of the damping coefficient was
965 < also varied from the undamped case ($\alpha = 0$) to a heavily damped
966 < case ($\alpha = 0.3$ \AA$^{-1}$) for all of the real space methods.  A
967 < sample was also run using the multipolar Ewald sum with the same
968 < real-space cutoff.
962 > temperature of 300K.  After equilibration in the canonical (NVT)
963 > ensemble using a Nos\'e-Hoover thermostat, this liquid-phase system
964 > was run for 1 ns in the microcanonical (NVE) ensemble under the Ewald,
965 > Hard, SP, GSF, and TSF methods with a cutoff radius of 12\AA.  The
966 > value of the damping coefficient was also varied from the undamped
967 > case ($\alpha = 0$) to a heavily damped case ($\alpha = 0.3$
968 > \AA$^{-1}$) for all of the real space methods.  A sample was also run
969 > using the multipolar Ewald sum with the same real-space cutoff.
970  
971   In figure~\ref{fig:energyDrift} we show the both the linear drift in
972   energy over time, $\delta E_1$, and the standard deviation of energy
# Line 988 | Line 990 | $k$-space cutoff values.
990    \includegraphics[width=\textwidth]{newDrift_12.eps}
991   \label{fig:energyDrift}        
992   \caption{Analysis of the energy conservation of the real-space
993 <  electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in
994 <  energy over time (in kcal / mol / particle / ns) and $\delta
995 <  \mathrm{E}_0$ is the standard deviation of energy fluctuations
996 <  around this drift (in kcal / mol / particle).  All simulations were
997 <  of a 2000-molecule simulation of SSDQ water with 48 ionic charges at
993 >  methods. $\delta \mathrm{E}_1$ is the linear drift in energy over
994 >  time (in kcal / mol / particle / ns) and $\delta \mathrm{E}_0$ is
995 >  the standard deviation of energy fluctuations around this drift (in
996 >  kcal / mol / particle).  Points that appear below the dashed grey
997 >  (Ewald) lines exhibit better energy conservation than commonly-used
998 >  parameters for Ewald-based electrostatics.  All simulations were of
999 >  a 2000-molecule simulation of SSDQ water with 48 ionic charges at
1000    300 K starting from the same initial configuration. All runs
1001    utilized the same real-space cutoff, $r_c = 12$\AA.}
1002   \end{figure}
1003  
1004 < \subsection{Reproduction of Structural Features\label{sec:structure}}
1005 < One of the best tests of modified interaction potentials is the
1006 < fidelity with which they can reproduce structural features in a
1007 < liquid.  One commonly-utilized measure of structural ordering is the
1008 < pair distribution function, $g(r)$, which measures local density
1009 < deviations in relation to the bulk density.  In the electrostatic
1010 < approaches studied here, the short-range repulsion from the
1011 < Lennard-Jones potential is identical for the various electrostatic
1012 < methods, and since short range repulsion determines much of the local
1013 < liquid ordering, one would not expect to see any differences in
1014 < $g(r)$.  Indeed, the pair distributions are essentially identical for
1015 < all of the electrostatic methods studied (for each of the different
1016 < systems under investigation).  Interested readers may consult the
1017 < supplementary information for plots of these pair distribution
1018 < functions.
1004 > \subsection{Reproduction of Structural \& Dynamical Features\label{sec:structure}}
1005 > The most important test of the modified interaction potentials is the
1006 > fidelity with which they can reproduce structural features and
1007 > dynamical properties in a liquid.  One commonly-utilized measure of
1008 > structural ordering is the pair distribution function, $g(r)$, which
1009 > measures local density deviations in relation to the bulk density.  In
1010 > the electrostatic approaches studied here, the short-range repulsion
1011 > from the Lennard-Jones potential is identical for the various
1012 > electrostatic methods, and since short range repulsion determines much
1013 > of the local liquid ordering, one would not expect to see many
1014 > differences in $g(r)$.  Indeed, the pair distributions are essentially
1015 > identical for all of the electrostatic methods studied (for each of
1016 > the different systems under investigation).  An example of this
1017 > agreement for the SSDQ water/ion system is shown in
1018 > Fig. \ref{fig:gofr}.
1019  
1016 A direct measure of the structural features that is a more
1017 enlightening test of the modified electrostatic methods is the average
1018 value of the electrostatic energy $\langle U_\mathrm{elect} \rangle$
1019 which is obtained by sampling the liquid-state configurations
1020 experienced by a liquid evolving entirely under the influence of the
1021 methods being investigated.  In figure \ref{fig:Uelect} we show how
1022 $\langle U_\mathrm{elect} \rangle$ for varies with the damping parameter,
1023 $\alpha$, for each of the methods.
1024
1020   \begin{figure}
1021    \centering
1022 <  \includegraphics[width=\textwidth]{averagePotentialEnergy_r9_12.eps}
1023 < \label{fig:Uelect}        
1024 < \caption{The average electrostatic potential energy,
1025 <  $\langle U_\mathrm{elect} \rangle$ for the SSDQ water with ions as a function
1026 <  of the damping parameter, $\alpha$, for each of the real-space
1027 <  electrostatic methods. Top panel: simulations run with a real-space
1033 <  cutoff, $r_c = 9$\AA.  Bottom panel: the same quantity, but with a
1034 <  larger cutoff, $r_c = 12$\AA.}
1022 >  \includegraphics[width=\textwidth]{gofr_ssdqc.eps}
1023 > \label{fig:gofr}        
1024 > \caption{The pair distribution functions, $g(r)$, for the SSDQ
1025 >  water/ion system obtained using the different real-space methods are
1026 >  essentially identical with the result from the Ewald
1027 >  treatment.}
1028   \end{figure}
1029  
1030 < It is clear that moderate damping is important for converging the mean
1031 < potential energy values, particularly for the two shifted force
1032 < methods (GSF and TSF).  A value of $\alpha \approx 0.18$ \AA$^{-1}$ is
1033 < sufficient to converge the SP and GSF energies with a cutoff of 12
1034 < \AA, while shorter cutoffs require more dramatic damping ($\alpha
1042 < \approx 0.36$ \AA$^{-1}$ for $r_c = 9$ \AA).  It is also clear from
1043 < fig. \ref{fig:Uelect} that it is possible to overdamp the real-space
1044 < electrostatic methods, causing the estimate of the energy to drop
1045 < below the Ewald results.
1030 > There is a very slight overstructuring of the first solvation shell
1031 > when using when using TSF at lower values of the damping coefficient
1032 > ($\alpha \le 0.1$) or when using undamped GSF.  With moderate damping,
1033 > GSF and SP produce pair distributions that are identical (within
1034 > numerical noise) to their Ewald counterparts.
1035  
1036 + A structural property that is a more demanding test of modified
1037 + electrostatics is the mean value of the electrostatic energy $\langle
1038 + U_\mathrm{elect} \rangle / N$ which is obtained by sampling the
1039 + liquid-state configurations experienced by a liquid evolving entirely
1040 + under the influence of each of the methods.  In table \ref{tab:Props}
1041 + we demonstrate how $\langle U_\mathrm{elect} \rangle / N$ varies with
1042 + the damping parameter, $\alpha$, for each of the methods.
1043 +
1044 + As in the crystals studied in the first paper, damping is important
1045 + for converging the mean electrostatic energy values, particularly for
1046 + the two shifted force methods (GSF and TSF).  A value of $\alpha
1047 + \approx 0.2$ \AA$^{-1}$ is sufficient to converge the SP and GSF
1048 + energies with a cutoff of 12 \AA, while shorter cutoffs require more
1049 + dramatic damping ($\alpha \approx 0.3$ \AA$^{-1}$ for $r_c = 9$ \AA).
1050 + Overdamping the real-space electrostatic methods occurs with $\alpha >
1051 + 0.4$, causing the estimate of the energy to drop below the Ewald
1052 + results.
1053 +
1054   These ``optimal'' values of the damping coefficient are slightly
1055   larger than what were observed for DSF electrostatics for purely
1056   point-charge systems, although a value of $\alpha=0.18$ \AA$^{-1}$ for
1057   $r_c = 12$\AA appears to be an excellent compromise for mixed charge
1058   multipole systems.
1059  
1053 \subsection{Reproduction of Dynamic Properties\label{sec:structure}}
1060   To test the fidelity of the electrostatic methods at reproducing
1061   dynamics in a multipolar liquid, it is also useful to look at
1062   transport properties, particularly the diffusion constant,
# Line 1060 | Line 1066 | the multipoles.  For the soft dipolar fluid, and the S
1066   \label{eq:diff}
1067   \end{equation}
1068   which measures long-time behavior and is sensitive to the forces on
1069 < the multipoles.  For the soft dipolar fluid, and the SSDQ liquid
1069 > the multipoles.  For the soft dipolar fluid and the SSDQ liquid
1070   systems, the self-diffusion constants (D) were calculated from linear
1071 < fits to the long-time portion of the mean square displacement
1072 < ($\langle r^{2}(t) \rangle$).\cite{Allen87}
1071 > fits to the long-time portion of the mean square displacement,
1072 > $\langle r^{2}(t) \rangle$.\cite{Allen87}
1073  
1074   In addition to translational diffusion, orientational relaxation times
1075   were calculated for comparisons with the Ewald simulations and with
1076 < experiments. These values were determined from the same 1~ns $NVE$
1077 < trajectories used for translational diffusion by calculating the
1078 < orientational time correlation function,
1076 > experiments. These values were determined from the same 1~ns
1077 > microcanonical trajectories used for translational diffusion by
1078 > calculating the orientational time correlation function,
1079   \begin{equation}
1080 < C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\gamma(t)
1081 <                \cdot\hat{\mathbf{u}}_i^\gamma(0)\right]\right\rangle,
1080 > C_l^\gamma(t) = \left\langle P_l\left[\hat{\mathbf{A}}_\gamma(t)
1081 >                \cdot\hat{\mathbf{A}}_\gamma(0)\right]\right\rangle,
1082   \label{eq:OrientCorr}
1083   \end{equation}
1084   where $P_l$ is the Legendre polynomial of order $l$ and
1085 < $\hat{\mathbf{u}}_i^\alpha$ is the unit vector of molecule $i$ along
1086 < axis $\gamma$.  The body-fixed reference frame used for our
1087 < orientational correlation functions has the $z$-axis running along the
1088 < dipoles, and for the SSDQ water model, the $y$-axis connects the two
1089 < implied hydrogen atoms.
1085 > $\hat{\mathbf{A}}_\gamma$ is the space-frame unit vector for body axis
1086 > $\gamma$ on a molecule..  Th body-fixed reference frame used for our
1087 > models has the $z$-axis running along the dipoles, and for the SSDQ
1088 > water model, the $y$-axis connects the two implied hydrogen atom
1089 > positions.  From the orientation autocorrelation functions, we can
1090 > obtain time constants for rotational relaxation either by fitting an
1091 > exponential function or by integrating the entire correlation
1092 > function.  In a good water model, these decay times would be
1093 > comparable to water orientational relaxation times from nuclear
1094 > magnetic resonance (NMR). The relaxation constant obtained from
1095 > $C_2^y(t)$ is normally of experimental interest because it describes
1096 > the relaxation of the principle axis connecting the hydrogen
1097 > atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion
1098 > of the dipole-dipole relaxation from a proton NMR signal and should
1099 > provide an estimate of the NMR relaxation time constant.\cite{Impey82}
1100  
1085 From the orientation autocorrelation functions, we can obtain time
1086 constants for rotational relaxation either by fitting an exponential
1087 function or by integrating the entire correlation function.  These
1088 decay times are directly comparable to water orientational relaxation
1089 times from nuclear magnetic resonance (NMR). The relaxation constant
1090 obtained from $C_2^y(t)$ is normally of experimental interest because
1091 it describes the relaxation of the principle axis connecting the
1092 hydrogen atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular
1093 portion of the dipole-dipole relaxation from a proton NMR signal and
1094 should provide an estimate of the NMR relaxation time
1095 constant.\cite{Impey82}
1096
1101   Results for the diffusion constants and orientational relaxation times
1102 < are shown in figure \ref{fig:dynamics}. From this data, it is apparent
1102 > are shown in figure \ref{tab:Props}. From this data, it is apparent
1103   that the values for both $D$ and $\tau_2$ using the Ewald sum are
1104 < reproduced with high fidelity by the GSF method.
1104 > reproduced with reasonable fidelity by the GSF method.
1105  
1106 < The $\tau_2$ results in \ref{fig:dynamics} show a much greater
1107 < difference between the real-space and the Ewald results.
1106 > The $\tau_2$ results in \ref{tab:Props} show a much greater difference
1107 > between the real-space and the Ewald results.
1108  
1109 + \begin{table}
1110 + \label{tab:Props}
1111 + \caption{Comparison of the structural and dynamic properties for the
1112 +  soft dipolar liquid test for all of the real-space methods.}
1113 + \begin{tabular}{l|c|cccc|cccc|cccc}
1114 +         & Ewald & \multicolumn{4}{c|}{SP} & \multicolumn{4}{c|}{GSF} & \multicolumn{4}{c|}{TSF} \\
1115 + $\alpha$ (\AA$^{-1}$) & &      
1116 + 0.0 & 0.1 & 0.2 & 0.3 &
1117 + 0.0 & 0.1 & 0.2 & 0.3 &
1118 + 0.0 & 0.1 & 0.2 & 0.3 \\ \cline{2-6}\cline{6-10}\cline{10-14}
1119  
1120 + $\langle U_\mathrm{elect} \rangle /N$ &&&&&&&&&&&&&\\
1121 + D ($10^{-4}~\mathrm{cm}^2/\mathrm{s}$)&
1122 + 470.2(6) &
1123 + 416.6(5) &
1124 + 379.6(5) &
1125 + 438.6(5) &
1126 + 476.0(6) &
1127 + 412.8(5) &
1128 + 421.1(5) &
1129 + 400.5(5) &
1130 + 437.5(6) &
1131 + 434.6(5) &
1132 + 411.4(5) &
1133 + 545.3(7) &
1134 + 459.6(6) \\
1135 + $\tau_2$ (fs) &
1136 + 1.136 &
1137 + 1.041 &
1138 + 1.064 &
1139 + 1.109 &
1140 + 1.211 &
1141 + 1.119 &
1142 + 1.039 &
1143 + 1.058 &
1144 + 1.21  &
1145 + 1.15  &
1146 + 1.172 &
1147 + 1.153 &
1148 + 1.125 \\
1149 + \end{tabular}
1150 + \end{table}
1151 +
1152 +
1153   \section{CONCLUSION}
1154   In the first paper in this series, we generalized the
1155   charge-neutralized electrostatic energy originally developed by Wolf

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines