Wilfried Meyer

Pavel Rosmus

Weber Thomas

Prof Werner

Thomas Weber

Vibrational motion of the radical cation C + 3 in its degenerate electronic ground state

Keywords: electronic structure, molecular vibration, C3

come

Vibrational motion of the radical cation C3+ in its

degenerate electronic ground state.

I. INTRODUCTION

The cation C + 3 has been intensively investigated since it was suggested that it may play a crucial role in the process of forming larger homonuclear carbon clusters in interstellar space via hydrogenation steps [1]. The room temperature reaction rate of the hydrogen abstraction reaction C + 3 + H 2 →C 3 H + + H was found abnormally low as compared to the rates of other C + n ion clusters [2,3] which has been tentatively linked to a cyclic C + 3 structure. Recent ion trap studies revealed a surprisingly steep increase of the reaction rate down to about 50 K where it suddenly levels off [4]. At this low temperatures even the radiativ association reaction C + 3 + H 2 → C 3 H + 2 + hν becomes competitive. So far there is no well-based explanation of this inverse temperature dependence, speculation envisions either a long-living precursor complex or a structural peculiarity of C + 3 such as very facile bentto-linear isomerization.

Apart from these reaction rates, there are indeed only very few experimental data available. Coulomb explosion experiments [5] have been interpreted as indicating a cyclic structure in contrast to older theoretical predictions of a linear structure. While early measurements of the first IP of C 3 left a range from 11 to 13 eV [6,7], the recently measured photoionization efficiency curve of C 3 [8] shows a clear step for vertical ionization at 11.70 ±0.05 eV and in particular a rather long tail to lower energies wich is believed to indicate unresolved vibrational states of a bent C + 3 . However, a firm experimental determination of its structure and a detailed spectroscopic characterization of C + 3 are still missing. Sunil et al [9] reported the first careful theoretical investigation of ionization potentials of C 3 and noted a particularly strong role of electron correlation for the properties of (linear) C + 3 . It can be traced to rather low-lying nonbonding 1π g orbitals which are not occupied in the Hartree-Fock (HF) configurations of the 2 Σ + u and 2 Σ + g states, 4σ 2 g 3σ u 1π 4 u and 4σ g 3σ 2 u 1π 4 u , respectively. They give rise to large left-right correlation via strong admixture of pair excitations into 1π 2 g or σ u,g 1π g . This results in a reordering of the electronic statess by moving the 2 Σ + u and 2 Σ + g states below the 2 Π u state, the lowest at the HF level. It turned out that SD-CI still gave a wrong ordering of states and the Moller-Plesset perturbation series MP n showed large oscillations of the correlation energies in particular from triple excitations. Only multi reference configuration treatments -MC-SCF as well as MR-CI -produced the correct ordering and a first IP of 11.5 eV in good agreement with the recent experimental value quoted above.

The observation of a complex electronic structure of C + 3 and the suggestion of a cyclic structure [5] triggered a rush of high-level ab initio calculations in which all possible methods and a large number of basis sets were explored, with the main purpose to ascertain the bentto-linear isomerization energy (which we henceforth abbreviate by ∆E e ). Raghavachari [10] used the quadratic configuration interaction method (QCISD(T)) with a relatively small basis set and was the first to predict a bent 2 B 2 global minimum with apex angle of 67 • . This minimum is only 3 kcal/mol lower than the 2 Σ + u minimum of linear C + 3 but separated from it by a barrier which rises about 4 kcal/mol above the linear minimum. An imaginary force constant for the asymmetric stretch mode indicated symmetry breaking even at the level of a restricted QCISD(T). Applying a similar basis and the same method, but based on a UHF reference function, Martin et al [11] proposed 1.9 kcal/mol as best estimate for ∆E e and stressed that this value would be further reduced by a zero point energy (ZPE) difference of about 0.7 kcal/mol. Grev et al. [12] tried the single reference CI method, exploring substitution levels up to quadruples and various active orbital spaces to conclude that the ∆E e should be about 4 kcal/mol, however with an uncertainty of ± 4 kcal/mol. They found a small imaginary bending frequency for the 2 Σ + u state indicating a direct decent to the 2 B 2 minimum and also report symmetry breaking for the HF wavefunctions of the 2 Σ + u and 2 B 2 states which prevails for SD-CI.

The newly established singles and doubles coupled cluster method with perturbative triples (CCSD(T)) was applied by Scuseria [13] with a focus on basis set saturation. Using atomic natural orbitals in symmetryrestricted calculations, ∆E e was found to increase monotonically with basis set size and reached 6.8 kcal/mol for the largest set. The performance of the CCSD(T) was investigated by Watts et al. [14] trying out various schemes for handling triple excitations, including the full CCSDT. They encountered variations up to 5 kcal/mol and concluded that the generally accepted CCSD(T) approximation artificially favors the linear structure by 2.0 kcal/mol. The results for the antisymmetric stretch frequency were indicative for the problems of singlereference methods with C + 3 : It changed from 1431i for CCSD(T) to 451i for CCSDT but was real for other variants.

Taylor et al. [15] stressed the strong multi-reference character of correlation in the lower electronic states of C + 3 as the primary cause of the problems. They performed a series of CASSCF/MRCI calculations with a full-valence active orbital space, state-averaged orbital optimization, reference configuration selection with varying thresholds for their coefficients and optionally applying a multireference analogue of the Davidson correction. Their best calculation yielded a ∆E e value of 5.18 kcal/mol. A critical evaluation of various energy increments and an estimate for the deficiency of their basis set data led them to propose an energy difference of 5.2 kcal/mol as best estimate, with a confidence interval from 4.2 to 6.7 kcal/mol. This nicely includes the estimate of ref. [12] at the lower end and the best calculated value of ref. [13] at the upper end and it has not been challenged since. It should not be overlooked, however, that they encountered imaginary frequencies for antisymmetric stretch of the linear species even at their full-valence CASSCF level.

The problems encountered with all the methods tried seem to have so far prevented any investigation of larger parts of the ground state energy surface. Only harmonic frequencies have been reported in some of the work discussed above, as they come along with geometry optimization. However, for bending and antisymmetric stretch vibrational modes they are of very limited use in view of the obvious strong anharmonicities. Symmetric stretch frequencies for the linear form vary between 1154 cm -1 from CASSCF [15] and 1261 cm -1 from SDCI [12], not very different from the experimental C 3 frequency of 1225 cm -1 . For the bent form, symmetric stretch frequencies reported from well correlated treatments vary between 1590 and 1724 cm -1 . The only attempt to provide a global potential energy surface has recently been made by Wang et al. [16]. An analytical expression was derived from a straight application of the many-body expansion theory [17] of polyatomic interaction potentials. The parameters were taken from potential curves of the dissociation products C 2 and C + 2 and from the harmonic force constants pertaining to the 2 B 2 global minimum. Since no reference was made to the peculiarities of the electronic structure of C + 3 this analytical surface can not expected to be really useful.

For the sake of completeness we mention some recent theoretical results which have been generated in the context of experimental work. McAnoy et al. [18] performed single-point CCSD(T) calculations with an aug-cc-pVTZ basis, albeit for B3LYP optimized geometries. While the ∆E e of 5.7 kcal/mol does fit into the picture of previous work, a rather high barrier of 2.9 kcal/mol is suggested. Fura et al. [19] applied a quite similar but unrestricted scheme to a larger number of ion clusters and reported a rather small ∆E e of only 1 kcal/molr. Nicolas et al. [8] presented CASSCF/MR-CI potential curves along the bending coordinate while fixing the CC bond length to a value favouring the linear form. Thus, the suggested ∆E e of 2.2 kcal/mol as well as a flat, suspiciously looking hump may be artifactual.

All these variations in the results for ∆E e and the transition state barrier indicate that great care has to be taken if a potential surface is to be generated that can be trusted. It is the purpose of the present paper to establish a global potential energy surface (PES) which is reliable up to some lower vibrational levels of linear C + 3 , i.e. the conformation in which it is most likely formed. Both sheets of the adiabatic PES are obtained from MR-CI calculations in numerical form on a grid that covers all conformations which are accessible within the given energy range (Section II). These data are used to generate an analytical PES suitable for integration with global vibrational basis functions (Section III). A transformation based on a J-T parametrization of the potential in the cyclic region and on the dipole moment matrix in the linear regian provides a PES in diabatic form for efficient account of vibronic coupling (Section IV). Vibrational states are finally obtained from variational calculations wich use hyperspherical coordinates for proper treatment of symmetry requirements and boundary conditions at D 3h and D ∞h geometries. They are characterized in terms of their distribution over the cylic and linear domains and in terms of approximate vibrational quantum numbers (Section V).

With respect to analytical fits and vibrational states we follow the technical procedures previously applied to degenerate states of Na 3 [20] and Li 3 [21][22][23]. Measured vibrational band positions of exited E" states could be reproduced to within 2 cm -1 and excellent agreement 3 with rotational resolved spectra was found. Of course, electronic structure calculations were facilitated by the small number of valence electrons and vibrational motion turned out to be confined to the lower sheets of the E" states. In the more recent analysis of vibrational structures in the 2 E ← 2 A 1 absorption spectrum of B 3 we arrived at a satisfactory simulation of a spectrum taken in a Ne matrix, but it was difficult to draw conclusions about the accuracy achieved for vibrational levels [24]. In the case of C + 3 the low-lying conical intersection and the presence of the three local linear minima will substantially complicate the vibrational analysis. Since there is not yet a convincing suggestion of how C + 3 will eventually be seen spectroscopically we restrict the present investigation to a vibrational analysis of the ground state. It is hoped that this might help understand the unexpected features of its reactivity.

II. AB INITIO ELECTRONIC STRUCTURE CALCULATIONS

From the discussion of the previous investigations of C + 3 it seems clear that a global PES requires a computational procedure which is based on a multiconfiguration reference wavefunction. Lacking a multistate, multi-configuration coupled cluster program, we have decided to use the multi-reference configuration interaction method (MR-CI) as implemented in the MOL-PRO program package [START_REF] Werner | MOLPRO: a package of ab initio programs[END_REF]. Its particular strength is the efficiency derived from the concept of internally contracted configurations, i.e. singles and doubles substitution operators are applied to the reference function as a whole. This was first formulated by Meyer [START_REF] Meyer | Methods of Electronic Stucture Theory[END_REF] but developed into an fully operative algorithm and efficiently implemented by Werner and Reinsch [START_REF] Werner | [END_REF].

The reference function is obtained from CAS-SCF calculations. The definition of the active orbital space is not trivial in the present case and no perfect solution has been found in test calculations along the minimum energy path (MEP) between the D 3h and C 2v conformations. The reason is the C-C bond breaking upon opening the apex angle from 60 • to 180 • which is -as well known for multiple bonds in diatomic molecules, e.g. N 2 [28] connected with a series of avoided crossings of Rydbergtype with valence-type orbitals and, consequently, crossings of the corresponding configurations. This is seen in Fig. 1 which shows CAS-SCF potential curves of the lower excited states of C + 3 . The optimized orbitals for a preselected set of states therefore undergo rather sudden changes which may affect also the PES of lower-lying states. The resulting humps may be rather large at the CAS-SCF level (e.g. 100-200 cm -1 for the lower component of the ground state). They are effectively reduced at the MR-CI level but may still be at a disturbing size of 10-20 cm -1 . Since extending the set of states did not improve the situation significantly, we restricted this set to the two components of the ground state which is suf- α

E x ' -A 1 -Σ g E y ' -B 2 -Σ u E x " -B 1 -Π u E y " -A 2 -Π u FIG. 1:
Cuts of PESs of lower electronic states of C + 3 along the ground state MEP, from CAS-SCF with basis AVTZ ficient to avoid the unphysical symmetry breaking. The avoided crossing of the upper component with the next higher state appears at about 75 • of the apex angle, i.e. in a region where the upper component is not significant in itself for the low-energy vibrations. Thus, the effect of the crossing could be diminished by introducing an angle-dependent weighting of the two components which reaches a 1:1 ratio only above 100 • . After some experimenting we found as best compromise a set of 13 active orbitals (more precisely 11a', 3a" orbitals in C s symmetry or 6a 1 , 5b 2 , 2b 1 , 1a 2 orbitals in C 2v symmetry), which is one orbital in excess of the full-valence set used previously. In order to reduce computing expense, an occupation restriction was imposed to the configurations of the reference function in the MR-CI: only two electrons were allowed to occupy simultaneously the four uppermost active orbitals. This restriction has an insignificant effect on the surface even though it may introduce an additional risk for small bumps.

Standard correlation consistent basis sets are used throughout [29]. For the full set of grid points, the cc-pVQZ basis was used [START_REF] Werner | MOLPRO: a package of ab initio programs[END_REF]. For a significantly smaller subset of points, the VTZ and V5Z basis sets were applied to perform complete basis set extrapolation (CBS) [30]. Since dynamic electron correlation is significantly larger for the more compact cyclic structures than for linear structures, basis set deficiencies are also more serious there. Thus, CBS extrapolation increases ∆E e by 13%, up from 2125 cm -1 to 2401 cm -1 . The truncation of the configuration space inherent to a SD-MR-CI expansion has the opposite trend since configuration mixing is stronger for linear structures. Thus, approximately account for higher substitutions by Davidson's correction as implemented in the MR-CI of MOLPRO, lowers ∆E e by 164 cm -1 to our final value of 2238 cm -1 . Adding up both extrapolations leads to a bending potential curve which is quite similar to the uncorrected one. Still, the small barrier for bending changes from 359 cm -1 to 289 method accounts better for dynamical correlation then uncorrected MR-CI but underestimates the effects from near-degenerate configuration mixing. Both differences add up, yielding after CBS extrapolation an ∆E e isomerization energy of 2660 cm -1 and in particular a strongly reduced linear-to-cyclic barrier of only 40 cm -1 . We note again that the CCSD(T) treatment converges nicely for C 2v geometries (as those of Fig. II) provided the occupied orbitals obey C 2v symmetry. But it is difficult to devise a proper one-determinant wavefunction for general C s geometries: convergence is not stable and may end up in unphysically distorted wavefunctions.

The much discussed question of the stability of the linear minima with respect to asymmetric stretch deformation is elucidated in Fig. 3 and seen to be somewhat academic. The harmonic force constant turns out to be very small indeed: It is weakly positive at equilibrium bond distances but soon changes sign with symmetric bond elongation, i.e. we find a valley-to-ridge transition. The effective force constant related to the fundamental frequency for asymmetric stretch is then largely determined from the positiv quartic force constant.

The molecular constants derived from the MR-CI PES are collected in Tab. II. The structural parameters agree quite well with the better calculations of the literature. Our value for the ∆E e of 2237 cm -1 = 6.4 kcal/mol is within the confidence range proposed by Taylor et al. [15] but close to its upper end. We are certain that the

III. ANALYTICAL REPRESENTATION OF POTENTIAL ENERGY SURFACES

MR-CI calculations as just described have been performed for a grid of non-redundant conformations defined by varying the bond lengths in steps of ∆r = 0.05 Bohr, starting from r = 2.56 Bohr, which is close to r at the D 3h geometry with minimum energy, and stopping only after the relative energy exceeded 4000 cm -1 above the D 3h minimum. With some extra points for bending from linear conformations, this amounted to 250 grid points.

As outlined in section V, the variational calculation of vibrational states involves basis functions that are global in the sense that they span the complete space of the cyclic hyperspherical coordinates θ and φ , e.g. harmonic functions of φ for pseudorotation. Numerical integration of the potential then requires substantial extension beyond the ranges covered by the grid of the electronic structure calculations. Therefore, we like to have the potentials in an analytical form which is well-behaved everywhere. Moreover, since finite basis sets are used which do not allow a vibrational wavefunction to be exactly zero over finite ranges in coordinate space, we have to make sure that the potential does nowhere assume very large values, including the ranges where two C atoms come very close to each other.

In previous applications to A 3 molecules (with JT distorted potentials) we found it convenient to base an analytical expression on the three bond distances since they are most directly related to bonding forces. In order to further adapt the coordinates to the typical shape of bonding potentials, we introduced 'Morse coordinates' defined as r = (1 -exp(-β(r -r o )))/β [31]. This transforms a Morse potential into a parabola with a finite range for positive r and usually allows to represent a bonding potential by a low-order polynomial of r. As compared to functions of the apex angle, the steep C + 3 potentials around 60 • are widened and the shallow potentials above 130 • are compressed. Thus, β has to be chosen with care so that bending forces can be accounted for also in the linear domain. For the present case we found β = 0.90 suitable.

The coordinates ri are combined to form symmetryadapted coordinates in the usual form:

ra = (r 1 + r2 + r3 )/ √ 3, (3.1) rx = (2r 3 -r1 -r2 )/ √ 6 , ry = (r 1 -r2 )/ √ 2 (3.2) re = (r 2 x + r2 y )/ √ 2, φ = arctan(r y /r x ) (3.3)
Symmetry is solely carried by functions of φ : cos(k φ), sin(k φ) form a pair of E species if (k mod 3) = 0, otherwise they are A 1 + A 2 species. Since the adiabatic PES belong to the A 1 representation, the following linear ansatz for the potential is appropriate:

V = ijk v ijk ri a rj e cos(3k φ) (3.4) i, j, k = 0, 1, 2, .. ; i + j ≤ n ; k ≤ k max (j) (3.5)
re is by definition non-negative and is zero only for D 3h geometries. At these points, the two components of an E state PES are smooth continuations of each other if taken as functions of rx and ry due to the fact that they are eigenfunctions of a regular Hamiltonian. Thus, the ansatz above can describe both components with a single set of v ijk if the range of re is simply extended to negative values, which then correspond to the upper component PES. In the close neighborhood of re = 0, v ijk should vanish for 3k > j but we only enforce this for j = 1. Admittedly, such a one-center polynomial expansion, which worked very well for our previous cases Li 3 and Na 3 , is somewhat under stress in the present case where linear geometries matter which have little to do with the conical intersection. In order to ensure sufficient flexibility, we adopted the pragmatic choice of k ≤ 0, 1, 4 for j ≤ 1, j = 2, j ≥ 3, respectively. The linear ansatz has the obvious advantage that it can easily handle a large number of parameters. It has been preferred over a non-linear fit in terms of the elements of a 2x2 diabatic potential matrix because the upper potential has less quality then the lower. In the present case we have chosen n=8 and i, j ≤ 7. The total number of terms then amounts to 132. It should be noted that powers of up to 7 in a polynomial expansion are acceptable here because the relevant ranges of ra and re are effectively limited to little more than 1 by the Morse transformation, and extended extrapolation is regularized as described below.

For a sensible least-squares fit the input energies should be weighted according to their relevance. We have chosen the weight function w(e) = (1 + e/e 0 ) -2 , where e = E -E min , i.e. squares of errors of energy differences are minimized for e < e 0 but squares of relative errors for e > e 0 . Above a certain threshold, e > e max , the weight is further reduced by the linear cut-off factor (2-e/e max ) and it is set to zero above 2e max . We have chosen e 0 = 3000 cm -1 and e max = 5000 cm -1 , respectively. Around α = 75 • a small ridge in the error surface was observed (see discussion in previous chapter) and the weights were reduced to half in that region. The final result comes with a rms error of 4.5 cm -1 for fitting the weighted MR-CI energies while the MR-CI-Q energies are fitted with an rms error of 2.9. (The Davidson correction reduces the ridge just mentioned.) Largest deviations amount to less then 15 cm -1 . The quality of the fit can also be judged by inspection of Fig. 3 for the most critical region around linear C + 3 . We need to deal with the fact that there are large regions of conformation space for which the PES is not controlled by calculated energies and for which the polynomial fit may go wild. Extending the grid does not help since the range of the fit can not easily be extended without loss of quality in the relevant region. Therefore, the analytical fit needs to be regularized in the regions which are beyond the grid but are accessed by (our) vibrational basis functions. The transformed coordinates ra and re are quite convenient for simple regularizations [20-22, 24, 31]: It is usually possible to determine an ellipsoidal cut-off boundary with respect to a reference point rao , reo which encloses all regions where the PES is low and well represented by the fit. Beyond this boundary the PES may be defined with reference to its values at the boundary or close to it. E.g., if the boundary is defined by requiring where ra = ra -rao and re = re -reo , the external potential may simply be derived by radial scaling of the coordinates ra , re whith t = rcut /D(r a , re ):

E(r a , re , φ) = E(r a t + rao , re t + reo , φ)(2 -t 2 ) (3.7)
This scheme has been chosen here. In addition, the potential is not allowed to exceed an energy E max , chosen here as 14000 and 19000 cm -1 for lower and upper surface, respectively. It has been made sure that the calculated vibrational energies are not affected by reasonable changes of all of these parameters.

IV. VIBRONIC COUPLING

In the region of the Jahn-Teller well there is significant vibronic coupling. In adiabatic framework the electronic coupling factor < Φ a 2 |∂/r∂φ|Φ a 1 > has a pole at the conical intersection, r → 0, which has to be removed by the vibrational factor < χ a 2 |∂/r∂φ|χ a 1 >. For practical reasons, among them program limitations, we prefer a diabatic framework in which the vibronic coupling turns into simple potential coupling. Diabatic electronic states are the specific combinations of the two adiabatic states which have minimal dependense on the angle φ. Applying a 2x2 rotation matrix the two adiabatic PESs V ad -, V ad + from the MR-CI calculations turn into the diabatic 2x2 matrix PES

2V d = V ad + + V ad -+ V ad + -V ad - cos δ sin δ sin δ -cos δ
where δ is twice the angle for rotating the wavefunctions.

The diabatic wavefunctions may be chosen to be identical to the adiabatic ones, δ = 0, at the C 2v geometries linked to φ = 0. In the cyclic region around the conical intersectioni, the angle δ can be determined from the familiar expansion of the Hamiltonian in a power series of the polar deformation coordinate r [32,33]. Since the diabatic wavefunctions Φ d x , Φ d y form a pair of E symmetry, their Hamiltonian matrix elements generate the representation A 1 ⊕ E and may be approximated up to cubic terms in r by only in the linear J-T region but generally at the angles φ = nπ/3, n=0,1,.., i.e. for all C 2v geometries. This has also to be required for a consistent continuation of δ(φ) beyond the J-T well.

H xx + H yy = a 2 r 2 + a 3 r 3 cos(3φ) (4.
At the other ends of the bending potential valleys, for linear C + 3 the separation of the surfaces is 6032 cm -1 , much larger than the vibrational energies considered there. Still, vibronic coupling is not negligible for asymmetric stretch motion due to a strong coupling between two diabatic states correlating to 2 Σ + u and 2 Σ + g , respectively, which is obvious from the shape of the asymmetric stretch potential, Fig. 3. As asymmetric stretch is always described by the angle φ, diabatization is again a matter of defining δ(φ). The pertinent quantity is the derivative (∂δ/∂φ) 0 since only small ranges of φ around the C 2v geometries are probed. (For linear structures the This procedure seems validated by the facts that µ 2 12 + ∆µ 2 is nearly independent of φ and that the resulting diabatic potential has nearly quadratic diagonal elements with a nearly linear coupling element. In Fig. 5 abatization procedure is illustrated for linear geometries but the picture changes little for apex angles down to about 120 • . Remarkable large values of δ/φ ≈ 17 are obtained which are somewhat difficult to reconsole with the requirement that δ → φ for φ → nπ/3. After some experimenting we have settled at an exponential approximation for π/3 -δ. This modeling, indicated by dashed curves in Fig. 5, does not affect the vibrational levels, however. The function δ(φ) is also shown in Fig. 4. The effect of the vibronic coupling is quite substantial for the zero-point energy of the asymmetric stretch mode: In the linear 1-D model, see Fig. 5, this energy is increased by about 13% or 85 cm -1 . After averaging over symmetric stretch and bending motions, there remains a shift of 65 cm -1 for low vibrational states that populate linear structures, as compared to the states that populate cyclic structures. Note that the level v 3 = 1 is hardly changed so that the fundamental transition is reduced by 85 cm -1 .

V. CALCULATION OF VIBRATIONAL STATES

We employ hyperspherical coordinates for ease of taking into account symmetry and proper boundary conditions for cyclic motions. For three equal masses they may be defined by expressing the squares of interatomic distances as [34] (k=1,2,3)

r 2 k = 3 -1/2 ρ 2 (1 + sin θ cos(φ + k )) ; k = k2π/3. (5.11)
where the hyperspherical radius ρ determines the size of the molecule, sin θ gives the strength of deformation and φ describes the pseudorotation motion around points of D 3h symmetry. With this definition θ = 0, π correspond to equilateral geometries (D 3h ) and θ = π/2 to linear ones (D ∞h , C ∞v ), φ = 0, ±2π/3 gives isosceles triangles (C 2v ).

For convenient visualization of the actual motion with changing φ it is useful to note the relations R = 3 -3/4 ρ cos(θ/2) ; r = 3 -3/4 ρ sin(θ/2) (5.12)

where R is the distance of the atoms from the center of mass in the D 3h reference geometry and r is the polar radius for the motion around this reference point, see Fig. 6. We shall use (r,φ) as coordinates for 2-D plots work. This latter version is applied here since the upper component of the E ground state has to be taken into account if vibrational states are investigated that reach out to linear conformations. The particular strength of our scheme is the tailoring of the vibrational basis functions to fit the given potential. The full-dimensional basis set is generated from lowerdimensional sets by a sequence of contractions based on the eigenfunctions of appropriate parts of the kinetic energy operator and corresponding effective potentials. The effective potentials are defined as lower-dimensional minima of the full-dimensional potential surface and ensure that the basis functions cover adequately all energetically accessible regions of the coordinate space while effectively excluding the regions with high potentials. The primitive 1-D functions are taken to be Morse functions for the radius ρ, Jacoby polynomials for the angle θ and harmonic functions cos(mφ), sin(mφ) with integer or half-integer modulus m for the angle φ . Two-dimensional contracted functions may either be defined for the θ -φ subspace or the ρ -θ subspace. Primitive functions as well as contracted functions and their derivatives are represented on a grid and Gauss integration is used for the Hamiltonian. In the present calculations we have used the following sets of primitive and contracted (numbers in parenthesis) functions: ρ : 28( 16), θ : 60(54), φ : 132(78), ρ -θ (69). The dimension of the Hamiltonian matrices were 8: Contour map of the X-PES in the r/φ polar plane for ρ =const., i.e. including the MEP of equ. 5.13. The polar radius r is given in equ.5.12 and shown in Fig. 6. Increments are 250 cm -1 up to 1250 cm -1 and 500 cm -1 beyond. Blue line at zero energy = lowest energy of conical intersection; red lines below, green lines above zero energy, respectively. The vibronic wavefunctions are symmetry-adapted products of electronic and vibrational wavefunctions. Since the coordinates ρ and θ are not affected by any symmetry operation, vibrational symmetry is determined only by the factor-functions of the coordinate φ : The pairs (cos(mφ), (-) l sin(mφ)) are e symmetry species if l = (2m modulo 3) = 3, otherwise they generate a 1 and a 2 species, respectively [35,36]. (The phase factor (-) l is required to ensure identical representation matrices of C 3 rotations; l is defined to allow for half-integer values of m required in the adiabatic frameA.) As customary, small case letters are used for the symmetry species of vibrational functions. With the diabatic electronic wavefunctions Φ x , Φ y one finds the following four species of vibronic wavefunctions:

Ψ A1 = Φ Ex χ ex + Φ Ey χ ey Ψ Ex = Φ Ex (χ ex -χ a1 ) + Φ Ey (χ ey + χ a2 ) Ψ Ey = -Φ Ex (χ ey -χ a2 ) + Φ Ey (χ ex + χ a1 ) Ψ A2 = Φ Ex χ ey -Φ Ey χ ex (5.14)
These symmetry species provide the only exact classification of vibronic states. The calculated vibronic energy levels are collected in Tab. V and are displayed in Fig. V. As the shape of the bending potential suggests, we find three groups of vibrational states. The first group, labeled by index b in Tab. V, starts lowest in energy and is confined to strongly bent regions around the conical intersection and the J-T minima. The second group, labeled by l, should appears above energies of 1870 cm -1 and is localized in the linear domain. The third group, labeled by g, is global in the sense that the barrier at about 2200 cm -1 is overcome and both regions are visited. A further characterization of the vibrational states may be attempted in terms of approximate quantum numbers. For group b states they may be derived by reference to the limiting case of a very high pseudorotation barrier, for which the states are superpositions of the three equivalent states of vibrators localized in the three deep wells. The states can then be labeled by the usual quantum numbers for small-amplitude vibration of a bent molecule. This situation is signaled by near degeneracy of pairs of E and A 1 states. The other limiting case, a very low barrier, lends itself to quantum numbers of a free internal rotation and is signaled by near degeneracy of pairs of A 1 and A 2 states. Comparison of the low-energy A 1 -E splittings of 201 and 285 cm -1 , respectively, with the first A 2 -A 1 splitting of 605 cm -1 indicates that the first option is more adequate. Although the pseudorotation barrier of only 376 cm -1 seems rather low, a quite effective separation of the three wells is caused by the fact that the saddle point region is rather narrow in bending direction (see Fig. II). Indeed, the lowest A 2 state, which heavily populates the saddle point region, lies as much as 806 cm -1 above the ground state. Thus, group b states are assigned with vibrational quantum numbers (v 1 , v 2 , v 3 ) for symmetric stretch, bend and asymmetric stretch modes of a C 2v vibrator, respectively. A particular set of such quantum numbers shoul appear twice, that is in E and A 1 symmetry for v 3 even but in A 2 and E symmetry for v 3 odd. Within such a pair, the energy order E' < A' 1 and E' > A' 2 , respectively, is usually observed. However, this assignment becomes less and less stringent ar higher energies as different local vibrators interact across the pseudorotaton barrier. From the energies given in Table II one derives frequencies for fundamental transitions of 1596-1620 , 681-765 and 1245 cm -1 for symmetric stretch, bent and antisymmetric stretch modes, respectively. In the case of group "l" states one finds E +A pairs which are exactly degenerate because the three components around φ = k , k=1-3, correspond to conformations with different center atoms which do not interact at all. Vibrational quantum numbers as used above apply even better to these states but it should be noted that states with v 2 odd do not appear in a rotation-free calculation. We observe only a short progression, v 2 = 0, 2, (4), of lowenergy bending vibrations with a fundamental frequency of only 73 cm -1 (level distance 146 cm -1 ), which is due to the broad barrier of 289 cm -1 at at α = 110 • . The levels v 2 = 4 at 272 and 261 cm -1 , respectively, show already a small splitting from interactions with group b vibrations. The first excitations of symmetric and antisymmetric stretch modes appear with energies 1195 and 1133 cm -1 , respectively, above the lowest linear level. (The latter value differs by -25 cm -1 from the result of the linear 1-D model for vibronic coupling, which is due to vibrational averaging.) Short bending progres-sions start with these excitations.

Finally, group g states spread rather evenly over cyclic and linear domains. The vibrational motions appear as 1-D bending vibrations within the three energy valleys, with rather regular interconnections in the cyclic region. (Quantum mechanics appears to effectively suppresses chaotic vibrational motion in the energy range under consideration here.) The states are simply denoted by the number of nodes along the minimum energy path (not double-counting the nodes as in the linear case). The average spacing between consecutive levels is around 110 cm -1 , but interactions with the pure cyclic and/or pure linear modes of similar energy cause significant variations. The assignment of quantum numbers becomes more and more questionable with increasing energy but the grouping according to domains turned out rather definite. The assignment has been carried out by applying small diagnostic step potentials and monitoring the resulting shifts of the energy levels. Of course, contour plots of the vibrational wavefunctions have been illuminating and very helpful in this regard.

Selected low-energy vibrational wavefunctions are displayed in figures 9 -11. For an are adequately representation of global motions the wavefunctions are shown for the (r,φ; ρ 0 ) surface as used for the potential in Fig. V. Since it appeared difficult to produce comprehensive pictures of the diabatic wavefunctions, we have transformed them to adiabatic ones. Only the component at the lower surface is shown since there is little population of the very narrow upper adiabatic surface. Note that the nonadiabatic wavefunctions are double-valued with a sign change for the second full period of φ. The plots show only the period -π to π and indicate by a bold black line the cut where the wavefunction continues smoothly on the other sheet. As to the symmetry assignment please note that we do not invoke the double group notation. Thus, adiabatic electronic states are considered of A symmetry and the symmetry species of adiabatic vibrational states are identical to that of the vibronic states. The lowest vibrational state then has E symmetry and shows only a single node over the 4π period (j psr = 1/2) while the lowest A 1 and A 2 states have three nodes (j psr = 3/2 ). The latter resides pretty much on top of the pseudorotation barrier. Plots for the ρ -θ plane have been instructive for the assignment of vibrational quantum numbers but are not shown here.

The vibrational analysis provides a vibronic value for the isomerization energy of ∆E ev = 1668 cm -1 ( 4.75 kcal/mol). The significant difference of 570 cm -1 to the purely electronic value of ∆E e = 2238 cm -1 results from the reduction of the ZPE from 1626 cm -1 for the cyclic ground state to the vibrational energy of the lowest linear state of only 1195 cm -1 , which is due to much weaker bending and antisymmetric stretch potentials. We have presented a first attempt to generate a reliable global potential energy surface for the electronic ground state of C + 3 . The internally contracted MR-CI method proved to be adequate to deal with the strong multiconfigurational character of the electronic structure. A large one-electron basis and beyond full-valence active orbital space as well as extrapolations to complete basis and complete configuration space are believed to support a rather uniform accuracy of about ±1 kcal/mol for the lower sheet of the potential surface. The electronic isomerization energy ∆E e is calculated as 2238 cm -1 (6.4 kcal/mol), cyclic and linear structures appear separated by a transition state barrier of only 290 cm -1 .

The variational calculations of vibronic states take full account of strong vibronic coupling in the cyclic (J-T) and the linear domains of the vibrational motion and are converged to better than 1 cm -1 . The complex manyfold of low-energy (< 3000 cm -1 ) states could be classified according to local motions in the cyclic or the linear domains or global motions connecting all low-energy parts of the PES. A more detailed characterization has been attempted in terms of quantum numbers of local vibrators. The vibrational analysis shows that The isomerization energy is reduced by vibrational effects to the vibronic value of 1626 cm -1 (4.75 kcal/mol), due to a rather low ZPE in the linear domain.

The present work is focused on the electronic ground state because of the substantial problems we encountered whith the accurate determination of the upper state wavefunction beyond the J-T region. Work is in progress to derive reliable transition moments for simulation of spectra which may be helpful for identifying C + 3 spectroscopically.

In a final remark we like to address again the low reaction rate for C + 3 + H 2 →C 3 H + + H and the puzzle of its inverse temperature dependence mentioned in the introduction. No new options seem to emerge from the knowledge of potential surface and vibrational manyfold of C + 3 : Only one vibrational level, the A 1 level at 202 cm -1 , is subject to temperature dependent population but its structure hardly differs from the ground vibrational state. More information is required on the ionic collision complex. Some exploratory calculations have revealed the following situation: Multipole Coulomb and induction forces form a van der Waals complex with several local minima, the deepest of which is bound by about 3500 cm -1 . On the one hand, it may be stabilized with respect to the entrance channel by internal vibrational redistribution. On the other hand, it is separated from the reactive channel by a barrier which rises up to about the educt energy -our calculations are not yet conclusive for the sign of the small difference in energy. This barrier may well explain the low reaction rate at room temperature and maybe even its inverse temperature dependence because the decay channel to the educts is effectively closed at low energies by rotational energy redistribution. Work is in progress to map out the van der Waals potential well and accurately locate the barrier to the reaction channel.

cm - 1 FIG. 2 :

 12 FIG.2: Cut of ground state PES along the MEP of the lower sheet from fits of MR-CI energies without and with Davidson correction (-Q) or/and CBS extrapolation. Also included is the CCSD(T) PES (cut along an approximate MEP)

D 2 (

 2 r a , re ) = d a r2 a + d e r2 e + 2d ae ra re =! r2 cut (3.6)

8 )HFIG. 4 :

 84 FIG. 4: Rotation angle δ as function of φ. Black curves: J-T region , α 60-15-105 • ; Red curves: linear domain, α = 105-15-180 •

  below 4 • for v 3 = 0). The rotation angle has been determined from the matrix elements of the dipole moment operator under the assumption that their changes are solely due to changing mixtures of the diabatic states upon asymmetric stretch. We are led to the relation tan(δ) = ∆µ/µ 12 (4.10) ∆µ = |µ 11 (φ) -µ 22 (φ) -(µ 11 (0) -µ 11 (0))|/2

FIG. 5 :

 5 FIG. 5: Transformation to diabatic frame for asymmetric stretch deformation of linear C + 3 . Upper panel: Dipole moment matrix elements and rotation angle. Lower panel: Adiabatic and diabatic potentials; vibrational levels from uncoupled adiabatic and coupled diabatic treatment.

FIG. 6 : 3 FFIG. 7 :

 637 FIG.6: Pseudorotation motion by varying φ with constant ρ, θ or R, r, respectively. Shown is the situation for φ = π/3

  FIG.8: Contour map of the X-PES in the r/φ polar plane for ρ =const., i.e. including the MEP of equ. 5.13. The polar radius r is given in equ.5.12 and shown in Fig.6. Increments are 250 cm -1 up to 1250 cm -1 and 500 cm -1 beyond. Blue line at zero energy = lowest energy of conical intersection; red lines below, green lines above zero energy, respectively.

  symmetries and 3700 for the E symmetry, respectively.

FIG. 9 :

 9 FIG. 9: Vibronic energy levels for XC + 3 . States with bent (b), linear (l) and global (g) domains indicated in red, green and blue, respectively. Tentative assignment by local vibrator quantum numbers (vs, v b , vas); related E-A pairs linked by dotted lines.

12 FIG. 10 :FIG. 11 : 13 FIG. 12 :

 1210111312 FIG. 10: Vibronic states of A 1 symmetry. 201 cm -1 966 cm -1 1571 cm -1

TABLE I :

 I Characteristic points of the X 2 E' PES. and bent structures is real, but it is difficult to qualify the accuracy of the 289 cm -1 obtained, given the large absolute amounts of CBS and Davidson extrapolations of together about 4700 cm -1 .

					5
	barrier between linear
	F o r
	P
	e e r	Methods and Basis sets: a: this work, MR-CI/cc-pVQZ; b: this work, MR-CI-Q/CBS; c: [10] QCISD(T)/6-31G*; d: [12] CISD(TQ)/TZ2P; e: [13] CCSD(T)/5s4p3d2f1g; f: [15] MR-CI-Q/5s3p2d1f. Numbers in brackets: ZPE from harmonic frequencies.
	E + 113 E h a -0.42535 b -0.44094 c -0.40392 d -0.29780 e -0.43724 f -0.32590 i C2v min e R 2 B2 v e C2v a sadd b w	Er cm -1 0 1.322 r1,2 Å 0 1.321 0 1.333 0 1.296 0 1.324 0 1.324 357 1.391 376 1.391	α • 66.9 1692 ZPE ω1(a1) cm -1 cm -1 1785 67.1 1626 1780 67.2 (1811) 1684 71.0 (1759) 1724 67.8 68.0 1675 1612 55.9 55.8 1760
		D 3h min C2v trans b a b a c f D ∞h a min b	1278 1.358 1334 1.358 2487 1.308 104.7 60.0 2526 1.308 110.1 2432 n O l 2026 2120 1.303 180.0 1251 2237 1.303 1195 y	1792 1465 1190 1188
		2 Σ + u	c		1051 1.318
			d		1400 1.283	1261
			e		2383 1.307
			f		1372 1.314	1154
		2 Σ + g	b	Tv	6035

TABLE II :

 II Calculated vibrational energy levels (in cm -1 ) and local vibrator assignments.

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


Johannes M Dieterich

Bernd Hartke
email: hartke@phc.uni-kiel.de

OGOLEM: Global Cluster Structure Optimization for Arbitrary Mixtures of Flexible Molecules A Multiscaling, Object-Oriented Approach

Keywords: global optimization, genetic algorithms, cluster structures, Lennard-Jones, Kanamycin global optimization, genetic algorithms, cluster structures, evolutionary computation, Lennard-Jones clusters, Kanamycin

In practical applications, global cluster structure optimization has so far been limited largely to homogeneous clusters of atoms or small molecules, with little or no choice in the calculation of interparticle forces. We eliminate these limitations by presenting a new program suite OGOLEM that is universal by design, both in cluster composition (including arbitrarily heterogeneous clusters of complicated molecules) and in its interfaces to force calculation backends. This is demonstrated by exemplary applications in two novel fields: strongly heterogeneous Lennard-Jones clusters (ternary, quaternary, quinary), and mixed clusters of the aminoglycoside Kanamycin A with sodium cations.

Introduction

Atomic and molecular clusters are important systems in the physical sciences, both for fundamental research [START_REF] Wales | Energy Landscapes[END_REF] and for nanotechnology applications [START_REF] Bittner | [END_REF]. Structures and properties of sufficiently small clusters differ drastically from the bulk, violating standard chemical intuition. Unfortunately, however, the number of local minima on the potential energy surface for a cluster of n particles increases at least exponentially [3,4] with n. Therefore, unbiased and non-deterministic global structure optimization techniques are required.

To motivate our present contribution to this area, we provide a brief (and possibly biased) sketch of the recent history in this area. One of the first general algorithms employed for this task was simulated annealing (SA) [5]. As a hopefully less compute-intensive approach, genetic algorithms (GA) were first applied to the global cluster structure optimization problem in a simple form [6]. After several algorithmic improvements, most notably the phenotype approach [7] as well as directed mutations and structural niches [8], it could be shown that the resulting evolutionary algorithm (EA) scaled only cubically in the standard benchmark application of homogeneous Lennard-Jones (LJ) clusters [8]. For several years, basin-hopping [START_REF] Li | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Wales | [END_REF][11] and this phenotype EA were the only unbiased algorithms able to find all known LJ n minima up to n ≈ 150. More recent developments [12][13][14] have pushed this limit further up to n = 561, by introducing dynamic constructions of intermediate grids for particle placement as a further key ingredient and by smaller variations of the basic techniques. By now the Cambridge Cluster Database [15] lists LJ structures up to n = 1000, but inclusion of a-priori information on the preferred grids (icosahedral and decahedral) is still needed for the upper half of this size range [16]. Also, the ultimate question at which cluster size the inevitable transition from icosahedral or decahedral structures to bulk-like grids (fcc or hcp) occurs is still open, although explicit investigations of LJ cluster phase change behavior has progressed to impressive sizes [17].

In the wake of this algorithmic development, both EA and basin-hopping methods have been applied to many kinds of homogeneous atomic clusters throughout the periodic system. The situation is different for heterogeneous clusters. There has been some limited work on binary LJ clusters [18][19][20] and binary Morse clusters [21] recently. Employing the Gupta potential, the Johnston group has studied several bimetallic cluster systems [22], including cluster sizes that are special for the LJ system, e.g., n = 38 [23] and n = 98 [24]. Literature on ternary clusters is extremely limited; in fact, we could locate only one significant paper [25] which apparently is starting off this line of research. (To avoid misunderstandings, we would like to emphasize that there are also related but traditionally separate research areas, for example inorganic chemistry style mixed metal clusters with bulky outer ligands, as those featuring in Refs. [26,27]; such systems are not of central interest here since for their successful treatment a proper level of ab initio electronic structure theory is more important than the performance of a global structure optimization scheme.) Clearly, one additional challenge of heterogeneous clusters is provided by the additional degrees of freedom due to exchange of atom types within otherwise identical structures, and the tight interplay of this with varying structural preferences. Another possibly less obvious aspect is the necessity to make the transition from experimental global optimization codes tuned to deal with one or very few given test systems to a general production code with the ability to handle arbitrary atom types with arbitrary interaction models.

Of course, the LJ potential is not an exact potential for anything, including rare gases. Likewise, other empirical potentials always are only crude approximations. Therefore, frequent attempts have been made to perform global cluster structure optimization directly at some better level of theory (as arbitrarily selected examples, see Refs. [28,29]). Such attempts have to be viewed with caution, since this causes energy (and gradient) evaluations to become more expensive by several orders of magnitude, which seriously compromises the convergence properties of any global optimization scheme. It is widely recognized by now that the global search should be done with a cheap model potential (or at most with a cheap semiempirical method), with the results then being refined at a suitable higher level of theory. An automatic scheme has been proposed to generate a system-specifically adjusted model potential on the fly [30], but this scheme has seen only a few applications [31,32] so far. In any case, this predicament indicates a pressing need for any global optimization strategy to be linked up easily with a broad spectrum of energy/gradient calculation methods.

Already the second application of GA methods to global cluster structure optimization was on benzene clusters [33]; followed by a few other studies of this system [34][35][36]. However, the vast majority of molecular cluster optimizations performed by EA methods were concerned with water clusters, in most cases with homogeneous neutral ones [4] (see the comprehensive overview in Ref. [37]), in rarer cases with water clusters containing a single different entity, e.g., an atomic ion [38][39][40][41]. Applications to clusters of other molecules, in particular to heterogeneous clusters of more complicated, flexible molecules, seem to be lacking entirely, with one ex- ception: prediction of crystal structures, either ab initio (see the review [42] and references cited therein) or using powder diffraction data [43,44]. In these cases, molecules more complicated than water or benzene have been studied, but in most cases with only very few dihedral angles (if any) that were actually allowed to vary. In addition, it can be argued that crystal structure prediction is somewhat simpler than cluster structure prediction, due to the necessary limitation to the known 230 space groups which allows to treat part of the problem by exhaustive enumeration. Nevertheless, again, the employment of suitably taylored model potentials is seen as a major asset [45].

From this historical sketch, the remaining challenges for global cluster structure optimization can be identified as: (i) heterogeneity with ≥ 3 species, (ii) flexible molecules with many internal degrees of freedom, and (iii) easy access to a broad variety of methods to calculate interparticle forces, from a variety of cheap model potentials to more expensive semiempirical and ab initio methods, if possible also with an option to tune a given model potential to the present system under study.

We have developed a new EA program suite OGOLEM that fulfills all of these requirements. It uses the experience gained from our previous, specialized implementations [8,39] but is written from scratch and extended to cover points (i)-(iii). Design and performance details are described in section 2. Section 3 provides application examples that highlight features (i)-(iii), namely heterogeneous atomic LJ clusters of more than three species and heterogeneous clusters of larger, more complicated molecules. Coordinate files for all clusters shown there can be obtained from the authors upon request.

Program Design

Nowadays, the most important feature of any program is scalability. This, in the advent of multicore processors, includes both SMP as well as MPP parallelization. Due to their inheret parallelizability, genetic algorithms are a well-suiting algorithmic choice.

To fulfill this kind of multiscalability, the OGOLEM framwork provides both a MPI frontend for MPP parallelization and a threading frontend for SMP parallelization omitting unecessary MPI overhead. In both cases, a pool algorithm was used to eliminate serial bottlenecks [46]. As Fig. 1 shows, the OGOLEM framework therefore provides excellent scaling both on SMP systems as well as in an MPP context. Our pool algorithm in combination with Java TM s build-in threading mechanisms provides a straightforward SMP parallelization. Therefore, the OGOLEM framework is well-equipped for the even further growing importance of multicore processing. Directly after scalability comes usability. This can be considered both from the developers point of view, wanting a framework to allow for rapid development cycles and containing clear yet flexible restraints, and from the users perspective, needing easy input and flexibility in algorithmic and methodological choice. The OGOLEM framework satsifies the developers needs through a rigorous objectoriented programming style and well-defined but flexible interfaces. The user is provided with an intuitive and easy input format and can choose between a variety of methods and algorithms for her problem. These algorithms include

• a straightforward genotype operator, • a phenotype operator, • a novel packing algorithm, • combinations of the above.

The genotype-Operator provides an unbiased 1-point-crossover, as described in Ref. [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF] but acting on real-number strings. The phenotype operator cuts a randomly rotated cluster using either a plane through the center of mass or, optionally, a plane randomly displaced from that one (similar to the operator described in Ref. [8]).

The OGOLEM framework also provides a pairwise collision detection engine and, as a novel development, a dissociation detection engine. The latter employs the efficient Warshall algorithm [START_REF] Warshall | [END_REF] to build a reachability matrix. These two ingredients are of crucial importance in the whole global optimization process, in particular for clusters containing larger molecules of non-trivial shape: Both collision detection and dissociation detection are activated before local optimization, discarding structures containing clashes or being dissociated. Since local optimization is responsible for ≫ 90% of the overall computer time [8], this erases much of the computational effort otherwise wasted into local optimizations that typically show both a very slow convergence behavior and non-competitive final energies.

The development of our novel packing algorithm was inspired by the design of new scheduling algorithms in the field of computational siences targeted at solving the NP-hard problem of finding an optimal dynamic packing of tasks. The main idea there is to order the tasks to be executed by their relative size (where the size definition is crucial and might contain priorities). The biggest tasks are then packed first and the smallest ones last, into the holes left by the bigger ones. This is similar to the classical knapsack problem [49], and it can be expected that optimal packing is a decisive ingredient also in global cluster structure optimization, quite independent of the interparticle potential(s). Therefore, we make use of the same idea and order the building blocks of the cluster by their size and pack them as desribed above. Performance benefits arise in particular in systems of very different building block sizes, as demonstrated in section 3.2. As a side effect, the combination of packing with classical genetic algorithms helps to avoid premature convergence of the latter, since new structures with a fresh genome are, with a certain probability, build from scratch, potentially bringing new genes to the pool.

Available methods include interfaces to • Orca and • Tinker and additional internal force fields. This includes force fields for homoatomic and heteroatomic Lennard-Jones clusters and for ionic clusters with the additional ability of system-specific global reparametrization. With some standard force fields, optimization of internal degrees of freedom could conceivably lead to force field definition gaps (bond dissociation, change of atom type) or to clashes. The former problem can be avoided by user input, which necessarily has to include bond definitions and lower/upper bond length limits; for the latter problem, re-usage of the collision detection engine is the obvious solution.

• AMBER, • DFTB+, • MNDO, • MOLPRO, • MOPAC, • NAMD,
The choice of method is made simplistic yet transparent for the user by specifiying program, method and, if applicable, further information like basis set or solvation model in a single input line. To give some examples, this starts with choosing AMBER as a backend, of course requiring an AMBER-style force field being available for this system, which can be used for bigger systems of biological relevance as shown later in this paper. Then there is, for smaller systems where a better description is needed, the possibility to choose MNDO's semiempirical methods. With MNDO's ability of using COSMO-style solvation for the system, this extends the possible input by an easy addition of the cosmo keyword to enable automatic water solvation. For even finer needs, e.g., the Orca interface provides combinations of quantum chemistry methods, like different DFT-functionals as well as ab-initio methods, and basis sets.

Through this selection of different program packages and methods therein, the user is likely to find a method for her system that provides the best compromise between accuracy and computational expense. Additionally, the user can (but does not need to) control every part of the global optimization in detail.

Completed are these features with platform-independence. Employing a pure Java TM approach, we can ensure platform-independence in principle, ranging from embedded systems to supercomputers, regardless of architecture and operating system. The only limitations arise from (i) the ratio between computational power and problem size concerning the architecture, e.g., MP2/aug-cc-pVTZ optimizations on embedded systems and (ii) the availability of, often binary, program packages on the specific platform, e.g., Microsoft Windows TM or the FreeBSD TM operating system.

We distribute the OGOLEM framework under a classical four-clause BSD-license to provide a high degree of freedom.

Applications

Mixed Lennard-Jones Clusters

The pairwise Lennard-Jones potential

v ij = 4ǫ ij σ ij r ij 12 - σ ij r ij 6 (1)
provides a tolerably good description for clusters of homoatomic rare-gas atoms and can easily be extended to heteroatomic clusters using the Lorentz-Berthelot

ǫ ij = √ ǫ i ǫ j (2)
σ ij = σ i + σ j 2 (3)
to obtain the mixed parameters ǫ ij and σ ij from their homoatomic equivalents. These rules have proven their applicatibility, e.g., in review [50].

As mentioned in section 1, one of the major problems with combinations of various atomic types is the exchange of atom types. To target this, we use an explicit XChange-Operator, providing a significant speedup in finding the global minimum over non-explicit exchange algorithms. For all computations in this section, we combine the XChange-Operator with an unbiased genotype-algorithm.

For homoatomic Lennard-Jones (LJ) clusters, there are, with few exceptions, well-defined and recurring structural units. One interesting question therefore is whether changes in composition will yield structural transitions at a given cluster size or not. As an example we looked into clusters composed of 19 LJ atoms. The homotamic LJ cluster has a double icosahedral shape in the global minimum, as depicted in Fig. 2. In the following three heteroatomic cases, we keep the total Figure 2. An Ar 19 cluster having a double icosahedral shape. This and all following cluster structure figures were generated using Jmol [START_REF]Jmol: an open-source Java viewer for chemical structures in 3D[END_REF] and POV-Ray [START_REF]POV-Ray -The Persistence of Vision Raytracer www[END_REF]. number of 19 atoms but start to modify the cluster composition. In all calculations, we rank our pool by energy, meaning that "rank 0" is our best guess for the global minimum, "rank 1" our best guess for the next-higher local minimum, etc. To all our experience, the low-energy end of such a ranked pool provides a very reliable impression of the true energy spectrum. Of course, the high-energy end will not be representative at all, due to the small size of the pool compared to the vast search space. Thus, looking at ranks 0-10, we are very confident that we are not missing any important structures within this energy range.

Our introductory example for mixed LJ clusters is a simple substitution of two argon atoms with xenon atoms, i.e., Ar 17 Xe 2 . This exchange modifies the composition only slightly, and the difference in pair potential characteristics (well depth and minimum-energy distance) between Ar and Xe is not large. Since there is a large energy difference between the double icosahedron and other cluster shapes in the homoatomic case, it is not surprising that this type of exchange does not have a big impact on the shape of the cluster: As presented in Fig. 3, the global minimum of Ar 17 Xe 2 still is a double icosahedron, only slightly distorted with respect to the homoatomic case, and retaining much of its symmetry. Obviously, instead of substituting the two interior argon atoms (the centers of the two icosahedra), the two xenon atoms are on the cluster surface. This can be rationalized by remembering the fact that if an icosahedron is built from perfect tetrahedra, small gaps remain on the outer surface. In other words, the larger xenon atoms are more easily accomodated on the surface than in the interior. As in nanoalloy clusters [22], another important issue is atom type mixing vs. segregation. Here, interestingly, the two xenon atoms stay next to each other, probably because this allows for the strain induced by the size difference to be compensated by bending the long axis of the double icosahedron which only has a negligible impact on most of the nearest-neighbor pair distances. This is confirmed by comparing to the lowest local minimum (rank 1 in Fig. 3): There, the two xenon atoms are separated but otherwise in very similar positions and surroundings. For the rank 0 structure, the long axis is bent by 5.6 • ; for the rank 1 structure, this value is 3.2 • , implying that a larger portion of the size-difference strain has to be compensated by other distortions, leading to a greater percentage of pair distances being slightly away from their optimal values, which in turn leads to a worse overall energy. It has to be pointed out, however, that these effects are very small, since the total energy difference is only 0.5 kJ/mol, or half a percent of the overall binding energy. This is a good indicator for the performance of the OGOLEM suite, and also a reminder of the importance of atom type exchange operators, as already noticed in the literature [22].

Of course, we can expect the energy landscape to get more complicated with every inclusion of more atom types. Considering the raw amount of possible compositions and their combinatorial explosion, already for ternary clusters any attempt at a complete treatment would be beyond the scope of this article. For this reason, we restrict the discussion to a few examplary results, trying to give explanations why certain cluster structures are superior to others and which computational challenges are important in these systems.

Looking into the quaternary LJ cluster Ar 2 Kr 5 Xe 5 Ne 7 , we can spot three different types of structures, within the best ten structures of our genetic pool and less than 0.5 kJ/mol apart. These three structures are depicted in Fig. 4; all of them show essentially the same distorted double icosahedral structure. The atom type distribution in the global minimum structure appears to be "messy" at first sight; it can, however, be rationalized a posteriori: One of the two icosahedra in the double icosahedral structure has the composition Ar 2 Kr 5 Xe 5 Ne 1 , with the biggest Xe atoms on one end, a middle layer of middle-sized Kr atoms, and the smallest Ar and Ne atoms at the other end. In effect, this brings all but one Xe atom into the central plane of the double icosahedron; there, they can participate in filling in the tetrahedral-packing gaps in both icosahedra. Then, the remaining six Ne atoms have two extreme options for forming the second icosahedron: They could either arrange around the single Ne atom end of the first icosahedron, or around the Xe atom end. The first option has the advantage that the size discrepancies at the border between the two icosahedra are minimized; the second option allows for energetically more favorable Ne-Xe contacts. It turns out that this second option wins, but at the expense of the system being unable to close this Ne icosahedron. Instead, the Ne atom that should sit at the far end of the long axis moves into a second geometrical layer, capping a Xe-Xe-Kr triangular face (i.e., in effect, five Ne-Ne interactions at distances much longer than the optimal one are replaced by one Ne-Xe and one Ne-Kr interaction at optimal distances).

The other two structures shown in Fig. 4 turn out to be very small variations of this rank 0 structure just described. From the rank 0 to the rank 7 structure, one additional Ne atom "jumps" from the Ne icosahedron into the second layer; no other changes occur, and the total energy decreases only by 0.13 kJ/mol (0.16%). Despite of these minor differences, disturbingly, six other structures are ranked inbetween these two, sometimes even with larger structural changes. As an example for these cases, we show the rank 2 structure. It differs from the rank 0 one in two respects: (i) the second-layer Ne atom moves to a different triangular face (but of the same composition, Xe-Xe-Kr, so the effect of this hop is very small), and (ii) one Ar-Kr pair switches its identity, leading to segregation of the two formerly neighboring Ar atoms, and to several changes in the numbers of nearestneighbor contacts per atom type pairs. Nevertheless, the overall energy change is even smaller than for the single-atom hop from rank 0 to rank 7. As a side note, the overall energy favours the most symmetric structure out of these three closelying possibilities. This discussion further documents the considerably increased difficulty of the mixed LJ problem and the utmost importance of unbiased but efficient search tools.

As an example representing the class of quinary LJ clusters we picked the Ar 5 Kr 5 Xe 5 Ne 2 He 2 cluster. As before, Fig. 5 displays three typical structures from the top ten in our pool. In comparison to the quaternary case, some patterns start to emerge. Again, all of these three quinary structures are almost identical, with only very few atom hops or atom type exchanges, as discussed below. Furthermore, in spite of the fairly different atomic composition, there are many similarities between these quinary structures and the quaternary ones discussed above. First of all, the general cluster framework is the same: It is again a double icosahedron, with one long-axis cap missing due to an atom having hopped into a second layer. Also, in both composition cases, the largest atom type (Xe) is arranged preferentially on the middle equator, joining the two icosahedra. The remaining atom types are distributed such that one icosahedron preferentially contains larger atoms, the other one smaller atoms (despite the fact that more symmetrical distributions could be realized, in particular in this quinary case). Finally, the atom hopping into the second layer is of the smallest type (here a He atom, above a Ne atom) and sits on a Xe-Xe-Kr triangular face in all cases. Also the differences between the quinary structures shown in Fig. 5 are qualitatively very similar to those discussed for the quaternary case. From the rank 0 structure to the rank 2 structure, an argon gets exchanged with a krypton, splitting the three adjacent argon atoms in the larger-atom icosahedron into a 2+1 segretation. This is accompanied by the outer He atom hopping to a different Xe-Xe-Kr triangular face. From the rank 2 to the rank 3 structure, the geometrical change is even smaller: It merely consists of turning one of the icosahedra relative to the other one around the long axis by 72 • (one fifth of a full turn around this five-fold axis; i.e., in the homoatomic case, these two clusters would be identical). For this particular atom type distribution, this rotation operation does not change the count of nearest-neighbor interactions for the different atom-type pairs. Hence, the overall energy diffence arises only from non-nearest neighbor interactions. This is in line with the fact that this energy change is very small (0.0015 kJ/mol or 0.0018%). Again, disturbingly, the seemingly larger change upon going from the rank 0 to the rank 2 structure has an even smaller overall effect on the energy.

Nevertheless, as Doye, Miller and Wales [START_REF] Doye | [END_REF] have shown, for 19 LJ atoms (in the homogeneous case) there is a particularly large energy difference between the (double) icosahedral global minimum and the best clusters of the other structural prototypes (decahedral and fcc). Therefore, as already mentioned above, it is not surprising that we have not found a fundamental structural change (away from the double icosahedron) even in the quaternary and quinary cases shown above. At other total atom counts, however, the difference between the structural prototypes is minimal, in particular at those famous sizes where the global minimum structure of small LJ clusters is not of the dominating icosahedral type anymore. The first such size is 38 atoms [START_REF] Doye | [END_REF]. Therefore, for LJ 38 , we expect to see major structural changes already for comparatively small changes in cluster composition. As a reminder, for the homogeneous case, the global minimum structure is a fcc packed one (see Fig. 6), despite almost all of the search space being dominated by icosahedral structures. For this reason, the fcc global minimum of LJ 38 is nontrivial to find already in the homogeneous case. Therefore, we cannot yet present final results in this first report. Nevertheless, our preliminary results do show an obvious tendency towards icosahedral structures upon going to mixed LJ clusters with a total of 38 atoms. As a typical example, we show our best candidate for the global minimum of Ar 30 Xe 4 Kr 4 in Fig. 7. As to be expected from the fairly intricate discussion for the various LJ 19 cases above, this structure is difficult to analyze. Nevertheless, the presence of (approximately) five-fold axes is strikingly obvious, ruling out the possibility that this could be a cut-out from an fcc crystal.

Actually, at least one pattern recognized above for LJ 19 of various compositions can be found here again: The largest atoms (Xe) preferentially sit at the boundary between (complete or partial) icosahedral motifs. Further discussion is referred to future publications of this ongoing work.

Kanamycin A Aggregation

In recent experimental studies in the institute for dermatology at the university of Kiel [54], larger aggregates of compounds related to Kanamycin A (KA) with certain physiological cations have been implicated as pathogen-associated agents, leading to a non-inflammatory innate immune response of the human skin [55]. Apparently, already for the parent compound KA, with the exception of a few isolated hints [56][57][58], nothing is known about molecular aggregates, neither experimentally nor theoretically. Therefore, we are currently investigating clusters of this compound class. Full details will be pushlished elsewhere [59]. Here we show first results on small clusters of the parent compound KA, to demonstrate the capabilites of the OGOLEM suite in dealing with more complicated species containing internal degrees of freedom with non-trivial influence on the outer molecular shape and properties. One of the main questions with this system is the role of physiological ions in the aggregation process. As a prelude to future studies, we therefore look into the differences arising from different KA:cation ratios in small clusters. All calculations shown here were carried out with NAMD [60] and AMBERs GAFF [START_REF] Case | AMBER 8[END_REF][START_REF] Wang | [END_REF] force fields as a backend, and with the packing operator in 1:1 mixture with a classical genotype algorithm. As a starting point, a pre-optimized KA geometry using Orca's B3LYP/TZVP [START_REF]Orca: an ab initio, DFT and semiempirical electronic structure package[END_REF][START_REF] Schaefer | [END_REF] was taken (see Fig. 8). Qualitatively, there are no structural differences between this DFT structure and the GAFF-optimized ones occuring during our global optimization runs. This gives some indication that this force field description is not totally inadequate. We do not claim, however, that this is the optimal or even a reliable description. A thorough inquiry into this question will be presented elsewhere [59]. For our present purposes, a possibly reasonable description is fully adequate. All claims made in this section should be taken as valid only within the particular force field description used here.

As an introductory example, we show the aggregate of two KAs with a single sodium cation. This is not only the smallest cluster possible of the described kind but there are also experimental indications supporting the hypothesis that this may be an important building block for clusters of bigger size. As can be seen in Fig 9, the two KAs encapsulate the sodium cation in a manner similar to two chelate ligands, coordinating it with eight oxygen atoms within 3.2 Å distance. In the KA-Na-KA closest contact area, not all Na-O distances are equal, instead they range from 2.4 Å to 3.11 Å. Since the next nearest possible coordination site is much farther away (4.2 Å), it does seem to be sensible to classify this is an eightfold coordination, in the form of a distorted tetragonal prism. While this is usually perceived as over-coordination for the comparatively small Na cation in aqueous biochemistry [65], it has been observed in zeolites [66] and inorganic compounds [67]. Of course, as indicated above, it could also be a force field artifact. This will be checked in ongoing studies. Looking into larger clusters, we would like to present KA tetramers both with a single sodium cation and with two sodium cations, showing differences arising from the addition of a single sodium cation to the system. As can be seen from Fig. 10 and the close-up provided there, the rank 0 structure for KA 4 Na + again coordinates the sodium cation with eight oxygens within 3.2 Å distance and with the same structural pattern. This similarity with the KA dimer could be a first support for the suspicion that KA 2 Na + could be a building block in larger KAcation clusters. The KA 4 Na + system turns out to be very stable: In all our global optimization runs, not a single dissociated structure was found. This changes drastically upon addition of a second sodium cation, generating KA 4 Na 2+ 2 . With two sodium cations (and neither counterions nor explicit or implicit solvent), the KA tetramer has a rather high tendency to dissociate. In fact, without employing our new packing algorithm as additional operator, it is very hard to find any non-dissociated structures at all. Even with packing, enabled the fraction of dissociated structures in a typical genetic pool after global optimization typically is at or above 50%, as shown in Fig. 11. Nevertheless, with packing as essential ingredient, we can effi- ciently find structures for KA 4 Na 2+ 2 that are not only fully associated but also very compact. One of them is our best candidate for the global minimum structure of this system (Fig. 12a). Even in this global minimum candidate structure, there are some indications that four KA molecules have difficulties in holding two positively charged sodium cations at a distance of only 9.6 Å. If we assume that the eightfold coordination described for KA 2 Na + above is ideal, this coordination should also be possible to attain in KA 4 Na 2+ 2 , since this cluster simply is (KA 2 Na + ) 2 . Clearly, however, the coordination surroundings of both sodium cations are different in our rank 0 structure: Both are only sixfold cordinated. Apparently, the less-than-optimal sodium coordination is partially compensated by favorable KA-KA non-bonding contacts. As support for this claim we also show another type of fully associated structure in Fig. 12b. There, KA 4 formation is obviously very loose, and at the same time the sodium coordination clearly is deficient. As a result, the total energy is even worse than for many dissociated structures (shown in the following paragraph). Apparently, there are two basic modes for KA 4 Na 2+ 2 to dissociate. One is loss of an isolated sodium cation, leaving a stable KA 4 Na + behind. The other one is dissociation into two KA 2 Na + building blocks. Typical candidates for both classes are shown in Fig. 13. It is not surprising but typical that the KA 4 Na + + Na + case has a better energy than the KA 2 Na + + KA 2 Na + case. Association of KA with itself happens via multiple hydrogen bonds, which overcompensates the loss of 5-8 Na-0 coordinative contacts. In fact, KA 4 Na + + Na + occur frequently in our pools, while KA 2 Na + + KA 2 Na + cases are hard to find. Preliminary results for KA 4 Na 3+ 3 indicate that despite our packing algorithm we are unable to find non-dissociated clusters for this system. An obvious but still speculative rationalization for this finding could be that four KA molecules in compact formation simply do not provide enough room and a sufficiently high number of coordination contacts to keep three sodium cations within a pair distance of 9 Å, which is the minimum sodium atom pair distance observed in KA

Summary and outlook

After identifying particle type heterogeneity, larger flexible molecules, and variability on calculating interparticle forces as the remaining challenges in global cluster structure optimization, we have presented our OGOLEM program that is designed to meet all these challenges, resulting in a general and user-friendly framework for diverse practical applications. As examples demonstrating these abilities we have first proceeded beyond the standard benchmark of homogeneous LJ clusters and the beginning research on binary systems, by making first steps into the area of ternary, quaternary and quinary LJ clusters. Depending on cluster size, we find both cases were the qualitative cluster structure is (almost) unchanged compared to the homogeneous case and others were even small changes in atom type composition lead to drastic structural transitions. Many of the observed effects occur within very narrow energy windows, re-emphasizing the need for particle type exchange operators, incorporated into OGOLEM.

As second demonstration system, we have examined small clusters of KA with sodium cations, as an example for heterogeneous clusters containing particles of very different kind, including a large, flexible species. Here it turned out that another new ingredient, a packing algorithm, is vital for an efficient treatment.

With these exemplary applications, we have demonstrated the ability of OGOLEM to treat complicated systems in global cluster structure optimization. This opens up a plethora of possible future applications. As already indicated in section 3.2, the KA aggregation results shown here only are the very first steps into a complex project. Further ongoing applications include Sb/Ge vanadates featuring curious arrangements of vanadate polyhedra, with possible future use as molecular magnets, and the influence of small structural variations on micelle building propensities in glycosylated lipids. Work on these and other systems is in progress in our lab.

  Global Cluster Structure Optimization for Arbitrary Mixtures of Flexible Molecules. A Multiscaling, Object-Oriented Approach.

Figure 1 .

 1 Figure 1. Scaling shown as relative speedup. Ideal linear scaling shown as line. Deviations arise from the heavy use of random numbers.

Figure 3 .

 3 Figure 3. Minimum structures of the Ar 17 Xe 2 cluster having a slightly distorted double icosahedral shape.

Figure 4 .

 4 Figure 4. Geometries of the Ar 2 Kr 5 Xe 5 Ne 7 cluster, each displayed in two different views (top, bottom).

Figure 5 .

 5 Figure 5. Geometries of the Ar 5 Kr 5 Xe 5 Ne 2 He 2 LJ cluster.

Figure 6 .

 6 Figure 6. The global minimum of the Ar 38 cluster, a fcc packing.

Figure 7 .

 7 Figure 7. The global minimum of the Ar 30 Xe 4 Kr 4 cluster, clearly not fcc.

Figure 8 .

 8 Figure 8. The starting point, a B3LYP/TZVP optimized monomer structure.

Figure 9 .

 9 Figure 9. The aggregate KA 2 Na + and a close-up of the coordinative region, indicating its tetragonal-prismatic form.

Figure 10 .

 10 Figure 10. The aggregate KA 4 Na + (Energy: 237.305 kJ/mol) and a close-up of the coordinative region.

Figure 11 .

 11 Figure 11. Dissociation ratios as a function of the binned structures of the genetic pool.

( a )Figure 12 .

 a12 Figure 12. Fully associated structures of the KA 4 Na 2+ 2 cluster.

( a )Figure 13 .

 a13 Figure 13. Dissociated structures of the KA 4 Na 2+ 2 system.

4

  Na 2+ 2 .

	October 20, 2009 20:13 Molecular Physics ogolem˙first
	14	REFERENCES
	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

Acknowledgements

On this occasion, both authors would like to thank Prof. Dr. Hans-Joachim Werner for his generous support. Both authors are also grateful to Prof. Dr. Jens-Michael Schröder for his enthusiastic initiative in starting a collaboration on KA aggregates. A large computer time grant at the HLRN supercomputer center (Berlin/Hannover) is gratefully acknowledged.