Alexey ; Alekseyev

Vidisha Rai-Constapel
email: vidishal@yahoo.com

Heinz-Peter Liebermann

Aleksey B Alekseyev

Robert J Buenker

Configuration interaction study including the effects of spin-orbit coupling for the electronic states of the LiX molecules (X=C, Si, Ge, Sn)

Keywords: MRD-CI, Spin-Orbit Interaction, Transition moment, Radiative lifetime MRD-CI, spin-orbit interaction, transition moment, radiative lifetime

Ab initio multireference configuration interaction calculations including spin-orbit coupling effects have been carried out for four LiX molecules (X=C, Si, Ge and Sn). Potential energy curves of the ground and low-lying excited states have been obtained in each case as well as the corresponding spectroscopic constants. Transition moments have also been computed in order to give estimates of the radiative lifetimes of the excited states for each system. Trends in a variety of quantities such as Te values, spin-orbit splittings, equilibrium bond lengths and vibrational frequencies for this series of molecules are discussed in detail and comparison with the corresponding data reported earlier for the PbLi system is also made.

Introduction

The complementary importance of theoretical chemistry to experimental work has now long been accepted. As a very recent example, the interpretation of the complex experimental spectra of PbLi received major help from the ab initio relativistic configuration interaction (CI) calculations [1]. PbLi has a number of bound states in the low-energy excitation region (within 14000 cm -1 above the ground state). Further complexity is a natural consequence of the strong spin-orbit (SO) coupling characteristic for the PbLi system. Ab initio calculations of spectroscopic constants (T e , ω e and r e ) and transition moments significantly simplified identification of the emission spectra of the hitherto unknown free radical PbLi [1].

The lighter diatomics of Li and Group IVA elements C-Sn are equally important in their chemical interest. The lightest molecule of this series LiC is known for its importance in radiocarbon dating [2] and as anode in Lithium-ion batteries [START_REF] Nazri | Lithium Batteries Science and Technology[END_REF]. Meanwhile LiSi is of special interest for the Li-ion batteries due to the large theoretical specific charge capacity of silicon [START_REF] Kubota | [END_REF]. LiSi and LiGe compounds are the so-called Zintl phases [5], which are known to have good semiconductor properties, whereas LiSn and PbLi belong to the class of intermetallic phases [5] possessing ceramic-like properties. In all these cases it is of interest to accurately determine the properties (e.g., bond lengths, vibrational frequencies etc.) of the single entity Li-X. Though DFT methods are being extensively applied to such cases (see, e.g., Refs. [START_REF] Kubota | [END_REF]6]), it is well known that not always they provide accurate and reliable data. Therefore, it is important to compare the DFT results with theoretical data obtained by more accurate ab initio approaches, like multireference CI used in the present study.

The low-energy electronic spectra of diatomic molecules containing the lithium atom are generally marked by a much higher density of states than for the corresponding hydrides. The reason for this distinction lies primarily in the well-known fact that hydrogen can form strong σ bonds with other atoms, whereas Li cannot. As a result, the energy gap between the highest-occupied σ molecular orbital (MO) and the lowest π is much smaller in lithium compounds, greatly increasing the number of their low-energy electronic configurations in comparison with the corresponding hydrides. Indeed, the latter relationship already can be seen from the fact that the ground states of the CLi and PbLi systems are not of 2 Π symmetry, with a σ 2 π configuration well known for CH and PbH, but rather prefer a σπ 2 occupation that causes the most stable state to be 4 Σ -in both Li compounds.

The SO interaction adds to the complexity in the spectra of the lithiumcontaining systems. This aspect is of little significance for LiC, but the effects of SO coupling are expected to grow steadily as heavier atoms replace carbon. For example, there are 19 bound electronic states that are found below 3.0 eV in the electronic spectrum of PbLi. Moreover, the composition of the resulting Ω states changes dramatically with the bond distance. The SO effects also have a strong influence on the radiative stability of the excited states of such systems. The purpose of the present study is to investigate trends in the spectra of the series of Li-Group IVA diatomic molecules by employing multireference single-and double-excitation configuration interaction (MRD-CI) calculations in conjunction with relativistic effective core potentials (RECP) to obtain potential energy curves and other properties of their ground and excited states. New calculations are reported for LiC, LiSi, LiGe and LiSn and comparison is also made with the results of earlier theoretical work dealing with PbLi.

Technical Details of the Theoretical Treatment

A common (12s7p4d)/[8s5p4d] basis set for Li [7] has been employed in the present study. The lithium 1s shell has been described using the RECP from Ref. [8]. The (12s7p4d1f )/[6s5p4d1f ] basis set for carbon has been taken from Ref. [7]. A RECP of the same type [8] has been employed to describe the 1s shell of the carbon atom. For Si a basis set of the (15s9p2d1f )/[5s4p2d1f ] type has been used which has been augmented with the diffuse s(0.033), p(0.0237), d(0.0556) and f (0.125) functions [9] along with a RECP from Ref. [8] to describe the 1s, 2s and 2p subshells. A semicore RECP [10] to describe the argon core and the augmented [6s5p3d2f ] basis set [11] have been used for germanium. The [7s6p4d2f ] basis set and the semicore relativistic core potential from Ref. [12] have been used to describe the tin atom. The original Sn basis set of Ref. [12] has been augmented by one s(0.01), p(0.01), d(0.08) and f (0.2) functions.

Self-consistent field (SCF) calculations has been performed for the neutral molecules in their ground 4 Σ -(σ 2 π 2 σ * ) states. The SCF calculations provide the symmetry-adapted molecular orbitals for constructing the reference spaces to be used in the multireference configuration interaction calculations [13,14]. Details of the CI calculations without the SO interaction are given for LiGe in [START_REF] Davidson | The World of Quantum Chemistry[END_REF][START_REF] Knowles | [END_REF], which takes into account the effects of higher than double excitations with respect to the reference configurations, has been applied to the extrapolated (T = 0) energies of each root to estimate the energies at the full CI level of treatment. There are no unbound states in the low-energy region of the spectra of these systems and hence no bound-free contribution is reported.

The variational wave functions are used to compute properties such as electric dipole moments and radiative lifetimes. These functions are also employed in the second step of the Λ -S contracted spin-orbit CI (LSC-SO-CI) procedure [17,18] to form a matrix representation of the full relativistic Hamiltonian including SO coupling. Finally, the calculated potential energy curves for the Ω states have been employed to solve the one-dimensional nuclear motion Schrödinger equations by means of the Numerov-Cooley method [19,20] in order to obtain the corresponding vibrational energies and wave functions. Further, dipole transition moments µ have been computed and used for calculation of the lifetimes of the various excited states according to the formula for the Einstein coefficient (s -1 ):

A = 2.026 × 10 -6 |µ| 2 ∆E 3 , (1)
where µ is given in units of ea 0 and ∆E is the corresponding transition energy in cm -1 . The lifetime τ is then obtained as

τ = i A i -1 , (2)
where the sum is over the lower-state vibrational levels.

Potential Energy Results Without Spin-Orbit Coupling

Potential energy curves for the 19 Λ-S states dissociating to the first three atomic limits of Li( 2 S) + X( 3 P ), Li( 2 S) + X( 1 D) and Li( 2 P o ) + X( 3 P ), with X = C, Si, Ge and Sn, are shown in Figs. 1234. In the actual calculations higher-lying states have also been included in order to take account of their effect on the more stable states. For purposes of clarity, however, only the lowest three atomic limits are shown in the figures. The ground state dissociation energy without including the SO interaction has been calculated to be energy calculation for the ground state of LiSn carried out within the framework of the DFT approach [6] has given 0.993 eV, which is 0.5 eV smaller than the present ab initio value. This indicates that the DFT method employed in Ref. [6] requires significant improvement of accuracy. No experimental or theoretical results are available for comparison in case of the LiGe radical.

The ground state symmetry is found to be 4 Σ -(σ 2 π 2 σ * ) in all cases. As mentioned in Ref. [21], the ion-pair ground state of LiC should correlate with Li + C -, but it converges to the lowest Li( 2 S) + C( 3 P ) neutral atom asymptote due to a strongly avoided crossing with the 2 4 Σ -state. The same is also true for the other molecules of this series. Three other states, namely 2 Σ -and 2,4 Π, correlate to the lowest dissociation limit. The 2 Π state is the first excited state resulting from the π → σ * electronic excitation. The 2 Σ -state is also bound and arises from the ground state electronic configuration, while 4 Π is repulsive in character and arises from the π → σ * * electronic transition. The potential curves of the states dissociating to this limit are rather similar (Figs. 1234) in all four cases, except for the fact that the equilibrium distances are progressively shifted towards larger values and the potential wells get shallower, as to be expected.

The states converging to the second atomic limit of Li( 2 S) + X( 1 D) are the 2 Σ + (σ 2 π 2 σ * ), 2 2 Π (σ 2 π 3 ) and 2 ∆ (σ 2 π 2 σ * ) states. The 2 2 Π state is bound in all four molecules, but the potential well becomes shallower as one moves downward from LiC to LiSn. The other two states corresponding to the second asymptote are also bound. The third atomic limit, which has also been included in the calculations, is Li( 2 P o ) + X( 3 P ). A total of 12 Λ -S states dissociate to this limit. Of these states, only 3 2 Π, 2 2 Σ + , 2 2 Σ -and 4 2 Π show strongly bound character.

In this case the second, third and fourth 2 Π states deserve a closer look. The avoided crossings of these three states can be clearly seen in the potential energy diagrams. An analysis of the leading electronic configurations of the 2 Π states at the Λ -S level makes the picture even clearer. In the case of LiSi, LiGe and LiSn, 2 2 Π (σ 2 π 3 ) and 3 2 Π (σ 2 πσ * σ * * ) exchange character at about 5.2 -5.4 a 0 . The dipole moment curves of the two states cross each other in this interval. In the case of LiC the 2 2 Π and 3 2 Π states have a minimal energy splitting at 4.70 a 0 . According to Ricca et al. [21] this crossing occurs at about 4.72 a 0 . At the larger internuclear distances of approximately 6.2 -6.4 a 0 , an avoided crossing between 3 2 Π and 4 2 Π is observed for LiSi, LiGe and LiSn. The corresponding electronic configurations [3 2 Π (σ 2 π 3 ) and 4 2 Π (σ 2 πσ * σ * * )] are exchanged around this point. The 2 2 Π state dissociates to the Li( 2 S) + X( 1 D) limit, whereas 3 2 Π and 4 2 Π converge to Li( 2 P o ) + X( 3 P ). For LiC the picture is somewhat more complicated. The third and fourth 2 Π states lie notably higher in energy than in the other three systems and are more heavily mixed with the higher-lying 2 Π states. Therefore the well-defined exchange in character observed in other cases does not occur in LiC, in which, starting at ≃ 6.00 a 0 , contributions from the other electronic configurations become important. Ricca et al. [21] calculated this avoided crossing in LiC to occur at 6.24 a 0 , while an estimate of 6.4 a 0 is obtained in the present study. respectively, for all states that correlate with the first two Λ-S atomic asymptotes. Altogether, there are nine states that correlate with the lowest Λ -S dissociation limit, five of Ω = 1/2 symmetry, three of 3/2 and one of 5/2 symmetry. The lowest seven of these Ω states are found to have bound potentials and correspond to the 4 Σ -, 2 Π, 2 ∆/ 4 Π and 2 Σ -states. The computed spectroscopic constants for these states in all four molecules are given in Tables 2 (LiC) and 3 (others). Five additional bound Ω states are found for the lightest of the four molecules, LiC, which correspond to the 2 Σ + , 2 2 Π and 4 Π states in the Franck-Condon region und dissociate to the Li( 2 S) + X( 1 D) limit.

The SO mixing of the various Λ-S states slowly increases as one goes along the series. The Λ-S composition of states for the heaviest of the four molecules studied, LiSn, is shown in Table 4 for a number of key bond distances. The LiSn Ω=1/2 ground state changes its character substantially on the way to the dissociation limit, becoming a nearly equal mixture of X 4 Σ -, a 2 Π and 1 4 Π at R = 9.0 a 0 . By contrast, the lowest 3/2 state continues to have a nearly pure composition over a large range of bond distances until the a 2 Π state starts to make a significant contribution beyond R = 5.5 a 0 . The splitting between these two states at their nearly equal equilibrium geometries is only 53 cm -1 , whereas it is much smaller for the three lighter systems (Tables 23). In the latter three cases, the lower of the two Ω states has 3/2 symmetry (inverted order). The opposite ordering found for LiSn is consistent with the results of the earlier SO-CI study on PbLi reported by Setzer et al. [1]. The value of the splitting is much larger in the latter case (2086.6 cm -1 obs. vs. 1600-1800 cm -1 calc.), so this is a good measure of the increasing importance of the SO effects along the series. The main reason for the small SO splitting for this state in the four lighter molecules is the fact that the effect in question is only of second order because it is a Σ state. One can give a fairly simple explanation for the different ordering of the 1/2 and 3/2 components of the X 4 Σ -ground state in the light and heavy members of the LiX series. The X 4 Σ -splitting is mainly caused by its SO interaction with the a 2 Π state. Due to the higher Wigner-Eckart factor such interaction is stronger for the 3/2 component and this leads to a stronger downward shift of 4 Σ - 3/2 in the light LiX (X=C, Si, Ge) molecules. One can also note that this effect is very weak, smaller than 1 cm -1 in LiC and LiSi, and only 4 cm -1 in LiGe. In the heavier LiSn and LiPb systems, however, the first-order splitting of the a 2 Π state is so large and this state lies so close to X 4 Σ -that the SO coupling of the lower 2 Π 1/2 component with 4 Σ - 1/2 begins to dominate and that makes 4 Σ - 1/2 the ground state. It is also worthwhile to note that at large internuclear separations the 4 Σ - 1/2 component becomes the lowest state in all LiX systems, since it is the only state which converges to the lowest Li( 2 S 1/2 ) + X( 3 P 0 ) dissociation limit. One can clearly see this for LiGe and LiSn in Figs. 5 and6, respectively.

The computed SO splitting for the a 2 Π state increases in a much more regular fashion due to its primarily first-order nature. It is 48 cm -1 in LiC (Table 2), 117 cm -1 in LiSi, 695 cm -1 in LiGe and 1669 cm -1 in LiSn (Table 3). The results in Table 4 for the latter system show that the corresponding two Ω states are composed predominantly of the a 2 Π state at their approximate equilibrium bond distance of 5.25 a 0 . The SO selection rules also keep the 2 ∆ splitting to relatively small values (only 99 cm -1 for LiSn, for example).

There is a general trend toward lower T e values for each state as the mass of the Group IVA atom increases. The greatest change is found for the a 1,2 2 Π 1/2,3/2 states, for which a decrease in excess of 6000 cm -1 is computed in going from carbon to silicon. Further decrease when Si is replaced by still heavier atoms is much more moderate (Table 3). The same is also true for the other low-lying states. For example, the T e value of the c 2 Σ -state is only 1055 cm -1 lower in LiSn than in LiSi. The changes are even smaller for the corresponding 2 ∆ states, but again, there are no exceptions to the general rule that these excitation energies decrease along the series. Similar remarks hold for the bond distances and vibrational frequencies for this group of systems. There is an increase in the R e value of 0.5 Å for the ground state of LiSi compared to LiC, and nearly the same shift is computed for each of the other low-lying states. The corresponding increment in replacing Si by Ge is much smaller, falling in the 0.03-0.05 Å range, but the changes at the next step are notably larger. For example, the R e value for the ground state of LiSn is 0.212 Å larger than for LiGe and it is 0.228 Å larger in the case of the c 2 Σ -state.

The situation for the vibrational frequencies is similar, with their values always falling along the series. The drop is 238 cm -1 at the first step for the ground state (Tables 23), and then only 44 and 50 cm -1 in going to LiGe and LiSn, respectively. It is obvious from the potential energy curves that the frequency for the a 2 Π state is always smaller than for X 4 Σ -. The corresponding values for the b 2 ∆ components always fall between these two pairs of frequencies, whereas the calculations indicate that the values for the c 2 Σ -state are always the smallest among this group of states for any of the LiX systems treated.

Transition Moments and Radiative Lifetimes

The SO-CI wave functions have been used to compute electronic transition moments among the various pairs of states. The results for selected transitions for LiSi, LiGe and LiSn are shown in Figs. 789. Since the values obtained for LiC are too small, as can also be inferred from the calculated lifetimes (Table 5), the corresponding figures are not presented here. These data are averaged over pairs of vibrational functions to obtain radiative lifetimes for the various upper states as described in detail in our previous work [1]. The computed lifetimes of the v ′ =0 levels for each of the five lowest Ω electronic states of all four LiX systems are given in Tables 5678, respectively. Both partial and total lifetimes are reported. They are obtained in the usual way by summing over the Einstein coefficients for transitions to all lower states and inverting [1]. Only transitions to the two components of X 4 Σ -are considered in obtaining the results shown in the lifetime tables. Results for partial lifetimes are given separately for parallel (∆Ω = 0) and perpendicular (∆Ω = 1) transitions to the X 1 and X 2 states (listed as τ 1(2) in the tables).

All of the above transitions are spin-forbidden. Their probabilities therefore depend strongly on the level of the SO mixing for the various Ω states involved (magnetic dipole and other mechanisms have not been considered but are expected to be of only minor significance). As a result, the lifetimes of the LiC states are quite long and are found to be of the order of a few seconds (Table 5). The corresponding values for the LiSi states are 20-50 times shorter, but still fall in the 0.1-0.4 s range (Table 6). The shortest lifetime among this group of the LiSi low-lying states is computed for c 2 Σ - 1/2 (41 ms). An even larger factor is noted in going to the LiGe system. Its lifetimes are computed (Table 7) to be 40-100 times shorter than those for LiSi. Again the shortest τ tot value is found for the c 2 Σ - 1/2 state, but this value is still only of the order of 1.0 ms and thus would be quite difficult to observe with currently available techniques [1]. The calculations show that the perpendicular c → X 1 transitions are primarily responsible for depopulation of the c state in all LiX molecules.

The lifetimes of the LiSn states (Table 8) are 5-15 times shorter than those of the corresponding LiGe states. Comparison with the lifetimes computed for PbLi in our earlier work [1] is difficult, because there is so much more SO-mixing in the latter system. For example, the lowest 4 Π state makes the dominant contribution to the third-lowest Ω=1/2 state in PbLi near the ground state equilibrium bond length, whereas in LiSn the corresponding state (c 2 Σ - 1/2 ) has virtually no 4 Π character until well beyond its R e value (see Table 4). The latter has the shortest lifetime among the low-lying LiSn states but it is still relatively long being 0.1 ms. The next smallest value occurs for the b 1 state and is about twice as large as for the c state.

Conclusion

The increasing importance of the SO interaction and other relativistic effects in the series of diatomic molecules including Li and a Group IVA atom has been studied by employing the multireference CI techniques in conjunction with effective core potentials. There is a significant decrease in the T e values of the excited states of LiSi by 4000-6000 cm -1 relative to LiC. An increase in the corresponding equilibrium bond distances of 0.5 Å for both ground and excited states is found, along with a substantial decrease in their vibrational frequencies. The trend towards decreased T e values and weaker bonds continues along the series, but the increments are considerably smaller in going from LiSi to LiGe and LiSn, respectively.

The electronic ground state of all four systems is composed primarily of the X 4 Σ - state and thus has only a small SO splitting due to the second-order character of this effect. The largest value is found for LiSn, for which the computed splitting is 53 cm -1 . The corresponding value for PbLi is much larger (2086.6 cm -1 ) by virtue of the large admixture of 2 Π character near the equilibrium geometry of the heavier molecule. By contrast, the splitting for the a 2 Π state itself is already 48 cm -1 for LiC, in which case the effect is of first order. This value grows steadily through the LiX series, reaching 1669 cm -1 for LiSn. The SO splitting for the b 2 ∆ state is only slightly larger than for X 4 Σ -.

Transition moments and radiative lifetimes have been computed for the low-lying electronic states of all four systems. As expected, the τ values become smaller as the mass of the Group IVA atom increases. The computed values for the lightest system are all greater than 1 s, whereas the shortest lifetime found for LiSn is 0.1 ms (for c 2 Σ - 1/2 ). Thus, there is a decrease of four orders of magnitude in these quantities over the entire series. In each case the c 2 Σ - 1/2 state has the shortest lifetime among the low-lying states. The effect is due primarily to the perpendicular transition to the X 1 state for the three lightest LiX systems, but in LiSn transitions to both X 1 and X 2 states are equally strong.

Acknowledgement

The authors wish to express their sincere gratitude to Prof. E. H. Fink for numerous discussions during the course of the present study. The financial support of the Deutsche Forshungsgemeinschaft is also hereby gratefully acknowledged.

F

State

R/a 0 X, 4 Σ -a 2 Π b 2 ∆ c 2 Σ -1 2 Σ + 1 4 Π X 1 4 Σ - 1/2

4 .

 4 Effects of the Spin-Orbit Interaction Diagonalization of the full spin-orbit Hamiltonian in the basis of the Λ-S states has been carried out to obtain potential energy curves for the corresponding Ω states of each of the four LiX molecules. Because of the small SO splittings in the light LiC and LiSi systems, their Ω potential energy curves are nearly indistinguishable from their Λ-S counterparts shown in Figs. 1-2 and therefore are not presented here. The corresponding Ω-state potentials for LiGe and LiSn are shown in Figs. 5 and 6,

  and I.T. Consultant

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Computed MRD-CI potential energy curves for the low-lying electronic states of LiC at the Λ -S level of treatment.

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

 4567 Figure 4. Computed MRD-CI potential energy curves for the low-lying electronic states of LiSn at the Λ -S level of treatment.

Figure 8 .Figure 9 .

 89 Figure 8. Computed dipole transition moments [(a) -parallel, (b) -perpendicular] between thea 1,2 2 Π 1/2,3/2 , b 1,2 2 ∆ 3/2,5/2 , c 2 Σ - 1/2 and X 1,2 4 Σ - 3/2,1/2 states of the LiGe molecule as functions of internuclear distance.

Table 1

 1 ground state equilibrium distance R e = 4.53 a 0 . The treatment has been similar for all molecules considered and hence individual tables for the other three systems are not presented. All calculations have been carried out in C 2v symmetry. A total of four roots have been calculated for 4 A 2 (five roots for LiSi), 2 A 1 and 4 B 1,2 symmetries. For the 4 A 1 , 2 A 2 and 2 B 1,2 symmetries the number of roots calculated are two, six and seven (eight in LiSi), respectively. The configuration selection threshold T is 0.02 × 10 -6 E h for LiSn and LiGe, 0.05 × 10 -6 E h for LiSi and twice this value for LiC. The 3d and 4d electrons of Ge and Sn, respectively, have been frozen as an inner core in the CI procedure. Thus, only five valence electrons have been active in all cases. The c 2 i values lie uniformly in the 92-94% range. The leading electronic configurations for each root are also listed in Table1. The multireferenceDavidson correction

	at the

Table 1 .

 1 Technical details of the MRD-CI calculations of LiGe. a N ref and Nroot refer to the number of reference configurations and calculated roots, respectively. The number of selected SAFs and the Σ c 2 p values over reference configurations are given for R = 4.53 a 0 .

	o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

a SAFTOT designates the total number of generated, SAFSEL the number of selected symmetry-adapted functions (SAFs),

Table 2 .

 2 Spectroscopic constants of LiC (transition energies Te, bond lengths Re and vibrational frequencies ωe) calculated in the present study (Calc.), and in References[21] (ICMRCI) and[22] (MP2). The results reported have been obtained without inclusion of SO coupling.

				Molecular Physics				Page 11 of 27
	State	T e /cm -1 Calc. ICMRCI a	R e / Å Calc. ICMRCI MP2 a	ω e /cm -1 Calc. ICMRCI MP2
	X 1 X 2 a 1 2 Π 1/2 4 Σ -3/2 4 Σ -1/2 a 2 2 Π 3/2 b 1 2 ∆ 5/2 b 2 2 ∆ 3/2 c 2 Σ -1/2 d 2 Σ + 1/2 e 1 2 Π 3/2 e 2 2 Π 1/2	0 0 10735 10783 12392 12409 13663 17347 19812 19826	0 10148 11339 12720 16563 19069	1.879 1.879 2.051 2.057 1.893 1.904 1.944 1.901 1.870 1.869	1.884 2.046 1.918 1.951 1.924 1.905	1.884 2.047 1.771	693 693 522 519 678 672 584 653 707 709	691 523 656 572 643 659	703 847

a URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 3 .

 3 Calculated (this work) spectroscopic properties of LiSi, LiGe and LiSn (transition energies Te, bond lengths Re and vibrational frequencies ωe). The ground state of LiSn is the 4 Σ - 1/2 state. b Boldyrev et al.[22] have calculated the Re and and ωe values for the X 4 Σ -and a 2 Π states of LiSi to be 2.355 Å, 472 cm -1 and 2.545 Å, 464 cm -1 , respectively.

	Page
	State	T e /cm -1 LiSi LiGe LiSn a	LiSi	R e / Å LiGe LiSn	ω e /cm -1 LiSi LiGe LiSn
	X 1 X 2 a 1 2 Π 1/2 4 Σ -3/2 4 Σ -1/2 a 2 2 Π 3/2 b 1 2 ∆ 3/2 b 2 2 ∆ 5/2 c 2 Σ -1/2	0 0.2 4430 3550 2214 0 53 4 0 4547 4245 3883 8403 8077 7620 8502 8148 7719 9299 8997 8244	2.384 b 2.412 2.617 2.384 2.412 2.624 2.535 b 2.585 2.762 2.534 2.578 2.753 2.401 2.442 2.637 2.401 2.445 2.644 2.484 2.531 2.759	455 b 411 455 411 402 b 364 402 365 429 391 435 390 362 322	367 361 347 343 353 348 279

a

12 of 27 URL: http://mc.manuscriptcentral.com/tandf/tmph

	February 2, 2011 11:0 Molecular Physics MS˙Mol˙Phys
	Molecular Physics
	Molecular Physics	11
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 4 .

 4 Composition of the seven lowest Ω states of LiSn (c 2 , %) at various bond distances R. Entries are only made for contributions with c 2 ≥ 1.0 %.

Table 5 .

 5 Contributions from the higher-lying roots are included. Partial and total radiative lifetimes (in s) for transitions from the v ′ = 0 level of the LiC excited states to the X1 4 Σ - 3/2 (τ1) and X2 4 Σ - 1/2 (τ2) states.

	4.30	96.4	1.8		1.4
	4.95 a 5.00 5.25 5.70 9.00 4.30 4.95 a 5.00 5.25 5.70 9.00 2 Π 1/2 4.30 X 2 4 Σ -3/2 a 1 4.95 5.00 5.25 a 5.70 9.00 a 2 2 Π 3/2 4.30 4.95 5.00 5.25 a 5.70 9.00 b 1 2 ∆ 3/2 4.30 4.95 5.00 a 5.25 5.70 9.00 b 2 2 ∆ 5/2 4.30 4.95 5.00 a 5.25 5.70 9.00 c 2 Σ -4.30 1/2 4.95 5.00 5.25 a 5.70 9.00	91.5 90.9 86.9 77.5 28.6 97.3 94.8 94.5 93.2 90.9 70.4 2.2 7.1 7.7 11.6 20.8 30.3 1.5 4.4 4.6 6.0 8.4 7.4 18.6 27.5 1.5	6.4 7.0 10.8 20.0 40.0 2.4 b 4.4 4.6 6.0 8.1 19.0 95.4 91.8 91.2 87.6 78.4 45.4 89.2 90.8 90.8 90.0 88.1 3.6 13.8 86.1 9.1 4.5 4.4 3.8 3.4 6.6 b 92.7 6.4 b 93.1 4.8 b 94.2 76.6 b 4.4 99.7 99.3 99.1 98.7 1.6 96.5 1.6 1.6 96.4 14.8 b	7.9 1.0 5.7 24.1 94.0 12.0 92.0 93.2 93.3 93.7 44.9	1.6 1.7 1.7 1.8 2.6 b 1.4 3.7 6.6 5.7 5.1 5.0 4.5 1.0	20.6 10.2 18.2 64.6 51.8 37.8

a An approximate equilibrium distance for this state. b

Table 6 .

 6 Partial and total radiative lifetimes (in s) for transitions from the v ′ = 0 level of the LiSi excited states to the X1 4 Σ - 3/2 (τ1) and X2 4 Σ - 1/2 (τ2) states.

	τ 2

Table 7 .

 7 Partial and total radiative lifetimes (in s) for transitions from the v ′ = 0 level of the LiGe excited states to the X1 4 Σ - 3/2 (τ1) and X2 4 Σ - 1/2 (τ2) states. Numbers in parantheses indicate powers of ten.

	Page 16 of 27
	State	τ 1	τ ⊥	τ 2 τ	τ	τ tot
	a 1 a 2 b 1 b 2 c 2 Σ -2 Π 1/2 2 Π 3/2 2 ∆ 3/2 2 ∆ 5/2 1/2	2.5 4.1 (-3) 113 (-3) 2.4 (-3) 1.2 (-3)	162 (-3) 132 (-3) 4.5 (-3) -2.8 (-3)	11.4 (-3) ---337 (-3)	10.7 (-3) ---2.8 (-3)	10.7 (-3) 3.9 (-3) 4.3 (-3) 2.4 (-3) 1.1 (-3)

a

Table 8 .

 8 Partial and total radiative lifetimes (in s) for transitions from the v ′ = 0 level of the LiSn excited states to the X1 4 Σ - 1/2 (τ1) and X2 4 Σ - 3/2 (τ2) states.

	State X 2 4 Σ -3/2 a 1 2 Π 1/2 a 2 2 Π 3/2 b 1 2 ∆ 3/2 b 2 2 ∆ 5/2 c 2 Σ -1/2	τ ⊥ 10.6 (+4) a 119 (-3) 242 (-3) 305 (-6) -352 (-6)	τ 1 τ -1.78 (-3) ---34.2 (-3)	τ -1.75 (-3) 242 (-3) 305 (-6) -348 (-6)	τ 2 -319 (-3) 680 (-6) 5.98 (-3) 355 (-6) 240 (-6)	τ tot 10.6 (+4) 1.74 (-3) 678 (-6) 290 (-6) 355 (-6) 142 (-6)

a Numbers in parantheses indicate powers of ten.

Taylor & Francis and I.T. Consultant

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

February 2, 2011 11:0


G Weigert

T Henlich

A Klemmt

Modelling and Optimization of Assembly Processes

Keywords: SCHEDULING, OPTIMIZATION, MODELLING, SIMULATION, ASSEMBLY scheduling, optimization, modelling, simulation, assembly planning

The following paper is a case study for scheduling of assembly processes. Two optimization methods, mixed integer programming (MIP) solver and heuristic simulation-based optimization are compared to each other. It can be shown, up to which model complexity the solver methods have the advantage over the simulation-based methods concerning a possible online application. Both, the mathematical formulation as well as the discrete event simulation (DES) system bases on a meta-model -a special Petri net. It is shown that bipartite graphs, i.e. Petri nets, avoid ambiguities in the case of alternative assembly sequences and have advantages in comparison to the usual simple precedence graphs. The results of optimization experiments are described by the example of assembling a numeric controlled (NC) milling table. The model is capable of reflecting multiple parallel assembly processes, alternative assembly steps and can predict characteristics like utilisation, stock or adherence to schedules.

Introduction

In mechanical and plant engineering as well as in other industrial fields such as aircraft construction and shipbuilding, assembly processes are very common. In many cases they are the last part of the manufacturing chain and contribute a great deal to the economic value added. This is one of the reasons for the big interest in the optimization of such processes. The concept of optimization means, in this context to find an optimal schedule for the assembly process; it is not subject of this work to optimize other technological aspects of assembly, for instance motion sequences. In addition our goal is the operational planning and real-time control of the process. This implicates some further important requirements:

1. The control variables as well as the objective variables are determined by the industrial user. 2. The required data must be transferred to the planning system without manual intervention and remarkable delay. Finally this also requires an automatic generation of simulation models. 3. The optimization run itself must not take a long time period to complete, because of the intended online application.

Because of the great importance of assembly processes, their optimization is widely discussed in the literature. As long ago as in 1989 you can find a comprehensive literature-review and analysis of the design, balancing and scheduling of assembly systems [START_REF] Ghosh | A comprehensive literature-review and analysis of the design, balancing and scheduling of assembly systems[END_REF]. Gosh and Gagnon distinguish between single or multi/mixed models respectively between stochastic or deterministic models. In this meaning our problem can be allocated as a multi/mixed deterministic model. They also stated that "…the criteria reported in the literature for designing and operating assembly line systems have typically been singular. … Therefore, the further use of multi-criteria objective functions in modelling assembly systems holds the promise of added research interest and of more closely reflecting the set of real-world objectives.". From this it follows that an investigation of optimization methods should consider multi-objective functions today. Among the optimization algorithms the Genetic Algorithm (GA) has been an important issue for years. Chaudhry and Luo give a general literature overview about the application of GA in production and operations management [START_REF] Chaudhry | Application of genetic algorithms in production and operations management: A review[END_REF]. They highlight the suitability and performance of GA for solving such problems, but deplore a lack of availability of user friendly commercial software packages. So the GA applications are still concentrated in only a few areas of production and operations management. Guo and Wong use indeed GA for solving an assembly line balancing problem with work-sharing and workstation revisiting [START_REF] Guo | A genetic algorithm based optimization model for solving the flexible assembly line balancing problem with work-sharing and workstation Revisiting[END_REF], but they do not compare this algorithm with other heuristics or MIP.

The following paper is a case study which investigates the possibility of an online optimization of assembling a numeric controlled (NC) milling table. Thus two main aspects are discussed in the following: on the one hand a suitable approach for modelling a general assembly process must be found and on the other hand efficient algorithms for optimization are needed. So the article is divided into the two parts, modelling (section 2) and optimization (section 3). In section 3 two different approaches are investigated and compared to each other:

1. Heuristic simulation-based optimization 2. Mixed integer programming (MIP)

While the first approach needs a classical simulation model, the second method is based on a mathematical formulation of the problem. The specific requirements, the advantages and disadvantages as well as the limits of both methods are discussed by means of the results. Both of these methods are fundamentally different and also based on a different approach of modelling. It should be interesting to see which method is suitable for an online application to schedule real processes. Because we are experienced in MIP and in simulation-based scheduling as well (see also [START_REF] Klemmt | Simulation-based optimization vs. mathematical programming: A hybrid approach for optimizing scheduling problems[END_REF]), an objective comparison should be possible as well as a connection of these methods should be possible.

Both approaches are brought together in a common meta-model as a database, which is described in the next section.

The paper is structured as follows: Section 2 introduces the example of an assembly process and distinguishes the importance of assembly alternatives for controlling this process. The advantages of a bipartite graph model are explained there, it is used as a meta-model for describing this processes. In section 3 two basic kinds of optimization methods are explained: the simulation-based optimization and a mixed integer programming (MIP) approach. Section 4 includes an optimization study, which compares several heuristics with the mathematical approach. The comparison should consider the later online application. Section 5 gives a short summery and an outlook as well.

A Suitable meta-model for assembly processes

First of all an assembly process is described by a graph-oriented meta-model, where i.e. the several assembly steps are represented by edges and the results of these steps represented by nodes of this graph, as shown in Figure 1 and 2. The advantage of this meta-model is its clearness for the operator. Using a computer for i.e. simulation, optimization or both, this graph must be mapped into a database. If this database is a relational one, it consists of some tables (Figure 3). One table could contain all nodes, another table contains the edges etc. -this is the database model of the assembly process. The structure of such a database is in turn based on a meta-model. Requirement for every operational planning is a database model describing the assembly process in all its characteristics in a sufficiently exact way and in real-time. For many practical applications the widely-known precedence graph has proved itself [START_REF] Prenting | The precedence diagram: a tool for analysis in assembly line balancing[END_REF], despite its obvious disadvantages. Among all the meta-models which are able to describe the static as well as the dynamic characteristics of the assembly process in an intuitive and formally correct way, bipartite graphs (i.e. Petri nets) have shown their suitability [START_REF] Weigert | Using Petri Nets for DES modeling and optimization of assembly processes[END_REF]. Unlike the well-known precedence graphs, which are commonly used for description of assembly processes in practice, they are unambiguous (Figure 1). In the following the example of the assembly process is described in detail before both graph models are compared to each other.

The example

The chosen process, the simplified assembling of the milling The structure of the assembly process is shown as a set of precedence graphs. It must be pointed out that the pre-products A through E are delivered at specified dates (delayed supply). There are n M mechanics available for assembly. (In the example n M = 2 was used). They can complete assembly steps of the same order in parallel or work at different orders. Thus the organizational order of the assembly steps is open. (see Figure 1, A3: AB can be completed before, after, or in parallel with CD).

Only the precedence relation (also known as technological order) of the respective assembly steps matter. In total there are six orders to be scheduled in this example, with the following objectives:

1. minimal total tardiness T in h 2. minimal average lead time (or cycle time), L in h 3. maximal capacity utilisation, U where 0 ≤ U ≤ 1 (this objective function is only observed, but not included in optimization)

Depending on the alternative chosen there are 3 or 4 assembly steps needed, to get the complete result.

The results of most of the assembly steps are intermediate products (e.g. ACD, BE, etc.).

Figure 1

Precedence graphs

Precedence graphs describe the predecessor-successor relationships between the respective assembly steps and are thus in mathematical terms the representation of an order relation on the set of assembly steps. In event-node form every node of the precedence diagram stands for an assembly step or an article, the directed edges always point from the directly preceding step to the one immediately following it. In this way the AND relations which are typical for assembly processes can be described. By assigning weights to edges in addition, also bill-of-materials relationships (as in a Gozinto graph) can be displayd. Precedence graphs are still widely used in practice [START_REF] Niu | A hierarchical approach to generating precedence graphs for assembly planning[END_REF], even though they have some inherent disadvantages. One of them is the inability to represent alternatives. Alternative assembly paths directly modify the order relation of the assembly process. The example of a machine tool producer described in the following shows this. It describes a subset of the assembly problem introduced in section 2.1. Each of the technological alternatives has its own respective order relation (O 1 -O 5 ) and thus its own precedence graph, as shown in Figure 1 for two of the alternatives. A "superposition" of the precedence graphs leads to ambiguity because AND and OR operations can not be distinguished. This ambiguity would have to be resolved in a yet to be determined manner, e. g. by introducing special labels for the assembly steps. Such a model extension would make the formal handling of the meta-model more difficult, especially in a computing context. That's why other means of description must be investigated

Bipartite graphs

Assembly alternatives can not be ignored because they make a great contribution to the flexibility of the assembly process and offer a substantial potential for optimization which needs to be exploited. So, the precedence graph can give an overview of the structure of the process but it is unsuitable as a metamodel for operational planning because of the missing OR relation. 1990). The practicality of AND/OR graphs in optimization was proven [START_REF] Lambert | Generation of assembly graphs by systematic analysis of assembly structures[END_REF].

Another type of bipartite graph is the Petri net which has places and transitions for nodes. Petri nets provide all requirements to model assembly processes in a structurally unambiguous way (Zhou and Venkatesh 1999, [START_REF] Mînzu | Linguistic Model and a Control Problem for Assembly Workstation[END_REF]. This even includes disassembly processes which can be interspersed with assembly steps (e.g. In the case, where an auxiliary fixture is temporarily used) as well as synchronous or asynchronous parallel completion of steps.

As an extension, Timed Petri nets make it possible to describe assembly durations in the metamodel [START_REF] König | Petri-Netze in der Steuerungstechnik[END_REF]. The type of meta-model used in this work is a Petri net with timed transitions and forward reservation of places [START_REF] Weigert | Using Petri Nets for DES modeling and optimization of assembly processes[END_REF]. It was chosen because it is best fitted to represent the assembly situation in a general manner. For instance it allows us to model some restraints of the assembly process in a consistent way, e.g. the maximum capacity of a station.

Figure 2 The Petri net in Figure 2 shows the simplified assembly process of the NC table including all of the 5 alternatives. In contrast to the precedence graph, there are obviously no ambiguities, because the AND and OR operations can be distinguished by the node type (transition or place) they are pointing to. So the Petri net was modelled in a database from which both, the simulation model for heuristic optimization, as well as the mathematical formulation for the CPLEX solver (the well-known state of the art library, developed by ILOG) were directly derived in a fully automated way. Petri nets are used as the foundation for the relational database so that the entity-relationship diagram in Figure 3 is just another shape of it. The transitions (work steps) and places (articles) have a table each, and the two types of arcs (input and output arcs) also have a table (requirements and results), which makes a total of four tables. The relationship between these entities is implemented with foreign key restraints. Figure 3 shows the entity-relationship diagram of the database. Transformation into simulation models is done via the simulator's database interface (e. g. ODBC or similar SQL based APIs). For both optimization methods described below, basically the same set of objective functions and control variables was used. Differences exist only in method-specific control variables or restrictions. Because many assembly processes are purely order-based, the compliance with agreed delivery dates is especially important. Thus one objective function used here is the total tardiness of all orders which needs to be minimized. The other objective function is the lead time which is also important in assembly processes. The primary control variables are alternative assembly processes which can be classified as a technological control variable. Restrictions are also taken into account in the respective models (e.g. stock level, availability of resources, setup time etc.)

The problems to be solved are usually NP hard, and their solution complexity is further increased by adding assembly alternatives: in a problem with n orders where each can be assembled in m alternative ways there exist m n assembly variants. So it can be shown, for instance, that some of these problems include the problem J2||∑C i , which is NP hard [START_REF] Brucker | Scheduling Algorithms[END_REF].

Simulation-based optimization

For optimization it is first necessary to create a simulation model in a discrete event simulation (DES) system. This simulation model includes both, control variables on the one side and target variables on the other side. A complete simulation run is started with fixed values for the control variables. It results in values for the determined target variables. A special algorithm evaluates the target variables and looks for new values of control variables. So we get a closed optimization cycle as shown in Figure 5. The advantage of this simulation-based optimization method is that the simulation model can be arbitrarily detailed created, inclusive diverse dispatch rules. So, also very complex manufacturing systems can be principally optimized. The model generating process described here was designed from the start with online and real-time capabilities in mind, so it runs fully automated. For storage a relational data-base with SQL interface is used, from which all simulation data is read (Petri nets, orders, availability of resources) and to which the results are written back.

Modelling of places and transitions of the Petri net and of the mechanics is done by using the basis building blocks of the DES system (stations) and the corresponding connections (technological chains, branches). In this manner arbitrary structures of Petri nets can be modelled [START_REF] Weigert | Using Petri Nets for DES modeling and optimization of assembly processes[END_REF]. The orders and their pre-products, intermediate and final products are represented by the simulator's element "job". Figure 4 shows the described model for the simcron MODELLER. Although it is not a Petri net simulator original, you can recognize the structure of the Petri net in Figure 2. The 5 pre-products (A, B, C, D and E) as well as the 6 possible intermediate products (CD, AC, AB, ACD, BE, ABCD) and the final product (ABCDE) -the finished NC milling table -are represented as round elements. The square elements on the bottom of the model are queues, which include the several orders. The elements below the queues connect these elements to each other for an executable simulation model. The figure is a screen shot of one instant state, where several colours (grayscale) represent different states of stations (busy, waiting or empty).

Figure 4 A special characteristic of the model presented here is the use of alternatives (branches in terms of the DES system) as control variables for optimization: Whereas branches have traditionally been used for dispatching purposes when using them as control variables the pass and therefore the choice of a certain alternative is fixed during a single simulation run and will be reset only by the optimization algorithm for the next run.

In order to increase the portability of the model generating method the use of script code to be executed during the simulation run was kept to an absolute minimum. This also results (at least theoretically) in a performance gain, because all operations run in the (usually much faster) simulation core and not in the (slower) script language interpreter. However, it turned out that this performance gain is smaller than originally estimated, because it is partially compensated by the increased model complexity caused by avoiding script code.

The optimization was done with 4 different algorithms: Blind Search, Old Bachelor II [START_REF] Hu | Old Bachelor Acceptance: A New Class of Non-Monotone Threshold Accepting Methods[END_REF]Tsao 1995, Beier 2007), Greedy Search, and a Genetic Algorithm. The algorithm determines the value of the control variables (the alternatives for each order) and sends it to the simulation model. Depending on the target value of the complete simulation run, a new control value is chosen randomly. The three algorithms Blind Search, Old Bachelor and Greedy Search distinguish in the probability of acceptance. The Blind Search accepts every target value and Greedy Search accepts only target values, which are better than all values before. The Old Bachelor algorithm is a special method of threshold acceptance. This means, the acceptance of a target value depends on a dynamic threshold, which becomes lower with an increasing cycle number, but can become higher again, if the last cycle was not accepted. The Genetic Algorithm deals with a number of simulation models with different values of their control variables (individuals of a population). New individuals are created by crossover and mutation. The best solutions of each population will survive, but all of the other solutions are discarded. The performance of these algorithms does not differ essentially, so a more detailed description of them is not necessary. Also the influences of the used population size, mutation type and crossover parameters as well as the threshold values of the other heuristics should be neglected here, because the given time limit is very tight. In addition, the goal of these investigations is not a comparison between different simulation-based heuristics but a comparison between simulation-based heuristics and several methods of mixed integer programming.

Figure 5

MIP-Solver

The goal is the creation of a general mathematical model in form of a Mixed Integer Program (MIP) for optimization of assembly graphs with alternatives and limited resources (resources include articles and mechanics). The problem will be solved by using the CPLEX solver. In order to be able to formulate the problem in a mathematical way, the following data are defined or extracted from the Petri net: Let S be the set of all places of the Petri net and n S their number. Let T be the set of all pre-products (articles) to be assembled into a final product. n T defines their number. B op specifies the supply date of the p-th pre-product of type o ∈ T. In total there are n J orders with due dates d i (i = 1, …, n J ) to be assembled by n M mechanics. The set of all assembly alternatives described by the Petri net is defined by A (see example: A = {A1, A2, …, A5}). Let SA l be the set of all places in alternative l ∈ A. P lm describes the necessary assembling time of place m ∈ S in alternative l ∈ A. Z ∈ S defines the destination node (final assembly). Let the order relation of the places in alternative l ∈ A be specified by O l (see example: O A4 ={{AC,ACD}, {ACD,ABCDE}, {BE,ABCDE}}). This results in n Ol , the number of order relations in alternative l ∈ A. Furthermore the set of places which directly access the pre-product o ∈ T, identified as T O , can be derived. R1 kl and R2 kl identify the first and second place of the k-th order relation (k = 1, …, n Ol ) in alternative l ∈ A (example: R1 2,A4 = ACD and R2 2,A4 = ABCDE).

For creating a mathematical model the following unknowns are to be defined: x il …Order i is produced in alternative l

s im …Start date of article i at place m (i = 1, …, n J ; m ∈ S) t i …Tardiness of order i (i = 1, …, n J ) v io
(i = 1, …, n J ; l ∈ A)
y imq …Order i is assembled at place m by mechanic q (i = 1, …, n J ; m ∈ S; q = 1, …, n M ) z imjn …Order i (place m) is assembled before order j (place n) (i, j = 1, …, n J ; i ≤ j ; i ≠ j ∨ m ≠ n ; m, n ∈ S);

s im ∈ ℝ + ; t i ∈ ℝ + ; v io ∈ ℤ + ; w iop ∈ {0, 1}
x il ∈ {0, 1} ; y imq ∈ {0, 1} ; z imjn ∈ {0, 1}

The described optimization problem can be formulated as follows: Let K be a sufficiently large constant. Equation ( 1) describes the objective function to be minimized. Restriction (2) makes sure that every order can be assembled in exactly one alternative. Equation (3) assures compliance with the order relation of assembly. Equation ( 4) means that if order i is assembled on route l no mechanic can carry out assembly at a place which is not part of that alternative. Restrictions ( 5) and ( 6) request the assignment of exactly one mechanic to one assembly. Equations ( 7) and ( 8) make sure that one mechanic will carry out one assembly step at a time. Restrictions ( 9) and ( 10) as well as ( 15) and ( 16) define which pre-product is used by order i. The equations ( 11), ( 12) and ( 14) convert this information unambiguously into the Boolean unknown w iop , which is used in equation ( 13) to implement the delayed supply of the individual pre-products. Equation ( 17) constricts the objective function to be minimized.

with

(1)

il l A x ∈ ∑ 1 = 1,..., J i n = (2) ( ) , 1 , 1 1 kl kl i R l R il s P K x + - - , 2 kl i R s ≤ 1,..., ; ; 1,..., J Ol i n l A k n ∈ ∈ =
(3)

il imq Kx y + K ≤ 1,..., ; 1,..., ; \ J M l i n q n l A m S SA = = ∈ ∈ (4) 1 M n il imq q Kx y = -∑ 1 K ≤ - (5) 1 M n il imq q Kx y = + ∑ 1 K ≤ +        1,..., ; J l i n l A m SA = ∈ ∈ (6)
( )

3 imq jnq imjn jn im K y y z s s - - - + - lm il l A P x ∈ ≥ ∑ (7) (
)

2 imq jnq imjn in jm K y y z s s - - + + - ln jl l A P x ∈ ≥ ∑        , 1,..., ; , ; 1,..., J M i j n i j i j m n m n S q n = ≤ ≠ ∨ ≠ ∈ = (8) imjn io jo K z v v + - 1 ≥ (9) (
)

1 imjn jo io K z v v - + - 1 ≥      , 1,..., ;
; ,

J O i j n i j o T m n T = < ∈ ∈ (10)
( )

1 iop io K w v - + p ≥ (11) (
)

1 iop K w p - + io v ≥      1,..., ;
1,...,

J J i n o T p n = ∈ = (12)
( ) International Journal of Production Research, Vol. X, No. Y, December 2010, xxx-yyy Due to the large number of unknowns and restrictions the computing time for the exact solution is very high even with small problems. It has to be marked, however, that in praxis there are at least some restrictions which considerably simplify the optimization problem under investigation. These include for instance:

1 iop im K w s - + op B ≥ 1,..., ; 1,..., ; J J O i n o T p n m T = ∈ = ∈ (13) 1 J n iop p w = ∑ 1 = (14) io v 1 ≥ (15) io v J n ≤        1,..., ; J i n o T = ∈ (16) iZ il lZ i l A s x p d ∈ + ⋅ - ∑ i t < 1,..., J i n = (17)
• A mechanic finishes an order alone (no order is handed from one mechanic to another).

• A mechanic can only start a new order after having finished the previous one (because the workspace is limited). • There is a fixed allocation of pre-products to final products.

In these cases an optimal solution can be found even for larger problems. Another approach for reducing problem complexity, without simplifications in model constraints, is the implementation of a time-window approach, identified as "MIP cyclic" in the following. Here, a decomposition of the optimization problem into smaller sub-problems (grouping of orders) and their sequential solution is performed. Then the following objective function is used for every sub-problem:

1 1 min J J n n i iZ il lZ i i l A K t s x P = = ∈   + + ⋅ →       ∑ ∑ ∑ (18)
instead of equation ( 1). Apart from minimization of tardiness this delivers at the same time the earliest possible completion of all orders. This leads to an initial situation for the respective next step of the timewindow approach, which is in each case beneficial. So, this method represents a type of hybrid between simulation and solver.

For further complexity reduction, some restrictions and variables, e. g. equations ( 7) and ( 8), and z imjn respectively, could be made redundant by defining suitable heuristics. Also, the omission of restriction (4) can result in a performance gain as well.

Results

In this section the results of the different optimization approaches to the described assembly problem are presented for two different manufacturing objectives. The first objective is the already introduced reduction of tardiness costs, which play an important role for all kinds of due-date dominated manufacturing problems (like supply chains). The second objective is the reduction of average lead time. Here the delay times of jobs in the manufacturing line are minimized which is important in cases where, due to the large size of the products, assembly space is limited, e.g. in plant engineering or the airspace industry.

Tardiness optimization:

The exact MIP approach was able to find the optimal solution for the NC table problem (for n J = 6) for objective total tardiness (1). This allows a comparison to the quality of the results from the simulation-based optimization. On average it was able to find the exact optimum after approximately 5 minutes of computing time (on an ordinary PC, 2 GHz, 2 GB RAM). A Gantt chart representing the optimal assembly plan is shown in Figure 6. It contains all relevant information such as supply dates of the respective pre-products and delivery dates of the orders. The assembly steps of one particular order are printed in the same colour (grayscale). The frame colour (gray-scale) of the box identifies the chosen assembly alternative. As can be seen, in the optimal solution four different alternatives (A2, A3, 3 times A4, and A5) are used. This shows clearly what big potential for optimization the flexibility of assembly processes with alternatives has to offer. Now in order to further compare simulation-based and solver-based methods, the problem dimension (number of orders) was increased gradually. As criterion for online capability an upper limit of 1 minute for the optimization time was chosen. Table 2 compares the results of the different methods with the problem size.

As can be seen, the exact MIP approach finds the optimum for 5 and 6 orders. However, for n J > 6 the quality of the solution found after 1 minute decreases rapidly. The reason for this is the enormously fast growing number of precedence variables. Already at n J = 10 finding a valid solution becomes difficult. Here the simulation-based optimization shows its advantages because every simulation run yields a valid solution.

Although this method exceeds the chosen limit of 1 minute too, the results for bigger problem dimensions are still clearly better as with the exact MIP approach.

Table 2 A combination of simulation and solver in the form of the cyclic MIP approach consistently yielded the best results. The orders were grouped into "batches" of up to 3 orders per optimization run. Thus the sub-problems remained always exactly solvable independent of the total number of orders to be scheduled.

Average lead time optimization:

To adapt the MIP model concerning the objective average lead time equation ( 1) has to be substituted by ( 19). This also holds for the cyclic MIP approach. Thereby r i described the start of the first possible processing of order i. The definition of variable t i as well as the constraints set ( 17) can be omitted.

1 1 min J n iZ il lZ i J i l A s x P r n = ∈   + ⋅ - →       ∑ ∑ (19)
As in the case of tardiness optimization the exact MIP approach was able to find the optimal solution for the given problem up to 6 orders (n J = 6). In the Table 3 the results for objective average lead time concerning the best MIP approach (MIP-cyclic) and the best simulation-based optimization approach (Genetic Algorithm) are displayed.

Table 3 The reason why the Old Bachelor, Greedy Search, and Genetic Algorithm are not significantly better than Blind Search is the absence of a meaningful distance measure for the control variable "choice of alternatives" (value 1 meaning A1 was chosen, etc.). For that reason in the second optimization series (average lead time) only the Genetic Algorithm was used. Table 4 shows the (order specific) correlation coefficients for the objective functions of total tardiness and lead time for the 6 job problem. As can be seen, they are significantly smaller than 1, thus eliminating the advantages of the three algorithms over Blind Search.

Table 4 This is visually confirmed by a plot of the objective function (total tardiness, z-axis) over the control variables of jobs 1 and 2 in Figure 7 which shows a relatively scattered distribution.

Figure 8 shows the average lead time over the same control variables, while Figure 9 shows all three investigated objective functions in one diagram. The main objective of the presented work was the evaluation of the suitability of different approaches for planning and controlling assembly processes. In addition, it should be shown that assembly alternatives include a big optimization potential. A practical application on NC table assembly of milling machines is currently in preparation. It requires the expansion of the algorithms and of the underlying meta-model, which permits the processing of more than 16 steps per alternative and more than 20 pre-products. This problem dimension excludes the use of MIP to reach an exact solution. First results have shown that both the cyclic MIP approach as well as the simulation-based optimization is suitable for practical adoption. It remains to be explored if and how a combination of these two methods is possible. The potential of this approach exists. One possibility could be to use a MIP solver for subsets of the model and to transfer the solution into heuristic simulation-based systems. Also the usage of improved "intelligent" search methods as well as multicriterial optimization will be continued to be investigated.

  Figure3

  …Position of order i for article o (i = 1, …, n J ; o ∈ T) w iop …Order i requires the p-th delivery of article o (i = 1, …, n J ; o ∈ T; p = 1, …, n J )

  Figure6

  Figure 7,8,9

  Figure 1: Ambiguity in modelling assembly alternatives

Figure 1 :Figure 2 :Figure 4 :Figure 5 :Figure 6 :

 12456 Figure 1: Ambiguity in modelling assembly alternatives Figure 2: Assembly process with 5 assembly alternatives as a Petri net

Figure 7 :Figure 8 :

 78 Figure 7: Scatter plot of objective function total tardiness Figure 8: Scatter plot of objective function average lead time Figure 9: Combined scatter plot of objective functions

Table 1

 1

	gear (A), table top (B), gearbox (C), oil pump (D) and measuring system (E). Assembly is possible in 5
	alternatives (A1-A5) including the intermediate products resulting from the respective steps:
	• Alternative A1: (((A AND C) AND D) AND B) AND E
	• Alternative A2: ((A AND C AND D) AND B) AND E
	• Alternative A3: ((A AND B) AND (C AND D)) AND E
	• Alternative A4: ((A AND C) AND D) AND (B AND E)
	• Alternative A5: (A AND C AND D) AND (B AND E)
	So one can assemble the spur gear (A) and the gearbox (C) at first, subsequently one assemble the oil
	pump (D), then assemble the table top (B) and at last assemble the measuring system (E). This is the
	alternative A1 of the list. Another possibility is to start similarly to alternative A1, but assemble the
	table top (B) and the measuring system (E) separately. At last, both intermediate products, the gear with
	the oil pump (ACD) and the top table with the measuring system (BE) are assembled. This is the F alternative A4 of the list. Both alternatives as well the other three lead to the same final product. The alternatives have a central importance here, because they are the control variables used for o optimization. The duration of the alternatives and their respective steps is shown in Table 1. The left r row of the table shows all the possible assembly steps and the right row shows the duration for this
	step, measured in hour and minute. Step A+B means that A and B are assembled in this step (duration 3 hour 20 minutes), AC+D means that D is assembled to AC (duration 2 hour 30 minutes) while A and C P were assembled in an assembly step before (duration 3 hour 20 minutes). The latter assembly operation e can also be done by one single step A+C+D (duration 4 hour 10 minutes). In the centre of the table the 5 possible alternatives are marked (look at the list ahead). The bottom line of the table finally shows the e total duration for every alternative which is the sum of the durations of every involved single assembly r step (i.e. A1: 18:20=3:20+2:30+6:40+5:50). The objective variables are freely selectable and will be
	discussed in the following sections. The durations for the assembly represent typical values from NC R table assembly practice.
	e
	v i e
	w
	O n l
	y

table, is a part of the assembly process of a NC milling machine. The final product, the milling table (ABCDE), consists of 5 pre-products: spur

uk International Journal of Production Research

	International Journal of Production Research, Vol. X, No. Y, December 2010, xxx-yyy
	F o F o
	r r
	P P
	e e r e e r
	R e v i e w R e v i e w
	O n O n
	l y l y

manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

Table 1 :

 1 Duration of assembly steps and of alternatives A1-A5

	Page 23 of 27
	Assembly step	A1	A2	A3	A4	A5	duration (h:m)
	A+B						3:20
	A+C						3:20
	B+E						2:30
	C+D						5:00
	AC+D						2:30
	A+C+D						4:10
	ACD+B						6:40
	AB+CD						4:10
	ABCD+E						5:50
	ACD+BE						8:20
	Total (h:m)	18:20	16:40	18:20 16:40 15:00

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

G.

Weigert et al.

Table 2 :

 2 Comparison of different optimization methods with respect to problem size (number of orders) when optimization time is limited to 60 seconds, optimum

							Page 24 of 27
	Number of orders
		n J =5	6	7	8	9	10
	Total tardiness/h
	Blind Search	5	7.5 14.16	18.3 20.83 23.33
	Old Bachelor 2	5 7.75 13.92 18.67 22.67 25.83
	Greedy Search	5.83 8.33 13.33	17.5	20	22.5
	Genetic Algorithm	5 7.58 12.5 16.75	19.5	22
	Exact MIP	3.33 5.83 11.66 63.33 135.8 259.1
	Cyclic MIP	3.33 5.83 10.83	15	17.5	20

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

  International Journal of ProductionResearch, Vol. X, No. Y, December 2010, xxx-yyy

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

Table 3 :

 3 Comparison of two different optimization methods with respect to problem size (number of orders) when optimization time is limited to 60 seconds, optimum

	Page 25 of 27
	Number of orders
		n J =5	6	7	8	9	10
	Average lead time/h
	Genetic Algorithm	17.8	19.0 20.1	20.0	20.6	21.0
	Cyclic MIP	17.5	18.6 19.5	19.4	20.1	20.5

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

G.

Weigert et al.

Table 4 :

 4 Correlation coefficients for total tardiness and average lead time

			Page 26 of 27
	Control variable for	Correlation coefficient	Correlation coefficient
	order	(tardiness)	(lead time)
	1	-0.02	-0.15
	2	-0.02	-0.30
	3	-0.002	0.16
	4	-0.009	-0.13
	5	0.04	-0.15
	6	0.07	-0.33

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

  International Journal of ProductionResearch, Vol. X, No. Y, December 2010, xxx-yyy

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

Table 1 :

 1 Duration of assembly steps and of alternatives A1-A5

Table 2 :

 2 Comparison of different optimization methods with respect to problem size (number of orders) when optimization time is limited to 60 seconds, optimum

Table 3 :

 3 Comparison of two different optimization methods with respect to problem size (number of orders) when optimization time is limited to 60 seconds, optimum

Table 4 :

 4 Correlation coefficients for total tardiness and average lead time

	Page 27 of 27

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.

uk International Journal of Production Research

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research


M.D Anita Hryncewicz-Gwóźdź

Tomasz Jagielski

Anita Dobrowolska

Jacek C Szepietowski

Eugeniusz Baran

Identification and differentiation of Trichophyton rubrum clinical isolates using PCR-RFLP and RAPD methods

Keywords: Trichophyton rubrum, dermatophytes, genotyping, PCR-RFLP, RAPD

Trichophyton rubrum represents the most frequently isolated causative agent of superficial dermatophyte infections. Several genotyping methods have recently been introduced to improve the delineation between pathogenic fungi, both at the species and strain level. The purpose of this study was to apply selected DNA fingerprinting methods to the identification and strain discrimination of T. rubrum clinical isolates. Fifty-seven isolates from as many tinea patients were subjected to species identification by PCR-restriction fragment length polymorphism (RFLP) analysis and strain differentiation using a randomly amplified polymorphic DNA (RAPD) method, with two primers designated 1 and 6. Using PCR-RFLP, 55 of the isolates studied were confirmed to be T. rubrum.

Among those, a total of 40 and 5 distinct profiles were obtained by RAPD with primer 1 and 6, respectively. The combination of profiles from both RAPD assays resulted in 47 genotypes and an overall genotypic diversity rate of 85.4%. A dendrogram analysis performed on the profiles generated by RAPD with primer 1, showed most of the isolates (87.3%) to be genetically related. The PCR-RFLP serves as a rapid and reliable method for the identification of T. rubrum species, while the RAPD analysis is rather a disadvantageous tool for T. rubrum strain typing.

Introduction

Dermatophytes comprise a highly specialized group of pathogenic fungi that infect keratinized tissues (skin, hair, and nails) of humans and animals, resulting in dermatophytoses, also referred to as tinea infections. Taxonomically, dermatophytes are classified into three genera: Epidermophyton, Microsporum, and Trichophyton, with the latter being the most complex, containing more than 15 species and several different variants within the species T. mentagrophytes [START_REF] Weitzman | The dermatophytes[END_REF]. The most frequently observed dermatophyte species worldwide is T. rubrum, whose infections usually manifest as tinea pedis and tinea unguium (onychomycosis). The prevalence of T. rubrum as a causative agent of onychomycosis is particularly high and exceeds 90% in Europe [START_REF] Tietz | Changes in fungal spectrum of dermatomycoses[END_REF].

The identification of T. rubrum, by means of conventional laboratory methods may not always be easy or straightforward, since T. rubrum exhibits substantial phenotypic variability. Contrastingly, a high degree of homogeneity of the T. rubrum genome, as revealed by using several anonymous molecular markers, significantly impedes discrimination at the strain level [START_REF] Gräser | Molecular markers reveal exclusively clonal reproduction in Trichophyton rubrum[END_REF][START_REF] Gräser | Microsatellite markers reveal geographic population differentiation in Trichophyton rubrum[END_REF]. Yet, as new molecular typing methods are becoming increasingly available, species determination and strain typing are still being improved.

The aim of this study was to apply some of the recently devised DNA fingerprinting methods to the identification and strain differentiation of T. rubrum.

Material and Methods

The study included 57 isolates of T. rubrum, recovered from as many dermatological patients from Lower Silesia, Poland (40), Kraków, Poland [START_REF] Gräser | Molecular taxonomy of the Trichophyton rubrum complex[END_REF], and Tübingen, Germany [START_REF] Jackson | Species identification and strain differentiation of dermatophyte fungi by analysis of ribosomal-DNA intergenic spacer regions[END_REF].

Clinically, onychomycosis showed the highest number of cases (40 patients; 70.2% of all patients), followed by tinea pedis (12; 21%), tinea corporis (2; 3.5%), tinea cruris (2; 3.5%) and tinea manuum (1; 1.8%).

The isolates were primarily identified as T. rubrum on the basis of their phenotypic characteristics, such as colony surface texture, reverse pigmentation, the ability of micro-and macroconidia formation, and urease activity, assessed by standard mycological procedures.

A rapid mini-preparation procedure was used to extract fungal genomic DNA [START_REF] Liu | Rapid mini-preparation of fungal DNA for PCR[END_REF]. One μL of purified DNA (ca 20 ng) was used for each PCR assay. Identification of T. rubrum was achieved by restriction fragment length polymorphism (RFLP) analysis of the internal transcribed spacer (ITS) region of rDNA, as previously described [START_REF] White | Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics[END_REF]. The amplified products were digested with MvnI, HinfI and where necessary with MvaI restriction endonucleases (Roche), at a 2 hours incubation at 37°C. The resulting restriction fragments were separated electrophoretically on 8% polyacrylamide gels and visualized under UV light after ethidum bromide (Et-Br) staining.

Two primers, designated 1 and 6, as devised by Beaza and Mendes-Giannini [START_REF] Baeza | Strain differentiation of Trichophyton rubrum by random amplification of polymorphic DNA (RAPD)[END_REF] were used in two separate RAPD assays. Amplification products were resolved by electrophoresis using 1.5% agarose gels, and photographed under UV after Et-Br staining. The gel images were analysed by GelDoc system and Quantity One (BioRad, USA) software. A dendrogram of the RAPD profiles obtained with primer 1 was constructed using Dice's coefficient of similarity and the unweighted pair-group arithmetic averaging (UPGMA) clustering method. An 80% genetic relatedness cutoff was used to define clusters.

Results

Of the 57 clinical isolates assessed in this study and recognized as T. rubrum by conventional identification methods, 55 (96.5%) were further confirmed as T. rubrum by means of PCR-RFLP analysis. Using MvnI, all the isolates, except one (841/05), yielded restriction patterns consistent to that of Trichophyton taxon. Using HinfI, 55 (96.5%) isolates produced fragments whose sizes coincided with those expected for T. rubrum. One isolate (899/05) generated a pattern corresponding to either T. mentagrophytes var. interdigitale or T. tonsurans. The isolate was eventually classified as T. tonsurans, based upon MvaI restriction analysis. One isolate (841/05) could not be assigned to any of the dermatophyte species, distinguishable by their unique RFLP profiles obtained with MvnI, HinfI, or MvaI enzymes.

Among the 55 T. rubrum isolates, a total of 40 distinct patterns (A-AN) were obtained by RAPD with primer 1. Two additional patterns were observed for non-T. rubrum isolates (designated as BB and CC for isolates no. 841/05 and 899/05, respectively). The four most prevalent patterns, designated as T, AJ, AE, and AN, contained 5, 4, 3, and 3 isolates respectively. Four patterns (C, AF, AG, and AL) contained two isolates each, and the remaining 34 patterns were represented by single isolates. The RAPD with primer 6 generated five different profiles in T. rubrum isolates (a-e), and two additional profiles in non-T. rubrum isolates (designated as "f" and "g" for isolates no. 841/05 and 899/05, respectively). Among the T. rubrum isolates, the most abundant RAPD profile was that designated as "b", found in 23 isolates, followed by profiles designated as "a" and "c", represented by 19 and 11 isolates, respectively. Patterns designated as "d" and "e" were each identified for one isolate.

The combination of the profiles obtained with both primers resulted in 47 different genotypes for the T. rubrum isolates. Six of those genotypes were common to 3 (type T-a and type AJ-a) or 2 (types: Ca, AE-b, AF-b, AN-a) isolates, whereas the remaining 41 genotypes were unique, that is represented by single isolates only.

The distribution of the genotypes varied among the three geographical regions, from which the isolates originated. Of the 40 RAPD genotypes produced with primer 1, each was restricted to only one geographical location. Twenty-six genotypes occurred in isolates from Lower Silesia, 8 genotypes were observed in isolates from Cracow, and 6in isolates from Tübingen, Germany. Regarding the genotypes obtained by RAPD with primer 6, genotypes "b" and "c" were found in isolates from all three regions studied, and genotype "a" was observed in the Lower-Silesian and Tübingen isolates, but not in the Kraków isolates. The genotypes "a", "b", and "c" were the most prevalent in single regions, i. e. in Lower Silesia (44.7% of isolates), Kraków (62.5%), and Tübingen (44.4%), respectively.

A dendrogram based on the RAPD profiles with primer 1, allowed the separation of the T. rubrum isolates into genetic similarity groups (clusters). A total of 48 isolates could be allocated into six main clusters (I-VI) with the similarity index between the isolates within the cluster being 80% or higher (Fig. 1). The largest cluster (VI) comprised 21 isolates, for which 11 different patterns were obtained.

Clusters IV and V consisted of 8 isolates each, being representative of either 4 (cluster IV) or 8 (cluster V) distinct patterns. The numbers of isolates (and corresponding patterns) belonging to the remaining clusters were 3 (2 patterns) for cluster I, 5 (5) for cluster II, and 3 (3) for cluster III.

Discussion

The traditional, culture-based methods of identifying dermatophytes are cumbersome, laborious, and often inconclusive, due to fungal phenotypic variability and pleomorphism [START_REF] Weitzman | The dermatophytes[END_REF][START_REF] Gräser | Molecular taxonomy of the Trichophyton rubrum complex[END_REF]. However, the advent of molecular biology tools has enabled the development of new molecular approaches to the diagnosis of dermatophyte infections. One such approach, introduced by Jackson et al., relies upon RFLP analysis of PCR amplified ITS regions of the rDNA gene complex [START_REF] Jackson | Species identification and strain differentiation of dermatophyte fungi by analysis of ribosomal-DNA intergenic spacer regions[END_REF]. This PCR-RFLP strategy was also employed in the present study. The only modification to the original procedure was the use of MvnI and HinfI instead of MvaI. Whereas digestion with MvnI allows discrimination between the three main dermatophyte genera (Trichophyton, Microsporum, and Epidermophyton) [START_REF] Dobrowolska | PCR-RFLP analysis of the dermatophyte isolates from patients in Central Poland[END_REF], the HinfI digestion results in species-specific restriction profiles [START_REF] Mochizuki | Rapid identification of Trichophyton tonsurans by PCR-RFLP analysis of ribosomal DNA regions[END_REF]. By using those two restriction enzymes, we wanted to verify their usefulness in the molecular identification of dermatophyte species, and T. rubrum in particular. The results from this study showed high concordance between conventional and molecular techniques for identifying T. rubrum. Only two isolates were identified as non-T. rubrum isolates: one was identified as T. tonsurans (upon restriction analysis with MvaI), and the other was concluded to be a non-dermatophyte fungus. A possible explanation for the latter may relate to the laboratory contamination of the specimen and/or culture. It is also to note that the application of the PCR-RFLP analysis, as in this study, does not allow to distinguish between T. rubrum and T.

soudanense [START_REF] Jackson | Species identification and strain differentiation of dermatophyte fungi by analysis of ribosomal-DNA intergenic spacer regions[END_REF]. The distinction between the two species is only possible with methods that target single nucleotide polymorphisms (SNPs) within the ITS regions of rDNA [START_REF] Kong | Rapid identification and differentiation of Trichophyton species, based on sequence polymorphisms of the ribosomal internal transcribed spacer regions, by rolling-circle amplification[END_REF]. Given however that T. rubrum differs substantially from T. soudanense in terms of phenotypic properties and geographical distribution (the former is cosmopolitan, while the latter is restricted mainly to the sub-Saharan part of Africa), the assignment to the T. rubrum species, provided by PCR-RFLP seems unambiguous. Differentiation of T. rubrum at the strain level has been attempted using a number of genotyping methods, though with unsatisfactory results [START_REF] Gräser | Molecular markers reveal exclusively clonal reproduction in Trichophyton rubrum[END_REF][START_REF] Liu | Use of arbitrarily primed polymerase chain reaction to differentiate Trichophyton dermatophytes[END_REF][START_REF] Zhong | Typing of common dermatophytes by random amplification of polymorphic DNA[END_REF]. However, Jackson et al., have recently reported on intraspecific variability within T. rubrum by PCR amplifying two tandemly repetitive subelements (TRSs), located in the non-transcribed spacer (NTS) region of the rRNA gene cluster [START_REF] Jackson | Strain identification of Trichophyton rubrum by specific amplification of subrepeat elements in the ribosomal DNA nontranscribed spacer[END_REF]. More recently, the RAPD analysis, performed by Baeza et al., with two decameric primers, designated 1 and 6, has been shown to produce a high degree of inter-strain polymorphism [START_REF] Baeza | Strain differentiation of Trichophyton rubrum by random amplification of polymorphic DNA (RAPD)[END_REF][START_REF] Baeza | Strain differentiation of Trichophyton rubrum by randomly amplified polymorphic DNA and analysis of rDNA nontranscribed spacer[END_REF].

The RAPD analysis with primers 1 and 6 was also applied in the present study, and this choice was motivated by a high discriminatory potential of the method, higher than that of the TRS typing system [START_REF] Baeza | Strain differentiation of Trichophyton rubrum by randomly amplified polymorphic DNA and analysis of rDNA nontranscribed spacer[END_REF]. The RAPD typing demonstrated a high genetic diversity of the T. rubrum population studied.

Based on the combined RAPD profiles, a total of 47 distinct genotypes were obtained. Hence, the overall genetic diversity rate (GDR), calculated as the number of different genotypes divided by the number of isolates, was 85.4%. It is noteworthy that most of the polymorphism was generated by RAPD with primer 1. It yielded 40 profiles, whereas RAPD with primer 6 resulted in only 5 profiles (GDRs of 72.7% and 9.1%, respectively). In the first study that used primers 1 and 6, among 10 clinical isolates of T. rubrum, 5 molecular patterns were observed for each primer [START_REF] Baeza | Strain differentiation of Trichophyton rubrum by random amplification of polymorphic DNA (RAPD)[END_REF]. In a subsequent study including 67 T. rubrum isolates, a total of 12 and 11 individual patterns were obtained by RAPD with primer 1 and 6, respectively (GDRs of 17.9% and 16.4%, respectively) [START_REF] Baeza | Strain differentiation of Trichophyton rubrum by randomly amplified polymorphic DNA and analysis of rDNA nontranscribed spacer[END_REF]. Another study that investigated the intraspecific diversity of T. rubrum isolates originating from Japan and China revealed a considerable tightness of the population structure. All 150 isolates tested were split into only 19 fingerprinting genotypes, based on the combined RAPD analyses with primer 1 and 6 (GDR of 12.7%) [START_REF] Yang | Differentiation of Trichophyton rubrum clinical isolates from Japanese and Chinese patients by randomly amplified polymorphic DNA and DNA sequence analysis of the non-transcribed spacer region of the rRNA gene[END_REF]. Much higher genotypic variability with the same primers was shown in a recent study of Santos et al. [START_REF] Santos | Trichophyton rubrum and Trichophyton interdigitale: genetic diversity among species and strains by random amplified polymorphic DNA method[END_REF]. Nineteen different molecular profiles were configured for 52 T. rubrum isolates, when each of the primers was used independently, resulting in a GDR of 36.5%. It was speculated by the authors that the greater genetic diversity revealed in their study might result from having used strains exclusively from patients with onychomycosis. This condition, with its frequent chronicity, has been associated with mixed infections by multiple T. rubrum genotypes [19]. Collectively, the results of the studies cited above differ in terms of genetic polymorphism achieved by RAPD with both of the primers. The polymorphism obtained in our study was exceptionally high, and this may relate to the specific, genetic structure of Polish (and German) isolates, different from that of the so far analysed T. rubrum populations.

The polymorphism obtained by RAPD with primer 1 was further investigated by a dendrogram analysis derived from the similarity coefficients between the RAPD patterns. Although the similarity coefficient value for the entire T. rubrum population studied was calculated at 52%, a vast majority of isolates (87.3%) were distributed into six clusters, within which all the isolates shared at least 80% similarity. This finding, considering that the similarity coefficient value of 0.8-0.99 is assumed to represent genetically related isolates [20], indicates that the clustered isolates might have originated from the same strain that had undergone microevolutionary changes [START_REF] Baeza | Strain differentiation of Trichophyton rubrum by randomly amplified polymorphic DNA and analysis of rDNA nontranscribed spacer[END_REF]. Thus, the high

  19. Santos DA, Araújo RA, Kohler LM, Machado-Pinto J, Hamdan JS, Cisalpino PS (2007) Molecular typing and antifungal susceptibility of Trichophyton rubrum isolates from patients with onychomycosis pre-and post-treatment. Int J Antimicrob Agents 29:563-569. 20. Marol S, Yücesoy M (2008) Molecular epidemiology of Candida species isolated from clinical specimens of intensive care unit patients. Mycoses 51:40-49. 21. Ohst T, de Hoog S, Presber W, Stavrakieva V, Gräser Y (2004) Origins of microsatellite diversity in the Trichophyton rubrum-T. violaceum clade (dermatophytes). J Clin Microbiol 42:4444-4448.

Figure legend Fig. 1

 legend1 Figure legend

Acknowledgements

The authors wish to thank Professor Martin Schaller (Mycological Laboratory, Universitaets-Hautklinic, Tübingen) and Professor Anna Macura (Mycological Department, Jagiellonian University, Kraków), for kindly providing the fungal strains.

polymorphism resolved by RAPD may relate to the reproducibility of the RAPD technique itself. Indeed, RAPD assays often suffer from poor reproducibility and variations in the amplification patterns between isolates may be caused by even the slightest changes in the PCR reaction conditions. This explains the reluctance on the part of researchers to use the RAPD method and their search for newer molecular tools with a better diagnostic performance. One such promising alternative method is multilocus microsatellite typing (MLMT) which has recently been developed by Gräser et al. [START_REF] Gräser | Microsatellite markers reveal geographic population differentiation in Trichophyton rubrum[END_REF].

Finally, the fact that every genetic cluster identified in RAPD with primer 1 was restricted to a single geographical locale, together with marked differences in the frequencies of the three most prevalent genotypes produced from RAPD with primer 6, may suggest that certain genotypes exhibit regional and/or geographical affinities. However, since this study was carried out on a relatively small sample of isolates, the geographical specificity of the genotypes would need to be confirmed by further research involving a larger population.

Interestingly, the geographical differentiation of T. rubrum populations has so far been revealed only by using the MLMT methodology [START_REF] Gräser | Microsatellite markers reveal geographic population differentiation in Trichophyton rubrum[END_REF]21].

In conclusion, the results from this study demonstrated the usefulness of the PCR-RFLP in the rapid identification of T. rubrum, yet insufficient suitability of the RAPD analysis for T. rubrum strain differentiation, due to its poor reproducibility.