Uncertainties in the Static Dielectric Constants computed from Molecular   Dynamics Simulations

Hern\'an R. S\'anchez

arXiv:1901.02127·physics.chem-ph·November 28, 2019

Uncertainties in the Static Dielectric Constants computed from Molecular Dynamics Simulations

Hern\'an R. S\'anchez

PDF

TL;DR

This study compares methods for calculating dielectric constants in molecular dynamics simulations, analyzing uncertainties, system size effects, and sampling frequency, with water as a test case.

Contribution

It introduces new uncertainty expressions for the external field method and compares them with the fluctuation method, highlighting differences and proposing alternatives.

Findings

01

Significant differences in system size dependence between methods

02

Best dielectric constant estimate for water: 70.46±0.31

03

Sampling frequency impacts uncertainty estimates

Abstract

In this work, the uncertainties in the dielectric constants of polar liquids, computed using molecular dynamics (MD) simulations, are compared for two different calculations schemes. Expressions of the uncertainty are derived for the external field method, and compared with those of the fluctuation method. Significant differences on their the system size dependence were found. In addition, alternative calculation procedures are proposed. The individual contributions of different parts of the system to the electric susceptibility, and their corresponding uncertainties, are also studied. Additionally, the effects of the sampling frequency on the uncertainties are analyzed. MD simulations of pure liquid water (SPC/E) at 298.15 K and 1 bar were performed in order to corroborate the obtained results. In such conditions, the best estimate of the dielectric constant obtained in this study is…

Figures4

Click any figure to enlarge with its caption.

Tables2

Table 1. Table 1: Results of MD simulations of dielectric constant. Columns 2, 5 and 7 correspond to the EFR, FR in two directions and the MVUE, respectively. Columns 3 and 4 contain the standard deviation (SD) for the EFR through the proposed and the numerical method, respectively. Columns 6, 7 and 9 contain the uncertainties corresponding to the FR and MVUE methods. See main text for details.

$E$ /(V/n)	$ε_{r, EFR}$	SD $_{EFR,Prop.}$	SD $_{EFR,Num.}$	$ε_{r, FR, x y}$	SD $_{FR x y}$	SD $_{FR x y, Num.}$	$ε_{r, mix}$	SD $_{mix}$
0.0005	55.551	58.693	54.455	69.877	1.032	1.178	69.873	1.031
0.0010	118.656	31.804	33.295	71.238	1.048	1.007	71.290	1.048
0.0050	57.762	6.262	5.996	71.541	1.069	1.069	71.151	1.053
0.0100	70.384	3.051	2.905	69.496	1.033	1.036	69.587	0.979
0.0500	70.471	0.604	0.636	68.891	1.018	0.997	70.060	0.519
0.1000	67.300	0.274	0.244	69.372	1.015	0.818	67.441	0.264
0.5000	40.045	0.016	0.018	39.949	0.435	0.468	40.045	0.016
1.0000	24.794	0.003	0.004	24.797	0.209	0.183	24.794	0.003
Molecules								—
256	75.040	4.481	4.359	68.994	1.002	0.960	69.282	0.977
512	70.384	3.051	2.905	69.496	1.033	1.036	69.587	0.979
1024	68.991	2.134	1.822	72.040	1.075	1.065	71.423	0.960
2048	69.166	1.479	1.441	69.611	1.029	0.915	69.465	0.845
4096	70.589	1.128	1.056	71.672	1.056	1.054	71.166	0.771
8192	69.761	0.778	0.714	70.843	1.050	1.019	70.144	0.625

Table 2. Table 2: Analysis of the partial susceptibilities. The first column on the left contains the percentage of molecules of the system considered (neglecting rounding errors). The second and third columns contain the estimated partial susceptibilities and their uncertainties computed with the proposed methods for the EFR. The fifth and sixth columns the same for the FR. The remaining columns contain the respective numerical estimations.

Percentage	$χ_{e, EFR}$	SD $_{EFR,Prop.}$	SD $_{EFR,Num.}$	$χ_{r, FR, x}$	SD $_{FR, x}$	SD $_{FR, x, Num.}$
10	7.189	0.315	0.261	7.048	0.200	0.188
25	17.384	0.608	0.507	17.749	0.431	0.402
50	34.257	1.130	0.941	35.240	0.771	0.703
75	50.646	1.632	1.382	52.792	1.113	1.030
90	61.087	1.939	1.662	63.414	1.337	1.231
100	67.991	2.134	1.822	70.401	1.516	1.389

Equations92

ε_{r} = 1 + \frac{Var [ M ]}{3 ε _{0} V k _{B} T}

ε_{r} = 1 + \frac{Var [ M ]}{3 ε _{0} V k _{B} T}

ε_{r} = 1 + \frac{E [ M \cdot u _{E} ^ ]}{ε _{0} V E \cdot u _{E} ^} = 1 + \frac{E [ M ]}{ε _{0} V E}

ε_{r} = 1 + \frac{E [ M \cdot u _{E} ^ ]}{ε _{0} V E \cdot u _{E} ^} = 1 + \frac{E [ M ]}{ε _{0} V E}

ε_{r}

ε_{r}

ε_{r}

ϕ (t) = \frac{E [( M ( t ^{'} ) - E [ M ( t )]) ( M ( t - t ^{'} )) - E [ M ( t )])}{Var [ M ( t )]}

ϕ (t) = \frac{E [( M ( t ^{'} ) - E [ M ( t )]) ( M ( t - t ^{'} )) - E [ M ( t )])}{Var [ M ( t )]}

r (k Δ t) = \frac{\sum _{i = 1}^{n - k} ( M ( i Δ t ) - M ˉ ) ( M (( i + k ) Δ t ) - M ˉ )}{\sum _{i = 1}^{n} ( M ( i Δ t ) - M ˉ ) ^{2}}

r (k Δ t) = \frac{\sum _{i = 1}^{n - k} ( M ( i Δ t ) - M ˉ ) ( M (( i + k ) Δ t ) - M ˉ )}{\sum _{i = 1}^{n} ( M ( i Δ t ) - M ˉ ) ^{2}}

f (t) = e^{- t / τ}

f (t) = e^{- t / τ}

n_{eff} := Var [M] / Var [\overset{ˉ}{M}]

n_{eff} := Var [M] / Var [\overset{ˉ}{M}]

n_{eff} = \frac{n}{1 + 2 \sum _{i = 1}^{n - 1} ( 1 - i / n ) ϕ ( i Δ t )}

n_{eff} = \frac{n}{1 + 2 \sum _{i = 1}^{n - 1} ( 1 - i / n ) ϕ ( i Δ t )}

n_{eff}

n_{eff}

= \frac{T}{2 \int _{0}^{T} e ^{- \frac{t}{τ}} d t} = \frac{T}{2 τ}

n_{eff} = n tanh \frac{T}{2 n τ}

n_{eff} = n tanh \frac{T}{2 n τ}

Var [ε_{r}] = \frac{1}{( ε _{0} V E ) ^{2}} \frac{Var [ M ]}{n _{eff}}

Var [ε_{r}] = \frac{1}{( ε _{0} V E ) ^{2}} \frac{Var [ M ]}{n _{eff}}

Var [ε_{r}] = \frac{1}{( ε _{0} V E ) ^{2}} \frac{2 τ _{E} Var [ M ]}{T}

Var [ε_{r}] = \frac{1}{( ε _{0} V E ) ^{2}} \frac{2 τ _{E} Var [ M ]}{T}

v_{eff} ≊ \frac{n}{1 + 2 \sum _{i = 1}^{n - 1} ϕ ( i Δ t ) ^{2}}

v_{eff} ≊ \frac{n}{1 + 2 \sum _{i = 1}^{n - 1} ϕ ( i Δ t ) ^{2}}

Var [ε_{r}] ≊ (\frac{Var ^{2} [ M _{x} ]}{v _{eff, x}} + \frac{Var ^{2} [ M _{y} ]}{v _{eff, y}} + \frac{Var ^{2} [ M _{z} ]}{v _{eff, z}}) \frac{2}{( d ε _{0} V k _{B} T ) ^{2}}

Var [ε_{r}] ≊ (\frac{Var ^{2} [ M _{x} ]}{v _{eff, x}} + \frac{Var ^{2} [ M _{y} ]}{v _{eff, y}} + \frac{Var ^{2} [ M _{z} ]}{v _{eff, z}}) \frac{2}{( d ε _{0} V k _{B} T ) ^{2}}

ν_{eff}

ν_{eff}

= \frac{T}{2 \int _{0}^{T} e ^{- \frac{2 t}{τ}} d t} = \frac{T}{τ}

Var [ε_{r}] ≊ \frac{2 τ}{T} \frac{Var ^{2} [ M ]}{( 3 ε _{0} V k _{B} T ) ^{2}}

Var [ε_{r}] ≊ \frac{2 τ}{T} \frac{Var ^{2} [ M ]}{( 3 ε _{0} V k _{B} T ) ^{2}}

\frac{Var [ ε _{r} ] _{FR}}{Var [ ε _{r} ] _{EFR}} \approx \frac{τ}{τ _{E}} \frac{E ^{2} Var ^{2} [ M ] _{F R}}{( d k _{B} T ) ^{2} Var [ M ] _{E F R}} \approx \frac{E ^{2} ε _{0} V ( ε _{r} - 1 )}{d k _{B} T}

\frac{Var [ ε _{r} ] _{FR}}{Var [ ε _{r} ] _{EFR}} \approx \frac{τ}{τ _{E}} \frac{E ^{2} Var ^{2} [ M ] _{F R}}{( d k _{B} T ) ^{2} Var [ M ] _{E F R}} \approx \frac{E ^{2} ε _{0} V ( ε _{r} - 1 )}{d k _{B} T}

Var [μ_{z}]

Var [μ_{z}]

= \frac{μ ^{2}}{3} - \frac{μ ^{2}}{15} (\frac{∣ μ ∣∣ E ∣}{k _{B} T})^{2} + \dots

Var [μ_{x}]

Var [μ_{x}]

= μ^{2} (\frac{k _{B} T}{∣ μ ∣∣ E ∣})^{2} [\frac{∣ μ ∣∣ E ∣}{k _{B} T} coth (\frac{∣ μ ∣∣ E ∣}{k _{B} T}) - 1] =

= \frac{μ ^{2}}{3} - \frac{μ ^{2}}{45} (\frac{∣ μ ∣∣ E ∣}{k _{B} T})^{2} + \dots

E [X]_{H_{1}} = \frac{\int _{Ω} X e ^{- β H_{1}} d q}{\int _{Ω} e ^{- β H_{1}} d q}

E [X]_{H_{1}} = \frac{\int _{Ω} X e ^{- β H_{1}} d q}{\int _{Ω} e ^{- β H_{1}} d q}

E [M]_{H_{1}} = \frac{\int _{Ω} M e ^{β M \cdot E} e ^{- β H_{0}} d q}{\int _{Ω} e ^{β M \cdot E} e ^{- β H_{0}} d q} = \frac{E [ M e ^{β M \cdot E} ] _{H_{0}}}{E [ e ^{β M \cdot E} ] _{H_{0}}}

E [M]_{H_{1}} = \frac{\int _{Ω} M e ^{β M \cdot E} e ^{- β H_{0}} d q}{\int _{Ω} e ^{β M \cdot E} e ^{- β H_{0}} d q} = \frac{E [ M e ^{β M \cdot E} ] _{H_{0}}}{E [ e ^{β M \cdot E} ] _{H_{0}}}

Var [ε_{r}]

Var [ε_{r}]

\overset{ˉ}{M}_{H_{1}}

\overset{ˉ}{M}_{H_{1}} = \frac{n ^{- 1} \sum _{i = 1}^{n} ( M _{i} - M ˉ ) e ^{β (M_{i} - \overset{ˉ}{M}) E}}{n ^{- 1} \sum _{i = 1}^{n} e ^{β (M_{i} - \overset{ˉ}{M}) E}}

\overset{ˉ}{M}_{H_{1}} = \frac{n ^{- 1} \sum _{i = 1}^{n} ( M _{i} - M ˉ ) e ^{β (M_{i} - \overset{ˉ}{M}) E}}{n ^{- 1} \sum _{i = 1}^{n} e ^{β (M_{i} - \overset{ˉ}{M}) E}}

n^{- 1} i = 1 \sum n e^{β (M_{i} - \overset{ˉ}{M}) E} = 1 + β^{2} E^{2} \overline{(M_{i} - \overset{ˉ}{M})^{2}} + O (M^{3}) \approx 1

n^{- 1} i = 1 \sum n e^{β (M_{i} - \overset{ˉ}{M}) E} = 1 + β^{2} E^{2} \overline{(M_{i} - \overset{ˉ}{M})^{2}} + O (M^{3}) \approx 1

n^{- 1} i = 1 \sum n (M_{i} - \overset{ˉ}{M}) e^{β (M_{i} - \overset{ˉ}{M}) E} =

n^{- 1} i = 1 \sum n (M_{i} - \overset{ˉ}{M}) e^{β (M_{i} - \overset{ˉ}{M}) E} =

\frac{e ^{- β \overset{ˉ}{M} E}}{n} [i = 1 \sum n M_{i} e^{β M_{i} E} - i = 1 \sum n \overset{ˉ}{M} e^{β M_{i} E}] \approx

\frac{1}{n} [i = 1 \sum n M_{i} e^{β M_{i} E} - \overset{ˉ}{M} n] = \frac{1}{n} i = 1 \sum n M_{i} (e^{β M_{i} E} - 1)

Var [f (M)] \approx (\frac{df ( E [ M ])}{d M})^{2} Var [M] + \frac{1}{2} (\frac{d ^{2} f ( E [ M ])}{d M ^{2}})^{2} Var^{2} [M]

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Uncertainties in the Static Dielectric Constants computed from Molecular Dynamics Simulations

Hernán R. Sánchez

[email protected]

Instituto de Física de Líquidos y Sistemas Biológicos. CONICET La Plata - UNLP, Argentina.

Abstract

In this work, the uncertainties in the dielectric constants of polar liquids, computed using molecular dynamics (MD) simulations, are compared for two different calculations schemes. Expressions of the uncertainty are derived for the external field method, and compared with those of the fluctuation method. Significant differences on their the system size dependence were found. In addition, alternative calculation procedures are proposed. The individual contributions of different parts of the system to the electric susceptibility, and their corresponding uncertainties, are also studied. Additionally, the effects of the sampling frequency on the uncertainties are analyzed. MD simulations of pure liquid water (SPC/E) at 298.15 K and 1 bar were performed in order to corroborate the obtained results. In such conditions, the best estimate of the dielectric constant obtained in this study is 70.46 $\pm$ 0.31.

I Introduction

The importance of the dielectric properties has been recognized for decadesvon1955dielectric . Their measurement is useful for many application in many fields, for example, in the characterization of heterogeneous systemsasami2002characterization . The influence of the relative (static) dielectric permittivity in the interaction between particles is significant for condensed phases. Due to this, in order to get accurate values for many properties, it must be accurately represented when performing computer simulations of physical systems.

For isotropic and homogeneous liquid systems, the dielectric constant can be computed from the total dipole moment. It in turn, can be obtained, among many alternatives, from molecular dynamics (MD) simulationsalder1959studies ; rapaport2004art . The expressions that link the dipole moment with the dielectric constant depend on the boundary conditions. Those used in this work are valid for the tin-foillevrel2008boundary and the Ewaldde1986computer boundary conditions, which are very common choices.

There are many routes of calculation of the static dielectric constant ( $\varepsilon_{r}$ ) of these kinds of systems heinz2001comparison . Two of them are the most used in the literature. In one, the dielectric constant is related to the fluctuations in time of the total dipole moment of the system ( $\boldsymbol{M}$ ), in absence of external fields, via the following equationneumann1983dipole ; kolafa2014static

[TABLE]

where $\varepsilon_{0}$ represents the vacuum permittivity, $k_{B}$ the Boltzmann’s constant, $V$ and $T$ system equilibrium volume and temperature, respectively. $\operatorname{Var}[\cdot]$ represent the variance operator. This calculation scheme will be referred here as the Fluctuation Route (FR).

The other calculation route requires the application of an external time-independent electric field ( $\mathbf{E}$ )riniker2011calculation . If the direction of the field is given by an unitary vector $\check{\mathbf{u}}_{E}$ , the dielectric constant can be obtained fromneumann1983dipole ; kolafa2014static

[TABLE]

where $\operatorname{E}[\cdot]$ represents the expected value operator, and $M$ and ${E}$ the euclidean norm of $\boldsymbol{M}$ and ${E}$ , respectively. Hereafter, this scheme of calculation will be referred as the External Field Route (EFR).

Molecular dynamics can be used for sampling the total dipole moment. As the number of simulation steps is finite, $\operatorname{Var}[\boldsymbol{M}]$ and $\operatorname{E}[M]$ can only be estimated. Unfortunately, obtaining converged values of the dielectric constant require uncommonly large simulation times. Because of this, it is much more difficult to calculate $\varepsilon_{r}$ than many commonly computed propertiescaleman2011force .

From a practical perspective, the dielectric constant cannot be obtained exactly, and commonly one of the following estimators is used,

[TABLE]

where the bar above a symbol denotes its sample arithmetic mean. Many times $\bar{\boldsymbol{M}}^{2}$ is not included in Eq. 3 neumann1983dipole .

In a previous work by the author an others, the analytical expressions for the variance of the dielectric constant computed through the FR (Eq. 3) were foundSANCHEZ2019546 . In this work, analytical expressions for the variance corresponding to the estimator of Eq. 4 (EFR) are derived, and both calculations routes are compared in terms of their corresponding uncertainties. These calculation schemes were compared before in terms of accuracy though numerical experiments, for example in Ref. kolafa2014static . However, the usage of analytical expressions provides clear advantages.

In the present study is also shown, that under usual circumstances, both calculation schemes can be combined for providing greater precision. This can be done by performing a MD simulation with an external electric field, computing the dielectric constant through the EFR, and using the information related to the remaining perpendicular directions in conjunction with the FR.

Afterwards, a new route of calculation is presented. It requires a simulation without an external field, and uses the equation for the dielectric constant corresponding to de EFR (Eq. 4). It can be considered an intermediate step between the EFR and FR, as it arises from the derivation of the FR from the EFR. The proposed scheme provides no added value in terms of accuracy, as its final values coincide with those from the FR. However, the use of the former is advantageous for the study of the contribution to the electric susceptibility of different parts of the system (partial susceptibilities) and their corresponding uncertainties.

The frequency at which the system configurations are stored is significant for the uncertainty in the computed dielectric constant. Although this has been studied before by means of numerical testsgereben2011accurate , it remains unclear in the specific literature. Because of this, the issue is addressed in this work.

II Theory

The largest part of this work is focused in the calculation of uncertainties. The total dipole moment of the system obtained with MD is auto-correlated, i.e, the latter correlates with a delayed copy of itselfgordon1968correlation . This has huge impact in the estimations of the uncertainties of the mean and sample variance, as is discussed in the sections II.1 and II.3.

The auto-correlation function ( $\phi(t)$ ) is a measure of the degree of correlation, of the signal with itself, as a function of the time delay. Its definition may depends on the context, in this work the following definition is usedbrockwell1991time

[TABLE]

The definition above reflects the strength of the statistical dependence at different times, as it is normalized. The auto-correlation function is significant for this study, because the former appears in the derived expressions for the uncertainty of the dielectric constant. Many estimators of the auto-correlation function exist. Their differences are not significant for this work, in which the most popular discrete estimator was used, and is defined by

[TABLE]

The dielectric permittivity can also be studied and modeled as a function of the frequency of the external field appliedhilfer2002h . Those representations in frequency domain have their counter-part in the time domainbottcher1978theory . The expression in time domain describe the dielectric relaxation process. It can be shown that, taking aside constants, the time domain expression coincide with the auto-correlation of the total dipole moment of the system. This is significant for this work because it allows to approximate the auto-correlation function with well known expressions of the dielectric relaxation in time domain. Among them, the Debye’s relaxation modelhilfer2002h will be used due to its simplicity. Its relaxation function which describes the decay of polarization isgarrappa2016models

[TABLE]

where $\tau$ represents the relaxation time.

II.1 Simulation length and sampling frequency

As mentioned above, computing the dielectric constant requires very large simulation times. This imply that huge file sizes will be obtained if the simulation coordinates are written to disk at every simulation step, or at least very frequently. In such cases, the subsequent processing of the full information may be computationally expensive. Fortunately, it is not necessary as the additional information provided for correlated steps is often less useful than the corresponding to uncorrelated steps. When observations are uncorrelated, the variance of the mean of a random variable $M$ satisfies $\operatorname{Var}[\bar{M}]=\operatorname{Var}[M]/n$ , where $n$ represents the number of observations. It is not true for autocorrelated variables, however, a similar equation is obtained by defining the effective number of observations asbayley1946effective

[TABLE]

which is given by

[TABLE]

For the cases taken into account in this paper, a good approximation is to consider the Debye’s model. This implies that $\phi(t)=e^{-\frac{t}{\tau}}$ . Then, if $\tau<<\mathcal{T}:=n\Delta t$ , we can approximate $1-i/n\approx 1$ . If both substitution are made in Eq. 9, and $\Delta t$ is very small

[TABLE]

An alternative and elegant derivation of the above equation was previously found, for auto-regressive processesbrockwell1991time , in the reference zikeba2010effective . In that work, it was also shown that for the process studied here

[TABLE]

This equation should be taken into account when deciding the sampling frequency. It is also significant in that it may serve to estimate the applicability of Eq. 10. The proposed derivation can be adapted to other relaxation processes which are not well represented by the Debye’s model. The Eq. 9 cannot be used in practice because $\phi(i\,\Delta t)$ is not known. However, it can be estimated from a small sample and the quotient $n_{\text{eff}}/n$ may help in choosing a sampling frequency in a more general cases.

II.2 Uncertainty in the dielectric constants through the external field route

Applying the variance operator to Eq. 4, and combining the result with Eq. 8, the following expression is obtained

[TABLE]

The Eq. 12 allows to evaluate the uncertainty in the dielectric constants computed though the EFR. The standard deviation turns to be inversely proportional to the external electric field, as stated previously in reference kolafa2014static . When the relaxation times are known, the Eq. 12 can be useful for the estimation of the simulation length required for a given target uncertainty, as it can be rewritten as

[TABLE]

where the subscript $E$ was included in $\tau_{E}$ for stressing the fact that, because of the external electric field, this is not the exact relaxation time of the substance. However, it is reasonably close to the true relaxation time for the present analysis. Notice that employing the Eq. 13 for the estimation of reasonable simulation times, requires the usage of the simulated (and not the experimental) relaxation times.

II.3 Uncertainty in the dielectric constants through the fluctuation route

In a previous work by the author an others2018arXiv180306421S , the uncertainty in the dielectric constant computed through the FR was found. An outline of the derivation of this uncertainty will be done, with some relevant remarks with regard to the individual treatment of each spatial direction. For simulations with typical length, the variance of the sample variance of an autocorrelated random variable, equals to two times the squared variance of the variable divided by the effective degrees of freedom, which can be asymptotically approximated by

[TABLE]

The number three in the Eqs. 1 and 3 arises by averaging over the three mutually orthogonal spatial directions. The coordinates of the dipole moment are not independent of each others. However, they are not correlated with each other in absence of external fields. Then, applying the variance operator to Eq. 3

[TABLE]

where $d$ represents the number of perpendicular directions considered, and $v_{\text{eff},i}$ the effective degrees of freedom computed for the $i$ direction. Predictably, as the three spatial directions are equivalent to each other, the standard deviation turns to be proportional to $1/\sqrt{d}$ in the long run.

In analogy to the treatment given for $n_{\text{eff}}$ ,

[TABLE]

and taking advantage of the three directions, as usual for this calculation route, the following equation is obtained for systems in which the Debye’s model is applicable

[TABLE]

II.4 Comparison of uncertainties

In terms of uncertainty, each calculation route has its advantages. They will be compared in terms of the standard deviation of the dielectric constant. As expected, uncertainties are proportional to the inverse of the square root of the total simulation time for both methods (Eqs. 13 and 18 ).

The uncertainty in the dielectric constant obtained through FR does not depends on the system volume. This very significant feature of the FR that has been missed in the literature, as far as the author is aware. The statement above is readily verified by replacing $V$ in Eq. 18, for example, according to Eq. 1.

The EFR returns uncertainties inversely proportional to the electric field and the square root of the volume. The relationship with the field strength and implications are well described in the reference kolafa2014static . The increment of the system size comes at expense of computational cost, both magnitudes are often more or less proportional to each other. Although increasing the electric field strength comes for free in terms of computational cost, it should be kept into the limits of validity of Eq. 2 when using this framework.

The uncertainties obtained with both methods can be compared by taking their quotient

[TABLE]

where it was supposed that $\tau_{E}\approx\tau$ and that, for each direction, the variance in absence of an external field, e.g. $\operatorname{Var}[{M}_{z}]_{FR}$ , equals the variance in the direction of the field of the remaining simulation ( $\operatorname{Var}[{M}]_{EFR}$ ). These approximations are valid only for weak fields.

The equation above highlights some of the expected qualitative relations between both routes. As $E$ grows, the mean dipole moment increases, and the relative importance of the uncertainty of the mean decreases. The temperature produces an opposite effect, but this is not significant here due to its change modifies the thermodynamic state for which the dielectric constant is to be computed. Nevertheless, this is significant in choosing a calculation route.

II.5 Fluctuation and external field routes combined

In the case of FR, the three mutually orthogonal spatial directions are used. Instead, only the electric field direction is considered in the EFR. In this work, the possibility of taking advantage of two discarded directions is considered. For them, the probability density function (PDF) remains symmetric. Intuitively, the PDF should be less affected by the field in those directions. Because of this, if the field is weak enough, they can be used for estimating the dielectric constant through the FR, and the this estimation can be combined with the one corresponding to the EFR for obtaining a better estimate.

The plausibility of the intuition above can be illustrated by the results of the following close related problem, which can be solved exactly: the variance of the dipole moment ( $\boldsymbol{\mu}=[\mu_{x},\mu_{y},\mu_{z}]$ ) of a rigid dipole in a bath of non-polar particles at constant temperature, which is under the influence of a constant electric field in the $z$ direction. Derivation details were omitted in order to conserve space. Solving for the $z$ direction,

[TABLE]

For $x$ and $y$ directions,

[TABLE]

It can be seen that for such idealized scenario, the change of the variance in the field direction, due to a weak electric field, is three times greater than the produced in the perpendicular directions.

II.6 An alternative calculation scheme: derivation

Below, a new approach for the calculation of the dielectric constant is presented. It can be briefly described as a three step process: 1. Perform a simulation with no external field. 2. Estimate the expected value that the dipole moment would have if an external field were applied. 3. Employ the equations for the calculation of the dielectric constant of the EFR. Although this method does not provides advantages in its direct application nor significant novelty, it simplifies the study of individual parts when the system is subdivided.

If a classical system evolves according to a hamiltonian $H_{1}$ , the expected value of an obervable $X$ satisfiestolman1979principles

[TABLE]

where $\boldsymbol{q}$ represent the state coordinates, $\Omega$ the domain of the hamiltonian, and $\beta^{-1}=k_{B}T$ . The subscript $H_{1}$ explicit which hamiltonian determines the evolution of the system. If $H_{1}$ is the hamiltonian of the system under the influence of the external field $\boldsymbol{E}$ , and $H_{0}$ is the hamiltonian of the system with no external field, then $H_{1}=H_{0}-\boldsymbol{M}\cdot\boldsymbol{E}$ griffiths2005introduction . Because of this, the expected value of the dipole moment under the influence of the electric field satisfies

[TABLE]

No approximation were made in the derivation of Eq. 26. However, large perturbations in the MD simulation leads to non representative samples.

II.7 An alternative calculation scheme: uncertainties

The uncertainties in the dielectric constant, computed using the procedure described in Section II.6, are analyzed bellow. The one-dimensional case is considered in order to facilitate the reading. As before, the uncertainties comes from the fact that the expectation can only be estimated. The problem can be stated as finding a computable expression for the variance

[TABLE]

As noted above, it is required to chose a very weak electric field in order to get a representative sample. In such case, the uncertainty in the computed mean value of the dipole moment will be significant when compared to the obtained estimate of $\operatorname{E}[\boldsymbol{M}]_{H_{1}}$ . A simple solution is to center the values of the sampled dipole moments by subtracting their mean ( $\bar{M}:=\bar{M}_{H_{0}}$ ). This change is small for large values of $n$ , as the distribution of $M$ is symmetric for the isotropic systems under study. The Eq. 28 is modified to

[TABLE]

The denominator of the above equation can be approximated to one by retaining only the first term of its expansion in Maclaurin series. This greatly simplifies the calculation of the variance.

[TABLE]

The approximation is exact up to first order, the remaining terms can be neglected by using a very weak field.

The numerator can be approximated to a convenient expression

[TABLE]

Using Taylor expansions for the variance of functions of random variables, the variance of a statistic $f(M)$ can be approximated usinghendebynonlinear

[TABLE]

then

[TABLE]

Here, the estimation of the variance of the mean must take into account that $M$ is auto-correlated, so

[TABLE]

Finally

[TABLE]

Under the suppositions considered for the equations 10 and 16, $n_{\text{eff}}=\nu_{\text{eff}}/2$ . By replacing $n_{\text{eff}}$ , the equation above transform into the one-dimentional case of Eq. 15.

The obtained result is not fortuitous, expanding in Maclaurin series the numerator of Eq. 29,

[TABLE]

By inserting this into the Eq. 4, the equation for the calculation of the dielectric constant through the one-dimentional case of FR is obtained.

The absence of the electric field in Eq. 37 is a consequence of the centering process. If this is not performed, the correct equation is

[TABLE]

which was derived in the same manner. Predictably, the additional term correspond to the right hand side of the Eq. 12.

II.8 Contribution of the individual components

The knowledge of the individual contributions to dielectric constant of the different parts of a system, may serve to increase the understanding of the interaction among them and the molecular mechanism that give rise to macroscopic properties. Their contribution should be understood in terms of the electric field that they generate, independently of how the other components influenced them. In this sense, is more natural to express the ideas in terms of the electric susceptibility ( $\chi_{e}=\varepsilon_{r}-1$ ), because it can be artificially decomposed in additive contributions, which is not possible for the dielectric constant.

Those contributions can be trivially obtained through the EFR. This is because of the linearity of the expectation operator. For example, if the system is subdivided in electrically neutral parts $A$ and $B$ ,

[TABLE]

where ${M}_{A}$ and ${M}_{B}$ represent the component of the dipole moment in the field direction of part $A$ and $B$ , respectively.

For the FR this partitioning may not be so obvious. As the variance operator is not linear, in general, $\operatorname{Var}[X_{A}+X_{B}]\neq\operatorname{Var}[X_{A}]+\operatorname{Var}[X_{B}]$ . The contribution of each component, to the variance of the total dipole moment, is obtained by summing over the respective row (or colum) of the covariance matrix. This may be understood in terms of the symmetry of the covariance operator. As expected, by summing over every component the expression for the variance of a sum of correlated variables is obtained.

The case of the proposed route of calculation can be treated as follows. Considering the linearity of the expectation operator in Eq. 26, and that the denominator in Eq. 28 is approximately one, the mean of the dipole moment of a component $K$ is given by

[TABLE]

The dipole moment of each component in the isotropic system can be centered individually, as the mean of the system is equal to the sum of the mean of its parts. Doing so helps to avoid some bias for individual components. In addition, this leads to the same expression obtained for the FR.

[TABLE]

Then, the partial electric susceptibility of the $K$ component is

[TABLE]

This expression coincides with the one expected for the FR.

II.9 Uncertainties of the individual contributions

In this section, the uncertainties in the partial susceptibilities are analyzed. For the EFR, they can obtained in a manner identical to the employed for the whole system. Then, the variance corresponding to a part of the system, for example, part $A$ is

[TABLE]

where $n_{\text{eff},A}$ is the effective number of observations computed in terms of ${M}_{A}$ .

Obtaining an equivalent expression for the FR is a cumbersome process. However, for simulations without an external electric field, the method proposed in Section II.6 can be used for the present purpose. This is because both methods returns exactly the same values as mentioned in the previous section. The advantage of this approach is that the problem is stated as obtaining the variance of sample means. If the centering process was performed over every part of the system, Eq. 42 can be expressed as

[TABLE]

which was derived using along the same lines that Eq. 33.

For the two parts system, expanding the exponential in Maclaurin series up to second order

[TABLE]

There is not prejudice in in neglecting the remaining terms of the expansion, as $E$ can be chosen arbitrarily small. Adding up every contribution and applying the expectation operator, the variance of the total dipole moment is obtained up to a constant, see Eq. 38. In other words, this verifies that the susceptibility of the whole system computed through the FR can be obtained from the individual contributions if they are individually centered.

The final expression for the uncertainty of the contribution of one part ( $K$ ) of the two parts system, can be obtained using the variable $\bar{M}_{K,H_{1}}$ in place of $\bar{M}_{A}$ in Eq. 45,

[TABLE]

where $n_{\text{eff},K,H_{1}}$ is the effective number of observations computed in terms of the variable $M_{K}^{2}+M_{K}M_{J}$ . This result can be generalized for an arbitrary number of parts,

[TABLE]

where $n_{\text{eff},K,H_{1}}$ must be evaluated for the variable $M_{K}\sum_{J}M_{J}$ .

III Numerical experiments

In order to verify the theoretical development described the previous sections, MD simulations of liquid water were carried out.

As a first step, the effects of different field strength were considered. Simulations with 512 molecules were performed using the followings values for this property: 0.0005, 0.001, 0.005, 0.01, 0.05, 0.5 and 1 V/nm. The main goal of this step is to test the validity of the equations obtained for the uncertainties, in dielectric constants of the whole system, through the EFR. It is expected for they to be satisfied within the application limit of Eq. 2. A second objective is to verify that the combination of both standard calculation routes may posses practical value. For this to be true, it is required that the FR can be applied for the directions perpendicular to the field, when the latter is large enough to allow the application of the EFR. For the systems studied, a reasonable compromise is found for 0.01 V/nm.

A second step is to analyze the volume dependence of the results. For this purpose, after finding that the electric field strengths of 0.01 V/nm allows the practical application of both routes, simulations with that field strength and the followings numbers of molecules were performed: 256, 1024, 2048, 4096 and 8192.

A third step is to analyze the calculation of the partial susceptibilities and their respective uncertainties. The system studied was again pure water. The advantage over a multicomponent system is that, for the studied system, the contributions of a given part should be more or less proportional to the number of molecules that it contains.

Finally, the relationship between the effective number of samples and the registration time step is briefly analyzed.

III.1 Simulation details

The MD simulations were performed in cubic boxes employing the Gromacs 2018.4 abraham2015gromacs program. The SPC/E modelberendsen1987missing was used for water. The simulation time used is 42 ns, of which the first 2 ns were used for equilibration.

The time step chosen is 2 fs, unless stated otherwise. Long range coulombic interactions were modeled with PMEdarden1993particle . The temperature was regulated at 298.15 K employing the Bussi-Donadio-Parrinello velocity rescaling algorithm bussi2007canonical . The simulation equilibrium pressure was set to 1 bar using the Berendsen’s barostatberendsen1984molecular . Notice that the applicability of the proposed methods is independent of most simulation parameters.

The posterior analysis was carried out using own Python3 van2014python routines based on Numpyvan2011numpy and MDTraj libraries mcgibbon2015mdtraj .

III.2 Results and discussion

In order to verify that, in terms of absolute value, the proposed methods provide reasonable uncertainties, they must be compared with any other procedure. A simple alternative used in this work is explained below.

Consider that the trajectory is evenly split in $p$ contiguous parts of length $\mathcal{T}/p$ . According to the equations 13 and 18, the standard deviation in the computed dielectric constant is proportional to $\sqrt{p}$ . For each $p$ in $\{p\in\mathbb{N}|1<p\leq 100\}$ , the dielectric constant of each part can be computed and the standard deviation, of the corresponding $p$ parts, calculated. An estimation of the standard deviation of the dielectric constant, corresponding to the the whole trajectory, is obtained by extrapolating to $p=1$ . For the extrapolation, a linear model without the constant term was used. The fitting procedure was performed by means of ordinary least squares.

The upper part of Table 1 contains results for the whole system obtained for different values of electric field strength. The estimates corresponding to the EFR can be found in the second column. As expected, for large fields there is a systematic error that decreases the estimate of the dielectric constant value. For these numerical experiments, it starts to be noticeable for electric field of about 0.1 V/nm. The error estimations obtained with the proposed method (Eq. 12) and the numerical method (NM) are found in the third and fourth columns. They are also plotted in Figure 1. This log-log representation present a linear behavior for low fields.

In order to compare how well the data obtained is represented by the linear model, the coefficient of determination ( $R^{2}$ ) of ordinary least square linear regressions of the logarithms of plotted data were calculated. The last two points were removed because no linear behavior is expected for those strong fields. The coefficients of determination obtained are 0.9975 and 0.9995 for the proposed and the numerical method, respectively. This mean that the linear model is a good representation of the data.

The fifth column of the top side of Table 1, contains the estimations of the dielectric constant obtained employing the FR for the $x$ and $y$ directions, in which the electric field is null. Their estimated uncertainties are found in the sixth and seventh columns, and were computed using Eq. 15 and the numerical method, respectively. It can be seen that, in the cases in which the dielectric constant grows linearly with the field, the FR is applicable for the directions in which the electric field is null. This behavior extends to almost 0.1 V/nm. In order to ensure the applicability of the FR, in the following, the value of 0.01 V/nm was used to study the volume dependence of the uncertainties.

The lower part of Table 1 contains the results of simulations performed with the mentioned field, and different numbers of molecules, ranging from 256 to 8192. Figure 2 is a log-log representation of the estimated standard deviations for the dielectric constant computed through the EFR, vs. the number of molecules in the system. For this case, the determination coefficient is 0.9978 for the NM, and 0.9918 for the proposed method (Eq. 12).

The standard deviations of the sixth and seventh columns allow to confirm that the volume of the simulation box does not affect the uncertainties in the FR framework. Although well conducted research was done comparing the FR and the EFRkolafa2014static ; riniker2011calculation or analyzing the behavior of the FR with system sizegereben2011accurate , this fact has not been made explicit nor suggested. However, the results reported in ref. gereben2011accurate clearly support that finding. For the present purposes, it is assumed that the system is large enough for considering that the effects of the periodic boundary conditions are negligible.

The eighth and ninth columns contain the final estimates of the mean and standard deviation, respectively. They were obtained by employing the minimum-variance unbiased estimator (MVUE)hartung2011statistical from the results of both calculation routes. The uncertainties employed are those obtained using the proposed methods. For the sake of completeness, the MVUE was employed in conjunction with the values of all those simulations containing more than 256 molecules, and in which the electric field strength was lower or equal to $0.01$ V/nm. The final result of this estimation is 70.46 $\pm$ 0.31. No systematic error due to system size is expected, as 256 water molecules suffices for reaching the thermodynamic limit in the studied casesnymand2000molecular .

In the interests of accuracy, some comments on the effects of the centering process will be made. Both the estimator of Eq. 3 and the one obtained by removing $\bar{\boldsymbol{M}}^{2}$ from it have the same expectation for the studied systems. This is because for them $\text{E}[{\boldsymbol{M}}]=\boldsymbol{0}$ . However, the estimator of Eq. 3 is biased. The above was expressed in other terms in Ref.kolafa2014static . The estimator of Eq. 3 has its advantages too; it is invariant under the centering process while the other estimator is not. This suggest that the estimator of Eq. 3 fits better in the framework of the present study. Despite it is biased, a simple modification allows to build an unbiased estimatorzikeba2010effective

[TABLE]

The inclusion of the quotient $n_{\text{eff}}/(n_{\text{eff}}-1)$ is not significant in practical terms because, in normal cases, it is very close to one. For example, a value close to 1.0007 is found for the systems simulated in this work. Then, the estimator of Eq. 3 is perfectly suitable for regular studies. The latter and its unbiased variant were interchangeably used for the derivation of Eq. 15.

The Table 2 contains the results for the partial susceptibilities, and their respective uncertainties. The system used for the analysis was the one containing 1024 molecules under the influence of an electric field of 0.01 V/nm. For each analyzed case, the system was subdivided into two parts, and one of them was used for the computations. The percentage of the molecules of the system belonging to the analyzed part is tabulated in the first column on the left. When it does not coincide with a natural number of molecules, the last one was truncated. It was verified that the sum of the contributions of both parts coincide with the susceptibility of the whole system. However, it is not true for the uncertainties. The remaining columns contain the calculated contributions and uncertainties. Please, see the caption of Table 2 for details.

For a given percentage, many different subsets of molecules could be chosen. In each case, only one subset was used. The results confirm that for both routes of calculation, the fraction of molecules considered approximately equals to the fraction of the susceptibility of the whole system in which those contribute.

Another simulation was performed in order to analyze the relationship between the effective number of samples and the registration time step. The simulated system contains 1024 water molecules in presence of an electric field of 0.01 V/m. All the parameters remain the same except that a registration time step of 0.04 ps was used. This implies that the number of samples employed is $10^{6}$ . Only one direction in which the electric field is null was considered. Many pairs $(n,n_{eff})$ were obtained after sub-sampling. Each sub-sample was obtained retaining only those samples whose corresponding ordinals are multiples of some integer. This elemental procedure retrieves samples which may have arisen by using different registration time steps on the same simulation. The pairs obtained are depicted in Figure 3. For comparison, the right hand side of Eq. 11 was plotted. The required relaxation time (12.6 ps) was computed using the whole sample. Predictably, the effective number of samples increases as the registration time step decreases, and there is not advantage in using registration time steps much shorter than the relaxation time. For similar systems, the authors of ref. gereben2011accurate found that registration time steps up to around 2 ps are adequate, but a value of 20 ps is not. This is consistent with the analysis and the results presented above. The Eq. 11 can be used for selecting a sampling frequency for the selected system, as can be seen by simple inspection of Figure 3. This is due to the relaxation process of the studied systems is well described by the Debye’s model.

IV Summary

Two routes for the calculation of the static dielectric constants were studied and compared. One of them is the FR, which consist in sampling the dipole moment of the system in absence of external electric fields, and then using the Eq. 3. The expression for the variance of this estimator, previously found by the author and others, was written considering an arbitrary number of directions (Eq. 15) and, assuming the applicability of the Debye’s model, as a function of the relaxation and simulation times (Eq. 18) . The other route of calculation is the EFR, which requires the application of a constant and external electric field on the system, and employing the estimator of Eq. 4. The variance for this estimator was derived (equations 12 and 13).

The standard deviation of the dielectric constant, obtained through the EFR, is inversely proportional to the electric field. In this work, it was found that the uncertainty does not depends on the volume for the FR, while the dependence is with the inverse of its square root for the EFR.

An alternative calculation scheme was presented. It consist in performing a simulation with no external field, and using Eq. 28 for the prediction of the mean dipole moment that should have an equivalent simulation, in which the electric field were applied. Then, the result of Eq. 28 is used in conjunction with Eq. 4 for the calculation of the dielectric constant. When these equations are used after centering the values of the dipole moment, by subtracting its mean, the uncertainty in the dielectric constant is given by Eq. 37. The equations above can be applied independently for each spatial direction in order to lowering the uncertainty. It was also shown that this method reduces to the FR. This fact was used to intuitively show how to compute the partial electric susceptibilities through the FR (Eq. 44). Also, their uncertainties were studied for the EFR and the proposed route. It was found that they are represented by Eq. 45 and Eq. 48, respectively.

The Eq. 9 is useful for choosing the sampling frequency provided that a small sample is previously available. When the Debye’s model is suitable, the Eq. 11 may be used for that purpose. The latter does not requires a previous simulation. One note of caution, however: the relaxation time included in Eq. 11, which refers to the one obtained through simulation, may differ significantly from the experimental one.

Te best estimate of the dielectric constant obtained in this work for pure liquid water (SPC/E) at 298.15 K and 1 bar is 70.46 $\pm$ 0.31.

Bibliography34

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) Arthur R von Hippel and SO Morgan. Dielectric materials and applications. Journal of The Electrochemical Society , 102(3):68C–68C, 1955.
2(2) Koji Asami. Characterization of heterogeneous systems by dielectric spectroscopy. Progress in polymer science , 27(8):1617–1659, 2002.
3(3) Berni J Alder and T‌ E Wainwright. Studies in molecular dynamics. i. general method. The Journal of Chemical Physics , 31(2):459–466, 1959.
4(4) Dennis C Rapaport and Dennis C Rapaport Rapaport. The art of molecular dynamics simulation . Cambridge university press, 2004.
5(5) L Levrel and AC Maggs. Boundary conditions in local electrostatics algorithms. The Journal of chemical physics , 128(21):214103, 2008.
6(6) SW De Leeuw, John W Perram, and ER Smith. Computer simulation of the static dielectric constant of systems with permanent electric dipoles. Annual review of physical chemistry , 37(1):245–270, 1986.
7(7) Tim N Heinz, Wilfred F van Gunsteren, and Philippe H Hünenberger. Comparison of four methods to compute the dielectric permittivity of liquids from molecular dynamics simulations. The Journal of chemical physics , 115(3):1125–1136, 2001.
8(8) Martin Neumann. Dipole moment fluctuation formulas in computer simulations of polar systems. Molecular Physics , 50(4):841–858, 1983.