Thermalized and mixed meanfield ADP potentials for magnesium hydrides
Supplementary Material
M.Molinos, M.Ortiz and P.Ariza
Contents
1 Euler Lagrange Equations Equilibrium Equations
1
1.1 Thermo-mechanical balance relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2 Chemical balance relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
2 Angular-dependent interatomic potential for a M-H system
2
2.1 Distance terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2.2 Energy density derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.3 Embedding energy term of the ADP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.4 Pair energy term of the ADP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.5 The Dipole contribution to the potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.6 The Quadrupole contribution to the potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
1 Euler Lagrange Equations Equilibrium Equations
1.1 Thermo-mechanical balance relations
The first condition enforces the quasi-static condition of the system:
∂
L
0
∂
̄
p
i
=
1
m
i
̄
p
i
= 0
,
therefore ̄
p
i
= 0
.
(1)
The second condition allow us to compute the mean position of the atomic sites
{
̄
q
}
:
∂
L
0
∂
̄
q
i
=
1
T
∂
⟨
V
⟩
0
∂
̄
q
i
=
1
T
⟨
∂V
∂q
i
⟩
0
= 0
.
(2)
Finally, the third condition is obtained deriving the functional
L
0
with respect the standard deviation of the atomic
positions ̄
σ
i
∂
L
0
∂
̄
σ
i
=
k
B
∂
log Ξ
0
∂
̄
σ
i
+
1
T
∂
⟨
V
⟩
0
∂
̄
σ
i
.
(3)
The first term in eq. (3) results in the following identiy:
k
B
∂
log Ξ
0
∂
̄
σ
i
=
3
k
B
̄
σ
i
.
(4)
Substituting the two preceding equations in eq. (3) we obtain the optimatily condition for the standard deviation of
the atomic positions:
3
k
B
̄
σ
i
+
1
T
1
̄
σ
3
i
⟨|
q
i
−
̄
q
i
|
2
V
⟩
0
−
3
̄
σ
i
∗
⟨
V
⟩
0
= 0
.
(5)
1.2 Chemical balance relation
Deriving at each site the free entropy
L
0
with respect the trial chemical multiplier ̄
γ
i
we get following equilibrium
equation
∂
L
0
∂
̄
γ
i
=
−
k
B
∂
log Ξ
0
∂
̄
γ
i
+
1
T
∂
⟨
V
⟩
0
∂
̄
γ
i
−
k
B
{
γ
}
T
∂
{
χ
}
∂
̄
γ
i
= 0
.
(6)
1
In order to obtain an analytical expression for
{
γ
}
at each site we proceed to further simplifications of eq. (6). Taking
derivatives with respect the grand canonical partition function, the first term can be simplified to
k
B
∂
log Ξ
0
∂
̄
γ
i
=
k
B
1
Ξ
0
∂
Ξ
0
∂
̄
γ
i
=
k
B
1
1 + e
̄
γ
i
+
γ
i
∂
e
̄
γ
i
+
γ
i
∂
̄
γ
i
=
k
B
e
̄
γ
i
+
γ
i
1 + e
̄
γ
i
+
γ
i
.
(7)
The second term can be simplified to
1
T
∂
⟨
V
⟩
0
∂
̄
γ
i
=
1
T
∂χ
i
∂
̄
γ
i
∂
⟨
V
⟩
0
∂χ
i
=
1
T
e
̄
γ
i
+
γ
i
(1 + e
̄
γ
i
+
γ
i
)
2
∂
⟨
V
⟩
0
∂χ
i
(8)
The third term can be simplified to
k
B
{
γ
}
T
∂
{
χ
}
∂
̄
γ
i
=
k
B
γ
i
∂χ
i
∂
̄
γ
i
=
k
B
γ
i
e
̄
γ
i
+
γ
i
(1 + e
̄
γ
i
+
γ
i
)
2
(9)
Introducing eqs. (7) to (9) into eq. (6)
∂
L
0
∂
̄
γ
i
=
k
B
e
̄
γ
i
+
γ
i
1 + e
̄
γ
i
+
γ
i
−
1 +
β
∂
⟨
V
⟩
0
∂χ
i
1
1 + e
̄
γ
i
+
γ
i
−
γ
i
1 + e
̄
γ
i
+
γ
i
=
k
B
e
̄
γ
i
+
γ
i
(1 + e
̄
γ
i
+
γ
i
)
2
−
(1 + e
̄
γ
i
+
γ
i
) +
β
∂
⟨
V
⟩
0
∂χ
i
−
γ
i
= 0
(10)
Therefore
−
(1 + e
̄
γ
i
+
γ
i
) +
β
∂
⟨
V
⟩
0
∂χ
i
−
γ
i
= 0
(11)
Factorizing out
γ
i
we get the following explicit equilibrium relation
γ
i
=
−
1
1
−
χ
i
+
β
∂
⟨
V
⟩
0
∂χ
i
=
−
1
1
−
χ
i
+
β
⟨
∂V
∂n
i
⟩
0
(12)
2 Angular-dependent interatomic potential for a M-H system
The Angular-Dependent interatomic Potential (ADP) is a class of many-body potentials developed by Mishin
et al.
as an extension of the Embedded-Atom Method (EAM) potentials.
⟨
V
(
{
q
}
)
⟩
0
=
X
i
⟨F
( ̄
ρ
i
)
⟩
0
+
1
2
X
i
X
j
̸
=
i
⟨
φ
i,j
⟩
0
+
1
2
X
i
⟨
μ
i
·
μ
i
⟩
0
+
1
2
X
i
⟨
λ
i
:
λ
i
⟩
0
−
1
6
X
i
⟨
(tr
λ
i
)
2
⟩
0
=
⟨
V
F
⟩
0
+
⟨
V
φ
⟩
0
+
⟨
V
μ
⟩
0
+
⟨
V
λ
⟩
0
.
(13)
For binary systems such as metal-hydrogen (M-H), and assuming that the molar fraction of the metal sites is fixed,
χ
= 1, and the molar fraction of the hydrogen sites may vary,
χ
∈
(0
,
1) the potential, and each individual term, is
fully developed in the subsequent subsections. In order to avoid unnecessary operations, each atomic site
i
interact
just with the nearest neighbour
j
.
B
i
(
r
c
) denotes the set of atomic sites within a distance
r
c
of site
i
, and
r
c
is the
cut-off distance of the ADP potential:
B
i
(
r
c
)
:
=
S
1
(
r
ij
,r
c
)
∩
(
I
M
∪
I
H
)
.
(14)
Therefore
B
i
(
r
c
) can be viewed as the set of atomistic positions which defines the neighborhood of the site
i
.
2.1 Distance terms
Let us define the vectorial space
V ∈
R
dim
. We will refer to the distance between any two atomic positions
i
and
j
as
r
ij
∈V
:
r
ij
=
q
i
−
q
j
∈V
,
(15)
with Euclidean norm:
r
ij
=
√
r
ij
·
r
ij
∈
R
.
(16)
2
The gradient of the norm
r
ij
can be computed as:
dr
ij
d
r
ij
=
r
ij
r
ij
.
(17)
Therefore
dr
ij
d
q
i
∗
:
dr
ij
d
q
i
∗
=
r
ij
r
ij
if
i
∗
=
i
dr
ij
d
q
i
∗
=
−
r
ij
r
ij
if
i
∗
≡
j
(18)
On the other hand, the hessian of norm
r
ij
has the following expression:
d
2
r
ij
d
r
2
ij
=
d
d
r
ij
r
ij
r
ij
=
1
r
ij
I
−
1
r
ij
r
ij
⊗
r
ij
r
2
ij
!
(19)
Therefore:
d
2
r
ij
d
q
2
i
=
d
d
q
i
dr
ij
d
q
i
=
d
d
q
i
r
ij
r
ij
=
d
d
r
ij
r
ij
r
ij
d
r
ij
d
q
i
=
d
2
r
ij
d
r
2
ij
(20)
The remaining terms are obtained similarly:
d
2
r
ij
d
q
i
d
q
j
=
−
d
2
r
ij
d
r
2
ij
,
d
2
r
ij
d
q
j
d
q
i
=
−
d
2
r
ij
d
r
ij
d
r
ij
,
and
d
2
r
ij
d
q
2
j
=
d
2
r
ij
d
r
ij
d
r
ij
(21)
For the sake of the clarity on the subsequent developments, we will not expand the expression
d
2
r
ij
d
r
2
ij
. Instead we will
leave the derivative.
2.2 Energy density derivatives
The phase-averaged pseudo-density function
⟨
̄
ρ
i
⟩
0
of binary system M-H is defined as:
⟨
̄
ρ
i
⟩
0
=
X
j
∈
I
M
j
̸
=
i
⟨
ρ
M
(
r
ij
)
⟩
0
+
X
j
∈
I
H
j
̸
=
i
χ
j
⟨
ρ
H
(
r
ij
)
⟩
0
(22)
Where
ρ
M
and
ρ
H
are experimentally-fitted 5th-order polynomial which depends on
r
ij
=
|
q
i
−
q
j
|
with first and
second derivatives
ρ
′
(
r
ij
) and
ρ
′′
(
r
ij
). The numerical evaluation of eq. (22) is detailed here below:
⟨
̄
ρ
i
⟩
0
≈
X
j
∈
I
M
∩B
i
(
r
c
)
j
̸
=
i
NQ
X
k
=1
ρ
M
(
r
ij,k
)
W
k
+
X
j
∈
I
H
∩B
i
(
r
c
)
j
̸
=
i
χ
j
NQ
X
k
=1
ρ
H
(
r
ij,k
)
W
k
(23)
Deriving
⟨
̄
ρ
i
⟩
0
with respect to
̄
q
i
∗
results in:
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
∗
=
X
j
∈
I
M
j
̸
=
i
⟨
ρ
′
M
(
r
ij
)
dr
ij
d
q
i
∗
⟩
0
+
X
j
∈
I
H
j
̸
=
i
χ
j
⟨
ρ
′
H
(
r
ij
)
dr
ij
d
q
i
∗
⟩
0
,
if
i
∗
=
i
(24)
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
∗
=
⟨
ρ
′
M
(
r
ij
)
dr
ij
d
q
i
∗
⟩
0
if
i
∗
̸
=
i,
and
i
∗
∈
I
M
(25)
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
∗
=
χ
j
⟨
ρ
′
H
(
r
ij
)
dr
ij
d
q
i
∗
⟩
0
if
i
∗
̸
=
i,
and
i
∗
∈
I
H
(26)
The numerical evaluation of eqs. (24) to (26) is detailed taking eq. (24) as representative:
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
∗
≈
X
j
∈
I
M
∩B
i
(
r
c
)
j
̸
=
i
NQ
X
k
=1
ρ
′
M
(
r
ij,k
)
dr
ij,k
d
q
i,k
∗
W
k
+
X
j
∈
I
H
∩B
i
(
r
c
)
j
̸
=
i
χ
j
NQ
X
k
=1
ρ
′
H
(
r
ij,k
)
dr
ij,k
d
q
i,k
∗
W
k
if
i
∗
=
i
(27)
3
The gradient of the distance norm,
dr
ij
d
q
i
∗
, can be obtained from eq. (18). Second order derivatives of
⟨
̄
ρ
i
⟩
0
with respect
to
̄
q
i
∗
:
∂
2
⟨
̄
ρ
i
⟩
0
∂
̄
q
2
i
=
X
j
∈
I
M
, j
̸
=
i
⟨
ρ
′′
M
(
r
ij
)
r
ij
⊗
r
ij
r
2
ij
⟩
0
+
⟨
ρ
′
M
(
r
ij
)
d
2
r
ij
d
r
2
ij
⟩
0
+
X
j
∈
I
H
, j
̸
=
i
χ
j
⟨
ρ
′′
H
(
r
ij
)
r
ij
⊗
r
ij
r
2
ij
⟩
0
+
χ
j
⟨
ρ
′
H
(
r
ij
)
d
2
r
ij
d
r
2
ij
⟩
0
(28)
Finally, we can compute the derivatives of
⟨
̄
ρ
i
⟩
0
with respect the molar fraction
χ
i
∗
∀
i
∗
∈
I
H
:
∂
⟨
̄
ρ
i
⟩
0
∂χ
i
∗
=
∂
∂χ
i
∗
X
j
∈
I
M
, j
̸
=
i
⟨
ρ
M
(
r
ij
)
⟩
0
+
X
j
∈
I
H
, j
̸
=
i
χ
j
⟨
ρ
H
(
r
ij
)
⟩
0
=
∂
∂χ
i
∗
X
j
∈
I
H
, j
̸
=
i
χ
j
⟨
ρ
H
(
r
ij
)
⟩
0
=
⟨
ρ
H
(
r
ij
)
⟩
0
,
if
i
∗
=
j
(29)
2.3 Embedding energy term of the ADP
The contribution of the embedding energy
F
to the ADP is given by:
⟨
V
F
⟩
0
=
X
i
∈
I
M
⟨F
M
( ̄
ρ
i
)
⟩
0
+
X
i
∈
I
H
χ
i
⟨F
H
( ̄
ρ
i
)
⟩
0
(30)
Where
F
M
and
F
H
are experimentally-fitted 5th-order polynomial which depends on the pseudo-density function ̄
ρ
i
.
For a three-dimensional material sample,
F
M
and
F
H
are 3(Q+1)-dimensional integrals, where Q denotes the number
of neighbor sites within the cut-off distance of the ADP potential, which in practice is typically of the order of 100.
Therefore, these integrals can not be evaluated analytically due to the complex form of the embedding energy function,
an thus It is noteworthy the evaluation of the
meanfield
embedding energy term,
⟨
V
F
⟩
0
, requires special attention due
to it non-linear relation with ̄
ρ
. Expanding the embedding energy
F
( ̄
ρ
i
) about the mean value of ̄
ρ
i
, we obtain:
F
( ̄
ρ
i
) =
F
(
⟨
̄
ρ
i
⟩
0
) +
F
′
(
⟨
̄
ρ
i
⟩
0
) ( ̄
ρ
i
−⟨
̄
ρ
i
⟩
0
) +
1
2
F
′′
(
⟨
̄
ρ
i
⟩
0
) ( ̄
ρ
i
−⟨
̄
ρ
i
⟩
0
)
2
+ h.o.t
(31)
Evaluating the mean of eq. (31) yields:
⟨F
( ̄
ρ
i
)
⟩
0
=
F
(
⟨
̄
ρ
i
⟩
0
) +
1
2
F
′′
(
⟨
̄
ρ
i
⟩
0
)
⟨
̄
ρ
2
i
⟩
0
−⟨
̄
ρ
i
⟩
2
0
+ h.o.t
(32)
If only the first term on the right-hand-side of eq. (32) is retained, we obtain a first-order approximation,
i.e.
:
⟨F
( ̄
ρ
i
)
⟩
0
≈ F
(
⟨
̄
ρ
i
⟩
0
)
(33)
This approximation relies on a point estimate to calculate the mean. In particular, if the embedding energy function
F
( ̄
ρ
i
) is convex, this first-order approximation represents a lower bound of the true value (
i.e.
Jensen
'
s inequality).
Following the same approximation, the derivatives of
⟨
V
F
⟩
0
with respect to
̄
q
i
∗
results in:
∂
⟨
V
F
⟩
0
∂
̄
q
i
∗
≈
∂
∂
̄
q
i
∗
X
i
∈
I
M
F
M
(
⟨
̄
ρ
i
⟩
0
) +
X
i
∈
I
H
χ
i
F
H
(
⟨
̄
ρ
i
⟩
0
)
!
=
X
i
∈
I
M
∂
F
M
(
⟨
̄
ρ
i
⟩
0
)
∂
̄
q
i
∗
+
X
i
∈
I
H
χ
i
∂
F
H
(
⟨
̄
ρ
i
⟩
0
)
∂
̄
q
i
∗
=
X
i
∈
I
M
F
′
M
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
∗
+
X
i
∈
I
H
χ
i
F
′
H
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
∗
.
(34)
Where
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
∗
can be evaluated using eqs. (24) to (26). The numerical evaluation of
∂
⟨
V
F
⟩
0
∂
̄
q
i
∗
is hence explicitly detailed
using the first summation of eq. (34) if
i
∗
=
i
and
i
∈
I
M
:
X
i
∈
I
M
F
′
M
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
∗
≈
X
i
∈
I
M
F
′
M
X
j
∈
I
M
∩B
i
(
r
c
)
j
̸
=
i
NQ
X
k
=1
ρ
M
(
r
ij,k
)
W
k
+
X
j
∈
I
H
∩B
i
(
r
c
)
j
̸
=
i
χ
j
NQ
X
k
=1
ρ
H
(
r
ij,k
)
W
k
...
X
j
∈
I
M
∩B
i
(
r
c
)
j
̸
=
i
NQ
X
k
=1
ρ
′
M
(
r
ij,k
)
dr
ij,k
d
q
i,k
∗
W
k
+
X
j
∈
I
H
∩B
i
(
r
c
)
j
̸
=
i
χ
j
NQ
X
k
=1
ρ
′
H
(
r
ij,k
)
dr
ij,k
d
q
i,k
∗
W
k
(35)
4
We want to compute the second order derivative of
⟨
V
F
i
⟩
0
with respect to
̄
q
i
:
∂
2
⟨
V
F
i
⟩
0
∂
̄
q
2
i
≈ F
′′
M
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
⊗
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
+
F
′
M
(
⟨
̄
ρ
i
⟩
0
)
∂
2
⟨
̄
ρ
i
⟩
0
∂
̄
q
2
i
if
i
∈
I
M
(36)
∂
2
⟨
V
F
i
⟩
0
∂
̄
q
2
i
≈
χ
i
F
′′
H
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
⊗
∂
⟨
̄
ρ
i
⟩
0
∂
̄
q
i
+
χ
i
F
′
H
(
⟨
̄
ρ
i
⟩
0
)
∂
2
⟨
̄
ρ
i
⟩
0
∂
̄
q
2
i
if
i
∈
I
H
(37)
The gradient of
⟨
V
F
⟩
0
with respect to ̄
σ
i
∗
results in::
∂
⟨
V
F
⟩
0
∂
̄
σ
i
∗
≈
∂
∂
̄
σ
i
∗
X
i
∈
I
M
F
M
(
⟨
̄
ρ
i
⟩
0
) +
X
i
∈
I
H
χ
i
F
H
(
⟨
̄
ρ
i
⟩
0
)
!
=
X
i
∈
I
M
∂
F
M
(
⟨
̄
ρ
i
⟩
0
)
∂
̄
σ
i
∗
+
X
i
∈
I
H
χ
i
∂
F
H
(
⟨
̄
ρ
i
⟩
0
)
∂
̄
σ
i
∗
=
X
i
∈
I
M
F
′
M
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂
̄
σ
i
∗
+
X
i
∈
I
H
χ
i
F
′
M
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂
̄
σ
i
∗
.
(38)
Finally, we can compute the derivatives of
⟨
V
F
⟩
0
with respect the molar fraction
χ
i
∗
∀
i
∗
∈
I
H
:
∂
⟨
V
F
⟩
0
∂χ
i
∗
≈
X
i
∈
I
M
∪
I
H
F
′
M
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂χ
i
∗
if
i
∈
I
M
and
i
∗
̸
=
i
χ
i
F
′
H
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂χ
i
∗
if
i
∈
I
H
and
i
∗
̸
=
i
F
H
(
⟨
̄
ρ
i
⟩
0
) +
χ
i
F
′
H
(
⟨
̄
ρ
i
⟩
0
)
∂
⟨
̄
ρ
i
⟩
0
∂χ
i
∗
if
i
∈
I
H
and
i
∗
≡
i
(39)
2.4 Pair energy term of the ADP
The contribution of the pair energy
φ
to the ADP is given by:
⟨
V
φ
⟩
0
=
1
2
X
i,j
∈
I
M
j
̸
=
i
⟨
φ
MM
(
r
ij
)
⟩
0
+
1
2
X
i
∈
I
M
,j
∈
I
H
j
̸
=
i
χ
j
⟨
φ
MH
(
r
ij
)
⟩
0
+
1
2
X
i
∈
I
H
,j
∈
I
M
j
̸
=
i
χ
i
⟨
φ
HM
(
r
ij
)
⟩
0
+
1
2
X
i,j
∈
I
H
j
̸
=
i
χ
i
χ
j
⟨
φ
HH
(
r
ij
)
⟩
0
(40)
Where
φ
MM
(
r
ij
),
φ
MH
,
φ
HH
are experimentally-fitted 5th-order polynomial which depends on
r
ij
, with first and
second derivatives
φ
′
(
r
ij
) and
φ
′′
(
r
ij
). Individual terms like
⟨
φ
MH
⟩
0
are evaluated using:
1
2
X
i
∈
I
M
,j
∈
I
H
j
̸
=
i
χ
j
⟨
φ
MH
(
r
ij
)
⟩
0
≈
1
2
X
i
∈
I
M
j
∈
I
H
∩B
i
(
r
c
)
j
̸
=
i
χ
j
NQ
X
k
=1
φ
MH
(
r
ij,k
)
W
k
(41)
Taking derivatives of
⟨
V
φ
⟩
0
with respect
̄
q
i
∗
:
∂
⟨
V
φ
⟩
0
∂
̄
q
i
∗
=
1
2
X
i,j
∈
I
M
j
̸
=
i
∂
⟨
φ
MM
(
r
ij
)
⟩
0
∂
̄
q
i
∗
+
1
2
X
i
∈
I
M
,j
∈
I
H
j
̸
=
i
χ
j
∂
⟨
φ
MH
(
r
ij
)
⟩
0
∂
̄
q
i
∗
+
1
2
X
i
∈
I
H
,j
∈
I
M
j
̸
=
i
χ
i
∂
⟨
φ
HM
(
r
ij
)
⟩
0
∂
̄
q
i
∗
+
1
2
X
i,j
∈
I
H
j
̸
=
i
χ
i
χ
j
∂
⟨
φ
HH
(
r
ij
)
⟩
0
∂
̄
q
i
∗
.
(42)
Particularising for
⟨
φ
MH
(
r
ij
)
⟩
0
we have the following first order derivatives:
∂
⟨
φ
MH
(
r
ij
)
⟩
0
∂
̄
q
i
∗
=
⟨
φ
′
MH
(
r
ij
)
r
ij
r
ij
⟩
0
,
if
i
∗
=
i,
(43)
∂
⟨
φ
MH
(
r
ij
)
⟩
0
∂
̄
q
i
∗
=
− ⟨
φ
′
MH
(
r
ij
)
r
ij
r
ij
⟩
0
,
if
i
∗
=
j.
(44)
5
eqs. (43) and (44) are evaluated numerically using:
∂
⟨
φ
MH
(
r
ij
)
⟩
0
∂
̄
q
i
∗
≈
NQ
X
k
=1
φ
′
MH
(
r
ij,k
)
r
ij,k
r
ij,k
W
k
,
if
i
∗
=
i,
and
j
∈
I
H
∩B
i
(
r
c
)
(45)
∂
⟨
φ
MH
(
r
ij
)
⟩
0
∂
̄
q
i
∗
≈ −
NQ
X
k
=1
φ
′
MH
(
r
ij,k
)
r
ij,k
r
ij,k
W
k
,
if
i
∗
=
j,
and
j
∈
I
H
∩B
i
(
r
c
)
.
(46)
Second order derivatives of
⟨
φ
MH
(
r
ij
)
⟩
0
:
∂
2
⟨
φ
MH
(
r
ij
)
⟩
0
∂
̄
q
2
i
∗
=
⟨
φ
′′
MH
(
r
ij
)
r
ij
⊗
r
ij
r
2
ij
⟩
0
+
⟨
φ
′
MH
(
r
ij
)
d
2
r
ij
d
r
2
ij
⟩
0
,
if
i
∗
=
i,
(47)
∂
2
⟨
φ
MH
(
r
ij
)
⟩
0
∂
̄
q
i
∗
∂
̄
q
j
∗
=
− ⟨
φ
′′
MH
(
r
ij
)
r
ij
⊗
r
ij
r
2
ij
⟩
0
− ⟨
φ
′
MH
(
r
ij
)
d
2
r
ij
d
r
2
ij
⟩
0
,
if
i
∗
=
i, j
∗
=
j,
(48)
∂
2
⟨
φ
MH
(
r
ij
)
⟩
0
∂
̄
q
j
∗
∂
̄
q
i
∗
=
− ⟨
φ
′′
MH
(
r
ij
)
r
ij
⊗
r
ij
r
2
ij
⟩
0
− ⟨
φ
′
MH
(
r
ij
)
d
2
r
ij
d
r
2
ij
⟩
0
,
if
i
∗
=
i, j
∗
=
j,
(49)
∂
2
⟨
φ
MH
(
r
ij
)
⟩
0
∂
q
2
i
∗
=
⟨
φ
′′
MH
(
r
ij
)
r
ij
⊗
r
ij
r
2
ij
⟩
0
+
⟨
φ
′
MH
(
r
ij
)
d
2
r
ij
d
r
2
ij
⟩
0
,
if
i
∗
=
j.
(50)
The gradient of
⟨
V
φ
⟩
0
with respect to ̄
σ
i
∗
results in::
∂
⟨
V
φ
⟩
0
∂
̄
σ
i
∗
=
∂
∂
̄
σ
i
∗
1
2
X
i,j
∈
I
M
j
̸
=
i
⟨
φ
MM
(
r
ij
)
⟩
0
+
1
2
X
i
∈
I
M
,j
∈
I
H
j
̸
=
i
χ
j
⟨
φ
MH
(
r
ij
)
⟩
0
+
1
2
X
i
∈
I
H
,j
∈
I
M
j
̸
=
i
χ
i
⟨
φ
HM
(
r
ij
)
⟩
0
+
1
2
X
i,j
∈
I
H
j
̸
=
i
χ
i
χ
j
⟨
φ
HH
(
r
ij
)
⟩
0
=
1
2
X
i,j
∈
I
M
j
̸
=
i
∂
⟨
φ
MM
(
r
ij
)
⟩
0
∂
̄
σ
i
∗
+
1
2
X
i
∈
I
M
,j
∈
I
H
j
̸
=
i
χ
j
∂
⟨
φ
MH
(
r
ij
)
⟩
0
∂
̄
σ
i
∗
+
1
2
X
i
∈
I
H
,j
∈
I
M
j
̸
=
i
χ
i
∂
⟨
φ
HM
(
r
ij
)
⟩
0
∂
̄
σ
i
∗
+
1
2
X
i,j
∈
I
H
j
̸
=
i
χ
i
χ
j
∂
⟨
φ
HH
(
r
ij
)
⟩
0
∂
̄
σ
i
∗
.
(51)
Finally, we can compute the derivatives of
⟨
V
φ
⟩
0
with respect the molar fraction
χ
i
∗
∀
i
∗
∈
I
H
:
∂
⟨
V
φ
⟩
0
∂χ
i
∗
=
1
2
X
i
∈
I
M
,j
∈
I
H
j
̸
=
i
∂χ
j
∂χ
i
∗
⟨
φ
MH
(
r
ij
)
⟩
0
+
1
2
X
i
∈
I
H
,j
∈
I
M
j
̸
=
i
∂χ
i
∂χ
i
∗
⟨
φ
HM
(
r
ij
)
⟩
0
+
1
2
X
i,j
∈
I
H
j
̸
=
i
∂χ
i
χ
j
∂χ
i
∗
⟨
φ
HH
(
r
ij
)
⟩
0
(52)
2.5 The Dipole contribution to the potential
The contribution of the dipole distortion
μ
to the ADP is given by:
⟨
V
μ
⟩
0
=
1
2
X
i
∈
I
M
∪
I
H
⟨
μ
i
·
μ
i
⟩
0
,
(53)
where
μ
i
=
X
i
∈
I
M
, j
∈
I
M
j
̸
=
i
u
MM
(
r
ij
)
r
ij
+
X
i
∈
I
M
, j
∈
I
H
j
̸
=
i
χ
j
u
MH
(
r
ij
)
r
ij
if
, i
∈
I
M
μ
i
=
X
i
∈
I
H
, j
∈
I
M
j
̸
=
i
χ
i
u
HM
(
r
ij
)
r
ij
+
X
i
∈
I
H
, j
∈
I
H
j
̸
=
i
χ
i
χ
j
u
HH
(
r
ij
)
r
ij
if
i
∈
I
H
(54)
6