6258
J. Phys. Chem. 1996, 100, 6258-6261
Molecular Models and Calculations in Microscopic Theory of Order-Disorder Structural Phase Transitions: Application to KH2PO4 and Related Compounds Alexandr A. Levin* and Sergey P. Dolin NS KurnakoV Institute of General and Inorganic Chemistry, Leninskii pr. 31, 117907 Moscow, Russia ReceiVed: September 20, 1995; In Final Form: January 15, 1996X
Applications of quantum chemistry approaches are considered to elucidate the mechanisms of structural phase transitions in H-bonded molecular crystals. The KH2PO4 (KDP) family and squaric acid (H2SQ) are studied as examples. Proton rearrangements within H-bonds OsH‚‚‚O a O‚‚‚HsO are treated in terms of variation of the oxygen atom states OOsH a OO‚‚‚H on the basis of the vibronic theory of ligand substitution effects in molecules and complexes. Proton positions within H-bonds are described using the pseudospin formalism. For the energy of crystal this approach results in the Ising model, where the Ising and Slater parameters are expressed in terms of the MO structure of molecular units of crystal. Quantum-chemical calculations are used to determine the MO structure of these units with dependence on the proton positions within the H-bonds. The values of the Slater parameters derived are in reasonable agreement with the experimental data, and the compositional trends for Tc behavior observed for the KDP-type compounds are explained. Applications of the direct quantum-chemical modeling to study of properties of H-bonded crystal are briefly described.
1. Introduction Although quantum chemistry is widely used at present to study various problems of the solid state, its applications to the theory of structural phase transitions (SPT) in crystals are nearly absent. The only exclusion is the vibronic theory of ferroelectricity,1,2 which is based on the second-order Jahn-Teller effect. However, this approach seems to be more successful rather for the description of the displacement-type SPT than for those of the order-disorder type. In this paper the quantum-chemical approach to study of microscopic mechanisms of the orderdisorder SPT in H-bonded molecular crystals3,4 is discussed in some details using the typical example of the KDP-family materials: MIH2AO4 and MID2AO4 (MI ) K, Rb, Cs, and NH4, ND4; A ) P, As). It is applicable to related systems such as squaric acid H2C4O4 (H2SQ), as well. 2. Theoretical Model Proton Redistribution as a “Substitution Effect”. It is well known that the paraelectric phase of KH2PO4 and its 13 isomorphs belong to the I4h2d (Z ) 4) space group (Figure 1) (CsH2PO4 and CsD2PO4 are quasi-1D systems). In these crystals the AO4 tetrahedral groups lie on the crystallographic 4h axes and have the effective S4 site symmetry. The adjacent AO4 tetrahedra are connected through the strong and nearly linear H-bonds (RO‚‚‚O ≈ 2.50 Å).5-8 Therefore, each oxygen atom of any AO4 tetrahedron can exist in one of the two possible states: OOsH or OO‚‚‚H. At the model level any proton transfer in such crystal may be considered as the “ligand substitution” OOsH a OO‚‚‚H at the vertices of AO4 tetrahedra. Such substitution has a cooperative character, i.e. the conversion OOsH w OO‚‚‚H at the vertex of one AO4 tetrahedron results in the opposite conversion OOsH W OO‚‚‚H at the vertex of a neighbor tetrahedron that is bonded to the first one with an H-bond: A
O
H• • •O
A
A
O• • •H
O
A
This treatment of the proton redistribution makes it possible to use the vibronic theory of heteroligand systems,9,10 which * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, March 1, 1996.
0022-3654/96/20100-6258$12.00/0
Figure 1. The crystal structure of the KDP (high-temperature phase).
describes the effects of ligand substitution ALn w ALn-kXY‚‚‚Z in a molecule or complex ALn on the lowest sheet of its adiabatic potential. In the case of closed electronic shells of ALn and ALn-kXY‚‚‚Z,9,10 the difference between the total energies of the systems ALn-kXY‚‚‚Z and ALn in their equilibrium configurations Qf and Q0 ) 0 can be obtained by perturbation theory (eq 1).
E(Qf) - E(0) = S0 - 2∑n*mSnm2/2∆nm -
4∑ν(1/Kν)[∑n*mSnmAνnm/2∆nm]2 (1)
Here Kν and Aνnm are the force constants and orbital vibronic constants, respectively, of the initial system ALn corresponding to the normal or symmetry coordinates Qν. Snm denotes matrix elements of the one-electron substitution operator Hs which represents the effect of ligand substitution in the molecular orbital (MO) basis set of the initial system ALn (φn are the occupied and φm are the unoccupied MO’s):
Aνnm ) 〈φn|(∂H/∂Qν)0|φm〉
Snm ) 〈φn|Hs|φm〉
(2)
2∆nm in eq 2 is the energy gap between MO’s φn and φm. © 1996 American Chemical Society
Order-Disorder Structure Phase Transitions
Jij )
{
J. Phys. Chem., Vol. 100, No. 15, 1996 6259
J| = -[ca2lt2∆R2/8∆][1 + 2(∆⊥ - ∆|)/
∆| + A⊥2/K⊥∆⊥ - 2A|2∆⊥/K|∆|]
J⊥ = [ca2lt2 ∆R2/8∆][1 + A⊥2/K⊥∆⊥]
∆R = ∆|2δP|/ca2lt2
Figure 2. The MO diagram (only σ-bonds) of the AO4 tetrahedron (A ) P, As).
The three terms in eq 1 represent, respectively, (i) the firstorder effect of ligand substitution on MO’s of ALn, (ii) the change in the electronic structure of the ALn system at its fixed geometry due to the admixing of unoccupied MO’s to occupied MO’s, and (iii) the change in the geometry of the initial ALn system. The Ising Model in Terms of MO’s. In accordance with the conventional Slater-Takagi concept11,12 the energy of the KDP-like crystal is determined by the sum of the energies of the four-proton configurations around the AO4 tetrahedra. One can write the energy of each of the AO4 tetrahedra in the form of eq 1 and sum over all tetrahedra in the lattice, taking into account the above-mentioned cooperative coupling between the states of oxygen atoms in adjacent tetrahedra. As a basis set we use the valence ns and np atomic orbitals (AO’s) of the A atoms and the σ AO’s of four oxygen atoms involved in the A-O bonds. Figure 2 shows the MO diagram for σ-orbitals of the AO4 tetrahedron. It refers to the unperturbed system with the local S4 symmetry corresponding to the protons located, hypothetically, at the central positions of the H-bonds (the unperturbed AO4 tetrahedron may be considered as the homoligand system ALn for these proton positions). The calculated e-b level splitting is ∼0.1 of the energy gap 2∆ between the highest occupied MO’s (HOMO’s) and the lowest unoccupied MO (LUMO) of the ideal AO4 tetrahedron with Td symmetry. We define the form of the substitution operator for the AO4 tetrahedra in the AO basis set by assuming that the transfer of the proton from the center of the ith H-bond Oi′‚‚‚H‚‚‚Oi′′ to a position corresponding to Oi′sH and Oi′′‚‚‚H bonds (or vice versa causes a change in the Coulomb (in terms of the Hu¨ckel theory) integral R for the σ(Oi′) AO by ∆Rσiz. The sign and the magnitude of ∆R are determined by calculations (see below). Here, σiz ) (1 is the pseudo-spin, which corresponds to two possible proton positions on the ith H-bond.5,13 Then, for any proton distribution, the energy of the crystal has the well-known Ising form (eq 3):
E ) -1/2∑i*jJijσizσjz
(3)
because the sum of the S0 terms over the lattice vanishes. Now, using the common frontier MO approximation for simplicity, and assuming that the substitution and vibronic operators in every AO4 tetrahedron mixes only HOMO and LUMO, we obtain the expressions for the Ising parameters in eqs 4 and 5.
(4)
∆R = ∆⊥2δP⊥/ca2lt2
(5)
Here, J| and J⊥ are the two types of the Jij parameters possible in KDP. They correspond to pairs of adjacent H-bonds lying in the same xy plane and in different xy planes, respectively. The values 2∆| and 2∆⊥ are the e-a* and b-a* energy gaps, respectively, and 2∆, ca and lt are the HOMO-LUMO gap and MO coefficients for the ideal AO4 tetrahedron, respectively. The value 2δP⊥ denotes the difference in the σ AO populations of the OO‚‚‚H and OOsH oxygen atoms in the case that the AO4 tetrahedron has two upper OO‚‚‚H atoms and two lower OOsH atoms. 2δP| denotes this difference in the case that the AO4 tetrahedron has two lateral OO‚‚‚H atoms and two other lateral OOsH atoms. (In both cases, the orbitals σ(O) involved in the A-O bonds are implied.) A|, A⊥, K|, and K⊥ are the orbital vibronic constants and force constants for the displacement of the central atom, with reference to the oxygen atoms in the xy plane and along the z axis, respectively. Now let us consider the expression for the total energies per formula unit of the ferroelectric (Ef) and antiferroelectric (Ea) phases of the KDP-family crystal. More exactly, we are dealing with the (H2AO4-)∞ molecular framework not taking into account the MI ions and long-range electrostatic forces. We derive from eq 3 the expressions in eq 6
Ef ) -2J⊥ - e0/2
Ea ) -2J⊥ + e0/2
(6)
where the Slater parameter e0 ) 4J| + 4J⊥ (and according to our quantum-chemical calculations, |e0| ≈ 0.1J⊥). 3. Results and Discussion Compositional Trends for KDP and Its Isomorphs. Equations 4-6 provide a qualitative explanation for some major trends in the behavior of the KDP-family materials. It is known that OsD bonds are shorter and O‚‚‚D are longer than the corresponding OsH and O‚‚‚H bonds in the KDP-type systems. Hence, J⊥(D) > J⊥(H) and |Ef| and |Ea| increase upon deuteration. This is in agreement with the experimental critical temperatures (Tc) behavior, which increase upon deuteration for both the ferroelectrics and antiferroelectrics of KDP-type materials.5,13 The effect of the P w As substitution can be explained as well. In molecules and crystals, energy gaps commonly decrease in the vertical groups of the Periodic Table with increasing atomic number. This trend is reproduced by the RHF/3G,3-21G calculations for some P(As)-containing clusters, e.g. [AO4]3-, [H4AO4]+, [H4AO4‚4H2O]+ (A ) P, As). The decrease in the HOMO-LUMO gap of the AO4 tetrahedron results in an increase in the values (∆⊥/∆|), (A|2/K|∆|), (A|2/ K|∆|)/(A⊥2/K⊥∆⊥), consequently, and the value |e0| increases (e0 < 0). Thus, the P w As substitution stabilizes Ea and destabilizes Ef, i.e. Tc decreases for ferroelectrics and increases for antiferroelectrics of the KDP type. This conclusion is in full agreement with the experimental data.5,13 Calculations and Estimation of the Ising (Slater) Parameters. Our electronic-structure calculations were carried out by the conventional MNDO method14 for the systems without H-bonds and by its modified version MNDO/H15,16 adapted to study H-bonded systems. In the latter only the core-core
6260 J. Phys. Chem., Vol. 100, No. 15, 1996 interaction function f[RYH] ) exp[-RRYH2] for the interaction between the H and Y ) N, O atoms in the H-bond XsH‚‚‚Y (or the functions f[RY(X)H] ) f[RY(X)H] - CY(X) exp[-(RY(X)H RX(Y)H)2] for the interacting H and Y(X) atoms in the H-bond X‚‚‚H‚‚‚Y with the proton transfer) have been implemented into the potential function EY(X)Hcore of the standard MNDO procedure.14 The recommended values for R0 and C0 are 2.0 Å-2 and 0.012, respectively.15,16 The molecular models with successively larger clusters are employed both to estimate the Ising parameters and to determine the equilibrium positions of the MI ions as a function of the proton ordering. Most of the clusters of interest are too large for ab initio treatments on the appropriate computational level (for example, HF/6-31G/MP2). By taking into account the experience of the MNDO/H procedure applications, it seems to be successful in calculating the relative values of heats of formation and atomic charges for the close atomic configurations distinguished, for example, by the proton ordering. For the estimation of the Ising parameters Jij involved in eq 4 or the Slater parameters (e0 ) 4J| + 4J⊥, w ) 2J| + 4J⊥) the electronic structure of some KDP structural fragments [[H4PO4]+ (I); [H4PO4 4OR]+, R ) H2, PH3, P(OH)3 (IIa-c); [H4PO4‚ 4OP(OH)3‚12OR]+, R ) H2, PH3 (IIIa,b)] was calculated primarily by the MNDO/H method. It should be noted that to yield the more correct estimations of δP| and δP⊥ it is essential to calculate these values for the oxygen atoms bonded rather with the bridge than with the terminal protons. Therefore, the pentamer IIc appears to be the most appropriate molecular model of the minimal size, and atomic charges on the oxygen atoms of the H-bonded central PO4 tetrahedron must be used to obtain δP| and δP⊥. Although four adjacent PO4 tetrahedra, surrounding the central one, are the constituents of the crystal, they play a role of the saturating groups in our modeling as well. One can expect that the nature of such groups has no decisive importance for the description of the OsH‚‚‚OR or O‚‚‚HsOR linkages. Therefore, the clusters IIa,b with R ) H2 and PH3 were studied too. These saturating groups have been used in the clusters III as more simple than R ) P(OH)3. The cluster I was examined for comparison, where the H-bonds were simulated by the “long” and “short” terminal O-H bonds with the experimental interatomic distances 1.07 and 1.42 Å, respectively. The use of the clusters IIIa,b, consisting of the cluster IIc with 12 terminal groups R ) H2, PH3 permit the understanding of the influence of the second nearest neighbors on the charge distribution within the central PO4 tetrahedron. Let us use the computational results for estimation of the Ising or the Slater parameters. For the sake of simplicity, the terms containing vibronic constants are neglected in eq 4, because, for large gap 2∆, the omitted terms are insignificant.10 The calculations of the systems II-III with the fixed (experimental) geometries exhibit the nearly monotonous decrease of e0 along with some increase of w: |e0| ≈ 210-420 K, w ≈ 840-1680 K (IIa) and |e0| ≈ 80-160 K, w ≈ 9001800 K (IIIb). These rough estimations of the Slater parameters are in reasonable agreement with the values w ≈ 800-1100 K and e0 ≈ 80-115 K, obtained by using the thermodynamic properties of KDP and DKDP.17 The estimations derived probably refer to other members of the KDP family [KDA (KH2AsO4), DKDA and RDP (RbH2PO4), DRDP] also because their electronic structure are similar to that of KDP and DKDP. The charge distributions derived in the calculations allow also the estimation of the dipole moment µ of the H2PO4- ion relative to the center of the PO4 tetrahedron (in the point-charge approximation). The calculated values [µcalc ≈ -4 - (-5 D)]
Levin and Dolin
Figure 3. The [MI(H4PO4+)6] cluster used in study the effect of MI ion on the proton ordering in the KDP-type compounds. The antiferroelectric proton ordering is shown (MI ) NH4).
are in reasonable agreement with the “experimental” one (µexp ≈ -6.7 D) obtained from the ferroelectric polarization data.18 The optimization bond lengths and bond angles of the PO4 tetrahedra as well as the mean length of the H-bonds (in finite clusters they are, of course, different) are close to the experimental diffraction data. The deflections do not exceed several hundredths and degrees, but the scatter of RO‚‚‚O is large (≈0.15 Å). In the above-mentioned analysis the MI ions are only considered as the structureless donors in relation to the (H2AO4-)∞ framework. Besides, the long-range interaction contributions to Ef and Ea (eq 6) are neglected too. The correct determination of a role of both factors is, of course, a difficult problem. An approximate estimation of the MI electronic and ionic polarization contributions to stabilization of the lowtemperature phases show that their sum increases in the series K f Rb f Cs. These contributions favor somewhat stabilization of Ef in comparison to Ea for the ferroelectrics of the KDP family. Influence of the MI Ions on the Proton Ordering. It should be noted that the direct quantum-chemical simulation is often a suitable way to investigate the SPT microscopic mechanisms as well.19 Particularly, the influence of the MI ions on the ferroelectric or antiferroelectric behavior of the KDP-family systems was examined using this approach.19,20 For this purpose the electronic structure of the clusters [MI(H2PO4-)6], [MI(H4PO4+)6], (Figure 3), [MI(H4PO4+)2(H2PO4-‚OH2)4], and [MI(H4PO4+)2(H2PO4-‚HOPH3+)4] was calculated (MNDO, MNDO/H). They simulate the nearest-neighboring environment of the MI ions in the KDP-type crystal structure, where every MI ion is coordinated by eight oxygen atoms forming two interpenetrating tetrahedral disphenoids (Figure 1). In the two latter clusters the second coordination sphere of the MI ion is taken into account as well. The calculations were carried out for different positions of protons, modeling their ordering in the ferroelectrics or antiferroelectrics of the KDP family. According to the results of our quantum-chemical calculations the alkali metal cations favor the ferroelectric ordering of the protons in the H-bonds, even when it is in the center of the disphenoids. On the contrary, NH4+ ions stabilize the antiferroelectric alternative, in particular, when the NsH‚‚‚O bonds (Figure 3) are taken into account in the frame of the MNDO/H procedure. The displacement of the MI ions along z or x (y) axes derived from the potential energy curves are not far from the experimental diffraction data: for the NH4+ ion they change in the range 0.20-0.25 Å (the experimental value is ≈0.10 Å)
Order-Disorder Structure Phase Transitions
Figure 4. The unit cell (b-axis projection) of the layered crystal structure of H2SQ (low-temperature phase). The energies of four-proton configurations are also shown.
and for the alkali metal ion they change in the range 0.030.09 Å (the experimental value is ≈0.04 Å).7 The results obtained allow the elucidation of the factors that determine the nature of the phase transitions in KDP-type crystals. Model and Calculations for H2SQ. As has been already mentioned, this approach can also be applied to other materials containing H-bonded polyhedra AOn and to more complex molecular units, for example, to crystalline squaric acid H2C4O4 (H2SQ). This material is a layered molecular crystal (Figure 4), which is usually referred to as the two-dimensional analogue of KDP. The structure of H2SQ is monoclinic with space group P21/m (Z ) 2) below Tc ) 371 K (or 516 K in D2SQ) and the SPT changes the symmetry of H2SQ (D2SQ) to tetragonal with space group I4/m.21 On the model level we shall consider the “C4-core” of the H2SQ molecule as a pseudoatom surrounded by four O(H) ligands. The valence orbitals of this pseudoatom are the σ and π MO’s of the C4-core containing σ (C) and π (C) AO’s, which participate in the C-O bonds. It is easy to obtain the expression for the total energy of the single layer in a form of eq 3. We omit the vibronic contributions for simplicity. Then, the independent Ising parameters result (in the frontier MO approximation):
Jcis = - cag2lbg2∆Rσ2/8∆bgag* - cau2cbu2∆Rπ2/8∆au*bu* Jtr = Jcis + cag2leu2∆Rσ2/4∆euag* + ceu2lbg2∆Rσ2/4∆bgeu* + cau2ceg2∆Rπ2/4∆au*eg* (7) Here Jcis and Jtr refer to the pseudospin coupling in the cis and trans positions within the monomer fragment, respectively; cag, lbg, etc., are the coefficients in the MO’s of the C4O4 fragment, and 2∆bgag*, etc., are the gaps between these MO’s. ∆Rσ and ∆Rπ are the variations of the σ(O) and π(O) AO energies depending upon the state the oxygen atom (OOsH or OO‚‚‚H). One can express ∆R in terms of the electronic populations of the σ(O) and π(O) AO’s (in analogy to eq 5). The quantum-chemical calculations for the appropriate structure fragments of H2SQ were also carried out. They permitted the estimation of the Ising parameters and the derivation of conclusions concerning the nature of the phase transition in H2SQ (D2SQ). In the MNDO/H calculations structural fragments of various types have been used: monomers [H4C4O4]2+, [H4C4O4‚4OH]2-, dimers [H7(C4O4)2‚6OH]3-, [H7(C4O4)2‚6H2O]3+, monolayer tetra- and pentamers [(H2C4O4)4], [(H2C4O4)5], double-layer tetrahedral tetramer [(H2C4O4)‚ 12H2O], and some other model clusters with the OH2-terminal group. The calculations yielded that the protons of each layer in the low-temperature phase could exhibit ferroelectric ordering and that in the polarizations the adjacent layers are opposite. In contrast to KDP, where e0/w ≈ 0,1, the corresponding values for H2SQ are: e0 ∼ w ≈ 500-1500 K (e0 ) 4Jtr + 4Jcis, w )
J. Phys. Chem., Vol. 100, No. 15, 1996 6261 2Jtr). This is in agreement with the results of the Monte Carlo simulation of the 16-vertex model22 (e0 ≈ 750 K, w ≈ 1100 K), but differs from the conclusion in ref 23 (e0 ≈ 0, w ≈ 1100 K), where the trans configurations (the second line in Figure 4) were excluded as forbidden by the classical structural chemistry. The transition from the low- to the high-temperature phase is due to proton disordering within each layer in the crystal. Thus, H2SQ as KDP is a system with competitive interactions (as Jtr ∼ |Jcis| as |J|| ∼ J⊥). However, in distinction from KDP the π-electron delocalization in H2SQ is large and results in stabilization trans configurations, the energies of which become comparative with those of the single-charged configurations. 4. Concluding Remarks In conclusion, we will highlight the physical meaning of the expressions for Jij in eqs 4 and 7. They result from eq 1 and correspond to pseudospin interaction due to the change in the internal energy of H-bonded molecular entities caused by the redistribution of protons. Thus, the mechanism of the indirect pseudospin interaction is considered, which is related to the reorganization in the electronic and geometrical structure of the non-hydrogen framework of crystal. The numerical estimations of the Ising parameters for the KDP-type crystals given above indicate that such a molecular mechanism is important and probably predominant. Acknowledgment. This work is supported by the Russian Fundamental Science Foundation, Project No. 93-03-18530. References and Notes (1) Bersuker, I. B. The Jahn-Teller Effect and Vibronic Interactions in Modern Chemistry: Plenum Press: New York, 1984. (2) Kristoffel, N.; Konsin, P. Phys. Stat. Sol. (B) 1988, 149, 11. (3) Levin, A. A.; Dolin, S. P. Dokladi RAS 1995, 341, 638 (in Russian). (4) Levin, A. A.; Dolin, S. P. Proceedings of XIIth Simposium on the Jahn-Teller Effect, Tartu, Estonia, 1994; p 45. (5) Lines, M. E.; Glass, A. M. Principles and Applications of Ferroelectrics and Related Materials; Clarendon Press: Oxford, 1977. (6) Matsubara, T. Proc. 6th Intern. Meeting on Ferroelectricity, Kobe 1985. Jpn. J. Appl. Phys. 1985, 27, Suppl. 24-2, 1. (7) Nelmes, R. J.; Tun, Z.; Kuhs, W. F. Ferroelectrics 1987, 71, 125. (8) Vaks, V. G.; Zinenko, V. I.; Schneider, V. E. Uspekhi Phys. Nauk 1983, 141, 629 (in Russian). (9) Levin, A. A. J. Coord. Chem. 1993, 19, 343 (in Russian). (10) Levin, A. A.; D’yachkov, P. N. Electronic Structure, Geometry, Isomerism and Transformation of Heteroligand Molecules; Nauka: Moscow, 1990 (in Russian). (11) Slater, J. C. J. Chem. Phys. 1941, 9, 16. (12) Takagi, U. J. Phys. Soc. Jpn. 1948, 3, 271. (13) Blinc, B.; Zeks, B. Soft Modes in Ferroelectrics and Antiferroelectrics: North-Holland: Amsterdam, 1974. (14) Dewar, M. J. S.; Thiel, W. J. Am. Chem. Soc. 1977, 99, 4899. (15) Burstein, K. Y.; Isaev, A. N. Theor. Chim. Acta. 1984, 64, 397. (16) Burstein, K. Y.; Isaev, A. N. J. Struct. Chem. 1984, 25, 25 (in Russian). (17) Vaks, V. G.; Zinenko, V. I. J. Exp. Theor. Phys. 1973, 65, 650 (in Russian). (18) Bystrov, D. S.; Popova, E. A. Solid State Phys. 1981, 23, 1461 (in Russian). (19) Levin, A. A.; Dolin, S. P. J. Chem. Phys. 1995, 14, 84 (in Russian). (20) Lebedey, V. L.; Dolin, S. P.; Levin, A. A. J. Inorg. Chem. 1995, 40, 1683 (in Russian). (21) McMahon, M. I.; Nelmes, R. J.; Kuhn, W. F.; Semmingsen, D. Z. Krystallogr. 1991, 195, 231. (22) Muser, H. E. Ferroelectrics 1984, 55, 275. (23) Zinenko, V. I. Phys. Stat. Sol. (b) 1976, 78, 721.
JP9527970
