# On the reaction field for interaction site models of polar systems

## 1 Introduction

The calculation of dielectric quantities by computer experiment requires an explicit consideration of effects associated with the truncation of long-range interactions. The concrete success in this direction has been achieved within the reaction field (RF) geometry \[1–5\]. As a result, computer adapted dielectric theories have been proposed \[6–10\]. In the framework of these theories, a bulk dielectric constant can be determined on the basis of a fluctuation formula via correlations obtained in simulations for finite samples. However, main attention in the previous investigations has been focused on polar systems with the point dipole interaction. As is now well established, the model of point dipoles can not reproduce adequately features of real polar liquids.

At the same time, attempts to apply the RF geometry for more realistic interaction site (IS) models have also been made \[11–13\]. However, acting within a semiphenomenological approach, it was not understood how to perform the truncation of intermolecular potentials. As a consequence, the molecular cut-off and the usual point dipole RF (PDRF) have been assumed. Obviously, such an approach includes effects connected with finiteness of the molecule inconsistently. Indeed, the interdipolar potential is replaced by site-site Coulomb interactions, whereas the RF is remained in its usual form. An additional complication for IS models consists in a spatial distribution of charges and this fact is not taken into account by the standard PDRF geometry.

In the present paper we propose two alternative approaches to remedy this situation. The first one follows from the usual fluctuation formula which is constructed, however, on the microscopic operator of polarization density for IS models. This leads to an ISRF geometry, where the cut-off radius is applied with respect to individual charges rather than to the molecule as a whole. Nevertheless, the molecular cut-off scheme can also be acceptable, but the reaction field together with the fluctuation formula need to be corrected. In the second approach a molecular RF (MRF) geometry is proposed and a new quadrupole term is identified. On the basis of a MCY water model we show that uncertainties of the dielectric quantities can be significant if the standard PDRF geometry is used in computer simulations.

## 2 Interaction site reaction field

We consider an isotropic, classical system of $`N`$ identical molecules enclosed in volume $`V`$. The microscopic electrostatic field created by the molecules at point $`𝒓V`$ is equal to

$$\widehat{𝑬}(𝒓)=\underset{i=1}{\overset{N}{}}\underset{a}{}q_\stackrel{}{a}\frac{𝒓𝒓_i^a}{|𝒓𝒓_i^a|^3}=\underset{V}{}𝑳(𝒓𝒓^{})\widehat{Q}(𝒓^{})d𝒓^{},$$
(1)   
where $`𝒓_i^a`$ denotes the position for charge $`q_\stackrel{}{a}`$ of $`i`$th molecule, $`\widehat{Q}(𝒓)=_{i,a}q_\stackrel{}{a}\delta (𝒓𝒓_i^a)`$ is the microscopic operator of charge density, $`𝑳(𝝆)=\mathbf{}1/\rho `$ and the summation extends over all molecules and charged sites. For the investigation of dielectric properties, it is more convenient to rewrite the electric field (1) in the polarization representation

$$\widehat{𝑬}(𝒓)=\underset{V}{}𝐓(𝒓𝒓^{})\widehat{𝑷}(𝒓^{})d𝒓^{}=\frac{4\pi }{3}\widehat{𝑷}(𝒓)+\underset{\rho +0}{lim}\underset{\rho <|𝒓𝒓^{}|}{\underset{V}{}}𝐓(𝒓𝒓^{})\widehat{𝑷}(𝒓^{})\mathrm{d}𝒓^{}.$$
(2)   
Here $`𝐓(𝝆)=\mathbf{}\mathbf{}1/\rho `$ is the dipole-dipole tensor, $`\widehat{𝑷}(𝒓)`$ denotes the microscopic operator of polarization density, defined as $`\mathbf{}\mathbf{}\widehat{𝑷}(𝒓)=\widehat{Q}(𝒓)`$, and the singularity $`lim_{\rho 0}𝐓(𝝆)=4\pi /3\delta (𝝆)𝐈`$ has been avoided, where $`𝐈`$ is the unit tensor of the second rank. The both charge (1) and polarization (2) representations are equivalent and applicable for infinite ($`N,V\mathrm{}`$) systems.

In simulations, which deal with finite samples, the sum (1) can not be calculated exactly taking into account an infinitely large number of terms. Therefore, we must restrict ourselves to a finite set of terms in (1) or to a finite range of the integration in (1) and (2) for which $`|𝒓𝒓^{}|R`$, where $`R`$ is a cut-off radius. Now the following problem appears. How to estimate the cut-off field caused by the integration over unaccessible region $`|𝒓𝒓^{}|>R`$? The solution of this problem has been found for the first time for systems with point dipoles in the RF geometry. The result for conducting boundary conditions is 

$$\widehat{𝑬}(𝒓)\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒓)=\frac{4\pi }{3}\widehat{𝑷}(𝒓)+\underset{\rho +0}{lim}\underset{\rho <|𝒓𝒓^{}|R}{\underset{V,\mathrm{tbc}}{}}\left(𝐓(𝒓𝒓^{})+\frac{𝐈}{R^3}\right)\widehat{𝑷}(𝒓^{})\mathrm{d}𝒓^{},$$
(3)   
where a cubic finite sample and toroidal boundary conditions (TBC) have been used, so that $`R\sqrt[3]{V}/2`$. The additional term $`𝐈/R^3`$ in the right-hand site of (3) describes the RF which is used for an approximation of the real cut-off field. For a pure spherical cut-off (SC) without the RF correction, we have $`\widehat{𝑬}_\stackrel{}{R}^{\mathrm{SC}}(𝒓)={\displaystyle \gamma (|𝒓𝒓^{}|)𝑳(𝒓𝒓^{})\widehat{Q}(𝒓^{})d𝒓^{}}`$, where $`\gamma (\rho )=1`$ if $`\rho R`$ and $`\gamma (\rho )=0`$ otherwise. Obviously, that $`lim_R\mathrm{}\widehat{𝑬}_\stackrel{}{R}^{\mathrm{SC}}(𝒓)=lim_R\mathrm{}\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒓)=\widehat{𝑬}(𝒓)`$.

Let us perform the spatial Fourier transform $`(𝒌)=d𝒓\text{e}^{\mathrm{i}𝒌\mathbf{}𝒓}(𝒓)`$ for arbitrary functions $``$. Then one obtains

$$\widehat{𝑬}_\stackrel{}{R}^{\mathrm{SC}}(𝒌)=𝑳(𝒌)\widehat{Q}(𝒌),\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒌)=\frac{4\pi }{3}\widehat{𝑷}(𝒌)+\left(𝐓(𝒌)+4\pi \frac{j_\stackrel{}{1}(kR)}{kR}𝐈\right)\widehat{𝑷}(𝒌),$$
(4)   
where

$$𝑳(𝒌)=4\pi \left(1j_\stackrel{}{0}(kR)\right)\frac{\mathrm{i}𝒌}{k^2},𝐓(𝒌)=\frac{4\pi }{3}\left(13\frac{j_\stackrel{}{1}(kR)}{kR}\right)\left(3\widehat{𝒌}\widehat{𝒌}𝐈\right),$$
(5)   
$`\widehat{Q}(𝒌)=_{i,a}q_\stackrel{}{a}\text{e}^{\mathrm{i}𝒌\mathbf{}𝒓_i^a}=\mathrm{i}𝒌\mathbf{}\widehat{𝑷}(𝒌)`$, $`𝒌=2\pi 𝒏/\sqrt[3]{V}`$ is one of the allowed wavevectors of the reciprocal lattice, $`𝒏`$ designates a vector with integer components, $`k=|𝒌|`$, $`\widehat{𝒌}=𝒌/k`$ and $`j_\stackrel{}{0}(z)=\mathrm{sin}(z)/z`$, $`j_\stackrel{}{1}(z)=\mathrm{cos}(z)/z+\mathrm{sin}(z)/z^2`$ are the spherical Bessel functions of zero and first order, respectively. In view of (5), the relations (4) transform into

$$\widehat{𝑬}_\stackrel{}{R}^{\mathrm{SC}}(𝒌)=4\pi \left(1j_\stackrel{}{0}(kR)\right)\widehat{𝑷}_\mathrm{L}(𝒌),\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒌)=4\pi \left(13\frac{j_\stackrel{}{1}(kR)}{kR}\right)\widehat{𝑷}_\mathrm{L}(𝒌),$$
(6)   
where $`\widehat{𝑷}_\mathrm{L}(𝒌)=\widehat{𝒌}\widehat{𝒌}\mathbf{}\widehat{𝑷}(𝒌)=\mathrm{i}𝒌\widehat{Q}(𝒌)/k^2`$ is the longitudinal component of the microscopic operator of polarization density.

It is easy to see from (6) that the both functions $`\widehat{𝑬}_\stackrel{}{R}^{\mathrm{SC}}(𝒌)`$ and $`\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒌)`$ tend to the same value $`\widehat{𝑬}(𝒌)=4\pi \widehat{𝑷}_\mathrm{L}(𝒌)`$ of the infinite system at $`R\mathrm{}`$ ($`k0`$). However, the results converge as $`R^1`$ for the pure SC scheme, while as $`R^2`$ in the RF geometry, i.e., more quickly, because a main part of the truncation effects is taken into account by the RF. This is very important in our case, where we hope to reproduce features of infinite systems on the basis of finite samples. That is why the pure truncation, which is standard for simple fluids with short-range potentials, is generally not recommended for polar systems with long-range nature of the dipolar interaction. The influence of the TBC and the difference between micro- and canonical ensembles are of order $`N^1R^3`$  and, therefore, they can be excluded from our consideration. It is worth mentioning that electrostatic fields are pure longitudinal. They can be defined via the longitudinal component of the microscopic operator of polarization density, that is confirmed by Eq. (6).

Let us enclose the system in an external electrostatic field $`𝑬_\stackrel{}{0}(𝒓)`$. The material relation between the macroscopic polarization $`𝑷_\mathrm{L}(𝒌)=\widehat{𝑷}_\mathrm{L}(𝒌)`$ in the weak external field and total macroscopic field is $`4\pi 𝑷_\mathrm{L}(𝒌)=\left(\epsilon _\stackrel{}{\mathrm{L}}(k)1\right)𝑬_\mathrm{L}(𝒌)`$, where $`\epsilon _\stackrel{}{\mathrm{L}}(k)`$ denotes the longitudinal wavevector-dependent dielectric constant. Applying the first-order perturbation theory with respect to $`𝑬_\stackrel{}{0}`$ yields for rigid molecules $`Vk_\mathrm{B}T𝑷_\mathrm{L}(𝒌)=\widehat{𝑷}_\mathrm{L}(𝒌)\mathbf{}\widehat{𝑷}_\mathrm{L}(𝒌)_\stackrel{}{0}𝑬_\stackrel{}{0}(𝒌)`$, where $`k_\mathrm{B}`$ and $`T`$ are Boltzmann’s constant and temperature, respectively, and $`\mathrm{}_\stackrel{}{0}`$ is the equilibrium average in the absence of the external field. Then, taking into account that $`𝑬_\mathrm{L}(𝒌)=𝑬_\stackrel{}{0}(𝒌)+\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒌)`$ and eliminating $`𝑬_\stackrel{}{0}(𝒌)`$, we obtain the fluctuation formula

$$\frac{\epsilon _\stackrel{}{\mathrm{L}}(k)1}{\epsilon _\stackrel{}{\mathrm{L}}(k)}=\frac{9yG_\stackrel{}{\mathrm{L}}(k)}{1+27yG_\stackrel{}{\mathrm{L}}(k)j_\stackrel{}{1}(kR)/(kR)}=9yg_\stackrel{}{\mathrm{L}}(k).$$
(7)   
Here $`G_\stackrel{}{\mathrm{L}}(k)=\widehat{𝑷}_\mathrm{L}(𝒌)\mathbf{}\widehat{𝑷}_\mathrm{L}(𝒌)_\stackrel{}{0}/N\mu ^2`$ is the longitudinal component of the finite-system wavevector-dependent Kirkwood factor, $`y=4\pi N\mu ^2/9Vk_\mathrm{B}T`$ and $`\mu =|𝝁_i|=|_aq_\stackrel{}{a}𝒓_i^a|`$ denotes the permanent magnitude of molecule’s dipole moment. It is necessary to note that we consider rigid IS molecules so that effects associated with molecular and electronic polarizabilities are not included in our investigation. In the case of $`R\mathrm{}`$, we have $`j_\stackrel{}{1}(kR)/(kR)0`$ and computer adapted formula (7) reduces to the well-known fluctuation formula for macroscopic systems in terms of the infinite-system Kirkwood factor $`g_\stackrel{}{\mathrm{L}}(k)=lim_R\mathrm{}G_\stackrel{}{\mathrm{L}}(k)`$.

As was mentioned earlier, the electric field $`\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}`$ in the form (3), (4) as well as the fluctuation formula (7) have been proposed for the first time to investigate polar systems of point dipoles . However, acting within a semiphenomenological framework, it was not understood how to perform the truncation of the intermolecular potential $`\phi _{\stackrel{}{ij}}`$ at attempts to extend this formula for IS models. As a result, the molecular cut-off $`r_{ij}=|𝒓_i𝒓_j|R`$, where $`𝒓_i`$ is the center of mass for $`i`$th molecule, and the usual PDRF have been suggested \[11–13\]:

$$\phi _{\stackrel{}{ij}}=\underset{a,b}{}\frac{q_\stackrel{}{a}q_\stackrel{}{b}}{|𝒓_i^a𝒓_j^b|}\frac{𝝁_i\mathbf{}𝝁_j}{R^3},r_{ij}R.$$
(8)   
It is essentially to emphasize that the fluctuation formula (7) takes into account finiteness of the system explicitly by the factor $`j_\stackrel{}{1}(kR)/(kR)`$. As a result, if the system size is sufficiently large (terms of order $`R^2`$ can be neglected), the bulk ($`N,V\mathrm{}`$) dielectric constant can be reproduced via the finite-system Kirkwood factor $`G_\stackrel{}{\mathrm{L}}(k)`$ which depends on $`R`$ in a characteristic way. However, to achieve this self-consistency in the evaluation of the bulk dielectric constant, the equilibrium averaging in $`G_\stackrel{}{\mathrm{L}}(k)`$ must be calculated for systems with the intermolecular potential which leads exactly to the microscopic electric field $`\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒓)`$ (3). As we shall below, the intermolecular potential (8) does not obey this condition.

To derive the exact intermolecular potential in the charge representation, we perform the inverse Fourier transform $`\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒓)=\frac{1}{(2\pi )^3}{\displaystyle d𝒌\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒌)\text{e}^{\mathrm{i}𝒌\mathbf{}𝒓}}`$ and obtain using (6)

$$\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒓)=\underset{i,a}{}q_\stackrel{}{a}\frac{𝒓𝒓_i^a}{|𝒓𝒓_i^a|^3}\left(1\frac{6}{\pi }\frac{|𝒓𝒓_i^a|^2}{R}\underset{0}{\overset{\mathrm{}}{}}j_\stackrel{}{1}(kR)j_\stackrel{}{1}(k|𝒓𝒓_i^a|)dk\right).$$
(9)   
Taking into account that $`\frac{6}{\pi }_0^{\mathrm{}}j_\stackrel{}{1}(kR)j_\stackrel{}{1}(k\rho )dk=\rho /R^2`$ if $`\rho R`$ and is equal to $`R/\rho ^2`$ if $`\rho >R`$, we have

$$\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒓)=\underset{i,a}{}q_\stackrel{}{a}\frac{𝒓𝒓_i^a}{|𝒓𝒓_i^a|^3}\left(1\frac{|𝒓𝒓_i^a|^3}{R^3}\right)\mathrm{if}|𝒓𝒓_i^a|R$$
(10)   
and $`\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒓)=0`$ otherwise, where the first term in the right-hand side is the Coulomb field, while the second contribution corresponds to the RF in the IS description.

In order to understand nature of this field, we consider a spherical cavity of radius $`R`$ with the center at point $`𝒓`$, embedded in an infinite conducting medium. Let us place a point charge $`q_\stackrel{}{a}`$ at point $`𝒓_i^a`$ in the cavity, so that $`|𝒓𝒓_i^a|R`$. The total electric field $`𝒆_i^a(𝒓)`$ at point $`𝒓`$ consists of the field due to the charge $`q_\stackrel{}{a}`$ and the field created by induced charges located on the surface of the cavity. According to the method of electrostatic images , this last field can be presented as the field of an imaginary charge $`q_\stackrel{}{a}^{}=q_\stackrel{}{a}R/|𝒓𝒓_i^a|`$ which is located at point $`𝒓_{}^{}{}_{i}{}^{a}=𝒓R^2(𝒓𝒓_i^a)/|𝒓𝒓_i^a|^2`$ outside the sphere. Then $`𝒆_i^a(𝒓)=q_\stackrel{}{a}(𝒓𝒓_i^a)/|𝒓𝒓_i^a|^3+q_\stackrel{}{a}^{}(𝒓𝒓_{}^{}{}_{i}{}^{a})/|𝒓𝒓_{}^{}{}_{i}{}^{a}|^3=q_\stackrel{}{a}(𝒓𝒓_i^a)(1/|𝒓𝒓_i^a|^31/R^3)`$ that is completely in line with the term of sum (10).

In the potential representation ($`\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒓)=\mathbf{}\mathrm{\Phi }(𝒓)`$), we obtain $`\mathrm{\Phi }(𝒓)=_{i,a}\varphi _i^a(𝒓)`$, where $`\varphi _i^a(𝒓)=q_\stackrel{}{a}(1/\rho _i^a+\frac{1}{2}\rho _{i}^{a}{}_{}{}^{2}/R^3+C)`$, $`\rho _i^a=|𝒓𝒓_i^a|`$ and $`C`$ is, in general, an arbitrary constant which for infinite systems is chosen as $`\varphi _i^a|_{\rho _i^a\mathrm{}}=0`$. In our case, according to the toroidal boundary conventional, $`\varphi _i^a|_{\rho _i^a=R}=0`$ whence $`C=3/2R^1`$. Then the intermolecular potential of interaction is $`\phi _{\stackrel{}{ij}}=_{a,b}q_\stackrel{}{b}\varphi _i^a(𝒓_j^b)=_{a,b}q_\stackrel{}{a}\varphi _j^b(𝒓_i^a)=_{a,b}\phi _{\stackrel{}{ij}}^{ab}`$, where

$$\phi _{\stackrel{}{ij}}^{ab}=\{\begin{array}{ccc}q_\stackrel{}{a}q_\stackrel{}{b}\left(\frac{1}{|𝒓_i^a𝒓_j^b|}+\frac{1}{2}\frac{|𝒓_i^a𝒓_j^b|^2}{R^3}\frac{3}{2R}\right)& ,\hfill & |𝒓_i^a𝒓_j^b|R\\ 0& ,\hfill & |𝒓_i^a𝒓_j^b|>R\end{array}$$
(11)   
and the site-site cut-off is performed.

It is easily seen from (11) that the ISRF part $`\frac{1}{2}_{a,b}q_\stackrel{}{a}q_\stackrel{}{b}|𝒓_i^a𝒓_j^b|^2/R^3`$ transforms into the usual form $`𝝁_i\mathbf{}𝝁_j/R^3`$ of point dipoles for $`r_{ij}Rd`$ only, where $`d=2\mathrm{max}|𝜹_i^a|`$ is the diameter of the molecule and $`𝜹_i^a=𝒓_i^a𝒓_i`$. In the case if the molecular rather than the site-site cut-off is applied to the potential (11), this transformation is valid for arbitrary $`r_{ij}R`$. Moreover, in the last case the constant $`C=3/2R^1`$ is canceled owing electroneutrality ($`_aq_\stackrel{}{a}=0`$) of the molecule and we recover the result (8) of previous work . However, the potential of interaction (11) corresponds completely to the conditions at which the fluctuation formula (7) is derived. Therefore, this potential, instead of (8), must be used in simulations to obtain a correct value for the dielectric constant.

## 3 Molecular reaction field

In the case of point dipoles, where $`d+0`$, $`q_\stackrel{}{a}\mathrm{}`$ provided $`\mu const`$, both (8) and (11) representations are identical and reduced to the well-known result

$$\phi _{\stackrel{}{ij}}=𝝁_i\mathbf{}𝐓(𝒓_{ij})\mathbf{}𝝁_j\frac{𝝁_i\mathbf{}𝝁_j}{R^3},r_{ij}R$$
(12)   
for the interdipolar interaction in the RF geometry. It is easy to see that in the case of IS models, the intermolecular potential (8) takes into account effects associated with finiteness of the molecule inconsistently. For example, the interdipolar potential is replaced by the real site-site Coulomb ones, whereas the reaction field is remained in its usual form of point dipoles. From this point of view a natural question of how to improve the RF within the molecular cut-off scheme arises. The simplest way to solve this problem lies in the following.

Let us consider the mentioned above spherical cavity, centered now at some fixed point $`𝒓_\stackrel{}{0}`$, in the infinite conducting medium. We place an $`i`$th molecule in such a way that all sites of the molecule would be located in the cavity. This condition is fulfilled providing $`|𝒓_i𝒓_\stackrel{}{0}|R_dRd/2`$. The potential of a molecular reaction field at point $`𝒓`$ belonging the cavity can be presented, according to the method of electrostatic images, as

$$\phi _i^{\mathrm{RF}}(𝒓)=\underset{a}{}\frac{q_\stackrel{}{a}^{}}{|𝝆𝝆_{}^{}{}_{i}{}^{a}|}=\underset{a}{}\frac{q_\stackrel{}{a}R/\rho _i^a}{\left|𝝆\left({\displaystyle \frac{R}{\rho _i^a}}\right)^2𝝆_i^a\right|}=\underset{a}{}\frac{q_\stackrel{}{a}}{\left|{\displaystyle \frac{\rho _i^a}{R}}𝝆{\displaystyle \frac{R}{\rho _i^a}}𝝆_i^a\right|},$$
(13)   
where $`𝝆=𝒓𝒓_\stackrel{}{0}`$ and $`𝝆_i^a=𝒓_i^a𝒓_\stackrel{}{0}`$. Differentiating (13) over $`𝒓`$ at point $`𝒓_\stackrel{}{0}`$ yields

$$\frac{\phi _i^{\mathrm{RF}}(𝒓)}{𝒓}|_{𝒓_0}=\frac{𝝁_i}{R^3},\frac{^2\phi _i^{\mathrm{RF}}(𝒓)}{𝒓𝒓}|_{𝒓_0}=\frac{𝐪_i^{𝒓_0}}{R^5},\frac{^3\phi _i^{\mathrm{RF}}(𝒓)}{𝒓𝒓𝒓}|_{𝒓_0}=\frac{𝐠_i^{𝒓_0}}{R^7},\mathrm{}$$
(14)   
Here $`𝝁_i=_aq_\stackrel{}{a}𝝆_i^a=_aq_\stackrel{}{a}𝜹_i^a`$ is the dipole moment of $`i`$th molecule, which does not depend on $`𝒓_\stackrel{}{0}`$ owing electroneutrality of the molecule, while $`𝐪_i^{𝒓_0}=_aq_\stackrel{}{a}(3𝝆_i^a𝝆_i^a\rho _{i}^{a}{}_{}{}^{2}𝐈)`$ and $`𝐠_i^{𝒓_0}`$ are the tensors of quadrupole and octupole moments, correspondingly, of $`i`$th molecule with respect to $`𝒓_\stackrel{}{0}`$. The third rank tensor $`𝐠_i^{𝒓_0}`$ has the following components $`𝐠_{i}^{𝒓_0}{}_{\alpha \beta \gamma }{}^{}=3_aq_\stackrel{}{a}\left(5𝝆_{i}^{a}{}_{\alpha }{}^{}𝝆_{i}^{a}{}_{\beta }{}^{}𝝆_{i}^{a}{}_{\gamma }{}^{}\rho _{i}^{a}{}_{}{}^{2}(𝝆_{i}^{a}{}_{\alpha }{}^{}\delta _{\beta \gamma }+𝝆_{i}^{a}{}_{\beta }{}^{}\delta _{\alpha \gamma }+𝝆_{i}^{a}{}_{\gamma }{}^{}\delta _{\alpha \beta })\right)`$. It is more convenient to present multipoles of higher order with respect to the molecular center of mass. For the tensor of quadrupole moment we obtain $`𝐪_i^{𝒓_0}=𝐪_i+𝐰_i`$, where $`𝐪_i=_aq_\stackrel{}{a}(3𝜹_i^a𝜹_i^a\delta _{i}^{a}{}_{}{}^{2}𝐈)`$ is the tensor of quadrupole moment of $`i`$th molecule with respect to its center of mass, $`𝐰_i=3(𝝁_i𝝆_i+𝝆_i𝝁_i)2𝝁_i\mathbf{}𝝆_i𝐈`$ and $`𝝆_i=𝒓_i𝒓_\stackrel{}{0}`$. It is necessary to underline that tensor $`𝐪_i`$ is split into dynamical $`𝝎_i=_aq_\stackrel{}{a}𝜹_i^a𝜹_i^a`$ and conservative $`_aq_\stackrel{}{a}\delta _{i}^{a}{}_{}{}^{2}𝐈`$ parts for rigid molecules.

Putting $`𝒓_\stackrel{}{0}=𝒓_j`$ and assuming $`dR`$, we obtain the energy of $`j`$th molecule in the MRF of $`i`$th molecule

$$\varphi _{ji}^{\mathrm{RF}}=𝝁_j\mathbf{}\frac{\phi _i^{\mathrm{RF}}(𝒓)}{𝒓}|_{𝒓_j}+\frac{1}{6}𝐪_j\mathbf{:}\frac{^2\phi _i^{\mathrm{RF}}(𝒓)}{𝒓𝒓}|_{𝒓_j}+\mathrm{}=\frac{𝝁_j\mathbf{}𝝁_i}{R^3}\frac{1}{6}\frac{𝐪_j\mathbf{:}𝐪_i^{𝒓_j}}{R^5}+\mathrm{},$$
(15)   
where multipoles of higher order have been neglected. Finally, using the RF potential $`\phi _{ij}^{\mathrm{RF}}=(\varphi _{ij}^{\mathrm{RF}}+\varphi _{ji}^{\mathrm{RF}})/2`$ yields the desired intermolecular potential

$$\phi _{\stackrel{}{ij}}=\{\begin{array}{ccc}\underset{a,b}{}\frac{q_\stackrel{}{a}q_\stackrel{}{b}}{|𝒓_i^a𝒓_j^b|}\frac{𝝁_i\mathbf{}𝝁_j}{R^3}\frac{𝐪_i\mathbf{:}𝐪_j3(𝐪_i\mathbf{:}𝝁_j𝒓_{ij}+𝐪_j\mathbf{:}𝝁_i𝒓_{ji})}{6R^5}& ,\hfill & r_{ij}R_d\\ 0& ,\hfill & r_{ij}>R_d\end{array}$$
(16)   
where equality $`𝐪\mathbf{:}𝐈=0`$ has been used.

The total reaction field, created by all molecules at point $`𝒓`$ near $`𝒓_\stackrel{}{0}`$ is

$$𝑬_{\stackrel{}{\mathrm{RF}}}(𝒓)=\underset{i}{\overset{\rho _iR_d}{}}\frac{\phi _i^{\mathrm{RF}}(𝒓)}{𝒓}=\frac{𝑴(R_d)}{R^3}+\frac{𝐐(R_d)+𝐖(R_d)}{R^5}𝝆+\mathrm{},$$
(17)   
where $`𝑴(R_d)=_i^{\rho _iR_d}𝝁_i`$ and $`𝐐(R_d)=_i^{\rho _iR_d}𝐪_i`$ denote the total dipole and own quadrupole moment, respectively, within the sphere of radius $`R_d`$ and $`𝐖(R_d)=_i^{\rho _iR_d}𝐰_i`$. In the case of point dipoles, we have $`R_dR`$, $`𝐪_i,𝐠_i,\mathrm{}0`$ and the MRF (17) transforms into $`𝑴(R)/R^3+𝐖(R)𝝆/R^5`$. This last formula shows that the reaction field of finite systems is inhomogeneous even for point dipoles. Only for macroscopic ($`R\mathrm{}`$) systems, we reproduce the well-known homogeneous reaction field $`𝑴(R)/R^3`$ introduced by Barker and Watts . For finite IS systems, additional higher multipole terms appear. This brings, for example, into existence of the new quadrupole-dipole and quadrupole-quadrupole interactions in the intermolecular potential (16). We note that the idea of using the higher multipole moments in the RF has been proposed for the first time by Friedman .

However, the modified intermolecular potential (16) still needs to be complemented by a self-consistent fluctuation formula as this has already been done in the preceding section by the fluctuation formula (7) for the potential of interaction in the site-site cut-off scheme (11). Unfortunately, it is not a simple matter to construct fluctuation formulas in the molecular cut-off approach. This problem will be considered in further studying.

The difference in the RF geometry between IS and PD models lies in the distinction for their microscopic operators of polarization density. For IS models

$$\widehat{𝑷}_\mathrm{L}(𝒌)=\frac{\mathrm{i}𝒌}{k^2}\underset{i=1}{\overset{N}{}}\text{e}^{\mathrm{i}𝒌\mathbf{}𝒓_i}\underset{a}{}q_\stackrel{}{A}\text{e}^{\mathrm{i}𝒌\mathbf{}𝜹_i^a}=\widehat{𝑴}_\mathrm{L}(𝒌)\frac{\mathrm{i}𝒌}{2}\widehat{𝒌}\widehat{𝒌}\mathbf{:}\underset{i=1}{\overset{N}{}}𝝎_i\text{e}^{\mathrm{i}𝒌\mathbf{}𝒓_i}+\mathrm{},$$
(18)   
where $`\widehat{𝑴}_\mathrm{L}(𝒌)=\widehat{𝒌}_{i=1}^N\widehat{𝒌}\mathbf{}𝝁_i\text{e}^{\mathrm{i}𝒌\mathbf{}𝒓_i}`$ is the microscopic operator of polarization density for point dipoles and an expansion over small parameter $`𝒌\mathbf{}𝜹_i^a`$ has been made . However, putting $`\widehat{𝑷}_\mathrm{L}(𝒌)\widehat{𝑴}_\mathrm{L}(𝒌)`$ in the microscopic electric field $`\widehat{𝑬}_\stackrel{}{R}^{\mathrm{RF}}(𝒌)`$ (6) at the very beginning and taking attempts to perform the inverse Fourier transform, we obtain that the corresponding integral is divergent in $`𝒌`$-space when $`k\mathrm{}`$. This divergence is involved by the specific nature of point dipoles for which the parameter $`𝒌\mathbf{}𝜹_i^a`$ becomes indeterminate in the limit $`k\mathrm{}`$ because of $`𝜹_i^a+0`$ and the expansion (18) fails. Therefore, we must manipulate with the full operator $`\widehat{𝑷}_\mathrm{L}(𝒌)`$ to obtain the interdipolar potential (12) consequently and let $`𝜹_i^a+0`$ at the end of the calculation only.

Since $`𝝁d`$ and $`𝐪d^2`$, the quadrupole contribution with respect to the dipole term is varied in (16) from of order $`(d/R)^2`$ at $`r_{ij}=0`$ to $`d/R`$ at $`r_{ij}=R_d`$. Therefore, as far as the usual intermolecular potential (8) is applied in simulations, the dielectric constant can not be reproduced with the precision better than $`d/R`$. It is evident that using the modified intermolecular potential (16) will lead to the uncertainties of order $`(d/R)^2`$. They decrease at increasing the size of the sample as $`R^2`$, i.e., with the same rate as those connected with the truncation of the potential. Effects of the octupole and higher order multipole contributions into the MRF are of order $`(d/R)^3`$ and can be ignored.

## 4 Applying the ISRF to a MCY water model

In the previous investigations \[11–13\], the standard PDRF geometry (8) has been applied to actual simulations of the MCY and TIP4P models. As a result, the static, frequency-dependent  and wavevector-dependent  dielectric constant has been determined. For these models $`d=1.837\mathrm{\AA }`$ and the cut-off radius $`R=9.856\mathrm{\AA }`$ has been used in the simulations. From the afore said in the preceding section, it is expected that the precision of these calculations can not exceed $`d/R20\%`$. We shall show now by actual calculations that this prediction indeed takes place.

As an example we apply the ISRF geometry (11) to the MCY potential . The calculations have been performed with the help of Monte Carlo (MC) simulations, details of which are similar to those reported earlier , at the density of $`\rho `$= 1.0 g/cm<sup>3</sup> and at the temperature of $`T=292`$ K, i.e., in the same thermodynamic point and yet with the same number $`N=256`$ of molecules and cut-off radius $`R=9.856\mathrm{\AA }`$ as considered in .

Our result of the calculation (7) for the longitudinal components of the wavevector-dependent infinite-system Kirkwood factor $`g_\stackrel{}{\mathrm{L}}(k)`$ and dielectric constant $`\epsilon _\stackrel{}{\mathrm{L}}(k)`$ obtained within the ISRF geometry is presented in Figs. 1 and 2, respectively, as the full circles connected by the solid curves. For the purpose of comparison, analogous calculations performed previously  within the PDRF are also included in these figures (the open circles connected by the dashed curves). It is obvious that oscillations observing in the shape of $`g_\stackrel{}{\mathrm{L}}(k)`$ and $`\epsilon _\stackrel{}{\mathrm{L}}(k)`$ obtained within the PDRF method are nonphysical and caused by the finite molecular size which is assumed to be zero in this approach. At the same time, the ISRF geometry gives the true, more smooth dependencies for the Kirkwood factor and dielectric constant because the influence of the finite molecular size is included here explicitly. As we can see from the figures, deviations of values for the wavevector-dependent dielectric quantities obtained using the PDRF from those evaluated within the ISRF geometry are significant. These deviations achieve maximal values about $`25\%`$ near $`k=3\mathrm{\AA }^1`$, where the Kirkwood factor has the first maximum. For great wavevector values $`(k>6\mathrm{\AA }^1)`$ the both geometries lead to identical results because the influence of boundary conditions is negligible in this range of $`k`$.

We remark that the wavevector-dependent quantities were calculated directly for the discrete set $`k=nk_{\mathrm{min}}`$ of grid points accessible in the simulations, where $`k_{\mathrm{min}}=0.319\mathrm{\AA }^1`$ and $`n`$ is an integer number. These quantities are marked in the figures by the symbols. To obtain intermediate values between the grid points we have used the cubic spline interpolation for the most smooth dependency, namely, for $`g_\stackrel{}{\mathrm{L}}(k)`$. Then values of $`\epsilon _\stackrel{}{\mathrm{L}}(k)`$ can be evaluated anywhere in the considered domain of $`k`$-space on the basis of the interpolation values of $`g_\stackrel{}{\mathrm{L}}(k)`$ via Eq. 7. In particular, the first singularity of $`\epsilon _\stackrel{}{\mathrm{L}}(k)`$ (see Fig. 2a) has been investigated in such a way.

## 5 Conclusion

Two alternative methods (ISRF and MRF) to overcome the difficulties associated with finiteness of the molecule with respect to the system size have been proposed for IS models of polar systems. It has been shown rigorously that the fluctuation formula, which is commonly used for the calculation of the dielectric constant in computer experiment, corresponds to the ISRF geometry with the site-site cut-off for Coulomb interaction potentials. The molecular cut-off scheme leads to the MRF geometry with an additional quadrupole term to the well-known PDRF.

It has been corroborated by actual calculations that the ISRF geometry exhibits to be much more efficient with respect to the usual PDRF method for the investigation of the dielectric properties of IS models. The modified MRF approach seem to be comparable in efficiency with the ISRF geometry. An application of the MRF to practical simulations we hope to perform in further studying.

Figure captions

Fig. 1. Longitudinal component of the wavevector-dependent Kirkwood factor for the MCY water. The results in the ISRF and PDRF geometries are plotted by the solid and dashed curves, respectively.

Fig. 2. Longitudinal component of the wavevector-dependent dielectric constant for the MCY water. Notations as for fig. 1. The vertical lines indicate positions of a singularity.