Computing the dielectric constant of liquid water at constant dielectric displacement Chao Zhang∗ and Michiel Sprik Department of Chemistry, University of Cambridge, Lensfield Rd, Cambridge CB2 1EW, United Kingdom (Dated: March 17, 2016) The static dielectric constant of liquid water is computed using classical force field based molec- ular dynamics simulation at fixed electric displacement D. The method to constrain the electric displacement is the finite temperature classical variant of the constant-D method developed by Stengel, Spaldin and Vanderbilt [Nat. Phys. 5, 304, (2009)]. There is also a modification of this scheme imposing fixed values of the macroscopic field E. The method is applied to the popular SPC/E model of liquid water. We compare four different estimates of the dielectric constant, two obtained from fluctuations of the polarization at D = 0 and E = 0 and two from the variation of polarization with finite D and E. It is found that all four estimates agree when properly converged. The computational effort to achieve convergence varies however, with constant D calculations being substantially more efficient. We attribute this difference to the much shorter relaxation time of longitudinal polarization compared to transverse polarization accelerating constant D calculations. I. INTRODUCTION The static dielectric constant of model polar liquids is usually computed from polarization fluctuations applying the linear response relations of Kirkwood-Fröhlich the- ory [1, 2]. These calculations are carried out in periodic molecular dynamics (MD) cells treating long range elec- trostatic interactions using Ewald summation[3]. This is an expensive calculation. Indeed, it is commonly ac- knowledged that for highly polarizable liquids such as water simulations at time-scales of nanoseconds are nec- essary to converge the fluctuation estimate for the dielec- tric constant.[3–26]. These time scales are accessible for classical force field based MD simulation. This is, how- ever, a major challenge for MD simulation with forces calculated “on-the-fly” using electronic structure calcu- lation methods, such as density functional theory [27–29]. The high costs of the Kirkwood-Fröhlich scheme is a strong motivation for the development for more efficient alternatives. Boundary conditions can have a drastic ef- fect on polarization fluctuations which has led to the search for optimal boundary conditions[7, 9, 19, 26]. The use of finite field methods is another option that has been investigated[26, 30, 31]. The rationale here is, of course, that converging an estimate of polarization should be quicker than converging its fluctuations. However it was found, in particular in the case of water, that the re- sponse of the polarization is non-linear for already mod- erate field strength requiring a careful extrapolation to zero field. As a result finite field calculations of the di- electric constant are in practice not that much cheaper than the computation from polarization fluctuations at zero field. A change of boundary conditions not only affects the size of polarization fluctuations but also the time scale. It has been shown that the standard Ewald summation ∗ cz302@cam.ac.uk method corresponds to constraining the average macro- scopic field to zero[3, 6]. The static dielectric constant under these conditions is dominated by fluctuations of the transverse polarization[32–34]. These are the slow modes. The relaxation time τL of longitudinal modes (in the k = 0 limit) is significantly faster compared to the relaxation time τT of transverse modes. The ratio according to Debye theory[2] is τT /τL = ε0/ε∞. For non- polarizable SPC models of liquid water (ε∞ = 1) this amounts almost to two orders of magnitude. The much faster relaxation of longitudinal fluctuations raises the question whether this can be exploited to ac- celerate the calculation of the static dielectric constant. In this paper we show that this can be achieved by chang- ing the boundary conditions for Ewald summation from zero macroscopic field (E = 0) to zero dielectric dis- placement (D = 0). The method to compute the to- tal energy of periodic supercells under fixed D has been developed by Vanderbilt and coworkers. The key refer- ence to this approach is the 2009 Nature Physics paper by Stengel, Spaldin and Vanderbilt[35] which will be re- ferred to as SSV (see also Ref. 36). The method is a recent spin-off of the modern theory of polarization de- veloped by Vanderbilt and Resta during the 90’s[37–39]. The modern theory of polarization caused a revolution in theoretical and computational solid state physics mak- ing it possible, for the first time, to investigate the elec- tric equation of state of ferroelectric systems. The ini- tial approach was to compute the total energy for fixed values of the polarization and to determine the electric field from the derivative[40, 41]. This was subsequently changed to a scheme using directly the macroscopic elec- tric field E or the electric displacement field D as the control variable, which has both computational and con- ceptual advantages[35, 36]. The SSV constant D and the related finite E method are easy to implement in a classical force field code. The method retains the regular “tin-foil” Ewald sum for the calculation of electrostatic energy extending it with an electric term which depends on the polarization P and 2 contains E or D as a parameter. The result is an ex- tended hamiltonian replacing the original hamiltonian in the MD simulation. The present paper is a feasibil- ity and validation study of this approach for a SPC/E model of liquid water [42]. We verify that the dielectric constant obtained from polarization fluctuations under D = 0 conditions agrees with the value estimated from standard E = 0 calculations. These results are then also compared to dielectric constant estimates computed from the change of the expectation value of P with applied E and D using the SSV hamiltonian in finite field mode. The approach in the present paper has multiple par- allels to similar work that has appeared in the litera- ture. Most of this work also involves macroscopic po- larization dependent energy terms extending the micro- scopic Coulomb interaction energy evaluated for an infi- nite periodic lattice of supercells using the Ewald summa- tion method[43]. These extensions are known as surface terms[3, 4, 44–50] or reaction fields[6–8] and are of a form similar to the D = 0 limit of the SSV constant D polar- ization coupling term. What is unique about SSV the- ory is that the extended Ewald Hamiltonians are derived strictly complying with the rules of dielectric thermody- namics as set out by Landau and Lifshitz[51]. The focus on thermodynamics has certain advantages as already pointed by Aragones et al.[52]. The paper starts there- fore with a fairly detailed outline of the finite temper- ature classical variant of the SSV approach (sections II and III) supplemented with three appendices with more formal theoretical considerations. Results are presented and discussed in section IV. We conclude with a sum- mary and outlook for future applications to interfaces. II. FINITE E AND D IN EXTENDED SYSTEMS A. Constant E and D hamiltonians and thermodynamics The theory behind the finite field method developed by SSV is summarized in the supporting information of Ref. 35. The central quantity is the electric enthalpy functional. The electric enthalpy of a system of volume Ω is written as F (E, v) = EKS(v)− Ω E ·P(v) (1) EKS(v) is the Kohn-Sham total energy with v denoting all microscopic degrees of freedom involved, ie the or- bital coefficients specifying the one-electron orbitals and the positions of the ions. EKS(v) is obtained for given v using the regular reciprocal space methods of compu- tational solid state physics explicitly excluding all k = 0 contributions. P(v) is the macroscopic polarization den- sity for the microscopic state specified by v. To compute the expectation value P of the polarization density the electronic structure v of the system is deter- mined by minimizing the electric enthalpy functional for fixed E F (E) = min v F (E, v)) = min v [EKS(v)− Ω E ·P(v)] (2) Taking the derivative gives the polarization dF dE = −ΩP (3) Eq. 3 can be regarded as the electric equation of state for a uniform insulator. Vanderbilts electric enthalpy scheme can be readily adapted to classical force field based MD. The KS total energy EKS(v) in Eq. 1 is replaced by the Hamiltonian H(v) of the system where v is now the set of momenta and positions of the particles. F (E, v) = HPBC(v)− Ω E ·P(v) (4) We have appended a subscript PBC to the Hamiltonian as a reminder that the electrostatic energies and forces are computed using standard Ewald summation as ap- plied also for the computation of the KS energy in Eq. 1. HPBC(v) can be formally written as the sum of a term Hsr describing the short range interactions and the re- ciprocal space representation of the electrostatic energy HPBC(v) = Hsr(v) + 2πΩ ∑ k6=0 ρ(k)2 k2 (5) where ρ(k) is the Fourier transform of the atomic charge distribution (the SPC charges). Ewald summation is a computationally more efficient method to calculate this energy carrying out the summation partially in real space. The result corresponds to zero average electric field and potential, the so called tinfoil boundary con- ditions (no surface terms)[3]. The finite electric field is introduced as a parameter in the second term of Eq. 4. The equivalent of the electric enthalpy of Eq. 2 is the free energy of the ensemble generated by the extended Hamiltonian F (E, v) of Eq. 4 F (E) = −kBT lnZE (6) ZE is the electric field dependent partition function ZE = ∫ dpNdrN exp[−β ( HPBC − Ω E ·P(rN ) ) ] (7) where HPBC is the Hamiltonian of the periodic MD sys- tem in Eq. 4. The coordinate and momentum arguments v = rN ,pN were suppressed. β = 1/kBT is the inverse temperature. The combinatorial prefactor 1/(h3NN !) has been omitted. The derivative of F (E) of Eq. 6 again gives the polarization according to Eq. 3. Recently, Stengel, Spaldin and Vanderbilt have modi- fied the constant E to a constant D method [35]. They introduced a new functional, the electric internal energy functional U (D, v). Transposed to classical MD, this functional is written as U (D, v) = HPBC(v) + Ω 8π (D− 4πP(v)) 2 (8) 3 The corresponding D dependent electric internal energy is again obtained from the partition function U (D) = −kBT lnZD (9) with ZD = ∫ dpNdrN exp [−β U(D, v)] (10) where as before v = rN ,pN . Note that U(D) is still a (Helmholtz) free energy with respect to temperature, similar to F (E). Evaluating the derivative with respect the control variable, D in this case, we recover the macro- scopic field dU dD = Ω 4π (D− 4πP) = Ω 4π E (11) The second identity follows from D = E + 4πP (12) which is the fundamental relation of Maxwell theory defining the dielectric displacement. That U of Eq. 9 is indeed the electric free energy be- comes evident when Eq. 11 is substituted in the Maxwell field expression for electrical work (see Landau and Lifshitz[51]) dW = Ω 4π E · dD = dU (13) The link to electrical work established in Eq. 13 is cru- cial. It is the ultimate justification for identifying E in the electrical enthalpy of Eq. 4 with the macroscopic field. As shown by SSV, the argument can be given a more for- mal thermodynamic basis in a Legendre transform frame- work. The derivation is repeated in Appendix A. B. Parallel plate capacitor and hybrid boundary conditions The SSV Hamiltonians of Eqs. 4 and 8 can be given more physical meaning when interpreted as a model of a macroscopic parallel plate capacitor[35, 36]. Fig. 1 shows a schematic picture of such a device as used in textbooks (see in particular Purcell[53]). We will assume that the normal to the electrodes is directed along the x axis. The electrodes are a distance l apart. The charge density on the left electrode is σm. Without dielectric material be- tween the plates the electric field E0 generated by this charge density is E0 = 4πσm corresponding to a poten- tial V0 = −E0l = −4πσml. In the presence of dielectric material the applied electric field E0 is screened by the induced polarization Px. The resulting macroscopic elec- tric field E can be written as E = 4π (σm + σp) = E0 − 4πPx (14) FIG. 1. Parallel plate capacitor at a) constant electric field Ex and b) at constant electric displacement Dx. σm is the surface charge density on the metal electrode. σp is the polarization surface charge density of the dielectric material (σm > −σp > 0 in the picture). The electric field E is determined by the net interface charge σm + σp (see text). where σp = −Px is the polarization charge density accu- mulating on the surface of the dielectric (see Fig. 1). Be- cause the polarization is aligned along the applied field, E0 and Px have the same sign (positive in the figure). σp and σm have opposite sign. The macroscopic electric field is determined by the net interface charge σm + σp and therefore E < E0. Similarly the potential V = −El is lower (in absolute value) than V0. This is how capaci- tors store charge[53]. Setting Ex = E,Ey = Ez = 0 in Eq. 4 we obtain the electric enthalpy hamiltonian Fx (E, v) = HPBC(v)− ΩEPx(v) (15) How can this relatively simple MD model without in- terfaces possibly represent the capacitor of Fig. 1? The idea is that interactions not affected by surface effects are described by the supercell hamiltonian HPBC. For a macroscopic capacitor these interactions are assumed to include everything except the coupling to the elec- tric field which is accounted for in the −E · P = EPx term. This term plays the role of a pair of “virtual” elec- trodes connected to a voltage source imposing a potential drop of ∆V = −EL over the length L of the MD cell in the x direction. Fx (E, v) is a microscopic Hamiltonian, not a free energy. All quantities, except Ex = E fluc- tuate in time. This applies to the polarization charge σp(t) = −Px(v(t)) but also to the charge on the virtual electrode σm(t). It is only the sum σm(t)+σp(t) = E/4π that is constant. The instantaneous electrode charge compensating the varying polarization charge is supplied by the voltage source. The implication is that for E = 0 Eq. 15 can be viewed as a capacitor in short circuit, con- sistent with the accepted view of the Ewald summation method. Next we introduce the displacement field. For the par- allel plate capacitor Dx = E0 = 4πσm. As explained above, Dx is interpreted as a property of the microscopic system of Eq. 15 and fluctuates in time because the elec- 4 tric field Ex = E is fixed. Applying the formalism of section II A, Dx, in turn, can be constrained to a value D in the dynamics driven by the Hamiltonian Ux (D, v) = HPBC(v) + Ω 8π (D − 4πPx(v)) 2 (16) The constraint on the displacement field applies only to the x direction. Eq. 16 is not a special case of Eq. 8. In the y and z directions the regular Ewald boundary con- ditions are maintained and hence Ey = 0 and Ez = 0. Because D is effectively the charge density on the virtual electrodes (D = 4πσm), fixing D in Eq. 16 corresponds to simulating an open circuit capacitor. Now it is the conju- gate variable Ex(t) that fluctuates. The time dependence of Ex is passed on to the potential ∆V = −ExL across the cell. As pointed by SSV Eq. 16 can be regarded as a hybrid form of Eq.8, obtained by a partial Legendre transform of the Hamiltonian of Eq. 4 (see further Ref. 35). In a plate capacitor D = E0 on average. However, it would be wrong to identify D(t) with E0 at every instant of time. The applied field, in the orientation of Fig. 1, is strictly along the x axis. While on average the y and z com- ponents of D vanish, instantaneous values can be finite. Dy(t) and Dz(t) are equal to the transverse polarization Py(t) and Pz(t)[34]. In fact, compared to the longitu- dinal polarization transverse polarization (Px(t) in our model capacitor) shows substantially larger fluctuations (see section III A). III. DIELECTRIC CONSTANT A. Dielectric constant from polarization fluctuations The constant E and D ensembles introduced in sec- tion II apply to different electrical boundary conditions. Fluctuations of the polarization P are therefore ex- pected to differ in magnitude and may occur on different timescales. However, following Kirkwood-Fröhlich the- ory [1, 2], it should be possible to obtain an estimate of the dielectric constant from polarization fluctuations ei- ther under constant E or D dynamics. Staying with the capacitor paradigm of section II B application of linear response to the system defined by the electrical enthalpy hamiltonian Eq 15 gives 〈Px〉 = βΩ ( 〈P 2 x 〉 − 〈Px〉2 ) E (17) where the second moment is evaluated at zero field (E = 0). The liquid is isotropic and we can set 〈P 2 z 〉− 〈Pz〉2 = (〈P2〉 − 〈P〉2)/3, which gains us some accuracy in the statistics. The response coefficient of Eq. 17 is the sus- ceptibility χ. Converting to the dielectric constant using 4πχ = ε− 1 we find ε = 1 + 4πβΩ 3 ( 〈P2〉E=0 − 〈P〉2E=0 ) (18) For clarity the thermal average brackets in Eq. 17 have been marked with a subscript indicating the condition under which the fluctuations have been obtained. Eq. 18 is identical to the standard fluctuation formula used in Ewald summation[3, 4, 8, 9]. We arrived at this established result without having to worry about how to relate the macroscopic field E to the applied field E0[7, 8, 33]. The electric field in the SSV extended Hamiltonian is directly equal to the macroscopic field. This is a key feature of the SSV scheme. We will return to this important point once more in section III C where we make a comparison to the methods commonly used in computational physical chemistry. Alternatively we can average Px over the ensemble de- fined by the open circuit Hamiltonian of Eq. 16. It is not difficult to show that linear response now leads to the relation 〈Px〉 = βΩ ( 〈P 2 x 〉 − 〈Px〉2 ) D (19) In Eq. 19 we recognize the definition of the polarizability α. Eq. 19 can be exploited to obtain another estimate of the dielectric constant via the relation 4πα = 1− 1/ε which should be consistent with the estimate from the susceptibility χ (Eq. 17). Similar to Eq. 18 we would like to retain the extra boost in accuracy provided by isotropy (which is even more critical here, see section IV). However, there is a complication. Eq. 19 resembles Eq. 17, but in contrast to a short circuit capacitor, the open circuit system is not isotropic (see the discussion in in section II B). Fluctu- ations in the y and z (transverse) direction are distinct from the longitudinal fluctuations in Px. Isotropy can be restored by imposing a D = 0 constraint also in the y and z direction which amounts to using the hamiltonian of Eq. 8 for D = 0. The corresponding estimate for ε obtained from the polarization fluctuations is written as ε = 1 1− 4πβΩ(〈P2〉D=0 − 〈P〉2D=0)/3 (20) However, it is not immediately clear what D = 0 ensem- ble represents. This is indeed an important question for the understanding of the SSV method and we will come back to it in section III C and appendices B and C. Finally, rearranging Eq. 18 and Eq. 20 leads to: ε = 〈P2〉E=0 − 〈P〉2E=0 〈P2〉D=0 − 〈P〉2D=0 (21) As shown in the Madden and Kivelson’s review [33](see also Ref. 34) polarization fluctuations are anisotropic in the k→ 0 limit, even if the dielectric tensor is isotropic. Transverse and longitudinal fluctuations differ by a factor ε (for non-polarizable polar molecules). The same ratio is found in Eq. 21 supporting our hypothesis (our argument is not a proof) that the fluctuations sampled at constant E and constant D can be identified with the k = 0 limit of transverse respectively longitudinal polarization. We 5 furthermore note that a relation similar to Eq. 21 can be derived in the framework of reaction field methods[6, 8]. This is discussed in Appendix C. B. Dielectric constant from finite field derivatives Simulations at finite field should enable us to deter- mine ε directly from the field derivative of polarization, which we then can compare to the fluctuation estimates at zero field. For sufficiently small fields these results should in principle agree but will differ in practice be- cause of different requirements on the accuracy of the sampling and possible finite system size effects. An as- sessment of these effects is the main purpose of our inves- tigation. Finite field calculations are also of interest as a test of the SSV scheme for additional reasons. While fluctuations vary with the electrical boundary conditions, the expectation value of P is a thermodynamic state vari- able which should, in the thermodynamic limit, be the same for a given thermodynamic state, whether obtained under constant E or constant D. For a sufficiently small electric field the dielectric con- stant can be estimated using the relation ε = 1 + 4π〈Px〉 E (22) with 〈Px〉 the expectation value of polarization obtained from a MD run using the constant Ex = E Hamiltonian of Eq. 15. A similar equation, valid in the linear regime, was already used to extrapolate to zero-field by Yeh and Berkowitz in their finite field calculation of ε of liquid water [30]. The corresponding finite D estimate for ε follows from an inverse relation ε = 1 1− 4π〈Px〉/D (23) where the expectation value of polarization is determined from an average over a trajectory generated by the con- stant Dx = E0 = D Hamiltonian of Eq. 16. C. Constant applied electric field E0 The form of the SSV electric enthalpy hamiltonian of Eq. 4 might be, at first, a surprise for readers familiar with the physical chemistry literature on polar liquids expecting to see the applied electric field E0 in the cou- pling term (see for example Refs. 8, 9 and 33). Eq. 4 is however a microscopic electric enthalpy intended for constant macroscopic field E, not constant applied elec- tric field E0. The difference Ep = E − E0 is the po- larization field, also referred to as the depolarising field in literature. Ep is the electric field generated by the polarization. Substituting E = E0 + Ep in Eq. 4 gives F (E, v) = HPBC(v)− Ω (E0(v) + Ep(v)) ·P(v) (24) Similar to the polarization, E0 in Eq. 24 depends on the microstate v. It is not constant all. The parallel plate capacitor of section II B is again the best example to understand Eq. 24. In this geometry Ep = −4πPx = 4πσp (see Eq. 14). The depolarizing field is determined by the polarization surface charge σp, while the applied field E0 = 4πσm is generated by the electrode charge σm. To keep the voltage constant, the voltage source adds or removes electrode charge compensating for the fluctuation in the polarization charge. Formu- lated in terms of fields, the applied electric field responds instantaneously to the fluctuations in the depolarizing field such that the sum E = E0−4πPx = constant. Sub- stituting in the enthalpy Hamiltonian for the capacitor (Eq. 15) we obtain Fx (E, v) = HPBC(v)− Ω(E0(v)− 4πPx(v))Px(v) (25) Even for a shortcircuited capacitor (V = E = 0), the applied field, while zero on average, will have finite in- stantaneous values, cancelling the fluctuations in Ep. To control E0 we must disconnect the voltage source (open circuit). The charge on the metal electrode is now fixed and therefore E0. The Hamiltonian to use for constant E0 is therefore not Eq. 15 but Eq. 16. Can Eq. 16 be rewritten in a form more recognizable to physical chemists? To answer this question we expand the coupling term (omitting volume for simplicity) 1 8π (D − 4πPx(v)) 2 = E2 0 8π − E0Px(v) + 2πPx(v)2 (26) where we have used that for the parallel plate capaci- tor D = E0. The first term of Eq. 26 gives the energy of the polarizing field. The second term is the sought for coupling of the polarization to the external field. To identify the last term we introduce the depolarizing field Ep = −4πPx to find that the third term is equal to −EpPx/2, the work done against the depolarizing field in the process of polarization. This energy has been shown to the equal to the surface term in Ewald summation carried out over layers[45, 50] instead of the more pop- ular summation over spherical shells[3]. The key point here is that the interaction of polarization with its own electric field is not included in the Ewald Hamiltonian HPBC. The coupling to the depolarizing field must be accounted for explicitly as part of the extended hamilto- nian. In conclusion Eq. 16 seems indeed equivalent to the hamiltonians used in theory of a polar liquids[33, 34] if we assume that the interaction is implicit in the Hamil- tonian. There are further similarities linking Eq. 26 to exist- ing methods. The surface term 2πP 2 x is well-known in the literature on simulation of liquid-solid interfaces[54– 56]. A benchmark in the field is the 1999 paper by Yeh- Berkowitz who introduced a correction term to decouple the electrostatic interactions between a slab of material and its periodic images[54]. The same term is used in sur- face science known there as the dipole correction[57, 58]. 6 These two corrections are identical and are moreover equal to the coupling term in Eq. 26 for E0 = 0. Note, however, that interface/surface modelling and the SSV scheme aim for different target systems. The ideal model system in computational surface science is an isolated slab suspended in vacuum. A similar setup is often used to model liquid-solid interfaces inserting a vacuum spacer in the solid. Periodic models include therefore vacuum layers of a width comparable to or larger than the parti- cle system. The Yeh-Berkowitz correction is intended as a computationally efficient replacement of the costly 2D Ewald sum method (it seems to be doing a very good job as shown in Ref. 54). SSV models, on the other hand, are designed to represent continuous condensed phase sys- tems. These systems can be homogeneous such as liquid water studied here. Vacuum layers opening up unwanted interfaces are avoided. Finally, moving on to the general SSV constant D Hamiltonian of section II A we set Dx = D,Dy = Dz = 0 in Eq. 8 to obtain U(D, v) = HPBC(v) + Ω 8π (D − 4πPx(v)) 2 (27) +2πΩ ( Py(v)2 + Pz(v)2 ) The coupling term is different from the one in Eq. 26 for the open circuit parallel plate capacitor. In Eq. 27 adds further quadratic terms for the polarization in the perpendicular y and z direction. However, while Eq. 27 may look unfamiliar, or even unphysical from the perspective of physical chemistry, it is supported by a thermodynamic foundation via Eq. 13 giving it a spe- cial status among other forms of constant applied field hamiltonians[26, 31, 52]. Setting D = 0 in Eq. 27 we obtain UD=0(v) = HPBC(v) + 2πΩP2 (28) This Hamiltonian is of special interest as it used to sam- ple the D = 0 fluctuations in Eq. 20. The self inter- action term in Eq. 28 resembles the (2πΩ/3)P2 surface term in Ewald summation for spherical vacuum bound- ary conditions[3] but is however a factor three larger. The difference can be traced back to the finite trans- verse polarization of a polarized sphere surrounded by vacuum. Further discussion is deferred to appendix B. Surprisingly, the Hamiltonian Eq. 28 is known in the framework of the reaction field method[6–8]. Eq. 28 can be reproduced by setting the dielectric constant of the embedding dielectric continuum to zero[47, 59](for de- tails see Appendix C). Such an Hamiltonian has in fact been used for simulation studies of the dielectric proper- ties of ionic solutions[13]. The Hamiltonian, while use- ful, was regarded as somewhat unphysical. It is not in the context of SSV theory. Moreover, Anthony Maggs has pointed out that an Hamiltonian of the form Eq. 28 can also be obtained from the time dependent Maxwell equations[59, 60]). His argument is summarized in Ap- pendix B. D. Model system and molecular dynamics The theory outlined above was verified by a classi- cal molecular dynamics simulation of liquid water at ambient conditions. The system consists of 706 water molecules in a fixed cubic box with length 27.7 Å. The interactions are described by the SPC/E model poten- tial [42]. The molecules are kept rigid using the SETTLE algorithm [61]. The MD integration time step is 2 fs. The Ewald summation is implemented using the Particle Mesh Ewald (PME) scheme[62]. Short-range cutoffs for the van der Waals and Coulomb interaction in the di- rect space are 10 Å. The temperature is controlled by a Nosé-Hoover chain thermostat [63] set at 298K. All sim- ulations are done with a modified version of GROMACS 4 package [64]. Two points need further comments for practice: one is the computation of the macroscopic polarization and the other is the force calculation in constant E and constant D simulations. Polar liquids such as water are extended systems as well. The total dipole moment of a periodic supercell depends, in principle, on how we decide to draw the boundaries [37–39] . If the bonds of a molecule are cut by the boundaries, the two halves of this molecule will end up on opposite sides of the MD cell resulting in a huge increase of the dipole moment. This problem can be ignored in practice because the molecular structure provides a natural gauge and it is automatically done for the rigid water used in simulations. Therefore, the macroscopic polarization (to be distinguished from the microscopic polarization) can be simply defined as the sum of the dipole moments of the molecules. For constant E simulation, the field-dependent force on atom i is qiE, where qi = ∂P ∂ri is simply the partical charge. For constant D simulation, the field-dependent force on atom i is qiD− 4πqiP. In this case, the force depends explicitly on the value and continuity of the macroscopic polarization. IV. RESULTS A. Structure and dynamics Theory predicts (Eq. 21) that a D = 0 constraint has the effect of suppressing polarization fluctuations com- pared to E = 0 conditions. The corresponding relaxation times are also faster. This is shown in Fig. 2 for the x component (the system is isotropic, so Py and Pz behave the same). The mean of Px vanishes for both the E = 0 and D = 0 time series (Fig. 2a)) but the amplitude of the E = 0 oscillations is significantly larger. For an analysis of the time dependence (Fig.2b) it is useful to recall the classical Debye theory of the relax- ation of polarization. The relaxation in Debye theory is exponential. Indeed, as Fig. 2b shows, the autocor- relation function of Px at E = 0 decays exponentially (with a hint of a slow oscillation). The relaxation time is 10.3 ps, which is close to the experimental Debye re- 7 FIG. 2. Simulation of bulk liquid water at E = 0 using the hamiltonian of Eq. 4 and D = 0 using the hamiltonian of Eq. 8: a) Time evolution of Px, the x component of the po- larization; b) Corresponding autocorrelation function defined as CPxPx = 〈Px(0)Px(t)〉/〈Px(0)Px(0)〉. The inset shows the short time behaviour of CPxPx for D = 0. laxation time τD for water. τD is the time for response to a sudden change in the electric field E. The other relaxation time defined in Debye theory is τL controlling the response to a change in D (or equivalently a jump in the charge of a solute). τL = τD/ε and τL and τD can be interpreted as a longitudinal respectively transverse relaxation time[33, 34]. In agreement with this picture, we find that switching from E = 0 to D = 0 accelerates the relaxation, but the time dependence in the open cir- cuit system appears to more complex. The short time behaviour (see inset) clearly shows oscillations. Estimat- ing an effective decay time from the time envelope, we obtain τL = 0.3 ps. The pronounced contrast in magnitude and time scale of polarization fluctuations raises the question whether this collective behaviour is reflected in the local molec- ular structure and dynamics. Fig. 3 shows the oxygen radial distribution function, molecular diffusion rate as measured by the mean square displacements and autocor- relation function of the molecular dipole moment charac- terizing molecular orientation. These properties are of- ten used as probe of the structure and dynamics in liquid water at the single molecular level. As can be seen from Fig. 3 the change in electrical boundary condition has little or no effect on the translational and orientational motion of the water molecules. FIG. 3. Simulations of bulk liquid water at E = 0 and D = 0. a) The oxygen-oxygen radial distribution functions gOO(r); b) The mean squared displacements (MSD) of water molecules; c) Autocorrelation function of the molecular dipole moment Cµµ. B. Static dielectric constant from fluctuations The formalism of sections III A and III B gives us four different estimates of ε, two from fluctuations, Eq. 18 and Eq. 20, and two from finite field derivatives, Eq. 22 and 23. We start with the fluctuation approach. The estimates of ε for SPC/E calculated from Eq. 18 and Eq. 20 are 71.4(8) and 76(4) in good mutual agree- ment and with literature values ranging from 67 to 81 [18, 21, 23, 24]. While the dielectric constant estimate should be the same, whether determined under E = 0 or D = 0 constraints, the polarization fluctuations under these conditions are very different. This is reflected in the r dependent Kirkwood G-factor GK(r), which is a orientational correlation function for (rigid) dipoles(see for example Ref. 8). It is defined as GK(r) = 1 +N(r) ∑ j, rij 22Å. The Kirkwood G-factor interpolates between local and global behaviour. The variation with r at short range (r < 6Å) is similar for E = 0 and D = 0. The two curves part for increasing values of r. Calculations of the static dielectric constant cannot be presented without error analysis. As demonstrated in 8 FIG. 4. Comparison of the distance dependence of the Kirk- wood G-factor GK(r) evaluated under E = 0 and D = 0 constraints Refs. 25 and 31 finite size effects are less of a concern for system sizes accessible to classical MD simulation. The MD cell used here containing 706 water molecules should be large enough for the purpose. The time scale needed to converge a second moment of the total dipole moment is a more critical issue. This is confirmed by the accu- mulating average of the normalized variance of the total dipole moment gK determined with regular Ewald sum- mation (E = 0) shown in Fig. 5a. Consistent with the literature we find that it takes at least several nanosec- onds to reduce the statistical uncertainty to a value below 1%. As can be seen from Fig. 5b, the same accuracy is reached within less than one nanosecond by changing the electrical boundary conditions to D = 0. Unfortunately, because of the troublesome inverse relation between fluc- tuations and dielectric constant (Eq. 20) the accuracy in the second moment must be proportionally higher, and much of the apparent gain in time scale is lost in practice. Returning to the question discussed in section III C how to compare the SSV method to methods used in physical chemistry, we note that the sensitivity of the Kirkwood G factor to a change of boundary condition has been studied in detail by Neumann for the Stockmayer fluid[8]. Neumann builds on the familiar cavity model of Kirkwood. His systems consist of a sphere containing the Stockmayer atoms (point dipoles with short-range Lennard-Jones pair interactions) embedded in a dielec- tric continuum. The dielectric constant of the continuum ε′ is varied from ε′ =∞ (conducting) to ε′ = 1 (vacuum). The distance dependent Kirkwood factors Neumann ob- tains for conducting and vacuum boundary conditions have a clear resemblance to our results of Fig. 4 for liq- uid water with E = 0 corresponding to ε′ =∞ and D = 0 to ε′ = 1. For conducting boundary conditions this was expected because, as mentioned, the expression relating the dielectric constant to the dipole fluctuations (Eq. 18) agree. However, the expression derived by Neumann for the ratio of the ε′ = ∞ and ε′ = 1 total dipole fluctu- ations is (ε + 2)/3, which (for large ε) is a factor three smaller than what we obtained in Eq. 21. Indeed, the ratio between the normalized variance of the total dipole FIG. 5. The accumulating average of the normalized variance of the total dipole moment gK of the MD cell with the length of the MD run. gK = (〈M2〉 − 〈M〉2)/(Nµ2), where N is the number of water molecules, µ is the (fixed) dipole moment of a single water molecule and M = ∑N µi. The shaded area is the margin for a 1% deviation from the final average. moment gK at E = 0 and D = 0 from our simulations gives 71 directly validating Eq. 21. This again raises the question about a possible geometric interpretation of the SSV D = 0 Hamiltonian. This will be discussed in detail in appendix B. C. Dielectric constant from field derivatives Finite E simulations necessarily involve a limited sub- set of state points. We selected five Ex = E values with increasing strengths, as listed in the first column of Ta- ble I. For these five values we carried out constant Ex simulations using the Hamiltonian of Eq. 15 determining for each of these runs the average of Px which is indicated in the Table. Next we used Eq. 12 to compute the D val- ues corresponding to the Px we had obtained. These values of D were then taking as the displacement field in constant Dx simulation using the Hamiltonian of Eq. 16. If the SSV constant D method works as promised, the resulting Px are the same as those of the constant Ex simulations from which the Px values were sampled. In- deed, as shown in the second and fourth columns of Ta- ble I, this is the case. Note that the values of D are an order of magnitude larger compared to E for the same value of Px, reflecting the efficient dielectric screening in liquid water. After this crucial consistency test, we computed the E and D derivative estimate of the dielectric constants using the method explained in section III B. The compar- 9 TABLE I. Simulation conditions at constant Ex or constant Dx and the corresponding observed 〈Px〉. Ex (V/Å) 〈Px〉 (10−3 e/Å2) Dx (V/Å) 〈Px〉 (10−3 e/Å2) 0.01 3.72(3) 0.684 3.724(2) 0.02 6.41(5) 1.180 6.418(2) 0.04 9.66(6) 1.788 9.660(1) 0.10 12.62(3) 2.385 12.632(1) 0.28 14.507(4) 2.907 14.512(1) FIG. 6. a) The static dielectric constant ε at constant E and constant D; b) The accumulating average of ε at Ez = 0.01 V/Å and Dz = 0.684 V/Å. ison of ε obtained form the polarizability as a function of D(Eq. 23) to ε computed from the susceptibility as a function of E(Eq. 22) is plotted in Fig. 6a. Regard- ing statistics and convergence, the small values of the field are the most critical and computer time consuming. Fig. 6b gives the running average for the state point cor- responding to our smallest electric field. Comparison to Fig. 5 confirms that the convergence for averages of the polarization is still faster than for the second moment even for small fields. This is of course as expected. It is encouraging to see that the convergence time of the dielectric constant under constant Dx turns out to be shorter than for constant Ex, even with the unfavourable inverse relation of Eq. 23 where the relative error δε/ε is proportional to ε. Using Eqs. 22 and 23 the way we did in Fig. 6 amounts to a global linear approximation. If the dielectric re- sponse of water was linear in the range of fields investi- gated the curves in Fig. 6a would be horizontal straight lines. Not surprisingly, they are not. The variation of ε with field strength is not even linear. Dielectric saturation for increasing E is known to follow the so- called Debye-Langevin equation (ε ∼ 1/E(coth(βµE) − FIG. 7. Polarization as a function of electric field E and dis- placement D determined from constant E respectively con- stant D molecular dynamics. 1/βµE)), derived for an independent dipole approxima- tion [65, 66]. Consistent with this simple picture the ε(E) dependence obtained from the constant E simulations is approximately exponential approaching an asymptotic value at about E = 0.2V/Å. The curvature in ε(D) is opposite to the curvature in ε(E). The non-linear effect in the ε(D) curve is much less pronounced. This is also evident from a direct comparison of the Px(D) to the Px(E) dependence(Fig. 7). This effect is, from technical point of view, perhaps the most encouraging observation made in this study, because it makes the extrapolation to zero field easier. V. SUMMARY AND OUTLOOK Including a term coupling polarization to electric fields is the first step in the derivation of fluctuation expres- sions for dielectric response coefficients of polar liquids in Kirkwood-Fröhlich theory [1, 2]. With the develop- ment of molecular dynamics methods, such hamiltoni- ans have also been used for finite field simulations. The electric field in the coupling term is traditionally the ap- plied electric field. Polar liquids are extended systems as are solids. The key innovation brought about by the modern theory of polarization in solids[37–39] is replac- ing the applied electric field by the macroscopic electric field which includes the internal field generated by the polarization[35]. This, from the perspective of physi- cal chemistry, rather bold move was made by Stengel, Spaldin and Vanderbilt(SSV) on the basis of a clear un- derstanding of what is included in Ewald summation of electrostatic interactions in periodic systems and what is not. This also gave the theory a firm foundation in the thermodynamics of dielectrics[51] which enabled SSV to transform their constant E to a constant D Hamiltonian. This study is a report on the application of the fi- nite temperature classical force field variant of the SSV scheme in a calculation of the dielectric constant of SPC/E liquid water. We started by rederiving the es- tablished fluctuation expression of the dielectric constant 10 for supercells which is generally credited to deLeeuw- Perram-Smith[3, 4] and Neumann[8, 9]. Coupling po- larization directly to the macroscopic electric field made this derivation a simple exercise in perturbation theory. A system under E = 0 constraints (equivalent to stan- dard Ewald summation) can be interpreted as a short cir- cuited capacitor (Fig. 1a). Using the new SSV constant D hamiltonian we also obtained the corresponding fluc- tuation formula for the dielectric constant under D = 0 conditions which can be compared to open circuit con- ditions (Fig. 1b). The theory was tested in E = 0 and D = 0 molecular dynamics simulation. Complementary finite E and D simulations were carried out to compare to the dielectric constant calculated directly from the field derivatives. The estimate of the dielectric constant ob- tained from E = 0 and D = 0 polarization fluctuations were found to be in good agreement with each other and with the estimates from finite field derivatives, validating the SSV method for classical force fields based MD. The motivation for this study was the possibility that application of constant D methods could reduce the costs of the computation of the static dielectric constant ε. This expectation was based on the theory of dielectric relaxation predicting that decay of longitudinal polar- ization is significantly faster compared to transverse po- larization. Arguing that polarization fluctuations under constant D are longitudinal, one can hope that fixing the dielectric displacement instead of the electric field will speed up the convergence of averages and second moments of polarization. This prediction was verified by the simulation. The significant gain in time scale turned out to be however of limited help for the computation of ε from fluctua- tions at D = 0 because of the more stringent demands on accuracy. The reason is that ε under constant elec- tric displacement must be computed from the inverse of a small number depending on the fluctuations(Eq. 20). The calculation of the dielectric constant from electric displacement derivatives suffers in principle from a sim- ilar problem(Eq. 23). Fortunately the statistics in this case is more favourable. A further advantage is that non- linear effects in the response to a change in displacement field are very modest compared to a change in the elec- tric field making the extrapolation to zero field easier. The constant D method may be therefore in the end the best option for the computation of dielectric response in DFT based MD, which is, for water, still effectively out of reach when using zero or constant E methods. We conclude with an outlook. The modern theory of polarization was developed to resolve the fundamental question of the treatment of polarization in solids. The founding fathers of the theory of polar liquids had no such problems, using a system consisting of point dipoles as their basic model.[1, 67] However, this simplicity is lost for realistic point charge models of polar molecules which led to rather confusing arguments about the contribu- tion of higher multipole moments to the polarization[22]. This issue is avoided in the modern of polarization by a strict focus on macroscopic polarization, at the expense of turning polarization into a multivalued quantity.[37– 39] The multivalued polarization could be ignored in the present application to liquid water. We simply held on to polarization as the sum of molecular dipole mo- ments viewing it as a special gauge appropriate for molecule systems. This will no longer work for appli- cations to interfaces between solids and electrolytic solu- tions. Liquid-solid interfaces are described in the frame- work of macroscopic Maxwell theory by dividing the sys- tem up in piece-wise uniform dielectric continua. This conventional approach is regarded by some as incompat- ible with microscopic theory[68]. To resolve these prob- lems it would be useful to extend the modern theory of polarization by reintroducing some form of local polar- ization. The question of local polarization has already been addressed in the context of modelling of solid-solid heterojunctions[36, 69–71]. It remains to be seen whether these concepts can be applied to the electrical double lay- ers formed at electrolyte-solid interfaces. If anything, the challenges in this area of research, and physical electro- chemistry in general, should be an inspiration for both physical chemists and solid state physicists and we are hopeful that progress can be made in the near future. ACKNOWLEDGMENTS Research fellowship (No. ZH 477/1-1) provided by German Research Foundation (DFG) for CZ is grate- fully acknowledged. CZ thanks for helpful discussions with P. Wirnsberger on the modified Green function ap- proach to Ewald summation. CZ and MS also thank R. M. Lynden-Bell for encouraging discussions. We are in particular grateful to A. Maggs for his explanation of an enlightening alternative view of polarization in peri- odic systems and to O. Steinhauser for his clarification of reaction field methods. Appendix A: Legendre transforms From the way they appear in the electric work relation Eq.13, the macroscopic field E and dielectric displace- ment D must be considered as thermodynamic conjugate variables[51]. This suggests that the electric enthalpy F (E) of Eq. 6 and internal energy U(D) of Eq. 9 are each others Legendre transform. However, while accord- ing to Eq. 11 the D derivative of U yields E, Eq. 3 is not consistent with a Legendre transform. The E deriva- tive of F is not reproducing −D. Following Landau and Lifshitz SSV adjust the definition of F (E) to[51] F̃ (E) = F (E)− Ω 8π E2 (A1) 11 Then with Eq. 12 dF̃ dE = − Ω 4π D (A2) as required. The term added to F (E) is a constant in the ensemble generated by the electric enthalpy Hamiltonian of Eq. 4 and will therefore not affect averages over the ensemble. Writing F̃ and U in thermodynamic form as the sum of an internal energy and entropic (−TS) term we find the expected Legendre transform relation F̃ (E) = U(D)− Ω 4π E ·D (A3) U can therefore be equated with the thermodynamic po- tential with respect to D and F̃ as the thermodynamic potential with respect to E referred to by Landau and Lifshitz as U respectively Ũ [51]. The negative sign of the ΩE2/8π term in Eq. A1 is consistent with this interpre- tation. The sign of field energy in electric enthalpy is op- posite (negative) to the sign in internal energy[51]. Sub- tracting ΩE2/8π in Eq. A1 therefore amounts to adding in the energy of the constant macroscopic field as ex- plained in the supporting information of Ref. 35 SSV also consider the Legendre transform of F (E) with respect to P using Eq. 3 leading to the thermodynamic potential E(P) = F (E) + ΩE ·P (A4) with the polarization derivative dE dP = ΩE (A5) The statistical mechanics generated by this potential cor- responds to an ensemble at fixed polarization P[35] which was used in the pioneering studies of the electrical equa- tion of state of ferroelectric materials[40, 41]. The function E(P) also has already a long history in physical chemistry[72] (going back to gas-phase chem- ical thermodynamics). This is presumably the reason why it was chosen by Aragones et al. as the starting point in their study of the electric field dependence of the phase diagram of ice[52]. They then change this into a constant E approach by applying the reverse Legen- dre transform giving them F (E) and the corresponding electric enthalpy Hamiltonian of Eq. 4. Finally to con- vert to a constant applied field method the macroscopic field E in Eq. 4 is replaced by E0. Aragones et al. base their approach on the thermodynamic theory of Landau and Lifshitz. Their finite field method is consistent with the theory in section II A which provides a further micro- scopic basis and a clarification of the distinction between the applied and macroscopic field (see section III C). Appendix B: Surface terms The issue of surface terms arose when it is was realized that the electrostatic interactions in a finite but large ionic crystal can be modelled by an intrinsic energy equal to the “tinfoil” Ewald energy of a supercell in the infinite crystal and a shape dependent extrinsic term[44]. The derivation familiar to physicists and chemist alike is due to deLeeuw, Perram and Smith (LPS)[3, 4]. The problem was revisited again and again generating an extensive literature from which we only quote a small subset[45– 50]. SSV add two field dependent coupling terms to the Ewald Hamiltonian. From Eq. 4 we have for constant E VE = −ΩE ·P (B1) and for constant D from Eq. 8 VD = Ω 8π (D− 4πP) 2 (B2) The question investigated in this appendix is whether VE and VD can be regarded as a surface terms. To this end we consider a piece of dielectric material with an homogeneous polarization density P representing the k = 0 component of the instantaneous polarization density in a liquid. Terminating the dielectric at an interface creates a polarization surface charge density σp σp = −n ·P (B3) where n is the normal to the bounding surface pointing inward. σp generates an electric field called the polariza- tion field Ep = −∇r ∫ dA σp(r ′) |r− r′| (B4) In the context of the physics of ferroelectricity Ep is usu- ally referred to as the depolarizing field. Ep is closely re- lated to the longitudinal component of the polarization, which will be defined below. The separation into a longitudinal component PL and transverse component PT is a rigorous mathematical re- sult of the Helmholtz theorem of vector calculus[73] P = PL + PT (B5) with PL and PT given PL(r) = −∇rφ(r) (B6) φ(r) = ∫ V dr′ ∇r′ ·P(r′) 4π|r− r′| + ∫ A dA n(r′) ·P(r′) 4π|r− r′| PT (r) = ∇r ∧Q(r) (B7) Q(r) = ∫ V dr′ ∇r′ ∧P(r′) 4π|r− r′| − ∫ A dA n(r′) ∧P(r′) 4π|r− r′| so that ∇∧PL = 0 and ∇·PT = 0. The reader is referred to Matyushov for a discussion of the role of transverse polarization in solvation and hydration[68, 74]. For homogeneous polarization, there are no volume contributions, only surface terms. The surface terms are in principle position dependent, even if their sum P is 12 homogeneous. However, because P is constant over the whole of the body, ∇ ·PL = ∇ · (P−PT ) = 0. Similarly ∇ ∧ PT = 0. PL and PT are cavity fields, satisfying vacuum electrostatic equations. A possible r dependence of PL and PT will be henceforward suppressed. Comparing Eqs. B3, B4 and B6 we obtain a general equation for the depolarizing field of a homogeneously polarized dielectric body of arbitrary shape identifying it with the longitudinal component of polarization. Ep = −4πPL (B8) Adding the applied field E0 gives the macroscopic field E = E0 + Ep (B9) and therefore with Eq. B8 E = E0 − 4πPL (B10) Next we reformulate Eqs. B9 and B10 replacing the ap- plied electric field E0 by the more fundamental displace- ment field D defined by the relation D = E + 4πP (B11) Substituting Eq. B5 and B10 we find D = E0 + 4πPT (B12) confirming that the longitudinal component DL of the electrostatic induction D can be treated as an applied electric field E0. However, Eq. B12 also states that elec- trostatic induction cannot be simply equated to the ap- plied field. Depending on the shape of the bounding sur- face D may contain a transverse residue. The E = 0 system is the original “tinfoil boundary geometry” of LPS. There is no surface term, and also VE of Eq. B1 is zero for zero field. This is straight forward. The difficulty is the D = 0 system. For zero displacement field the function VD of Eq. B2 becomes equal to VD(0) = 2πΩP2 (B13) Can this term can be interpreted as a LPS-type surface term? The answer to this question is negative. We base our argument on work by Kantorovitch[49], who showed that the surface term of an ellipsoid shaped cluster of supercell replicas can be written in the form Vs (E) = −Ω 2 Ep ·P (B14) where Ep is the depolarizing field given by Eq. B4. Re- cently Ballenegger showed that this relation is valid for a smooth surface of arbitrary shape[50]. The dependence on surface geometry is implicit in the polarization field Ep. Clearly VD(0) of Eq. B13 and Vs of Eq. B14 are equal if Ep = −4πP (B15) Comparing to Eq. B8 we see that for a geometry to satisfy Eq. B15 the polarization must be entirely longitudinal P = PL or PT = 0. This can be realized in special directions. For a isotropically fluctuating polarization it must hold for all directions. The dielectric response considered by LPS is isotropic. However, for a sphere PL = P/3 or equivalently PT = 2PL (see for example Ref. 73). The depolarizing field of a sphere fails Eq. B15. Do closed surfaces for which PT is strictly zero exist? As far as we are aware they don’t. There will always be some direction in which the depolarizing field is less than what the full polarization would give (Eq. B15). This also implies that there is no dielectric body for which the displacement field is equal to a finite applied electric field of arbitrary orientation (Eq. B12). We are unable to provide a mathematical proof for this statement which must therefore remain a hypothesis. This hypothesis is however strongly supported by the work of Maggs who pointed out that Eq. B15 is in fact satisfied for a periodic system without surfaces[59, 60]. While such a geometry cannot be realized experimentally, it can be considered as a theoretical and computational construction. The argument is briefly summarized below. For further justi- fication and discussion we refer to the original papers. Maggs starts from the Helmholtz theorem Eqs. B5-B7, but now applied to a general inhomogeneous electric field E(r). E(r) = −∇φ(r) +∇∧Q(r) + Ē (B16) The first term is the longitudinal field EL = −∇φ derived from a scalar potential φ. The formalism also allows for a transverse component ET = ∇ ∧ Q. The system is periodic, there are no boundaries generating surface terms. This is the crucial difference with LPS approach. Instead we have in Eq. B16 the r independent term Ē, which is an as yet unspecified uniform field. Using an infinite system from the start disregarding surfaces is also the basis of the derivation of the frequency and wavevector dependent dielectric and magnetic lin- ear response functions by Fulton[32] and Madden and Kivelson[33]. The crucial step proposed by Maggs is to link the uniform (k = 0) electric field Ē to the polariza- tion using the time dependent Maxwell equation ∂E(t) ∂t = −4πJ + c∇∧B (B17) where J is the current and B the magnetic field. c is the velocity of light. Eq. B17 is a microscopic Maxwell equation in Gaussian units. Integrating over time and space gives the uniform polarization Ē(t) = −4π Ω ∫ t t0 dt′ ∫ cell dr J(t′) ≡ −4πP(t) (B18) where we have used that the spatial integral of a curl over the unit cell of a periodic system vanishes[59]. We can imagine that the system was subjected to some action starting from an unpolarized reference state at time t0. 13 The adiabatic current this action creates establishes a final polarized state. This is also how polarization is defined in the modern theory of polarization. Eq. B18 relates the homogeneous polarization P to an internal electric field Ē. This the only electric field there is, hence with Eq. B11 we conclude that D = 0 which then allows us to equate Ē of Eq. B18 to the depolarizing field Ep satisfying Eq. B15. The electric energy of a volume Ω cut out of the system is ΩE2/8π. Replacing the filed by the polarization and adding the Ewald sum representing all k 6= 0 contributions to the energy leads to the D = 0 SSV Hamiltonian of Eq. 28. The time dependent scheme proposed by Maggs and the modern theory of polarization have much in common. Both schemes use the current to define polarization in a periodic system without specifying a surface. This ap- plies to D = 0 boundary conditions as well, which are normally associated with a finite body in vacuum. This may seem counterintuitive, but is consistent with the ge- ometry invariance of the constitutive relations[51]. Appendix C: Reaction fields In parallel to the Ewald summation based methods of the deLeeuw, Perram and Smith, an alternative approach based on reaction field methods was developed by Neu- mann and Steinhauser(NS)[6–8]. The central result is a general fluctuation expression for the static dielectric constant ε of a spherical body of polar fluid embedded in a dielectric continuum. ε− 1 ε+ 2 = 4π 3 ( βΩ〈P2〉 3 )[ 1− 3 4π ε− 1 ε+ 2 Tmod(εR) ] (C1) Following the notation of section III A, P is the polariza- tion computed from the total dipole moment of a system of volume Ω. The factor Tmod(εR) is the volume aver- age of the dipolar field tensor adapted (“modified”) for interaction with the embedding dielectric continuum of dielectric constant εR. Tmod(εR) is closely related to the Onsager reaction field of a dipole[2]. Tmod(εR) = 4π 3 2 (εR − 1) 2εR + 1 (C2) For the derivation of Eqs. C1 and C2 we refer to the orig- inal papers by NS[6, 8]. Substituting εR = ε in Eq. C1 recovers the Kirkwood-Fröhlich expression for the dielec- tric constant[8]. However, εR can also be different from ε. While this changes the magnitude of the polarization fluctuations, Eq. C1 combined with the reaction field fac- tor Eq. C2 still gives the correct relation to the dielectric constant. In particular, the limit εR → ∞ corresponds to a polar fluid in a spherical cavity in a metal. The reaction field factor Eq. C2 becomes Tmod(∞) = 4π 3 (C3) which inserted in Eq. C1 yields ε− 1 = 4πβΩ 3 〈P2〉εR=∞ (C4) Eq. C4 is identical to the relation of Eq. 18 between 〈P2〉E=0 and the dielectric constant (for simplicity here we have assumed that 〈P〉 = 0 as it should be in a con- verged MD run). This led NS to identify Ewald summa- tion with their εR =∞ reaction field limit[6–8] The opposite case of a sphere in vacuum is obtained by setting εR = 1, and hence Tmod(1) = 0 (C5) giving the fluctuation relation ε− 1 ε+ 2 = 4πβΩ 9 〈P2〉εR=1 (C6) As can already be expected from the argument of ap- pendix B, the vacuum equation C6 is not equal to the fluctuation relation Eq. 20 under D = 0 conditions. However, as noticed by Caillol and coworkers[13, 47] the NS reaction field approach is in fact capable of repro- ducing the D = 0 Hamiltonian Eq. 28 (see also Ref.59). This is achieved by setting the dielectric constant of the embedding medium to the “unphysical” value of εR = 0. Inserting in Eq. C2 gives Tmod(0) = −8π 3 (C7) which should be compared to Eq. C3 valid for E = 0 (note the change of sign). Substituting in Eq. C1 we find ε− 1 ε = 4πβΩ 3 〈P2〉εR=0 (C8) which is indeed the counterpart to Eq. 20. Similarly, dividing Eq. C4 by Eq. C8 gives the same ratio as Eq. 21 〈P2〉εR=∞ 〈P2〉εR=0 = ε (C9) confirming that D = 0 boundary conditions are con- tained in the NS reaction field formalism, all be it for an environment with the unphysical dielectric constant of zero. The interpretation of this rather surprising connection is better left to the experts in reaction field methods. As mentioned it was noticed later and is not discussed in the original papers on the NS method. We furthermore point out, that the thermodynamic perspective underlying the SSV method greatly simplifies the derivation of the zero field fluctuation relations avoiding the complication of the dipolar tensor. Moreover, it allows for natural exten- sion to finite displacement fields D as we have shown in the present paper. 14 [1] J. G. Kirkwood, J. Chem. Phys. 7, 911 (1939). [2] H. Fröhlich, Theory of Dielectrics: Dielectric Constant and Dielectric Loss, 2nd ed., Monographs on the physics and chemistry of materials (O.U.P., 1958). [3] S. W. De Leeuw, J. W. Perram, and E. R. Smith, Proc. R. Soc. A 373, 27 (1980). [4] S. W. De Leeuw, J. W. Perram, and E. R. Smith, Proc. R. Soc. A 373, 57 (1980). [5] S. W. De Leeuw, J. W. Perram, and E. R. Smith, Annu. Rev. Phys. Chem. 37, 245 (1986). [6] M. Neumann and O. Steinhauser, Chem. Phys. Lett. 95, 417 (1983). [7] M. Neumann and O. Steinhauser, Chem. Phys. Lett. 102, 508 (1983). [8] M. Neumann, Mol. Phys. 50, 841 (1983). [9] M. Neumann, O. Steinhauser, and G. S. Pawley, Mol. Phys. 52, 97 (1984). [10] M. Neumann, J. Chem. Phys. 82, 5663 (1985). [11] M. Neumann, J. Chem. Phys. 85, 1567 (1986). [12] G. N. Patey, D. Levesque, and J. J. Weis, Mol. Phys. 45, 733 (1982). [13] J.-M. Caillol, J. Chem. Phys. 91, 5555 (1989). [14] M. R. Reddy and M. Berkowitz, Chem. Phys. Lett. 155, 173 (1989). [15] K. Watanabe and M. L. Klein, Chem. Phys. 131, 157 (1989). [16] M. Sprik, J. Chem. Phys. 95, 6762 (1991). [17] S. W. Rick, S. J. Stuart, and B. J. Berne, J. Chem. Phys. 101, 6141 (1994). [18] D. van der Spoel, P. J. van Maaren, and H. J. C. Berend- sen, J. Chem. Phys. 108, 10220 (1998). [19] P. H. Hünenberger and W. F. van Gunsteren, J. Chem. Phys. 108, 6117 (1998). [20] G. Lamoureux, E. Harder, I. V. Vorobyov, B. Roux, and A. D. MacKerell Jr., Chem. Phys. Lett. 418, 245 (2006). [21] J. L. Aragones, L. G. MacDowell, and C. Vega, J. Phys. Chem. A 115, 5745 (2011). [22] K. K. Mandadapu, J. A. Templeton, and J. W. Lee, J. Chem. Phys. 139, 054115 (2013). [23] C. Zhang and G. Galli, J. Chem. Phys. 141, 084504 (2014). [24] D. Braun, S. Boresch, and O. Steinhauser, J. Chem. Phys. 140, 064107 (2014). [25] D. C. Elton and M.-V. Fernández-Serra, J. Chem. Phys. 140, 124504 (2014). [26] J. Kolafa and L. Viererblová, J. Chem. Theory Comput. 10, 1468 (2014). [27] V. Dubois, P. Umari, and A. Pasquarello, Chem. Phys. Lett. 390, 193 (2004). [28] M. Sharma, R. Resta, and R. Car, Phys. Rev. Lett. 98, 247401 (2007). [29] D. Pan, L. Spanu, B. Harrison, D. A. Sverjensky, and G. Galli, Proc. Natl. Acad. Sci. USA 110, 6646 (2013). [30] I.-C. Yeh and M. L. Berkowitz, J. Chem. Phys. 110, 7935 (1999). [31] S. Rinike, A.-P. E. Kunz, and W. F. van Gunsteren, J. Chem. Theory Comput. 7, 14691475 (2011). [32] R. L. Fulton, Mol. Phys. 29, 405 (1975). [33] P. Madden and D. Kivelson, Adv. Chem. Phys. 56, 467 (1984). [34] D. Kivelson and H. Friedman, J. Phys. Chem. 93, 7026 (1989). [35] M. Stengel, N. A. Spaldin, and D. Vanderbilt, Nat. Phys. 5, 304 (2009). [36] M. Stengel, D. Vanderbilt, and N. A. Spaldin, Phys. Rev. B 80, 224110 (2009). [37] R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993). [38] R. Resta, Rev. Mod. Phys. 66, 899 (1994). [39] R. Resta and D. Vanderbilt, in Topics in Applied Physics Volume 105: Physics of Ferroelectrics: a Modern Per- spective, edited by K. M. Rabe, C. H. Ahn, and J.-M. Triscone (Springer-Verlag, 2007) pp. 31–67. [40] N. Sai, K. M. Rabe, and D. Vanderbilt, Phys. Rev. B 66, 104108 (2002). [41] O. Dieguez and D. Vanderbilt, Phys. Rev. Lett. 96, 056401 (2006). [42] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987). [43] P. P. Ewald, Ann. Phys. 369, 253 (1921). [44] A. Redlack and J. Grindlay, J. Phys. Chem. Solids 36, 73 (1975). [45] E. R. Smith, Proc. R. Soc. A 375, 475 (1981). [46] E. R. Smith, J. Chem. Phys. 128, 174104 (2008). [47] J.-M. Caillol, J. Chem. Phys. 101, 6080 (1994). [48] L. M. Fraser, W. M. C. Foulkes, G. Rajagopal, R. J. Needs, S. D. Kenny, and A. J. Williamson, Phys. Rev. B 53, 1814 (1996). [49] L. N. Kantorovich and I. I. Tupitsyn, J. Phys. Condens. Matter 11, 6159 (1990). [50] V. Ballenegger, J. Chem. Phys. 140, 161102 (2014). [51] L. D. Landau and E. M. Lifshitz, Electrodynamics of continuous media. Course of theoretical physics volume 8 (Pergamon Press, Oxford, 1960). [52] J. L. Aragones, L. G. MacDowell, J. I. Siepmann, and C. Vega, Phys. Rev. Lett. 107, 155702 (2011). [53] E. M. Purcell, Electricity and magnetism (Cambridge University Press, Cambridge, 2011). [54] I.-C. Yeh and M. L. Berkowitz, J. Chem. Phys. 111, 3155 (1999). [55] V. Ballenegger and J.-P. Hansen, J. Chem. Phys. 122, 114711 (2005). [56] D. J. Bonthuis, S. Gekle, and R. R. Netz, Langmuir 28, 7679 (2012). [57] J. Neugebauer and M. Scheffler, Phys. Rev. B 46, 16067 (1992). [58] L. Bengtsson, Phys. Rev. B 59, 12301 (1999). [59] A. C. Maggs, J. Chem. Phys. 120, 3108 (2004). [60] A. C. Maggs and R. Everaers, Phys. Rev. Lett. 96, 230603 (2006). [61] S. Miyamoto and P. A. Kollman, J. Comput. Chem. 13, 952 (1992). [62] T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). [63] G. J. Martyna, M. L. Klein, and M. Tuckerman, J. Chem. Phys. 97, 2635 (1992). [64] B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008). [65] P. J. W. Debye, Polar Molecules (Chemical Catalog Company, Incorporated, 1929). [66] F. Booth, J. Chem. Phys. 19, 391 (1951). [67] L. Onsager, J. Am. Chem. Soc. 58, 1486 (1936). 15 [68] D. V. Matyushov, J. Chem. Phys. 140, 224506 (2014). [69] M. Stengel and D. Vanderbilt, Phys. Rev. B 80, 241103(R) (2009). [70] M. Stengel, Phys. Rev. Lett. 106, 136803 (2011). [71] F. Giustino and A. Pasquarello, Phys. Rev. B 71, 144104 (2005). [72] R. A. Alberty, Pure Appl. Chem. 73, 13491380 (2001). [73] G. Lehrer, Electromagnetic Field Theory for Engineers and Physicists (Springer, 2010). [74] D. V. Matyushov, J. Chem. Phys. 120, 7532 (2004).