XP-PCM Calculations of High Pressure Structural ... - ACS Publications

Mar 4, 2016 - ABSTRACT: The structure and the vibrational properties of the P4S3 crystal at high pressures are discussed by application of the XP-PCM...
0 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCA

XP-PCM Calculations of High Pressure Structural and Vibrational Properties of P4S3 Marco Pagliai,*,†,‡ Roberto Cammi,*,¶ Gianni Cardini,*,‡ and Vincenzo Schettino*,‡ †

Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy Dipartimento di Chimica, Università di Firenze, Via della Lastruccia 3, 50019 Sesto Fiorentino, Italy ¶ Dipartimento di Chimica, Università degli Studi di Parma, Parco Area delle Scienze 17/A, 43124 Parma, Italia ‡

S Supporting Information *

ABSTRACT: The structure and the vibrational properties of the P4S3 crystal at high pressures are discussed by application of the XP-PCM method. The vibrational assignment has been clarified. The structure and the electron distribution changes as a function of pressure are analyzed. The pressure effect on the vibrational frequencies is satisfactorily reproduced and discussed in terms of confinement and structure relaxation contributions.

■

INTRODUCTION In this paper, we present the results of a quantum chemical study on the effect of pressure on the structural and vibrational properties of the P4S3 (tetraphosphorus trisulfide or phosphorus sequisulfide) molecular system. The study exploits a computational methodology (XP-PCM) systematically developed by one of the present authors to study molecular properties and processes of systems dissolved in a dense medium at high pressure.1−3 The XP-PCM method is conceived as an extension of the polarizable continuum model used in the study of molecules in solution at standard pressure conditions.4 Therefore, the target molecule is treated at the quantum mechanical (QM) level and the environment is represented by a continuum hosting the solute into a molecularly shaped void cavity. The medium is characterized by its dielectric permittivity and averaged electronic charge density. The electronic distribution is treated as a uniform continuum defined outside the hosting cavity with density equal to the valence electron distribution of the bulk environment at the given pressure conditions. In the XP-PCM approach the primary action of the pressure is given by the increase of the Pauli repulsion between the solute and the medium5−8 modeled by shrinking the cavity volume with respect to the reference value at standard pressure conditions (p = 1 bar). The shrinking determines an increase of the electronic charge density outside the cavity and of the overlap with the charge distribution of the solvent, with a corresponding increase of the Pauli repulsion. This approach has been employed to study structural and spectroscopic properties of the P4S3 system, which is characterized by a cage molecule as shown in Figure 1. P4S3 is particularly suited for the XP-PCM approach as an almost spherical cage molecule with small directional © XXXX American Chemical Society

Figure 1. Molecular structure of P4S3 seen along and perpendicular to the C3 symmetry axis.

intermolecular interactions. In addition, the vibrational infrared and Raman spectra of the crystal have been studied in some details at pressures up to 86 kbars, showing that no phase transitions or chemical reactions occur in this range.9,10 It will be reported that the observed properties as a function of pressure can be nicely reproduced with the adopted computational approach.

■

THEORY A target molecular system dissolved in a dense medium at highpressure is described by the effective Schrödinger equation o ̂ (Ψ)]|Ψ⟩ = E|Ψ⟩ Ĥ |Ψ⟩ = [Ĥ + Vint

(1)

Special Issue: Piergiorgio Casavecchia and Antonio Lagana Festschrift Received: January 19, 2016 Revised: March 4, 2016

A

DOI: 10.1021/acs.jpca.6b00590 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A where Ho is the electronic Hamiltonian of the isolated molecule and V̂ int(Ψ) the molecule-environment interaction operator given by ̂ (Ψ) = Vê (Ψ) + Vr̂ Vint

the geometrical second derivatives of Ger. Furthermore, to give a rationale of the computed effect of the pressure on the equilibrium geometries and on the vibrational frequencies the XP-PCM model is supported by a specific theoretical tool.1 Theoretical Analysis: Direct and Indirect Effects. In this section the theoretical framework to analyze the computed effects of pressure on equilibrium geometries and vibrational frequencies is summarized. The potential energy function for the nuclei motion, Ger(Q,p), may be formally expressed as a function of the external pressure, p, as

(2)

with V̂ e(Ψ) and V̂ r the usual PCM electrostatic interaction term and the nonelectrostatic Pauli repulsion contribution, respectively, and |Ψ⟩ the wave function of the molecular solute.11 The electrostatic operator V̂ e(Ψ) is represented in terms of a set of polarization point charges, located at the cavity surface and depending on the permittivity of the dielectric medium, on the cavity boundary and on the total charge density of the molecular system. The Pauli repulsion operator V̂ r of eq 2 is represented by a repulsive step barrier potential located at the boundary of the cavity with a height V depending on the number density of the valence electron of the external medium at the given pressure. The basic energy functional to be minimized during the QM procedure has the thermodynamic status of a free-energy, Ger (that corresponds to the hentalpy in the conditions of the present calculations), for the whole molecule-environment system: o

Ger = ⟨Ψ|Ĥ +

1 Q̅ (Ψ) ·V̂ + Vr̂ |Ψ⟩+Vnñ 2

Ger(Q,p) = Ger(Q,0) + p ∑ ΓiQ i + p ∑ ΓiiQ i 2 + ... i

(5)

where Ger(Q,0) is the potential energy function for the nuclei at zero pressure (isolated molecules), as a function of the normal coordinates Q, and the second term is a mixed second-order contribution with ⎛ ∂ 2G ⎞ ⎛ ∂GQ i ⎞ er ⎟⎟ Γi = ⎜⎜ = ⎜ er ⎟ ⎝ ∂p∂Q i ⎠Q = 0 ⎝ ∂p ⎠Q = 0

(6)

representing the direct effect of pressure on the energy gradients component along the normal coordinate Qi.12 The third term of eq 5 is a mixed third-order contribution with

(3)

where Ṽ nn is the nuclei−nuclei interaction contribution in the presence of the external medium, Q(Ψ)·V̂̂ is the PCM electrostatic interaction term,4 and V̂ r is the non electrostatic Pauli repulsion term. The reference state for Ger is given by a hypothetical state comprised by the noninteracting electrons and nuclei of the solute and by the unperturbed medium at the chosen thermodynamic conditions.1 A canonical ensemble approach leads to the pressure p as the negative of the derivative of the free-energy functional, eq 3, with respect to the cavity volume Vc

⎛ ∂G ⎞ p = −⎜ er ⎟ ⎝ ∂Vc ⎠

i

⎛ ∂ 3G ⎞ ⎛ ∂k ⎞ er ⎟ Γii = ⎜⎜ = ⎜ i⎟ ⎟ 2 ⎝ ∂p∂ Q i ⎠Q = 0 ⎝ ∂p ⎠Q = 0

(7)

representing the direct effect of pressure on the force constant ki of the ith normal mode.13 As shown in ref 1, the variations of the equilibrium geometry induced by pressure correspond to atomic displacements along totally symmetric (TS) normal coordinates Qi given by Γ Q ieq(p) = − i p ki

(4)

The cavity is operatively defined in terms of overlapping spheres centered on the atomic nuclei with radius related by a uniform scaling factor (f) to the corresponding atomic van der Waals radii (RivdW). The cavity volume may be freely reduced by decreasing the scaling factor with respect to a reference value f 0 = 1.2. Potential Energy Surface: Equilibrium Geometries and Vibrational Frequencies at High Pressure. For the description of the vibrational frequencies of solvated molecules in a dense medium at high pressure a nonequilibrium solvation approach is adopted which accounts for the dynamical aspects of the solute−solvent interaction due to differences in the relaxation time of the molecular process under consideration (i.e., the vibrational motion of the nuclei) and to characteristic relaxation times of the rotational and translational motion of the molecules of the medium. As the latter are much slower, XP-PCM assumes that the shape of the molecular cavity does not follow the solute vibrational motion. Under this approximation, Ger acts as a potential energy functional of the vibrational motion of the nuclei and, therefore, the equilibrium geometries of the molecule under pressure correspond to the minimum of the potential energy surface and the vibrational harmonic frequencies are obtained from the Hessian matrix of

i ⊆ TS

(8)

where ki is the harmonic force constants of the ith normal mode. Equation 8 clearly shows that the effects of the pressure on the geometry of molecular systems is governed by the pressure coupling coefficients Γi and by the stiffness of the normal modes. This equation is in agreement with previous classical phenomenological formulation of the effect of pressure on the molecular geometry.14−16 The pressure coupling coefficients Γi is related to the direct effect of the pressure on the electronic charge distribution of the target molecule due to the increase of the Pauli repulsion with the environment. The change of the vibrational harmonic force constants with pressure is given by ⎛ kĩ (p) = ki + p⎜⎜Γii − ⎝

TS

∑ j

giij Γj ⎞ ⎟ kj ⎟⎠

(9)

where giij is the cubic anharmonic constant coupling the normal mode i and a totally symmetric normal mode j. Although the overall effect of pressure on the harmonic force constants can be traced back to the primary effect of pressure on the intermolecular Pauli repulsion and the ensuing change of the electron density, eq 9, as discussed in ref 1, separates the B

DOI: 10.1021/acs.jpca.6b00590 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A pressure effect on the harmonic force constants in two basic components: • a direct (or curvature) effect, represented by Γii, which relates the influence of pressure to changes of the electron density of the system; • an indirect (or relaxation) effect, involving the cubic anharmonicities giij, which is related to the pressure induced shift of the equilibrium geometry of the molecular system. Consequently, the pressure behavior of the vibrational frequencies can likewise be partitioned into a curvature and a relaxation contribution:1 ⎛ ∂ν(p) ⎞ ⎛ ∂ν(p) ⎞ ⎛ ∂ν(p) ⎞ =⎜ ⎜ ⎟ ⎟ +⎜ ⎟ ⎝ ∂p ⎠Q(p) ⎝ ∂p ⎠cur ⎝ ∂p ⎠rel

adopting a self-consistent procedure to reach the converged structure for a fixed cavity. After structure optimization a new geometric convergence was searched with a new cavity. The procedure has been repeated to self-consistence (from 2 to 4 iterations were necessary).

■

RESULTS AND DISCUSSION The optimized geometry and normal modes, along with the infrared and Raman intensities, of the P4S3 molecule have been computed with the four chosen functionals and with the 631+G(d) basis set. Detailed results of the calculations are reported in the Supporting Information. The calculated bond lengths of the isolated molecule are reported in Table 1 and are compared with the averaged

(10)

Table 1. Bond Lengths (Å) of Isolated Molecules for Different XC Functionals (P′ Labels the Apical P)

with ⎛ ∂ν(p) ⎞ ⎛ ∂ν(p) ⎞ ⎜ ⎟ =⎜ ⎟ ⎝ ∂p ⎠cur ⎝ ∂p ⎠Q(0)

(11)

⎛ TS ⎛ ∂ν(p) ⎞ Γj ⎞ ⎟ ⎜ ⎟ = ⎜⎜ −∑ giij̃ (0) kj(0) ⎟⎠ ⎝ ∂p ⎠rel ⎝ j Q(0)

(12)

XC

P−P

S−P

S−P′

B3LYP M06-2X PBE0 wB97XD exp29

2.276 2.243 2.246 2.249 2.244

2.141 2.120 2.116 2.120 2.094

2.149 2.121 2.120 2.123 2.091

experimental values from X-ray diffraction.29 The agreement is good and only slightly worse for the B3LYP functional, which also gives a less satisfactory prediction for the vibrational frequencies, particularly for the two highest frequencies modes (Supporting Information). Therefore, in the following we shall mostly refer to calculations with the M06-2X functional. Calculated bond angles reported in Table 2 are in excellent agreement with experiments and do not show any particular dependence on the XC functional.

where ν(p) denotes the vibrational frequencies evaluated at the equilibrium geometries and Q(0) is the equilibrium geometry of the isolated molecule. In essence, the curvature contribution, due to the variation of the second derivative of the potential energy surface, accounts for the confinement effect on the repulsive intermolecular interactions, evaluated at zero-pressure equilibrium geometry, whereas the relaxation contribution is related, through the anharmonic force constants, to the relaxation of the equilibrium geometry of the molecule. A similar partitioning, useful to rationalize observed frequency shifts with pressure, has previously been carried out with a semiclassical approach.17 The present XP-PCM method has the definite advantage of accurately estimating the two contributions for a general molecular system. From a computational point of view, the curvature contribution is easily evaluated by calculation of the vibrational frequencies as a function of pressure at fixed equilibrium geometry in vacuo. Instead, the evaluation of the relaxation contribution is, in general, more complex and computationally expensive, because it requires the calculation of the anharmonic vibrational constants, g̃iij(0), for the isolated molecules. However, this contribution can be obtained as the difference between eqs 10 and 11.

Table 2. Bond Angles (deg) of Isolated Molecule for Different XC Functionalsa XC

SP′S

PSP′

SPP

B3LYP M06-2X PBE0 wB97XD exp29

99.6 99.6 99.8 99.8 99.5

102.4 102.4 102.2 102.2 103.0

103.6 103.6 103.6 103.6 103.0

a

The PPP and SSS bond angles are 60° by symmetry and are not reported.

Vibrational Spectrum at 0 Pressure. In the C3v molecular point group the representation of the normal vibrations of P4S3 classifies as

■

COMPUTATIONAL DETAILS All the calculations have been performed at the DFT level adopting the 6-31+G(d)18 basis set19 with a modified version of the Gaussian 09 rev.C.01 suite of programs.20 Four different exchange correlation (xc) functionals (B3LYP,21,22 M06-2X,23 PBE0,24 and wB97XD25) have been used in the present calculations. The XP-PCM calculations have been performed using the same protocol adopted in a previous paper on fullerenes26 using cyclohexane as solvent and a van der Waals radius of 1.89 and 1.90 Å for S and P, respectively.27 The only difference in the protocol lies in the molecular geometries optimization that has been performed, with very tight convergence criteria,28

Γ = 4A1 ⊕ A 2 ⊕ 5E

(13)

All the modes are infrared and Raman active except for the inactive A2 vibration. The vibrational spectrum of P4S3 has been discussed on the basis of the infrared spectrum of the crystal9,30 and the Raman spectrum of the solution and of the crystal in polarized light.10,31 In principle, the vibrational spectrum of this model cage molecule should be quite simple. However, complications arise from various factors including: • The near frequency coincidence of various modes of the same or different symmetry species is already evidenced in the analysis of the experimental spectra9,10,31 and C

DOI: 10.1021/acs.jpca.6b00590 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A confirmed in early theoretical analysis32 and in the present calculations as well; • The complexity of the infrared and Raman spectrum is enhanced in the crystal by removal of the E modes degeneracy. Indeed, in the crystal sites the molecules only retain a symmetry plane parallel to the ab plane.29 Therefore, the A1 and A2 modes (this latter becoming active in the crystal) split in four Davydov components of symmetry (Ag, B3g, B1u, and B2u) and (B1g, B2g, Au, and B3u), respectively, in the D2h unit cell group. The degenerate modes will split into 8 components, one for each species of the unit cell group. The number of Davydov components is further doubled by the presence in the unit cell of two independent sets of sites; • The isotopic structure of the system with approximate 32 34 abundance of the two main species, P32 4 S3 and P4 S2 S, of 95 and 11% adds to the complexity and asymmetry of the spectrum. In Table 3 the calculated frequencies (in cm−1) of the normal modes of P4S3 (at the M06-2X level) are reported together with

the calculated infrared intensities and Raman relative activities. The calculated intensities of the present work are in overall agreement with those reported by Jensen et al.,32 except for the highest frequency A1 mode. No significant differences in the calculated intensities for the four functionals are found. An overall view of the calculated infrared and Raman spectrum is shown in Figure 2. A comparison with the experimental spectra shows that the calculated intensities can be a reliable guide for the vibrational assignment, suggesting in particular that the strongest Raman (443 cm−1) and infrared (435 cm−1) peaks correspond to two different normal modes. The lowest frequency band, weakly observed in the crystal spectra, can be confidently assigned as the inactive A2 mode. As to the A1 modes, the highest Raman intensity calculated mode is assigned with certainty to the 443 cm−1 Raman band, weakly observed also in the infrared spectrum. The highest frequency modes are calculated at 478 (E) and 483 (A1) cm−1: actually bands are observed at 485 cm−1 in the infrared spectrum at 80 K and at 490 cm−1 in the Raman spectrum. The highest E mode is resolved into a shoulder in the high pressure Raman spectrum, in good agreement with calculated intensities. The lowest A1 can be assigned at 342 cm−1, nearly overlapping the E mode at 340 cm−1, as expected from the calculations. The remaining A1 mode is then assigned as the 421 cm−1 shoulder of the very strong infrared band at 435 cm−1. Once the totally symmetric modes have been localized the assignment of the remaining degenerate modes is straightforward. The resulting vibrational assignment is summarized in Table 3. It can be appreciated that the agreement between calculated and observed frequencies is generally good. Cavity Volume as a Function of Pressure. The free energy has been fitted to the following equation33

Table 3. Calculated and Scaled (0.966) Normal Modes Frequencies (cm−1) and Infrared Intensities and Raman Activities (Arbitrary Units) of P4S3 at 0 P with the M06-2X XC-Functional and Proposed Vibrational Assignment calculated

scaled

symmetry

IR intensity

Raman activity

experimental

196 231 303 363 365 431 451 454 495 500

190 224 293 351 352 417 435 438 478 483

A2 E E E A1 A1 E A1 E A1

0 7 0.5 1.2 0.1 17 27 5 8 0.1

0 6 23 20 9 2 9 55 2 8

185a,b 223a,b 286a,b 342b 343a 421a,b 435a 443b 485a 488b

a

K ⎡⎛ ⎞ V ⎞ p⎛ 1 ⎟ Ger(V ) = G0 + K 0/K p·V (f ) ·⎢⎜ 0 ⎟ ·⎜⎜ ⎢⎣⎝ V (f ) ⎠ ⎝ (K p − 1) ⎟⎠

⎤ K 0·V0 + 1⎥ − ⎥⎦ (K p − 1)

From ref 9. bFrom ref 10.

(14)

Figure 2. Calculated infrared (top panel) and Raman (bottom panel) spectrum of P4S3 with the M06-2X functional. D

DOI: 10.1021/acs.jpca.6b00590 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A where K0 and Kp are the bulk modulus and its pressure derivative, f is the scaling factor of the van der Waals spheres associated with the atoms of the solute molecule, and V0 and V( f) are the initial and varying cavity volume. The relation between the f parameter and pressure is given by the underlying Murnaghan type equation of state34 P(V ) =

K ⎤ ⎡ K 0 ⎢⎛ V0 ⎞ p ⎜ ⎟ − 1⎥ ⎥⎦ K p ⎢⎣⎝ V (f ) ⎠

In the present calculations the f parameter has been allowed to vary in the range 1.0−1.2 and the corresponding pressure (in GPa) obtained by the various functionals ranges up to 21.2 GPa, with only inappreciable variations between the various functionals as shown in Figure 3 Figure 4. Relative cavity volume variation with the pressure in GPa. The reference cavity volume (V0) is for f 0.

Table 4. Bond Lengths (Å) with Pressure According to the M06-2X XC-Functionala

a

Figure 3. Pressure as a function of f for different functionals. Black, blue, green, and red correspond to B3LYP, PBE0, wB97XD, and M062X, respectively.

P (GPa)

P−P

S−P

S−P′

S−S

1.0 2.4 4.9 7.1 10.2 14.7 21.6 rs

2.240 2.237 2.232 2.225 2.219 2.212 2.200 −8.7

2.118 2.116 2.111 2.108 2.104 2.096 2.087 −7.1

2.120 2.118 2.114 2.109 2.105 2.098 2.089 −7.1

3.235 3.231 3.224 3.217 3.210 3.200 3.186 −7.4

rs is the angular coefficient of the linear regression multiplied by 104.

and S−P bonds and the S−S distance as a function of pressure. The relative bond contraction is reported in Figure 5. Although

The lower extrapolated pressure volume of the cavity (241.6 Å3) is significantly larger than the volume per molecule in the crystallographic unit cell (174 Å3).29 This depends on the cavity definition and on the consideration of a solvent accessible surface. A similar effect, but not so pronounced (720 against 650 Å3), has also been reported in XP-PCM calculations on fullerene C60.26 The calculated bulk modulus, 3.87 GPa, is rather small indicating that P4S3 is a highly compressible solid. By comparison, it can be mentioned that for α-phosphorus, a system with similar intermolecular interactions, a bulk modulus of 11.2 GPa has been reported.35 It is possible that the overestimation of the initial cavity volume somehow lowers the estimated bulk modulus even though in the case of fullerenes26 the method gave a bulk modulus in substantial agreement with experiments. The cavity volume variation with pressure, which should mimic the crystal equation of state, is shown in Figure 4. As expected from the low value of the bulk modulus, a substantial volume contraction V/V0 = 0.57 is found in the 0−21.3 GPa range. It can be noted that in the 0−8 GPa range the calculated volume contraction in P4S3 (0.65) nicely compares with experimental findings in α-phosphorus (0.67). Molecular and Electronic Structure at High Pressure. The molecular structure (bond lengths and angles) as a function of pressure is dictated by isotropic compression of the molecular cavity. Table 4 reports the length of the P−P, P−S,

Figure 5. Relative pressure variation of bond distances in P4S3: red, PP; green, SS; blue, SP; black, SP′.

the absolute variations are slightly different for the various bonds, it is evident that the percentage variation are equal within the numerical errors of the computational procedure. This is due to the isotropic compression of the molecule as is more evident from the variations of the bond angles which are always lower than 0.2°, as reported in Table 5. The pressure effect on the molecular structure of P4S3 can be more conveniently analyzed with reference to the normal coordinates Qi of the isolated molecule. In this reference frame, E

DOI: 10.1021/acs.jpca.6b00590 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 5. Variation of Internal Angles (deg) with the Pressure (M06-2X Exchange and Correlation Functional) P (GPa)

PSP′

SPP

SP′S

SSP

SSP′

1.0 2.4 4.9 7.1 10.2 14.7 21.6

102.47 102.49 102.53 102.50 102.48 102.50 102.47

103.59 103.59 103.59 103.61 103.63 103.63 103.66

99.48 99.46 99.41 99.42 99.41 99.39 99.38

76.41 76.41 76.41 76.39 76.37 76.37 76.34

40.26 40.27 40.30 40.29 40.30 40.30 40.31

shifts along the totally symmetric modes 1A1, 2A1, 3A1, 4A1, which preserve the molecular point group of P4S3, will be observed. The corresponding values of the pressure coefficients dQi/dp can be further analyzed according to previous eq 8 in terms of the pressure coupling constants Γi and of the harmonic force constants ki. The pressure coupling constants have been evaluated by a linear fit of the nuclear gradients as a function of pressure and at the equilibrium energy in gas phase. The pressure coefficients have been evaluated by a linear fit of the normal modes coordinates. The results are reported in Table 6.36 It can be seen that the dominant effect on the molecular structure is a shift along the

Figure 6. Position of the Boys orbitals centroids in P4S3 at 1 GPa.

These data show that there is a tendency with increasing pressure to concentrate the electron density toward the interior of the cage. Apart from this, it is clear that bond ibridization in this molecule is far from sp3. An overall view of the change in the electron distribution with pressure is shown in Figure 7 where the regions where the lone pair electron density reduces is shown in red and the regions in the interior of the cage were the electron density increases are shown in blue.

Table 6. Pressure Coefficients ∂Qi/∂p, Pressure Coupling Constant Γi, Harmonic Force Constants in the Gas Phase ki, and the Ratio Γi/ki for Totally Symmetric Normal Modes Qi

Γi (Ehao−1 GPa−1)

ki (Ehao−1 Å−1)

Γi/ki (Å GPa−1)

∂Qi/∂p (Å GPa−1)

1A1 2A1 3A1 4A1

0.00029 0.00026 −0.00177 0.00050

0.2826 0.3739 0.4394 0.5535

0.00100 0.00070 −0.00403 0.00089

0.00090 0.00064 −0.00388 0.00089

normal mode coordinate 3A1 as a consequence of the larger value of the pressure coupling constant Γ3A1. It is of some interest to consider the pressure effects on the electron distribution in the molecule. In Figure 6 the Boys centers at low pressure are shown.37 It is found that the Boys centers are generally slightly displaced from the line joining the bound atoms toward the center of the cage and this is enhanced at higher pressures. In particular it is found that • for all the S atoms the distance of the lone pairs reduces insignificantly (in the fourth decimal figure) at the highest pressure and their orientation angle reduces by 2° to ∼139.14°; • for all the S atoms the bonding center position reduces at high pressure from ∼0.944 to ∼0.933 Å, and the angle Boys center−S−Boys center increases from ∼93.58° to ∼93.95°; • for the three in-plane phosphorus atoms the lone pair distance increases from ∼0.406 to ∼0.416 Å, with a similar slightly smaller increase from ∼0.438 to ∼0.444 Å for the apical P atom; • in the basal plane the angle of the Boys centers through the P atom increases from ∼78.9° to ∼79.4°, and the angle between a planar Boys center the P atom and the P−S Boys center changes from ∼97.1° to ∼97.2°; • the angle through the S atom of the Boys centers of two connected S−P bonds changes from ∼95.1° to ∼95.3°.

Figure 7. Electronic density difference (Δρ) between P = 1.1 GPa and P = 21.5 GPa in P4S3 (isosurface Δρ = 2.5 × 10−4 au) with the B3LYP functional.

Pressure Effects on the Vibrational Modes. The pressure dependence of the frequencies of the Raman active modes in the crystal at room temperature has been reported in the range up to 8.6 GPa.10 For several of the modes a number of Davydov components is resolved as the pressure increases, most likely because of the frequency separation of the various components. In Figure 8 the calculated pressure dependence of F

DOI: 10.1021/acs.jpca.6b00590 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A the normal modes frequencies is reported and compared with experiments.

Table 7. Normal Modes Frequencies Variation with Pressure symm A2 E E E A1 A1 E A1 E A1

ν 195.5 227.2 293.4 339.2 352.7 398.7 412.8 430.7 466.3 471.9

∂ν ∂p

1.19 1.48 0.86 1.86 0.61 2.1 2.34 2.2 2.04 2.7

( ) ∂ν ∂p

cur

1.69 1.54 0.48 0.67 0.23 0.57 0.37 0.45 0.41 0.47

num

( ) ∂ν ∂p

rel

−0.50 −0.06 0.38 1.19 0.38 1.53 1.97 1.75 1.63 2.23

anal

( ) ∂ν ∂p

rel

−0.58 −0.03 0.40 0.98 0.53 1.50 1.41 1.70 2.19 2.33

Figure 8. Vibrational frequencies of the P4S3 crystal as a function of pressure. Dots, cross, etc. are experimental values for the various normal modes. Colors identify the isolated molecule symmetry: A1, A2, and E are in red, green, and blue, respectively. Full lines: computational results from the present work with the M06-2X functional.

It can be seen, as already noted above, that some of the modes are overlapped and the assignment of the specific Davydov components in the Raman spectrum can be difficult. The calculated trend is in general satisfactory agreement with experiments. The present calculations extend over a larger pressure range than experiments and could serve as a useful guide for experiments at higher pressures which could, in turn, constitute a more sensitive test of the method. The pressure coefficients (dν/dp) of the vibrational normal modes of P4S3 can be analyzed in terms of direct and indirect effects as described by eqs 10−12. The direct (or curvature) effect, (dν/dp)cur, is determined by the pressure coupling constant Γii of eq 7 and can be computed by linear fitting of the vibration frequencies as a function of pressure at the fixed gas phase geometry. The indirect (or relaxation) effect, (∂ν(p)/ ∂p)rel, is related to the anharmonicity (giij(0)) and to the pressure coupling constants Γi and harmonic force constants ki(0). The relaxation contribution has been computed both analytically from eq 12 and numerically as the difference between eqs 10 and 11. The linear pressure coefficients of the vibrational frequencies of the normal mode of P4S3 and their partitioning in curvature and relaxation contributions are reported in Table 7.38 These data unveil interesting aspects of the nature of the pressure effects on the vibrational modes of P4S3: • the curvature effect is positive for all the vibrational normal modes. • the relaxation effect is negative for the 1A2, 1E normal modes whereas it is positive for all other normal modes. This behavior is determined by the opposite sign of the anharmonic coupling constants (Supporting Information) • the curvature is the dominant pressure effect on the lower frequencies normal modes (1A2, 1E) • the relaxation is the dominant pressure effect on the normal modes (3E, 2A1, 3A1, 4E, 3A1, 5E, 4A1) A pictorial view of this trend is shown in Figure 9.

Figure 9. Calculated change of normal modes frequencies between normal and highest pressure: red, overall change; green, confinement contribution; blue, relaxation contribution.

■

CONCLUSIONS The essentials of the XP-PCM method to compute the pressure effects on the structural and vibrational properties of a molecular system embedded in a continuum have been reviewed. It has been shown that application of the approach to systems with substantially isotropic intermolecular interactions, like the P4S3 molecule, can give quite satisfactory results in comparison with experiments. In particular, the analysis of the pressure effects on the vibrational frequencies in terms of confinement and relaxation contributions, as a clue to the different pressure behavior of various normal modes, is of considerable interest. In this respect, and from a more general point of view as well, a further extension of the method is desirable to encompass anisotropic compression of the molecular cavity to cope with some of the extended anisotropic effects observed in molecular crystals at high pressures. Work in this direction is in progress.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b00590. Dependence of the pressure from f. Comparison of the cavity volume with f. Normal modes of isolated P4S3 and IR intensity and Raman activity. Absolute value of the projection of atomic forces on the eigenvectors of normal modes at f = 1.2. Projection of atomic forces on eigenvectors of normal modes at f = 1.2. Projection of atomic displacements on eigenvectors of normal modes at f = 1.2. Normal mode variation at fixed geometry (f = 1.2) with f. (PDF) G

DOI: 10.1021/acs.jpca.6b00590 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

■

(16) Thiery, M. M.; Besson, J. M.; Bribes, J. L. High-Pressure SolidPhases of Benzene. 2. Calculations of the Vibration Frequencies and Evolution of the Bonds in C6H6 and C6D6 up to 20 GPa. J. Chem. Phys. 1992, 96, 2633−2654. (17) Moroni, L.; Ceppatelli, M.; Gellini, C.; Salvi, P. R.; Bini, R. Excitation of Crystalline All Trans Retinal under Pressure. Phys. Chem. Chem. Phys. 2002, 4, 5761−5767. (18) Spitznagel, G. W.; Clark, T.; von Ragué Schleyer, P.; Hehre, W. J. An Evaluation of the Performance of Diffuse Function-Augmented Basis Sets for Second Row Elements, Na-Cl. J. Comput. Chem. 1987, 8, 1109−1116. (19) Adding a diffuse function has allowed to obtain a better agreement with the experimental data as reported in the text. (20) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian’09 Revision C.01; Gaussian Inc.: Wallingford, CT, 2010. (21) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (22) Lee, C.; Yang, W.; Parr, R. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron-Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (23) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (24) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0Model. J. Chem. Phys. 1999, 110, 6158−6169. (25) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (26) Pagliai, M.; Cardini, G.; Cammi, R. Vibrational Frequencies of Fullerenes C60 and C70 under Pressure Studied with a Quantum Chemical Model Including Spatial Confinement Effects. J. Phys. Chem. A 2014, 118, 5098−5111. (27) Alvarez, S. A Cartography of the van der Waals Territories. Dalton Trans. 2013, 42, 8617−8636. (28) Maximum force less than 2 × 10−6, RMS force less than 1 × 10−6, maximum atomic displacement less than 6 × 10−6, and RMS atomic displacement less than 4 × 10−6. An ultraf ine grid has been adopted for the evaluation of the numerical integrals. (29) Leung, Y. C.; Waser, J.; et al. The Crystal Structure of P4S3. Acta Crystallogr. 1957, 10, 574−582. (30) Gardner, M. Infrared and Raman spectra of some phosphorous sulphides. J. Chem. Soc., Dalton Trans. 1973, 691−696. (31) Burns, G. R.; Rollo, J. R.; Syme, R. W. G. Raman Spectra of Single Crystals of α-P4S3. J. Raman Spectrosc. 1988, 19, 345−351. (32) Jensen, J. O.; Zeroka, D.; Banerjee, A. Theoretical Studies of the Infrared and Raman Spectra of P4S3 and P4S7. J. Mol. Struct.: THEOCHEM 2000, 505, 31−43. (33) Fukuda, R.; Ehara, M.; Cammi, R. Electronic Excitation Spectra of Molecules in Solution Calculated Using the Symmetry-Adapted Cluster-Configuration Interaction Method in the Polarizable Continuum Model with Perturbative Approach. J. Chem. Phys. 2014, 140, 064114. (34) Murnaghan, F. D. The Compressibility of Media under Extreme Pressures. Proc. Natl. Acad. Sci. U. S. A. 1944, 30, 244−247. (35) Rissi, E. N.; Soignard, E.; McKierman, K. A.; Benmore, C. J.; Yarger, J. L. Pressure-Induced Crystallization of Amorphous Red Phosphorous. Solid State Commun. 2012, 152, 390−394. (36) The reported values refer to the PCM/B3LYP/6-31+G(d) level of calculation, only minor differences are found with different xc functionals. (37) The B3LYP functional along with the SDDALL basis set has been adopted to avoid localization problems due to core orbitals.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: *E-mail: *E-mail: *E-mail:

[email protected]. [email protected]. gianni.cardini@unifi.it. vincenzo.schettino@unifi.it.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by the Italian Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR) under the Contract FIRB - Futuro in Ricerca 2010 No. RBFR109ZHQ. Support and suggestion by Vincenzo Barone are acknowledged.

■

REFERENCES

(1) Cammi, R.; Cappelli, C.; Mennucci, B.; Tomasi, J. Calculation and Analysis of the Harmonic Vibrational Frequencies in Molecules at Extreme Pressure: Methodology and Diborane as a Test Case. J. Chem. Phys. 2012, 137, 154112. (2) Fukuda, R.; Ehara, M.; Cammi, R. Modeling Molecular Systems at Extreme Pressure by an Extension of the Polarizable Continuum Model (PCM) Based on the Symmetry-Adapted Cluster-Configuration Interaction (SAC-CI) Method: Confined Electronic Excited States of Furan as a Test Case. J. Chem. Theory Comput. 2015, 11, 2063−2076. (3) Cammi, R. A New Extension of the Polarizable Continuum Model: Toward a Quantum Chemical Description of Chemical Reactions at Extreme High Pressure. J. Comput. Chem. 2015, 36, 2246−2259. (4) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3094. (5) Schettino, V.; Bini, R. Molecules under Extreme Conditions: Chemical Reactions at High Pressure. Phys. Chem. Chem. Phys. 2003, 5, 1951−1965. (6) Bini, R.; Schettino, V. Materials under Extreme Conditions: Molecular Crystals at High Pressure; Imperial College Press: London, 2014. (7) Grochala, W.; Hoffmann, R.; Feng, J.; Ashcroft, N. W. The Chemical Imagination at Work in Very Tight Places. Angew. Chem., Int. Ed. 2007, 46, 3620−3642. (8) Labet, V.; Hoffmann, R.; Ashcroft, N. W. A Fresh Look at Dense Hydrogen under Pressure. II. Chemical and Physical Models Aiding our Understanding of Evolving H-H Separations. J. Chem. Phys. 2012, 136, 074502. (9) Chattopadhyay, T.; von Schnering, H. Infrared Absorption Studies on Crystalline α-P4S3. Phys. Status Solidi B 1981, 103, 637− 642. (10) Chattopadhyay, T.; Carlone, C.; Jayaraman, A.; von Schnering, H. Temperature and Pressure Dependence of the Raman Spectrum of Crystalline P4S3. Phys. Rev. B: Condens. Matter Mater. Phys. 1981, B23, 2471−2483. (11) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3094. (12) For symmetry reasons, Γi is different from zero only for totally symmetric (TS) normal mode coordinates, and therefore the summation of eq 5 runs only over these normal coordinates. (13) Γii may be different from zero for any vibrational normal mode. In the expansion (5) third-order terms pΓijQiQj describing the effects of pressure on normal modes mixing have been neglected. (14) Fujimura, Y.; Lin, S. H.; Eyring, H. A Theoretical Study of Resonance Raman Scattering from Molecules: High Pressure Effect. Proc. Natl. Acad. Sci. U. S. A. 1980, 77, 5032−5035. (15) Etters, R. D.; Helmy, A. Pressure Dependence of Intramolecular Mode Frequencies in Solid N2,O2, and CO2. Phys. Rev. B: Condens. Matter Mater. Phys. 1983, 27, 6439−6445. H

DOI: 10.1021/acs.jpca.6b00590 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (38) Linear dependence computed analytically and numerically with the B3LYP functional.

I

DOI: 10.1021/acs.jpca.6b00590 J. Phys. Chem. A XXXX, XXX, XXX−XXX