minor empirical modifications provide a reasonable way to do this
4.1 Model of the potential energy
4.1.2 Potential energy of anionic clusters
The contributions to the potential energy used by Yourshawet al.[3]were employed in a similar fashion. In detail, these included
• binary potentials as laid out in Table C.1 and Table C.4,
• Axilrod-Teller triple dipole interactions as described in the previous section on neutral clusters with the appropriate parameters from Table C.6,
• the interaction of the halide anion with the exchange quadrupole formed by pairs of closed-shell noble gas atoms using the parameters from Table C.5, and
• the interaction of the halide anion with those di- and quadrupole moments in-duced by dispersion.
For the latter, the expressions given by Yourshawet al.[3]were simplified to yield
Vdisp =qx
∑
i6=q
C6iiBi αi
ri−3
3
2~ri·
∑
j6=q,i
~rijr−ij8−1 8
∑
j6=q,i
3
"
~ri·~rij ririj
#2
−1
r−ij6
where~ri is the vector from the noble gas atom to the halide and~rij from the jth to the ith noble gas atom, analogous to the Axilrod-Teller potential in Figure 4.2. The corresponding particle parameters are given in Table C.6.
Finally, the quantitatively most important – except for the binary potentials – and simultaneously from a computational standpoint most costly contribution is induction nonadditivity. The treatment outlined by Yourshawet al.[3]can be summarized as fol-lows: The constituents of the cluster, i.e. the halide anion and the noble gas atoms, are characterized by point dipole and quadrupole moments, both of which are induced by the point charge of the halide as well as all the respective moments of the atoms.
Terms higher than linear in either of these moments are neglected. The working equa-tions for the determination of the induced moments are (originally labeled Eqs. (21)
4.1. Model of the potential energy through (24) in [3]):
E(αi) =
∑
j6=i
−Tα(ij)qj+Tαβ(ij)µ(βj)− 1
3Tαβγ(ij)Θ(βγj)
(4.7a) E(αβi) =
∑
j6=i
−Tαβ(ij)qj+Tαβγ(ij)µ(γj) −1
3Tαβγδ(ij) Θ(γδj)
(4.7b)
µ(αi) =αiEα(i) (4.7c)
Θ(αβi) =CiEαβ(i) (4.7d)
whereE(αi) and E(αβi) are the electric field and its gradient, respectively, at centeri. Also, qi, µ(αi), andΘ(αβi) are the charge, the α-component of the induced dipole moment and the αβ-component of the induced quadrupole moment of center i. Further, αi and Ci are the dipole and quadrupole polarizabilities of the ith atom. Finally, Tαβij···ν =
∇α∇β· · · ∇νR−ij1and Rij is the vector from thejth to theith atom in accordance with Buckingham’s notation[176].
In their work, Yourshaw et al.[3] continue to describe an iterative scheme for the determination of the induced multipoles. This process could be confirmed to converge satisfactorily and can be sped up in successive computations at slightly modified ge-ometries – i.e. in molecular dynamics calculations or geometry optimizations – by initiating the procedure with the multipole moments obtained for the last geometry.
However, this still remains a very time-consuming procedure and, as an additional disadvantage, analytical derivatives are not available.
A more direct scheme, which also remedies the lack of analytical derivatives, was found to be more efficient. This broadly follows the polytensor formulation proposed by Applequist[177, 178]and is essentially an extension along the lines of the treatment by Kong and Ponder[179, 180].
Due to their linear form and incorporating the inter-dependence of the multipole moments, Equations 4.7a through 4.7d can be cast into a single matrix-vector equation as follows (vectors bear single and matrices double underlining for clarity):
µ = α q + T µ
(4.8)
α is a diagonal matrix containing the dipole and quadrupole polarizabilities. µ is a vector whose elements are the symmetry-unique components of all dipole and quadrupole moments. The T-tensor elements coupling the charge – or, in the more
Chapter 4. Interactions in weakly bound noble gas–halide clusters
general case, charges – and the various multipoles constitute the elements of vector q. Finally, the matrix T contains all the remaining T-tensor components coupling the multipoles with each other. Details on the construction of these vectors and matrices shall be outlined below as some transformations are still to be considered.
At this point, it is obvious that transformations leading to µ =α−1 − T−1
q (4.9)
will provide the multipole components directly. With these known, Yourshawet al.[3]
continue to compute the energy of the charge-multipole interaction as Vind,total =
∑
j
∑
i6=j
qj
1
2Tα(ij)µ(αi)+1
6Tαβ(ij)Θ(αβi)
(4.10)
which is essentially equation (27) in Yourshaw et al.’s work[3]. Note the change of indices, which eventually leads to a change of sign in the dipole-related summand.
Introducing a new vector ˜q, which differs from the original q in the weights 1/2 for components of dipole moments and 1/6 for those of quadrupole moments (1/3 in the case of symmetry-equivalent quadrupole moments), equation Equation 4.10 can be translated into matrix notation:
Vind,total = q˜† µ = q˜†
α−1 − T−1
q (4.11)
As only real values occur, the dagger “†” is used only to denote transpose matrices and row vectors. Thus, the computation of the inductive contribution to the cluster energy can be viewed as a scalar product, or – when looking more closely – as a simple vector-matrix-vector operation.
The notation introduced here allows for making another simplification. Splitting off the weights of the ˜q-vector in a separate diagonal matrix W
˜
q† = q† W (4.12)
4.1. Model of the potential energy leads, when inserted into eq. Equation 4.11, to
Vind,total = q† W
α−1 − T−1
q
= q†
α−1 W−1 − T W−1−1
q
=−1
2q† D−1 q
(4.13)
Since W is diagonal, inversion of this matrix is straightforward. The effect of including W in D is to essentially lift all of the weighting factors in T mentioned above along with the factors arising with symmetry-equivalent elements, e.g. Θ(xyi) andΘ(yxi). All of these are now solely vested in the diagonal elements of D, which arise from α−1 W−1. The coupling matrix D itself is then real-symmetric.
At this point, the expressions for the individual elements of the D-matrix shall be given explicitly. The diagonal elements can directly be obtained as reciprocal polariz-abilities:
Dξξ =
1
αi for a dipole component
3
Ci for a symmetry-unique quadrupole component (e.g. xx)
3
2Ci for a symmetry-non-unique quadrupole component (e.g. xy) whereξ denotes a multipole component, i.e. either a cartesian dipole component µ(αi) or a component of the quadrupole tensor of an atomΘ(αβi). The off-diagonal elements are either zero if referring to the same nucleus or respective elements of the T-tensor where those bearing a column index of a dipole component carry a negative sign:
Dξξ0 =
0 for multipole components at the same center
−Tαβ(ij) ifξandξ0are dipole components
−Tαβγ(ij) ifξbelongs to a quadrupole andξ0to a dipole Tαβγ(ij) ifξbelongs to a dipole andξ0 to a quadrupole Tαβγδ(ij) ifξandξ0are quadrupole components
Here, nucleusibelongs toξand jtoξ0.
Chapter 4. Interactions in weakly bound noble gas–halide clusters
With eq. Equation 4.13 it is also clear how partial derivatives with respect to any given coordinateνcan be obtained analytically:
∂νVind,total =−µ˜† ·∂νq
−1
2µ˜† ∂νD
˜
µ (4.14)
where ˜µ = D−1 q = q D−1†
was used, the “effective multipole moment vector”.
It is, of course, advantageous to store this latter vector in memory after energy calcula-tions for computacalcula-tions of the gradient.
Finally performance aspects need to be discussed. In the original development of the underlying model Vesely pointed out the option of an algebraic treatment and explicitly mentioned the high symmetry of the matrix concerned.[181] Due to the large size of the matrix in the case at hand, the problem seemed intractable at the time, however. Presently and particularly for the problem of noble gas–halide clusters, both pathways are certainly well feasible.
In terms of theoretical scaling limits, both methods can be regarded – in worst case scenarios – as scaling as O(n3). For the algebraic method, this result holds for the least optimized algorithms[182], while in the case of the iterative method it is based on the assumption that the number of iterations scales proportionally with the number of nuclei. The latter is certainly not the case if a suitable updating scheme is imple-mented. Regarding the former, however, it is known[183] that algorithms with better scaling properties exist. For this particular problem, it seems to be generally possible to assume that D is positive-definite and that Cholesky-type methods may be used for better performance.[166, 184] Only at unphysically short interatomic distances do non-positive-definite matrices occur.
For systems of rather limited size, such as the halide–noble gas clusters, overhead considerations are probably of greater importance than asymptotic scaling behavior.
In this respect, the iterative scheme is inherently more prone to producing overhead, due to the up to multi-dimensional T-tensor elements. Recalculation at every step of the iteration will certainly be the worst alternative, but storage of and access to up to four- or five-dimensional arrays is also more costly than the two-dimensional matrices required for the algebraic solution.