Sebastian Höfener

Wim Klopper
email: klopper@kit.edu

Klopper

Analytical nuclear gradients of the explicitly correlated Møller-Plesset second-order energy

Keywords: electron correlation, second-order Møller-Plesset perturbation theory, explicitly correlated theory, analytical nuclear gradients electron correlation, second-order Møller-Plesset perturbation theory, explicitly correlated theory, analytical nuclear gradients

The analytical computation of nuclear gradients has been derived and implemented for the explicitly correlated second-order Møller-Plesset method (MP2-F12). The implementation has been accomplished in the Turbomole program package for ansatz MP2-F12/2 * A. A Slatertype geminal expanded in six Gaussian geminals (STG-6G), a complementary auxiliary basis set (CABS), and robust density fitting approximations are used. In addition, a second-order perturbation theory correction for single excitations into the complementary auxiliary basis set (CABS singles) is included to reduce the Hartree-Fock error. Smooth convergence towards the basis set limit is observed for a selection of molecules. For computations on dimers of weakly interacting molecules in small basis sets, explicitly correlated second-order Møller-Plesset theory outperforms conventional second-order Møller-Plesset theory because basis set superposition errors are largely avoided at the MP2-F12/2 * A level.

Introduction

The determination of molecular equilibrium geometries and other stationary points (e.g., transition-state geometries) is one of the most important tasks of modern ab initio quantum chemistry, for molecules in electronically excited as well as ground states. For many quantum-chemical methods, analytical nuclear gradients are available, which is a prerequisite for computationally efficient geometry optimizations. Concerning wave-function based methods that account for electron-correlation effects, analytical first derivatives have been derived and implemented for general coupled-cluster and configuration-interaction models, that is, basically for all relevant methods [1].

Altough these methods are general in the sense that both single-and multireference-type methods are included, which are able to handle wave functions truncated at arbitrary excitation level or which include excitations selected via a complete active space, the computation of analytical nuclear gradients thus far has been restricted to coupled-cluster and configuration-interaction wave functions that are expanded in terms of orbital products. The drawback of orbital-product expansions, however, is the very slow convergence of the computed results (energies, geometries, spectra, etc. [2][3][4]) to the limit of a complete basis set when the orbital basis set is enlarged. In the present article, we shall be concerned with explicitly correlated methods. These methods go beyond the orbital-product expansion by the inclusion of additional two-electron basis functions (geminals). By virtue of these geminals, it becomes possible to satisfy the electron-electron cusp condition and to accelerate the basis set convergence towards the limit of a complete basis set.

The simplest wave-function based method accounting for electron-correlation effects is the second-order Møller-Plesset (MP2) method. For the standard MP2 theory in terms of an orbital-product expansion, analytical nuclear gradients were derived and implemented by Pople and co-workers already more than 30 years ago [5]. In the present work, we derive and implement the computation of analytical nuclear gradients of the explicitly correlated variant of MP2 theory, which today is either denoted as MP2-R12 or MP2-F12 theory. Explicitly correlated secondorder Møller-Plesset theory has recently been reviewed in Ref. [6]. Here and in the following, we shall refer to explicitly correlated theory as MP2-R12 method if the geminal consists of an orbital product multiplied by r 12 , the distance between the two electrons in the two orbitals of the orbital product, and as MP2-F12 method if the geminal is constructed from a function f (r 12 ) multiplied by an orbital product. In particular, we shall in the present article be concerned with the MP2-F12 method, where f (r 12 ) is a Slater-type geminal (STG) of the form f (r 12 ) = γ -1 {1 -exp(-γr 12 )} , (1) expanded in a basis of Gaussian geminals [7][8][9][10].

The computation of certain analytical first derivatives of the MP2-R12 and MP2-F12 energies have been reported before [11][12][START_REF] Kordel | Analytische Gradienten für die MP2-R12-Methode[END_REF][START_REF] Höfener | [END_REF], but a number of aspects have not been accounted for thus far.

In our previous work [START_REF] Höfener | [END_REF], only the analytical first derivatives needed for the computation of first-order molecular properties had been implemented in the Turbomole program package [15]. In the present article, the implementation of Ref. [START_REF] Höfener | [END_REF] is extended to the computation of analytical nuclear gradients.

In Ref. [11], Kordel et al. report the implementation of analytical nuclear gradients in the Dalton program package [START_REF]a molecular electronic structure program[END_REF] at the level of MP2-R12 theory, that is,

with linear r 12 geminals. This implementation had been accomplished for an ansatz that is usually denoted as "standard approximation A" (cf. Ref. [6]). Furthermore, in Ref. [11], the auxiliary basis set (ABS) used for the resolution-of-the-identity (RI) approximation of explicitly correlated theory was identical with the basis set of atomic orbitals (AO) used to expand the Hartree-Fock molecular orbitals (MO). The latter basis is referred to as Hartree-Fock basis (HF basis) in the following.

In modern applications of explicitly correlated theory, not this HF basis set but an especially optimized auxiliary basis set is used for the RI approximation. In particular, today, researchers in the field almost exclusively use a complementary auxiliary basis set (CABS) for this approximation [START_REF] Valeev | [END_REF]. Furthermore, only the MP2-F12 method using STG's is currently of interest.

In the present work, analytical nuclear gradients are derived and implemented for the MP2-F12/2 * A-[T+V] method. Elsewhere, we could show that applying standard approximation A and the extended Brillouin condition (EBC, indicated by the asterisk: 2 * ; cf. Ref. [6]) is sufficiently accurate for obtaining molecular geometries close to those obtained at the basis set limit [START_REF] Höfener | [END_REF]. Therefore, we decided to restrict the implementation of analytical nuclear gradients to the MP2-F12/2 * A-[T+V] level. Extending the work towards the MP2-F12/2B-[T+V] level does not seem worth the effort [START_REF] Höfener | [END_REF].

Compared with our previous work on analytical first derivatives [11,[START_REF] Höfener | [END_REF], the new aspects of the MP2-F12/2 * A-[T+V] analytical nuclear gradients reported in the present article are the following:

• Use of an complementary auxiliary basis set (CABS) [START_REF] Valeev | [END_REF];

• Use of the [T+V] approximation [18,19]; • Use of density fitting [20];

• Inclusion of CABS singles [21,22].

The consequences of these new aspects will be described in detail in the remainder of this article. The CABS consists of basis functions that usually are centered at the positions of the atoms' nuclei in the molecule. Hence, the CABS changes when the nuclei are moved when computing the nuclear gradient. This gives rise to CABS related reorthonormalization terms, that is, terms related to the first derivative of the CABS overlap matrix (and of overlap matrix elements between the CABS and the HF basis set). These new terms are discussed in Section 2.3. The [T+V] approximation involves matrix elements of the kinetic energy (T) and electron-nuclear attraction (V) operators, and the first derivatives of these matrices thus enter the expression of the analytical MP2-F12/2 * A-[T+V] nuclear gradient, see Section 2.4. In the MP2-F12 implementation in Turbomole, all of the two-electron integrals are computed using density fitting techniques. Hence, the analytical computation of the first derivatives of these integrals requires taking care of changes in the density fitting approximation when nuclei are moved (Section 2.5). Finally, when the MP2-F12 energy is supplemented by a correction for single excitations into the CABS, then also the analytical nuclear gradient of the corresponding energy contribution must be considered. This is done in Section 2.6.

Apart from the above mentioned sections, the present article is organized as follows: In Section 2, we present the MP2-F12 theory and derive its nuclear gradient. Numerical issues are discussed in Section 3. Then, the performance of the MP2-F12 method is assessed and compared with conventional MP2 results (Section 4). Section 5 concludes the article. In explicitly correlated MP2-F12 theory, the first-order wave function is expanded not only in a basis of double excitations containing the product of two virtual orbitals (as in standard MP2 theory) but also in a basis of two-electron basis functions

F

Q12 f 12 |xy , (2)
where f 12 = f (r 12 ) is the correlation factor which depends on the interelectronic distance r 12 of electrons 1 and 2. In Eq. ( 2), |xy = ϕ x (1)ϕ y (2) denotes a pair of "F12 geminal" orbitals. The projection operator Q12 ensures that the new pairs are strongly orthogonal on the occupied Hartree-Fock orbitals. Throughout this article, we will use

Q12 = (1 -Ô1 )(1 -Ô2 ) -V1 V2 , (3)
known as (modified) ansatz 2 in F12 theory [START_REF] Valeev | [END_REF]23,24]. The projection operators onto the one-particle spaces are given by

Ô = m |m m| , (4)
V = a |a a| , (5) P
= p |p p| , (6)
V = a |a a | , (7)
P = p |p p | . ( 8
)
The geminal orbitals are denoted as x, y. Often, these are active (non-frozen) occupied orbitals, denoted as i, j, k, l. All occupied orbitals (frozen as well as active) are denoted as m, n, o. Since a large number of different orbital sets are used in F12 theory, the index conventions are summarized in Table 1.

In F12 methods, an auxiliary basis set (ABS) is used for the standard approximation, that is, an approximate resolution of the identity (RI) is inserted to avoid three-and four-electron integrals [23]. Typically, in the MO basis, this RI basis set consists of the orbital basis supplemented with a complementary auxiliary basis set (CABS). Using this CABS, the projection operator Q12 is approximated as

Q12 ≈ QCABS 12 = 1 -Ô1 V 2 -V 1 Ô2 -P1 P2 . ( 9
)
When computing molecular gradients, reorthonormalization terms occur since the unperturbed molecular orbitals, which are orthonormal at the reference geometry (x = 0), are no longer orthonormal at the new geometry (x = 0). (OMO) [25][26][27][START_REF] Helgaker | Calculation of geometrical derivatives in molecular electronic strucutre theory[END_REF]:

|p = r |r [S -1 2 ] rp , ( 10
)
where S is the overlap matrix at the new geometry. The transformation to the OMO basis ensures orthonormality of the orbitals, since

q|p = rs [S -1 2 ] qs S sr [S -1 2 ] rp = δ qp . ( 11
)
The final molecular orbitals at the new geometry are obtained by orbital rotations,

|p = e -κ |p , (12)
which are chosen such that κ = 0 at the reference geometry (x = 0). For the various projection operators, this means that we must write

Ô = m qr e -κ |q [S -1 2 ] qm [S -1 2 ] mr r|e κ , ( 13
) V = a qr e -κ |q [S -1 2 ] qa [S -1 2 ] ar r|e κ , (14)
P = qr |q [S -1 ] qr r| , (15)
V = P -P ,

= q r |q [S -1 ] q r r | . (16) P
Note that the projection operators P and P are invariant to orbital rotations. Furthermore, V is expressed as P -P because the complementary auxiliary basis set {a } must remain orthogonal on the orbital basis set {p} for all geometries x.

The MP2-F12/2 * A-[T+V] Lagrangian

For second-order Møller-Plesset perturbation theory, the inclusion of two-electron basis functions (geminals) in the F12 fashion is formulated using the Hylleraas functional, which consists of a conventional MP2 part, a pure F12 part, and a coupling term,

H MP2-F12 = H MP2 + H F12 + H coupling , (18)
where H MP2 contains the Hartree-Fock contribution. In the following, we will only discuss the F12 contributions to the total MP2-F12 gradient in detail. The conventional MP2 gradient is considered to be available.

In the framework of molecular gradients, a Lagrange functional L is constructed that is variational with respect to both orbital rotations and amplitudes. In ansatz 2 * , by virtue of the EBC, the coupling term H coupling in Eq. (18)

where T12 denotes the sum of the operators of kinetic energy of electrons 1 and 2.

The operator Pxy,vw symmetrizes an arbitrary four-index object E(x, y, v, w) as

Pxy,vw E(x, y, v, w) = 1 2 {E(x, y, v, w) + E(v, w, x, y)} . ( 20
)
g 12 = 1/r 12 is the operator of interelectronic repulsion. With c xy ij being the F12 amplitudes, let us introduce the quantities:

d xy ij = 2c xy ij -c yx ij , (21)
e vw xy = ij c xy ij d vw ij . ( 22
)
Inserting the projection operators into Eq. ( 19) yields different types of integrals. Except for integrals over the interelectronic distance, g 12 , integrals over the correlation factor f 12 have to be evaluated. Additionaly, integrals over the product f 12 g 12 have to be computed, as well as f 2 12 r 2 12 , arising from the kinetic energy operator acting on the correlation factor. These will be denoted f g and f 2 r 2 integrals, respectively. Furthermore, the integral over the commutator

t12 = [ T12 , f 12 ] (23)
needs to be computed. Although it can be calculated analytically, we introduce an additional RI approximation (using ABS) to avoid the calculation of these integrals [18]. This RI approximation is most easily invoked using the transformed orbitals

|p = P ĥ|p = r |r r | ĥ|p , (24)
where ĥ denotes the sum of kinetic energy and electron-nuclear attraction. The inclusion of the latter has turned out to increase numerical stability (using ĥ is denoted as [T+V] approximation) [START_REF] Bischoff | [END_REF]. Using these transformed orbitals, an integral over the commutator can be expressed as

pq|[ T12 , f 12 ]|rs ≈ Xpq,rs pq|f 12 |rs , ( 25
)
where the operator Xpq,rs carries out the appropriate transformations, Xpq,rs E(p, q, r, s) = E(p, q, r, s) + E(p, q, r, s) -E(p, q, ȓ, s) -E(p, q, r, s) . ( 26)

Here, E is an arbitrary four-index object.

As for all gradients, there are three types of contributions for MP2-F12 gradients: Differentiated one-and two-electron integrals, orbital response, and change of metric. In the framework of density fitting, differentiated two-electron integrals are not computed explicitly but rather split into their two-and three-index contributions (vide infra). Hence, analytical first derivatives with respect to x can be calculated as (cf. Ref. [30])

L [x] F12 = µ ν H ao µ ν h [x] µ ν - µ ν F eff,ao µ ν S [x] µ ν + µνP f g,f 2 r 2 f f ∆ P ,ao µν (µν|f 12 |P ) [x] - P Q g,f ,f g,f 2 r 2 f f γ P Q (P |f 12 |Q) [x] + µ ν P g,f f f ∆ P ,ao µ ν (µ ν |f 12 |P ) [x] , (27)
where the upper index [x] denotes the partial derivative according to Lagrange's method, and where F eff,ao is the effective Fock matrix taking care of the reorthonormalization of the united basis (HF + CA basis). The differentiated one-electron integrals arising from the [T+V] approximation are contracted with the density H ao . Greek letters indicate atomic orbitals while the superscript "ao" labels the contravariant quantities that transform as the MO's. In the following, we will derive expressions for the various contributions. The computation of the orbital rotation multipliers κ has been described elsewhere [START_REF] Höfener | [END_REF].

New terms due to the CABS

We are now in the position to take derivatives of the projection operators including the operator QCABS

12

. Let us first consider Ô, for which we write

∂ Ô ∂x x=0 = ∂ Ô ∂x x=0,reorth + ∂ Ô ∂x x=0,other . ( 28
)
The first term on the r.h.s. collects all terms in which the differentiated overlap matrix S [x] appears, where it is understood that only the overlap matrix elements in the AO basis are differentiated with respect to x, not the molecular orbitals. Since

∂S -1/2 ∂x x=0 = - 1 2 S [x] , (29)
∂ V ∂x x=0,reorth = - 1 2 ar |r S [x] ra a| + |a S [x] ar r| , (31)
∂ P ∂x x=0,reorth = - qr |q S [x] qr r| , (32)
∂ P ∂x x=0,reorth = - q r |q S [x] q r r | , (33)
∂ V ∂x x=0,reorth = - q r |q S [x] q r r | + qr |q S [x] qr r| = - a b |a S [x] a b b | - a r |a S [x] a r r| - qb |q S [x] qb b | . ( 34
)
Note that κ = 0 at the reference geometry (x=0). For the strong-orthogonality projection operator QCABS

12

, we obtain

∂ QCABS 12 ∂x x=0,reorth = pqr S [x] qr (|pq pr| + |qp rp|) + 1 2 mra S [x] rm |ra ma | + |a r a m| + 1 2 mra S [x] mr |ma ra | + |a m a r| + ma b S [x] a b |ma mb | + |a m b m| + ma r S [x] a r |ia mr| + |a m rm| + mqb S [x] qb |mq mb | + |qm b m| . ( 35
)
It is obvious how these derivatives contribute to the renormalization term

-µ ν F eff,ao µ ν S [x]
µ ν . The fourth term on the r.h.s. of Eq. ( 35), for example, yields the following contribution to the CABS-CABS block of F eff,ao µ ν : Of course, the [ T12 , f 12 ] integrals in Eq. ( 36) are evaluated using the [T+V] approximation. The other terms on the r.h.s. of Eq. ( 35) contribute accordingly to other parts of F eff,ao µ ν .

F eff,ao a b ← -4 m ij

New terms due to the [T+V] approximation

In the [T+V] approximation, the integrals over the commutator t12 = [ T12 , f 12 ] are computed using the matrix representation of the one-electron Hamiltonian ĥ in the united basis set {p } (HF + CA basis) and two-electron integrals over f 12 . Hence, apart from differentiated two-electron f 12 integrals, renormalization terms occur with respect to the united basis set {p }. Morover, differentiated integrals of the one-electron Hamiltonian must be accounted for.

Let us first consider the renormalization terms arising from the RI approximation in Xpq,vw pq|f 12 |vw , that is, from the projection operator P ,

∂ Xpq,vw pq|f 12 |vw ∂x x=0,reorth = - q r S [x]
q r h pq r q|f 12 |vw

+ h qq pr |f 12 |vw -h q v pq|f 12 |r w -h q w pq|f 12 |vr + other terms . ( 37
)
The "other terms" refer to renormalization terms arising from the orbitals p, q, v, w.

From this derivative of P , we obtain the following contribution to the matrix elements F eff,ao µ ν ,

F eff,
+ ij b m d g ij,b m b m|g 12 |ij [x] + ij pq d g ij,pq pq|g 12 |ij [x] + vw b m d t vw,b m b m|t 12 |vw [x] + vw pq d t vw,pq pq|t 12 |vw [x] . ( 49
)
The two-electron integrals over the Coulomb operator g 12 are evaluated using "standard" density fitting, whereas all other integrals are computed using robust density fitting. When computed by means of robust density fitting, the general structure of a differentiated integral is for all integrals x = f g, f 2 r 2 , f :

(µν|x 12 |κλ) [x] DF = P (µν|x 12 |P ) [x] g Γ P κλ + P g Γ P µν (P |x 12 |κλ) [x] + P [ x Γ P µν + x ΓP µν ](P |g 12 |κλ) [x] + P (µν|g 12 |P ) [x] [ x Γ P κλ + x ΓP κλ ] - P Q [ g Γ P µν x ΓQ κλ + g Γ P κλ x ΓQ µν + x Γ P µν g Γ Q κλ + g Γ Q µν x Γ P κλ ](P |g 12 |Q) [x] - P Q g Γ P µν (P |x 12 |Q) [x] g Γ Q κλ , (50)
where we have used the abbreviation for (e.g., for instance x = f ),

f Γ P κλ = Q (P |g 12 |Q) -1 (Q|f 12 |κλ) = Q (P |g 12 |Q) -1/2 F Q κλ , (51)
f ΓP κλ = - QR (P |g 12 |Q) -1/2 U f QR G R κλ , (52)
and

∂ ∂x M -1 = -M -1 M [x] M -1 . ( 53
)
Note that the differentiated two-electron integrals are not evaluated explictly. Rather, they are contracted on the fly with appropriate two-and three-index densities, respectively.

In the following two sections, we present the gradient expressions involving the differentiated two-and three-index integrals arising from the density fitting approximation. According to Eq. ( 27), the contributions from the differentiated two-index integrals have the structure

E [x] ← P Q (P |f 12 |Q) [x] • f γ P Q (54)
for various two-electron operators f. Collecting all terms from Eq. (49) yields

f g γ P Q = - xy ij d f g xy,ij g Γ P xi g Γ Q yj , ( 55
)
f 2 r 2 γ P Q = - xyvw d f 2 r 2 xy,vw g Γ P xv g Γ Q yw , ( 56
)
f γ P Q = + xy b m d f xy,b m g Γ P xb g Γ Q ym + pq d f xy,pq g Γ P xp g Γ Q yq + vw b m d t vw,b m ( g ΓP vb g Γ Q wm + g Γ P vb g ΓQ wm ) + pq d t vw,pq ( g ΓP vp g Γ Q wq + g Γ P vp g ΓQ wq ) . (57)
The corresponding expression for g γ P Q (i.e., the density to be contracted with the differentiated two-index Coulomb integrals) is

g γ P Q = - xyvw d f 2 r 2 xy,vw g Γ P xv f 2 r 2 ΓQ yw + g Γ P yw f 2 r 2 ΓQ xv + f 2 r 2 Γ P xv g Γ Q yw + g Γ Q xv f 2 r 2 Γ P yw - xy ij d f g xy,ij g Γ P xi f g ΓQ yj + g Γ P yj f g ΓQ xi + f g Γ P xi g Γ Q yj + g Γ Q xi f g Γ P yj - xy pq d f xy,pq g Γ P xp f ΓQ yq + g Γ P yq f ΓQ xp + f Γ P xp g Γ Q yq + g Γ Q xp f Γ P yq - xy b m d f xy,b m g Γ P xb f ΓQ ym + g Γ P ym f ΓQ xb + f Γ P xb g Γ Q ym + g Γ Q xb f Γ P ym - vw pq d t vw,pq Xvwpq g Γ P vp f ΓQ wq + g Γ P wq f ΓQ vp + f Γ P vp g Γ Q wq + g Γ Q vp f Γ P wq - vw b m d t vw,b m Xvwb m g Γ P vb f ΓQ wm + g Γ P wm f ΓQ vb + f Γ P vb g Γ Q wm + g Γ Q vb f Γ P wm - ij b m d g ij,b m g Γ P ib g Γ Q jm + pq d g ij,pq g Γ P ip g Γ Q jq . ( 58
2.5.2. Differentiated three-index integrals

In analogy to the previous section, three-index quantities contribute to the gradient as

E [x] ← µνP (µν|f 12 |P ) [x] • f ∆ P ,ao µν . ( 59
)
Note that for convenience, the three-index contributions are expressed using a back transformation,

f g ∆ P ,ao µν = xy ij d f g xy,ij C xµ C iν g Γ P yj + C yµ C jν g Γ P xi ( 60
)
f 2 r 2 ∆ P ,ao µν = xyvw d f 2 r 2 xy,vw C xµ C vν g Γ P yw + C yµ C wν g Γ P xv f ∆ P ,ao µ ν = - xy b m d f xy,b m C xµ C b ν g Γ P ym + C yµ C mν g Γ P xb - vw b m d t vw,b m Xvwb m C vµ C b
Note that for the operators f 12 and g 12 , the back transformation must be carried out in the united basis due to the [T+V] transformation. When computing the CABS singles energy, we correct for the Hartree-Fock basis set error by means of second-order perturbation theory,

g ∆ P ,ao µ ν = xy ij d f g xy,ij C xµ C iν [ f g Γ P yj + f g ΓP yj ] + C yµ C jν [ f g Γ P xi + f g ΓP xi ] + xy vw d f 2 r 2 xy,vw C xµ C vν [ f 2 r 2 Γ P yw + f 2 r 2 ΓP yw ] + C yµ C wν [ f 2 r 2 Γ P xv + f 2 r 2 ΓP xv ] + ijpq d g pq,ij C iµ C pν g Γ P ip + C jµ C qν g Γ P jq + ijb m d g b m,ij C iµ C b ν g Γ P ip + C jµ C mν g Γ P jm + xy pq d f xy,pq C xµ C pν [ f Γ P yq + f ΓP yq ] + C yµ C qν [ f Γ P xp + f ΓP xp ] + xy b m d f xy,b m C xµ C b ν [ f Γ P ym + f ΓP ym ] + C yµ C mν [ f Γ P xb + f ΓP xb ] + vw pq d t vw,pq Xvwpq C vµ C pν [ f Γ P wq + f ΓP wq ] + C wµ C qν [ f Γ P vp + f ΓP vp ] + vw b m d t vw,b m Xvwb m C vµ C b ν [ f Γ P wm + f ΓP wm ] + C wµ C mν [ f Γ P vb + f ΓP vb ] (63
H S = mn a b c a m a | Q(S mn F -F mn ) Q|b c b n + 2 m a c a m a | Q F |m , (64)
with Q = 1 -P . Note that this expression does not account for single excitations into the HF basis, since we work within the EBC (f a a ≈ 0). In the original implementation, when computing the CABS singles energy, we used an approximated RI-JK Fock matrix whenever one index is in the CABS but the exact Fock matrix for the occupied-occupied matrix elements F mn . Due to these matrix elements, which depend on exact four-index two-electron integrals, differentiated four-index integrals would have entered the expression for the analytical nuclear gradient. To avoid the calculation of differentiated four-index integrals, we here invoke the RI-JK approximation also when computing the occupied-occupied matrix elements F mn . Numerical tests showed that the change in energy due to invoking the RI-JK approximation for the occupied-occupied matrix elements is completely negligible compared with the remaining basis set errors.

Differentiated integrals

We obtain for the derivatives of the one-electron integrals

E [x] S ← mn a b c a m (δ mn h [x] a b -δ a b h [x] mn )c b n + 2 a c a m h [x] ma . ( 65
)
Furthermore, also the differentiated two-electron integrals are taken into account

E [x] S ← mn a b c a m (δ mn G [x] a b -δ a b G [x] mn )c b n + 2 a c a m G [x] ma , ( 66
)
where G p q denotes the sum of Coulomb and exchange contributions in the RI-JK approximation.

Contributions to the effective Fock matrix

For the effective Fock matrix, we obtain contributions from the outer orbitals

E [x] S ← - a m d a m c S [x] a c F c m + o S [x] om F a o + 2 r S [x] a r F rm - 1 2 a b d a b c S [x] a c F c b + c S [x] b c F c a + 4 o S [x] a o F ob + 1 2 mn d mn S [x] mr F rn + S [x] nr F rm . ( 67
)
with

d mn = a c a m c a n , d a b = m c a m c b m , d a m = d ma = c a m . ( 68
where Pro E(r, o) = E(r, o) + E(o, r).

Computational details

General remarks

We have used Dunning's orbital basis sets aug-cc-pVXZ with X = D, T, Q, 5 [31,32] together with the MP2 fitting basis sets C-aug-cc-pwCV{T,Q,5}Z for X = D, T, Q. For the aug-cc-pV5Z orbital basis, due to the lack of a C-aug-cc-pwCV6Z MP2 fitting basis, we used the C-aug-cc-pV6Z MP2 fitting basis [33,34]. Note that also for the standard MP2 calculations, we have used the C-aug-cc-pwCV(X+1)Z MP2 fitting basis, that is, one cardinal number higher than what is usually done.

For the RI-JK approximation, we have applied the auxiliary basis sets JK-augcc-pVXZ [35]. For the intrinsic RI approximation of F12 theory, we have used the CABS of Yousaf and Peterson [36].

The convergence criterion for the density-matrix norm was set to 10 -9 . The geometries were considered to be optimized when a change less than 10 -8 E h in the energy and less than 10 -4 E h /a 0 in the gradient was reached. Tightening the convergence criteria did not affect the significant digits reported. The frozen-core approximation was applied in all of the calculations.

For conventional methods using density fitting, the HF basis and the level of correlation treatment should be the factors that limit accuracy, whereas errors due to integral screening or density fitting should be negligible. With respect to the CABS, one has to distinguish between CABS as auxiliary basis for the MP2-F12 method and as basis for the computation of the CABS singles energy correction. Concerning the use as auxiliary basis, we note that the total basis for the RI approximation of MP2-F12 theory consists of the HF basis supplemented with the CABS. Hence the importance of the CABS decreases with increasing HF basis, because then the HF basis already becomes sufficiently large.

Density fitting error

As pointed out in Section 2, robust density fitting is applied to all two-electron integrals in our explicitly correlated methods. Since thus far no density fitting auxiliary basis sets have been optimized especially for the F12 methods, it is not a priori clear, which (non-optimized) auxiliary basis sets one should select. Typically, as a default, when using an aug-cc-pVXZ orbital basis in conjunction with F12 geminals, we use the C-aug-cc-pwCV(X+1)Z MP2 fitting basis for the density fitting approximation. This large MP2 fitting basis is selected to be on the safe side. It is taken with increased cardinal number, and furthermore, it is a fitting basis that would allow for correlating the core electrons. The effect of the MP2 fitting basis is investigated by performing calculations on C 2 H 4 and one of its dimer structures. Selected results are summarized in Table 2. It is observed for the C 2 H 4 molecule that different density fitting basis sets yield the same geometry when the explicitly correlated amplitudes are optimized (var-MP2-F12 approach). This is not true for its dimer. Using the var-MP2-F12 approach, the differences in the optimized dimer geometries are significant. When the C-augcc-pwCVTZ is used as MP2 fitting basis, the intermolecular distance is 4 pm longer than when the C-aug-cc-pwCVQZ or C-aug-cc-pwCV5Z MP2 fitting basis sets are used. Fortunately, we find no differences in the optimized dimer geometries when the explicitly correlated amplitudes are kept fixed (f ix-MP2-F12 approach) [37]. Hence, it seems that in particular the nondiagonal elements of the matrices B and V of F12 theory are sensitive to the MP2 fitting basis in this weakly interacting system. Further investigations of the requirements on the MP2 fitting basis set appear to be needed, but treating a weakly interacting system with the var-MP2-F12 approach is certainly a special and challenging case.

Results

Molecules

Table 3 shows the convergence of the geometries of selected molecules using the var -MP2-F12/2 * A-[T+V] method. The geometries exhibit a smooth and fast convergence to the basis set limit. Even for the small double-ζ basis set, reasonable results are obtained both in terms of energies and geometrical parameters (i.e., bond lengths and angles). See also Ref. [38]. The geometries are improved, that is, brought closer to the corresponding basis set limit geometries, by including the CABS singles energy correcions. By accounting for CABS singles in a given aug-cc-pVXZ basis, the optimized geometries are close to those obtained in the next larger aug-cc-pV(X+1)Z basis without CABS singles. Basis set convergence is practically reached at the quadruple-ζ level. Calculations with fixed F12 amplitudes yield very similar results.

In Table 3, the row denoted as "CABS" displays the results obtained with the larger C-aug-cc-pwCV(X+2)Z MP2 fitting basis. The effect of taking this larger fitting basis appears to be negligible.

Hydrogen-bonden complexes

We have first used the fix -MP2-F12/2 * A-[T+V] approach for the hydrogen-bonded complexes, that is, the approach with fixed F12 amplitudes. Results of geometry optimizations of HF..H 2 O and NH 3 ..H 2 O are summarized in Table 4. The geometries of the individual molecules and their dimers have been optimized independently, and hence, the reported energy differences are equilibrium binding energies (-D e ).

We observe that basis sets larger than double-ζ provide a smooth and rapid convergence to the basis set limit both for binding energies and intermolecular distances. At the aug-cc-pV5Z level, the inclusion of CABS singles has little effect, and we expect that we have approached the basis set limit closely at this level.

Next, we assess the performance of the var -MP2-F12/2 * A-[T+V] approach (Table 5). When using the var approach, in which all F12 amplitudes are optimized, the results may be slightly corrupted by the so-called "geminal basis set superposition error" (geminal BSSE) [39]. This error may occur because pairs of occupied orbitals located on one fragment may be correlated by geminals consisting of a Slater-type geminal muliplied by geminal orbitals that are located on the other fragment. To investigate this geminal BSSE, we decided to use a large MP2 fitting basis (aug-cc-pwCV(X+2)Z in place of aug-cc-pwCV(X+1)Z) to avoid density fitting errors. Indeed, when comparing the results of Tables 4 and5, we observe that in the aug-cc-pVDZ and aug-cc-pVTZ basis sets, the binding energies obtained from the var -MP2-F12 approach are more negative than those obtained from the fix -MP2-F12 approach. In fact, the latter approach yields results closer to the basis set limit. Considering that the fix -MP2-F12 approach is computationally less involved than the var -MP2-F12 approach, and considering that the var -MP2-F12 approach requires a larger MP2 fitting basis, it is clear that the fix -MP2-F12 approach is the method of choice for studies of weakly interacting complexes.

Increasing the CABS in the fix -MP2-F12 approach (Table 4) leads to increased intermolecular distances of 0.23 pm and 0.12 pm, respectively, for HF..H 2 O and NH 3 ..H 2 O in the double-ζ orbital basis (cf. row denoted "CABS"). Thus, the double-ζ CABS appears to be slightly too small. However, in the larger orbital basis sets, the effect of the CABS is much smaller, of the order of only 0.05 pm or less.

Comparison with standard MP2 results

In this section, we shall compare our MP2-F12 results with the standard MP2 approach. We shall consider the same two dimers as in the previous section as well as the water molecule.

Table 6 shows the convergence of the bond length and bond angle of the water molecule, as obtained from geometry optimizations. The correlation energy as well as the geometrical parameters converge smoothly to their basis set limits. Concerning the geometry, one may consider the geometry being converged with respect to the basis set at the aug-cc-pV5Z level. However, even at this level, the error in the second-order correlation energy is of the order of 7.5 mE h . When performing a single-point MP2-F12/aug-cc-pV5Z calculation at the standard MP2/augcc-pV5Z geometry, however, the basis set limit second-order correlation energy is obtained. In this respect, it seems legitimate to optimize equilibrium geometries at the standard MP2 level in conjunction with MP2-F12 computations of the final total energy.

A slightly different picture is observed for the dimers (Table 7). A fortuitous error compensation occurs at the standard MP2 level when the dimer geometries are optimized without correction for the basis set superpostion error (BSSE) by means of a function counterpoise (CP) correction. This error compensation occurs since the basis set truncation error of the orbital-product approach partially cancels the BSSE. It is most pronounced for the calculations in the aug-cc-pVDZ basis. Without CP correction, the intermolecular bond lengths are ca. 0.7 pm longer than the basis set limit value for both complexes. A geometry optimization including CP corrections, however, yields much larger bond lengths, ca. 5-6 pm longer. The interaction energies are too large in magnitude without CP correction while they are too small when the CP correction is applied.

Comparison with literature data

For the dimers studied here, only few calculations at the basis set limit are known, see Table 8. Concerning the HF..H 2 O complex, a good agreement is obtained when comparing our results with literature data. For example, in Ref. [40], an intermolecular distance of 171.44 pm was computed using a quintuple-ζ orbital basis, which may be compared with our standard MP2/aug-cc-pV5Z value of 171.43 pm (Ta- ). The fact that we use density fitting may explain the difference of 0.01 pm. In Ref. [40], the authors emphasize that a counterpoise correction has to be carried out for every HF basis. Moreover, it is unclear whether Q5 extrapolations or sixtuple-ζ results exhibit higher accuracy (both corrected for BSSE) [41]. Compared with the basis set limit computed in this work, the MP2/aug-cc-pV5Z intermolecular distance is about 0.6 pm too long when corrected for BSSE (0.1 pm too short without correcting). The large counterpoise corrections show that the standard MP2 results converge only slowly to the basis set limit. In contrast, the MP2-F12 method provides results much closer to the basis set limit (Table 4). We expect that our MP2-F12 results for the water..ammonia complex are similarly accurate.

Conclusions

In the present article, we have reported the implementation of analytical nuclear gradients of the explicitly correlated Møller-Plesset second-order energy. In particular, the implementation has been accomplished for the MP2-F12/2 * A-[T+V] model as implemented in the Turbomole program package. Equilibrium geometries close to the limit of a complete basis are obtained at the MP2-F12/2 * A-[T+V] level already in relatively small basis sets. Concerning weakly (hydrogen-bonded) complexes, it is strongly recommended to use the MP2-F12/2 * A-[T+V] method in the fashion, in which the geminal amplitudes are kept fixed (fix-MP2-F12 approach). The fix-MP2-F12 approach appears less amenable to shortcomings in the MP2 fitting basis for the density fitting approximation, and the fix-MP2-F12 approach is free of geminal BSSE.

Acknowledgements

The development and implementation of the MP2-F12 method, including the calculation of analytical nuclear gradients, have been supported by the a Basis set limit for binding energies using MP2-F12 at given geometry (basis: aug-cc-pV5Z).

b CABS was increased by one cardinal number.

c Used as basis set limit in the following. a Single-point MP2-F12 energy calculation (orbital basis: aug-cc-pV5Z).

b CABS was increased by one cardinal number.

c See Table 4. a Single-point MP2-F12/aug-cc-pV5Z energy.

b See Table 4.

1 .

 1 Differentiated two-index integrals

  mb xy|f 12 |mb mb |f 12 |r w + xy|f 12 |b m b m|f 12 |r w .(40)

	March 22, 2010 12:16 Molecular Physics Hoefener-Klopper March 22, 2010 12:16 Molecular Physics Hoefener-Klopper
						Molecular Physics	11
	H b r , and H r v , respectively: Using these intermediates, the gradient contribution from differentiated two-
	electron integrals can be written as
	E F12 ← [x]		d f g xy,ij xy|f 12 g 12 |ij [x] +	d f 2 r 2 xy,vw xy|f 2 12 r 2 12 |vw [x]
		ij xy			xy,vw
		+	H pr ← -2 xy b m d f xy,b m xy|f 12 |b m [x] + vw xy e vw xy Pxy,vw	xy pq	d f xy,pq xy|f 12 |pq [x]
		H mr ← -2	e vw xy Pxy,vw
						vw xy
		H b r ← -2	e vw xy Pxy,vw
						vw xy
	+	F o r H r v ← w xy F o r	e vw xy + e xy vw	pq	xy|f 12 |pq pq|f 12 |r w
			P P
					e e
						e r e r
						R R
						e e
		×{ d f g xy,ij = 2d xy ij d f 2 r 2 xy,vw = e xy vw	µ ν	ao v ← -2 vw xy i e w e vw xy Pxy,vw v i e w	(41) (42)
	. Next, we consider the contributions in which the differentiated one-electron Hamil-(38) tionan integrals occur, O n l y d f xy,b m = -4 ij d xy ij b m|g 12 |ij -vw (e xy vw + e vw xy ) b m|t 12 |vw (43) d f xy,pq = -2 ij d xy ij pq|g 12 |ij -1 2 vw (e xy vw + e vw xy ) pq|t 12 |vw (44) d g ij,b m = -4 xy O n l d xy ij b m|f 12 |xy (45) y
	∂ Xpq,vw pq|f 12 |vw ∂x d g ij,pq = -2	x=0,h [x] xy d xy ij pq|f 12 |xy = r	h [x] pr r q|f 12 |vw + h	[x] qr pr |f 12 |vw	(46)
		d t vw,b m = -	xy	-h [x] r v pq|f 12 |r w -h [x] r w pq|f 12 |vr xy ) xy|f 12 |b m vw + e vw (e xy	(39) (47)
		d t vw,pq = -	1 2 xy	(e xy vw + e vw xy ) xy|f 12 |pq	(48)

pq xy|f 12 |pq h pq r q|f 12 |vw -h q v pq|f 12 |r w + mb xy|f 12 |mb h mq r b |f 12 |vw + h b q mr |f 12 |vw -h q v mb |f 12 |r w -h q w mb |f 12 |vr This leads to the following contributions to the density matrix elements H pr , H mr , q xy|f 12 |pq r q|f 12 |vw , b xy|f 12 |mb r b |f 12 |vw , m xy|f 12 |mb mr |f 12 |vw , 2.5. New terms due to density fitting Let us introduce the following intermediate densities:

Table 1 .

 1 Summary of index conventions. , γ , . . . complementary virtual orbitals in formally complete basis a , b , c , . . . CABS representation of complementary virtual orbitals p , q , r , . . .

	Deutsche Forschungsgemeinschaft's Priority Programme SPP 1145 through grant
	No. KL 721/2. S. H. gratefully acknowledges support by the Deutsche Telekom
	Stiftung.

Table 2 .

 2 Investigations concerning the density fitting error using aug-cc-pVDZ as orbital basis set. CABS singles were not included. Different auxiliary basis sets of the type C-aug-cc-pwCVYZ have been used. ||∇|| denotes the norm of the gradient for the inital geometry in E h /a0. Convergence criteria are described in the text. REt-Et denotes the intermolecular distance in the dimer. Energies E are given in mE h , distances R in pm, and angles in degrees.

				DF basis:
			Y=T	Y=Q	Y=5
			( = X+1) ( = X+2) ( = X+3)
	C 2 H 4		var -MP2-F12/2 * A-[T+V]
	1st Iteration	||∇||	27.11	27.10	27.10
		E F12	-97.03	-97.02	-97.02
	Optimized	R C-C	133.77	133.77	133.77
		R C-H	108.69	108.69	108.69
		H-C-H	117.53	117.53	117.53
	(C 2 H 4 ) 2		var -MP2-F12/2 * A-[T+V]
	1st Iteration	||∇||	20.97	20.51	20.50
		E F12	-192.45	-192.20	-192.20
	Optimized	R C-C	133.80	133.82	133.82
		R C-H	108.75	108.75	108.76
		H-C-H	117.33	117.44	117.44
		R Et-Et	370.00	366.42	366.36
	(C 2 H 4 ) 2		fix -MP2-F12/2 * A-[T+V]
	1. Iteration	||∇||	20.21	20.21	20.21
		E F12	-181.08	-181.06	-181.06
	Optimized	R C-C	133.84	133.84	133.84
		R C-H	108.72	108.71	108.71
		H-C-H	117.43	117.43	117.43
		R Et-Et	369.69	369.69	369.69

Table 3 .

 3 Geometry optimization of different molecules using var -MP2-F12 and the [T+V] commutator approximation. Bond lengths R are given in pm, angles in degrees.

		CABS		H 2 O	NH 3	HF
	Basis	singles	R O-H	H-O-H	R N-H	H-N-H	R H-F
	aVDZ yes	96.02	104.32	101.08	106.85	92.01
		no	96.24	104.08	101.50	106.42	92.21
		CABS a	96.24	104.07	101.49	106.44	92.21
	aVTZ	yes	95.86	104.39	100.98	106.94	91.87
		no	95.97	104.39	101.03	106.94	92.05
		CABS a	95.97	104.39	101.03	106.94	92.05
	aVQZ	no	95.85	104.39	100.97	106.98	91.87
		yes	95.84	104.39	100.96	106.97	91.84
	aV5Z	no	95.84	104.39	100.96	106.99	91.85
		yes	95.84	104.39	100.95	106.99	91.85
	a CABS was increased by one cardinal number.

Table 4 .

 4 Geometry optimizations at the fix -MP2-F12 level without counterpoise correction. aVXZ is a shorthand notation for the aug-cc-pVXZ orbital basis set. O NH 3 ..H 2 O HF..H 2 O NH 3 ..H 2 O

	Orbital CABS	Bond length / pm	Binding energy / mE h
	basis HF..H 2 aVDZ singles yes 171.38	196.56	-13.35	-10.10
		limit a			-13.88	-10.40
		no	171.72	196.00	-13.58	-10.54
		CABS b	171.95	196.12	-13.44	-10.49
	aVTZ	yes	170.82	195.96	-13.86	-10.40
		limit a			-13.89	-10.37
		no	170.25	195.45	-14.05	-10.51
		CABS b	170.29	195.50	-14.02	-10.50
	aVQZ	yes	170.81	196.07	-13.91	-10.41
		no	170.72	196.01	-13.97	-10.44
	aV5Z c	yes	170.84	196.09	-13.90	-10.40
		no	170.82	196.08	-13.90	-10.40

Table 5 .

 5 Geometry optimizations at the var -MP2-F12 level without counterpoise correction. aVXZ is a shorthand notation for the aug-cc-pVXZ orbital basis set. Auxiliary basis sets of the type C-aug-cc-pwCV(X+2)Z were used for density fitting. O NH 3 ..H 2 O HF..H 2 O NH 3 ..H 2 O

	Orbital CABS	Bond length / pm	Binding energy / mE h
	basis HF..H 2 aVDZ singles yes 171.09	196.24	-13.63	-10.40
		limit a			-13.89	-10.41
		no	171.46	195.69	-13.85	-10.83
		CABS b	171.73	195.83	-13.75	-10.79
	aVTZ	yes	170.69	195.79	-13.98	-10.51
		limit a			-13.90	-10.41
		no	170.13	195.28	-14.17	-10.62
		CABS b	170.18	195.33	-14.13	-10.60
	Limit c		170.84	196.09	-13.90	-10.90

Table 6 .

 6 MP2 optimized geometry of H2O. At the MP2-F12/aug-cc-pV5Z level.

		R O-H / pm	H-O-H / •	E (2) / mE h
	aVDZ	96.59	103.87	-220.143
	aVTZ	96.14	104.11	-268.704
	aVQZ	95.89	104.27	-286.052
	aV5Z	95.84	104.34	-293.002
				-300.587 a
	aV6Z	95.84	104.37	-296.048
	Limit b	95.84	104.39	-300.586

a Single-point MP2-F12/aug-cc-pV5Z energy. b

Table 7 .

 7 MP2 geometry optimizations. O NH 3 ..H 2 O HF..H 2 O NH 3 ..H 2 O

	Basis	CP	Bond length / pm	Binding energy / mE h
	correction HF..H 2 aVDZ yes 176.07	201.98	-12.59	-9.31
		no	171.52	196.78	-14.39	-11.07
		limit a			-13.89	-10.44
	aVTZ	no	170.39	195.74	-14.28	-10.69
		limit a			-14.02	-10.53
	aVQZ	no	170.37	195.83	-14.17	-10.55
		limit a			-13.89	-10.40
	aV5Z	yes	171.43	196.74	-13.69	-10.26
		no	170.76	196.01	-14.02	-10.45
		limit a			-13.89	-10.40
	aV6Z	no	170.79	196.03	-13.95	-10.42
	Limit b		170.84	196.09	-13.90	-10.40

Table 8 .

 8 Comparison with literatur data. HF..H 2 O NH 3 ..H 2 O HF..H 2 O NH 3 ..H 2 O

	Method	Basis	Reference	Bond length / pm	Binding energy / mE h
	MP2	daV5Z	[40]	170.78
	MP2 (CP) a	aV5Z	[40]	171.44		-14.37
	MP2 (CP) a	Q5 b	[40]	171.02
	MP2 (CP) a	aVTZ	[42]				-9.94
	MP2	aV6Z		170.79		-13.95
	MP2-F12 c	aV5Z		170.84	196.09	-13.90	-10.40

a CP corrected results.

b Extrapolated results.

c See Table

4

. Page 22 of 22 URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics


Charles Brunet

Jean Guillaume Malherbe

Saïd Amokrane

Monte Carlo simulation of confined fluids of polarizable particles:an efficient iterative treatment of the local field in slab geometry using Ewald

Keywords: dipolar interaction, polarizability, confined fluid, simulation

published or not. The documents may come     Monte Carlo simulation of confined fluids of polarizable particles:an efficient iterative treatment of the local field in slab geometry using Ewald summation

Introduction

The study of dipolar fluids by simulation has become rather common since the mid-eighties (see for example [1][2][3][4]). However, relatively few studies dealt with the case of polarizable particles. One reason for this is a practical difficulty. Indeed, when the dipole moments are induced by the applied field, the polarization in the whole system must be updated after each particle move, in the Monte Carlo (MC) method (see for example [5]). The result is a critical slowing down of the simulation, in comparison with the situation of permanent dipoles. Special MC techniques were recently proposed to overcome this problem, such as the one proposed by Moučka et al. [6] that treats the polarization by multi-particle force-biased moves. Shortly before, Předota et al. [7] proposed an approximate method, the pair approximation for polarization interactions (PAPI) to treat the polarization by introducing a cutoff distance beyond which the local field should not not be affected significantly by a local MC move. While both methods result in a drastic increase in the simulation efficiency, they share the same ingredient to handle long range interactions, namely the reaction field method. The latter is not adapted for the slab geometry which occurs frequently in surface and interfacial systems because of its mean-field nature. As a result, one usually prefers to use in this case the Ewald summation method (see [START_REF] Holm | Efficient Methods for Long Range Interactions in Periodic Geometries Plus One Application[END_REF] and refs therein).

The specific situation of induced dipoles in slab geometry, raises thus two different problems whose combination has, to our knowledge, not been considered before. The first one is the treatment of the long-range dipolar interaction, for a given distribution of the dipole moments. The second one is the time needed to compute the induced dipoles. We thus show in this paper that it is possible to combine the slab adapted 3D-Ewald summation with the appropriate treatment of the local field to efficiently treat induced dipoles confined in slab geometry.

The practical motivation of this study is the development of experiments on con- fined colloidal particles. In this domain, permanent dipoles are usually found in ferromagnetic colloids and ferroelectric nematic ones (spherical ferroelectric particles that can be dispersed in solution have been elaborated more recently [9]).

Polarizable colloidal particles are on the other hand rather common, since polarizability of the materials forming the colloids, say the electronic one, is always present to some extent. Examples are polymethylmethacrylate (PMMA) or polystyrene Latex particles. Polarizability is also quite frequent in core-shell nanoparticles [10]). Furthermore, the consideration of an external field in confined geometry gives additional flexibility in the control of the behavior of a confined dipolar fluid [11][12][13].

Studying polarizable models in confinement is thus of much practical interest, besides the general question of the treatment of non-additive interactions [6]. We hence show in this work that the local field can be computed efficiently also in slab geometry. To this end, this paper is organized as follows: we first present the model of confined polarizable fluids that we will study. We then propose the method to treat the polarization for slab geometry by the Ewald summation. The efficiency of the proposed method is shown by simulating mixtures of polarizable and non-polarizable particles in the canonical and grand canonical ensembles. This will extend our previous results [14,15] relative to permanent dipoles to the case of induced ones. We consider polarizable particles confined between two parallel walls. In this model, the polarizability is accounted for by a point dipole µ i , taken proportional to the local electric field E loc at the position of the center of mass r i of particle i:

µ i = α p E loc (r i ). ( 1
)
α p is here taken as a scalar for simplicity. The local field involves two contributions :

E loc (r i ) = E ext (r i ) + E µ (r i ), (2)
where E ext the external field and E µ is the field due to the other induced dipoles :

E µ (r i ) = j =i T ij • µ j (3)
with

T ij = 1 r 3 ij 3r ij r ij r 2 ij -I ( 4
)
where I is the identity matrix. The induced dipoles are then given by the system of equations: Since the direct solution of eq. ( 5) is impractical [START_REF] Vesely | [END_REF], it is more convenient to compute the local field by an iterative scheme. Indeed, the local field E loc (r i ) inducing the dipole moment µ i at the position r i depends on all the dipole moments for the actual spatial configuration of the system E loc (r i ) = E loc (r i ; µ j ({r k })). Now, the difference with permanent dipoles is that one needs to evaluate the local field before computing the energy. A convenient expression of the local field can be derived starting from the potential energy for polarizable particles in an external field U written (see for example [START_REF] Böttcher | Theory of Electric polarization[END_REF][START_REF] Vesely | [END_REF]) in the form :

µ i = α p (E ext + j j =i T ij • µ j ) (5
U = U µµ + U ext + U pol , (6)
The first two terms in eq. ( 6) are the dipole-dipole and dipole-external field energies, as if the dipoles were permanent :

U µµ = - 1 2 i =j µ i • T ij • µ j (7)
U ext = - i µ i • E ext (8)
and U pol is the polarization energy [START_REF] Böttcher | Theory of Electric polarization[END_REF] or the work done against the internal forces, to form the dipoles in the local field :

U pol = 1 2 i µ i • E loc (r i ) (9)
When the relation eq. ( 1) between the induced dipoles and the local field is used, this expression becomes

U pol = 1 2α p i µ 2 i ( 10
)
Since this additional term involves only the local field any technique to handle long range interactions for permanent dipoles will also work for polarizable particles. The local field can then be obtained as

E loc = - ∂(U -U pol ) ∂µ i (11)
Once the induced dipoles are known, the total energy in eq. ( 6) can be conveniently evaluated with the help of eq. ( 7), eq. ( 8) and eq. ( 9) as:

U = - 1 2 i µ i • E ext (12)
The problem of the treatment of the periodic boundary conditions is however still present in the evaluation of the induced dipole moments µ i . This is discussed in the next section. In slab geometry, adapted 3D-Ewald summation techniques [18][19][20][21] are usually used and the local field is derived from the Ewald energy. In the Ewald sum in which the dipole-dipole term is written as U µµ = U (r) + U (k) + U (s) , each term in the potential energy gives a contribution to the local field computed from eq. (11) (see for example [START_REF] Holm | Efficient Methods for Long Range Interactions in Periodic Geometries Plus One Application[END_REF]22]). The short-range contribution U (r) may be written as :

F

U (r) = - 1 2 i =j µ i,s (T (r) ) st ij µ j,t (13)
with

(T (r) ) st ij = C(r ij , α)r s ij r t ij -δ st B(r ij , α) (14)
where s and t stand for the space coordinates {x, y, z} (Einstein's convention for repeated indexes is used); δ is the Kronecker symbol and the functions B and C are defined by :

B(r, α) = erfc (αr) /r 3 + 2α/π 1/2 exp(-α 2 r 2 )/r 2 (15)
C(r, α) = 3 erfc (αr) /r 5 + 2α/π 1/2 2α 2 + 3/r 2 exp(-

α 2 r 2 )/r 2 (16)
The corresponding contribution to the local field reads:

E (r) i,s = j j =i (T (r) ) st ij µ j,t (17)
The long-range contribution to the energy is given by :

U (k) = - 1 2 i,j µ i,s (T (k) ) st ij µ j,t (18)
with

(T (k) ) st ij = - 1 L 3 k =0 4π k 2 exp - k 2 4α 2 k s k t exp[ik • (r i -r j )] (19)
and the corresponding long-range contribution to the field in the s-direction is :

E (k) i,s = j (T (k) ) st ij µ j,t (20)
The self-energy term that eliminates unwanted self-interactions in the central box is :

U (s) = - 2α 3 3 √ π i µ i • µ i (21)
E (s) i,s = 4α 3 3 √ π µ i,s (22)
In the standard 3D-Ewald-scheme, an additional term must be added that depends on the dielectric medium surrounding the replicated boxes. For slab geometry, one actually incorporates a vacuum region between the central slab and its replicas as proposed by Yeh and Berkowitz [18]. This introduces a correction to the energy and thus to the field given by :

E (c) = 4πP z u z ( 23
)
where P is the total polarization in the slab. As a result, in Slab geometry, the component in the s-direction of the local field acting on the particle i is :

E local i,s = E ext i,s + E (r) i,s + E (k) i,s + E (s) i,s + δ sz E (c) s ( 24
)
As a verification, one can check that for permanent dipoles, the potential energy obtained with this local field as s) . For polarizable particles the energy is obtained from eq. ( 12) as recalled above. One known problem with the original method is that the vacuum region can be quite large. To reduce the gap width the so-called electrostatic layer correction (ELC) has been introduced, initially for point charges [20] and later for dipoles [21]. Furthermore the ELC method and its variants allow an a priori control of the errors, this being important to determine the optimum conditions for using the method (see for example [23]). Here, however, our purpose was to investigate the question of the local field, which is a distinguishing feature in comparison with permanent dipoles. Of course, all the improvements made in the latter case can and should be implemented, in a second stage to study induced ones.

U = U ext -1 2 i µ i • (E loc (r i ) -E ext ) is consistent with the standard Ewald Energy route U ext + U (r) + U (k) + U (

3

Iteration process for the local field with polarizable particles

The full update process

We now turn to the implementation of these considerations in Monte Carlo simulation. For polarizable particles with point dipoles, only the positions r i need to be sampled. But even a local displacement of one particle yields a change of the local field acting on all the other particles. After a move of particle i from

r old i to r new i , the new induced dipole moment is µ new i = α p E new loc (r new i
). This modifies the local field acting on the other particles and hence one has µ new j = α p E new loc (r j ). As a result, after each local change (particle move, or more generally insertion, etc.) one needs to update all the dipole moments in an iterative way to compute the local field. In an algorithm with full update, all the induced dipole moments are sequentially updated according to µ

(n) i = α p E loc ({µ (n) j < i }, {µ (n-1)
j ≥ i }) until one observes for all dipoles : where |∆µ i | is the magnitude of the variation of the dipole moment between iterations (n -1) and (n). With ǫ = 10 -4 , we checked that the energy is affected only by a term of the order of the machine single-precision error. This criterion is satisfied after approximately 5 iterations. The time required to evaluate the change in energy by this method after N particle position updates is then much larger than in simulations with permanent dipoles. This full update algorithm would require a considerable amount of computer time, for example in order to obtain very accurate results in a specific situation (see for example the benchmark result on water of Jedlovszky et al. [5]). Faster -albeit probably less accurate -variants are thus needed to explore the parameters space of a generic model in a reasonable time.

|∆µ i | |µ i | < ǫ (25)
One possibility is considered below.

Update process with a Cut-off

Following this objective, the pair approximation for polarization interaction was proposed by Předota and coworkers [7] in the context of bulk polarizable fluids. The idea of this method is to update only the dipole moments of the particles that are within a given distance R cut from a displaced particle (see figure 2 in ref. [7] for implementation details). Only interactions between pairs are taken into account during the update process. Because of the long-range nature of the dipolar interaction it is clear that introducing cutoff distances can give only an approximation of the true local field. However, when the cut-off distance is taken large enough and the local change not too important, the thermodynamics should remain correct in this approximate treatment.

In slab geometry, one cannot use the reaction field method both because of the loss of the spherical symmetry and the problem of the meaning of the dielectric constant at high confinement. Since the PAPI method uses the reaction field, it must be adapted to deal with such a geometry. As depicted in Fig. 1 the zone inside which convergence of the local field will be made is a cylinder of height H and radius R cut around the displaced particle. During the iteration process with the cut-off distance, only the dipole moments µ in inside the cylinder will be updated. Since the others, µ out , are kept unchanged, their contribution to the local field acting on the particles inside the cylinder is fixed. The iteration process to compute the new dipole moments µ in is as follows :

• A random trial move dr for particle i mv is chosen:

r new imv = r old imv + dr
• The ensemble D cut of the particles inside the cylinder is selected :

D cut = {i; (r i,x -r imv,x ) 2 + (r i,y -r imv,y ) 2 ≤ R cut }
• The contribution of the out-dipoles to the local field acting on these particles is evaluated from the non-updated values µ out (j ∈ D cut ) :

E out i,s = j ∈Dcut (T (r) ) st ij µ j,t + j ∈Dcut (T (k) ) st ij µ j,t + δ sz 4π j ∈Dcut µ j,z
In practice, E out i,s is more efficiently computed by subtracting the contribution of the actual dipoles located inside the cylinder from the stored values of the local fields computed during the previous Monte Carlo step ( not possible for the displaced particle itself).

• The dipoles moments inside the cylinder µ in i (i ∈ D cut ) are next modified :

µ in i = α p E loc i with E loc i,s = E ext s + E out i,s + j∈Dcut j =i (T (r) ) st ij µ j,t + j∈Dcut (T (k) ) st ij µ j,t + 4α 3 3 √ π µ i,s + δ sz 4π j∈Dcut µ j,z
until the criterion eq. ( 25) is verified.

Finally, the energy of the new configuration is computed with eq. ( 12) and the move is accepted or rejected according to the standard Metropolis algorithm. Note that this process is just a convenient way for solving eq. ( 5) ie computing, for a given set {r i } of particles coordinates, the distribution of dipole moments. The electrostatic energy depends thus only on {r i }, as with ordinary (pair additive) interactions. Since the sampling of the configuration space is done here by single particle translation moves (the first step above from to r new imv to r old imv ), the detailed balance condition is indeed satisfied by the standard acceptation rule. Now, the error in the change in energy after one particle move, due to the approximate treatment of the local field, should of course be minimized by taking a large cutoff, as shown below. We first tested this method on a model considered in our previous work [14] : a mixture of hard-spheres (diameter σ s ) and polarizable hard-spheres (diameter σ) confined between two parallel walls. The pore width is H = 3σ, and the size ratio σ/σ s = 2. The other parameters are the same as in part C of [14] but the permanent dipole µ * = 1.2 in the external field E * = 20 (defined for convenience as E * = E σ 3 /ǫ w from the actual field strength E and wall-small spheres interaction strength ǫ w ) is replaced by an induced one with α p = 0.06σ 3 .

For the Ewald summation we took a parameter α = 7/L were L is the lateral box size. The gap width is γL -H. As in reference [19], we took γ = 10. For the terms in reciprocal space we took n 2 = 80. No significant deviation was found [19] between the exact electrostatic energy for a lattice of permanent dipoles and the one obtained with these values of the Ewald parameters and gap width. We made a similar comparison for a slab of three layers of induced dipoles arranged on a square lattice. For α p = 0.06σ 3 and typical field strengths, we again found no significant deviation between the exact and the Ewald sum, with γ = 10. To estimate the sensitivity to this parameter, we considered the value γ = 1. We found no significant variation of the average quantities, compared to γ = 10. For example, when going from γ = 1 to γ = 10 we found a difference in the average energy per particles δU/U ≈ 3.10 -6 for E * = 20 and the densities considered below. For the slab width considered here, H = 3σ and γ = 1, the ratio γL/H is indeed already greater than three even for the smallest lateral width L = 10.25σ. Since optimizing this aspect of the method was not our objective here, we kept the somewhat rather large value γ = 10 for safety. Furthermore, when the dipoles are induced by an external field, there is a significant dipole-field contribution to the variation in electrostatic energy after one particle move. A very accurate treatment of the long-range dipole-dipole contributions seems then less pressing, especially in view of the fact that the largest error comes from the partial update itself, as discussed below.

To test the update method, we performed a converged run of N step = 40 Kcycles (one cycle involves N L + N s trial moves ) with a cut-off distance R cut = 3σ. The same run was repeated with full update. To begin with, N L = 129 polarizable particles and N s = 1032 small hard-spheres were used in both cases. While it takes 3 days on a 4-processors computer to complete the first simulation, the same run with full update takes about 15 days on the same machine.

As the time saving is important, the effect of the approximate treatment for the local field must be discussed. Fig. 2 shows the mean change in the magnitude of the dipole moment for the particles that are separated by a lateral distance greater than R from the displaced particle, for full and partial update with R cut = 3σ.

The direct effect of the cut-off distance is to ignore the variation of the dipole moment beyond R cut . The latter is more than 10 times smaller than the maximum value around R = 1. Inside the cylinder the variations differ slightly, but are on the whole quite consistent.

We now evaluate the impact of the cut-off method on the energy of the system. With the same parameters a relative difference |U Full -U Rcut=3σ |/U Full ≈ 8.10 -6 is found between the two methods, for each configuration. This difference does not increase during the simulation. This suggests that the equilibrium states sampled by the two methods are not too different. This is confirmed by Fig. 3 which shows the change in energy ∆U induced by the trial move of one particle, in the two methods, for two box sizes corresponding to the same state: the former box with N L = 129 and N s = 1032 and a bigger one having a twice larger surface, with N L = 258 and N s = 2064. Fig. 3a and Fig. 3c show that ∆U Rcut=3σ remains close to ∆U Full . The rare cases where the latter is zero but not the former show directly the magnitude of the error due to the cut-off, compared to k B T . Since the difference between the two quantities is always very small compared to k B T , the acceptance probability of the moves should remain nearly the same. Using ∆U Rcut=3σ , instead of the exact one, should thus not introduce a significant bias, as mentioned above. Fig. 3b and Fig. 3d show ∆U Rcut=3σ versus ∆U Full . In this example, the difference between the two quantities increases with the magnitude of ∆U Full . The points corresponding to ∆U Rcut=3σ seem to be distributed on a line with a slope slightly different from unity and not as a random noise about the line of unit slope. This observation could be used to improve the method by correcting ∆U Rcut=3σ (say by a linear fit that forces the equality with the exact result). But without this, we find that the cut-off induces a change in the Markov chain after approximately 250 trial moves in the smaller box. We repeated the same tests with R cut = σ. Then the Markov chain obviously deviates faster (after about 50 trials). The choice R cut = 3σ seems thus here a good compromise for reducing the computer time without significantly affecting the simulation dynamics.

We finally compare in Fig. 4a the density profiles of the dipoles for the system with N L = 129 in the two methods. We observe no significant difference due to the cut-off. To investigate the size effect in a converged simulation, we considered the bigger box with N L = 258 and N s = 2064. In the cut-off method (with R cut = 3σ) a simulation with 40 Kcycles, takes 9 days on the same 4-processors computer. This gives a scaling time factor of N 3/2 for the algorithm with cut-off. The corresponding density profile of the polarizable hard-spheres is shown in Fig. 4b. No significant size effect is found in comparison with the simulation with 129 dipoles. In the full update method now, we estimated that the same simulation would take about 6 months (on the basis of limited runs of about 100 cycles). Even so, a single check would be insufficient. We thus cannot discuss here the size effect in the method with full update, beyond the data shown as a function of the number of particle moves in Fig. 3a and Fig. 3c. From a more physical perspective now, we may note that the results are similar to those shown in [14] for permanent dipoles: without the field, the density profile of the large particles is structured at the scale of the small ones. In a normal electric field, this layering is replaced by one at the scale of the large polarizable particles.

Additional tests of the size effects and of different choices of R cut should of course be performed to assess more quantitatively the consequences of this partial update method (this last question is not specific to the slab geometry considered here). At the present stage, we found it more important to examine if the general idea remains applicable in an open system. This model is motivated by our recent study [15] in which we have shown how the composition of a fluid comprising dipolar particles in an open pore may be controlled by applying an electric field. For mixtures with non-additive interactions, a population inversion occurs in the pore beyond a certain threshold value of the applied field: starting from the situation of a pore in equilibrium with a bulk mixture having widely different concentrations of the two components, a jump in the adsorption of the minority component occurs when the field strength is varied. After the threshold, the pore becomes filled with this species, while the former majority component simultaneously desorbs.

In ref. [15] the dipoles were taken permanent. We relax here this constraint for the reasons explained in the introduction (polarizable colloids). We thus consider a mixture of hard-spheres (species (1)) and polarizable hard-spheres (species (2)) with non-additive diameters. For species (2) we took a polarizability α p = 0.04σ 3 (in order to model colloidal particles). A reduced field strength E * = E(σ 3 /kT ) 1/2 will be used. For simplicity, we took σ 2 = σ 1 = σ. For the non additivity parameter (defined from the cross diameter σ ij = 1/2(σ i + σ j )(1 + δ)) we took δ = 0.2. The mixture is confined in a pore of width H = 3σ. The bulk fluid that fixes the chemical potentials in the pore is the same mixture of non additive hard-spheres, but without polarization, because the electric field is applied only in the pore. The total reservoir density is ρ b = 0.53. The bulk mole fraction of the polarizable particles is x 2 = 0.02. The simulation in the pore is done in the Grand-canonical ensemble. The reduced density in the pore is ρ = N σ 3 /V , with N the average number of particles and V = S(H -σ).

For the translation moves, we first study the effect of the cut-off on the relative mean variation of the dipole moment after the move of one particle. Fig. 5 shows the result for R cut = 3σ and E * = 5. Although < |∆µ| |µ| > is larger than in Fig. 2, due to a different dynamics -the maximum displacement of the particles ( δr = 0.3σ) is 5 times bigger, we again find a negligible effect beyond R cut .

The consequences on the change in energy is shown in Fig. 6. The same features as in Fig. 3 to the change in energy after one particle move (with the same remark as with for Fig. 3a when DeltaU Full = 0) . The smaller energies are due to lower values of the field and polarizability. We finally measure that a change in the Markov chain occurs only after 100 trial moves, probably because of the larger accepted displacement.

Similar tests may be done for the deletion or the insertion of the dipoles. To anticipate the effect of larger changes than with translation moves only, we took a larger cut-of, R cut = 5σ. The iterations are then much longer : about one month on the same 4-processors machine, for comparable values of the parameter than in the canonical simulation. Needless to say, the simulation with full update remains equally problematic.

In Fig. 7, we present the population inversion in the pore that has been discussed in ref. [15], but now for polarizable particles. In this case, the inversion is still present but it occurs for a field that is 4 times stronger than for permanent dipole [15]. These results depend of course on the value of α p and the comparison with permanent dipole is not immediate, since the dipole moment changes with the applied field. An accurate treatment of the local field will of course be necessary to obtain accurate values (for example for the threshold field). The very question of the occurrence of the population inversion should however not depend on this, since it has been observed even without dipoles, by simply changing then the bulk density.

Conclusion

In this paper, we presented a method for studying confined polarizable fluids by Monte Carlo simulation that combines a slab-adapted Ewald method with an efficient strategy for computing the induced dipoles at each Monte Carlo step. Since we consider a slab geometry, the long-range interactions are handled by the method of Ewald sums adapted to this geometry rather than the reaction field method. To reduce the amount of computer time required by a full update of the dipole moments, we have proposed a simple update strategy that uses a cut-off distance, in the same spirit as the PAPI methods for bulk systems. We have shown that the partial update allows a large gain in computer time while preserving the consistency of the Markov chains over a significant number of Monte Carlo particle displacements. We suggested how this method could be improved to ensure a more consistent calculation of the energy when necessary. This method opens up the possibility to study confined polarizable fluids by MC simulation with a reasonable restitution time. As with bulk systems, however, the choice of the cut-off distance must be checked on a representative number of Monte Carlo steps before attempting a fully converged simulation. Other improvements specific to the slab geometry, such as the ELC method for reducing the gap width, should also be considered to reduce even further the length of the simulation.

  It contributes to the local field by :

4. 2

 2 Grand-canonical simulation of a non-additive mixture of pure and polarizable hard-spheres

  are observed. The error induced by the cut-off is still very small compared

Figure 1 .Figure 2 .

 12 Figures

Figure 3 .

 3 Figure 3.: Effect of the cut-off distance on the energy change ∆U/k B T induced by the trial moves. (a) and (b): N L = 129 and N s = 1032; (c) and (d) N L = 258 and N s = 2064

Figure 4 .Figure 5 .PageFigure 6 .Figure 7 .

 4567 Figure 4.: Density profiles of polarizable particles in an external normal electric field.

  this is of course

	April 13, 2010 18:37 Molecular Physics Brunet˙rev
	7
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics