Home
Search
Collections
Journals
About
Contact us
My IOPscience
Structural evolution of amino acid crystals under stress from a non-empirical density functional
This article has been downloaded from IOPscience. Please scroll down to see the full text article. 2012 J. Phys.: Condens. Matter 24 424209 (http://iopscience.iop.org/0953-8984/24/42/424209) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 128.179.149.88 The article was downloaded on 21/08/2013 at 15:34
Please note that terms and conditions apply.
IOP PUBLISHING
JOURNAL OF PHYSICS: CONDENSED MATTER
J. Phys.: Condens. Matter 24 (2012) 424209 (9pp)
doi:10.1088/0953-8984/24/42/424209
Structural evolution of amino acid crystals under stress from a non-empirical density functional 1,2,3 , Brian Kolb4 , T Thonhauser4 ¨ ¸ ukbenli ¨ Riccardo Sabatini1 , Emine Kuc and Stefano de Gironcoli1,2 1 2 3 4
International School for Advanced Studies (SISSA), Trieste, Italy CNR-IOM DEMOCRITOS Simulation Center, Trieste, Italy Ecole Polytechnique F´ed´erale de Lausanne (EPFL), Lausanne, Switzerland Department of Physics, Wake Forest University, Winston-Salem, NC 27109, USA
E-mail:
[email protected]
Received 23 March 2012, in final form 25 May 2012 Published 3 October 2012 Online at stacks.iop.org/JPhysCM/24/424209 Abstract Use of the non-local correlation functional vdW-DF (from ‘van der Waals density functional’; Dion M et al 2004 Phys. Rev. Lett. 92 246401) has become a popular approach for including van der Waals interactions within density functional theory. In this work, we extend the vdW-DF theory and derive the corresponding stress tensor in a fashion similar to the LDA and GGA approach, which allows for a straightforward implementation in any electronic structure code. We then apply our methodology to investigate the structural evolution of amino acid crystals of glycine and L-alanine under pressure up to 10 GPa—with and without van der Waals interactions—and find that for an accurate description of intermolecular interactions and phase transitions in these systems, the inclusion of van der Waals interactions is crucial. For glycine, calculations including the vdW-DF (vdW-DF-c09x) functional are found to systematically overestimate (underestimate) the crystal lattice parameters, yet the stability ordering of the different polymorphs is determined accurately, at variance with the GGA case. In the case of L-alanine, our vdW-DF results agree with recent experiments that question the phase transition reported for this crystal at 2.3 GPa, as the a and c cell parameters happen to become equal but no phase transition is observed. (Some figures may appear in colour only in the online journal)
1. Introduction
of the damped dispersion type (DFT-D) [9–13] have been widespread and often successful but cannot be deemed satisfactory from a conceptual point of view. Recently, the new scheme vdW-DF [14] (together with its variants) has been proposed to account for dispersion interactions in a seamless and unified way, fully within DFT formalism. This approach has been shown to yield reliable binding energies for a wide variety of systems with dispersion interactions [1–7]. Current implementations of vdW-DF, calculating the total energy, allow us to study several problems related to energetics. However, often a complete structural relaxation with van der Waals interactions is crucial. The common practice with vdW-aware functionals has been to use experimentally determined cell parameters, or—for simple
The van der Waals (vdW) interactions are of critical importance for a wide variety of physical and chemical phenomena and materials, ranging from the stability of molecular crystals [1] and simple dimers [2], to physisorbed molecules [3], water [4], DNA [5], and drug design [6]. They further play a crucial role in the search for novel materials in energy applications [7]. Traditional density functional theory with local or semilocal exchange–correlation functionals cannot account for van der Waals (vdW) interactions [8], leaving a wide range of problems beyond the scope of first-principles studies. Augmentation of local or semilocal DFT approximations with empirical vdW interaction energies 0953-8984/12/424209+09$33.00
1
c 2012 IOP Publishing Ltd Printed in the UK & the USA
J. Phys.: Condens. Matter 24 (2012) 424209
R Sabatini et al
Bravais lattices—to deduce structural parameters using equation of state fits. Therefore, as the electronic structure strongly depends on the structural parameters, many studies of materials under high pressure and phase transitions are still beyond the scope of vdW-DF calculations. One such area that requires both the inclusion of dispersion forces and calculation of the stress tensor is a first-principles study of organic materials under pressure. These studies are of great interest in various areas, from the industrially relevant issue of pressure treatment of food products to academically significant studies of extremophile bacteria that live in the Earth’s crust or in the depths of the oceans. High-pressure studies are also important while exploring new phases of molecular crystals and their stability, e.g. during a crystal structure prediction. In this paper, we address this issue by extending the theory of vdW-DF to the calculation of the stress tensor, which allows for a full relaxation of the crystal structure and, hence, pressure studies. We then apply our scheme to the case of glycine and alanine crystals under pressure, and test the adequacy of several variants of the vdW-DF functional while investigating the relative stability of known glycine polymorphs and a debated pressure-driven phase transition in alanine.
our discussion will be kept as general as possible, thereby maintaining applicability to any present and future functional in this class. The non-local correlation described above is typically a small effect, in magnitude; in most cases it is only on the order of 1–2% of the overall exchange–correlation energy. Nevertheless, its contribution is often vital and can have significant effects. For example, many simple molecular dimers do not bind without it [2]. And even for crystals, this small effect can result in completely different ground state crystal structures [20]. The contribution of the non-local correlation term to the Kohn–Sham effective potential [21] is obtained by taking the functional derivative of the corresponding energy with respect to the density for a general point Z ∂ (n8) 0 nl dr vc (r) = n(r0 ) ∂n Z X ∇α n 0 ∂ (n8) 0 − ∇α n(r ) dr . (2) |∇n| ∂|∇n| α Notice that this expression reduces to the well known one appropriate for a local or semilocal functional [22] if the explicit r − r0 dependence in the non-local kernel is reduced to a delta function, thus transforming the non-local functional into a local one: √ nn0 8(n, n0 , |∇n|, |∇n0 |, |r − r0 |) 2 p LDA/GGA √ 0 ( nn , |∇n| |∇n0 |) δ(r − r0 ). (3) −→ Fc
2. Theory Dispersion interactions originate from subtle correlations between charge fluctuations in different parts of an extended system [15–17], which cannot be accounted for by any local or semilocal functional that focuses on the density behavior in the neighborhood of a single point [8]. A fully ab initio description of these effects can be obtained from the detailed knowledge of the frequency-dependent electronic polarizability in the adiabatic coupling fluctuation-dissipation (ACFD) formulation [18] of DFT and much work is currently being devoted to this approach. Unfortunately, its computational cost remains unaffordable for even moderately large systems. The non-local correlation functional vdW-DF [14], proposed by Dion et al, is based on the ACFD approach, where a simplified LDA-like treatment of the electronic polarizability is assumed. This approach has proven to be a breakthrough in the field; although more complex than ordinary LDA/GGA functionals, it can still be applied to large systems. Despite its simplifications, it has been shown to successfully address vdW interactions in a number of otherwise problematic situations [19]. vdW-DF and all of the variants in this family add to a standard LDA/GGA term an explicitly non-local correlation term of the form ZZ Ecnl = 21 n(r)8(n, n0 , |∇n|, |∇n0 |, |r − r0 |)n(r0 ) dr dr0 ,
In order to evaluate the stress tensor, Nielsen and Martin’s procedure [23] is applied, where the crystal cell is homogeneously strained, r → (1 + ) r, and the wavefunctions and density are rigidly scaled. Thanks to the variational principle for the total energy, no further adjustment is required. After some tedious but straightforward manipulations, the non-local correlation contribution to the stress tensor can be expressed as Z 1 nl nl nl σc αβ = n vc dr δαβ E − c ZZ 1 ∂8 ∇α n ∇β n − nn0 dr dr0 ∂ |∇n| |∇n| ZZ 1 1 ∂8 Ecnl δαβ + nn0 + 2 ∂ |r − r0 | (r − r0 )α (r − r0 )β 0 × dr dr . (4) |r − r0 | Here, the first term has the same form as the LDA contribution to the stress tensor. The second term is similar to the additional term present for gradient corrected functionals and reduces to it in the special case where the non-local kernel 8 is forced to be local (see equation (3)). The third term in the above expression constitutes the real novelty introduced by the use of an explicitly non-local functional and can be shown to vanish if 8 is taken local. The vdW-DF kernel depends on the value of the density and its gradient at two points in space but only through a
(1) where n and n0 are shorthand forms for n(r) and n(r0 ), respectively, and different functional variants differ in the functional form of the non-local kernel 8. In the following, 2
J. Phys.: Condens. Matter 24 (2012) 424209
R Sabatini et al
one being the γ phase [29]. The evaporation of aqueous solutions yields the α phase. Crystallization of glycine from water/ethanol solutions, however, results in the β phase. At ambient conditions, the order of stability is γ > α > β [29, 30]. The γ phase transforms to the α phase upon heating, at around 170 ◦ C, depending on the thermal history of the sample [29, 31]. High humidity drives a phase transition from the α to the γ phase [29, 32]. The β phase, metastable in dry air, rapidly transforms to α or γ phases in the presence of moisture at room temperature [29, 33]. Under pressure, the α phase is found to be stable up to 23 GPa, and the β phase is shown to undergo a phase transition to the δ phase at 0.76 GPa. Two independent experimental studies showed that the γ phase transforms to a new glycine phase called , starting from 1.9 or 2.74 GPa. The experimental studies by Boldyreva et al [34] resulted in a proposed structure for the phase based on the data obtained at 7.85 GPa. An alternative model for this phase is given by the work of Dawson et al [35] at 4.3 GPa. In our study, we model the phase on the basis of the latter structure. All known phases of glycine crystals are monoclinic, with the exception of the trigonal γ form. L -alanine is the smallest naturally occurring chiral amino acid and has been studied extensively. It is found − to crystallize in the zwitterion form CH3 CH(NH+ 3 )COO , in an orthorhombic structure with space group P21 21 21 with more than one molecule in its asymmetric unit cell (Z 0 = 4) [36–38]. Previous high-pressure, Raman-scattering experiments on L-alanine have indicated a structural phase transformation at about 2.3 GPa [39] and further synchrotron XRD experiments confirmed this transition and identified the new phase as tetragonal in symmetry [40]. Moreover, the same XRD experiments also suggested another transition, to a monoclinic phase at about 9 GPa. Both of these findings were questioned by Tumanov and co-workers as their Raman, x-ray, and optical microscopy did not confirm the previously suggested phase transitions. Instead, with increase in pressure, cell parameters a and c are found to happen to equate at around 1.5–2 GPa, without changing the cell symmetry. A new phase, proposed to belong to a solvate of L-alanine, was observed at around 3.5 GPa [41].
combined function q(r) = q(n(r), |∇n(r)|). Exploiting this feature, Rom´an-P´erez and Soler [24] have introduced an important technical device for efficiently computing vdW-DF integrals. In their approach the kernel is computed on a two-dimensional grid of fixed q values and any needed value in the integral is then obtained by interpolation: X 8(q, q0 , |r − r0 |) ≈ Pi (q) 8(qi , qj , |r − r0 |) Pj (q0 ), (5) ij
where the interpolating functions Pi (q) are defined such that Pi (qj ) = δij on a sufficiently dense grid of qs. The double spatial integrals appearing in the vdW-DF functional can then be computed as a series of convolutions that are most efficiently evaluated in reciprocal space, as for instance in ZZ Ecnl = 21 nn0 8(q, q0 , |r − r0 |) dr dr0 =
XX ∗ θ (G) 8(qi , qj , |G|) θj (G) 2 ij G i
=
XX ∗ u (G) θj (G), 2 j G j
(6)
where is the crystal cell volume, 8(qi , qj , |G|) and θi (G) are the Fourier transforms of the real space valued functions 8(qi , qj , |r|) and θi (r) = n(r)Pi (q(n(r), |∇n(r)|)), respecP tively, and the auxiliary function uj (G) = i θi (G)8(qi , qj , |G|) is further introduced. With these definitions, the general expression for the stress tensor, equation (4), can efficiently be computed as Z 1 nl nl nl E − n vc dr δαβ σc αβ = c Z X ∂Pj ∂q ∇α n∇β n − uj (r) n dr ∂q ∂|∇n| |∇n| j 1 X X Gα Gβ ∂8 qi , qj , |G| ∗ − θi (G) θj (G). 2 ij G |G| ∂|G| (7) The above formula has been implemented in the PWscf code of the Q UANTUM ESPRESSO distribution [25] that efficiently solves Kohn–Sham self-consistent equations in a plane-wave pseudopotential formalism. Numerical tests on simple benzene and methane crystals confirmed that the stress tensor computed in this way agrees perfectly (well within 0.01 GPa) with the numerical derivative, the residual discrepancy to be attributed to discretization error in the latter.
4. Computational details Ab initio calculations are performed in the framework of density functional theory as implemented in the Q UANTUM ESPRESSO distribution [25]. The exchange–correlation functionals used in our study include two popular GGA functionals and the fully non-local vdW-DF correlation functional combined with two flavors of GGA exchange proposed in the literature. The two GGA functionals considered are the standard Perdew–Burke–Ernzerhof (PBE) [42] functional and the Zhang–Yang revised PBE (revPBE) [43] functional. In the work introducing vdW-DF [14] the exchange part of the revPBE functional was selected to be combined with the fully non-local correlation term, based on the observation that more standard GGAs predict substantial binding in rare gas dimers from exchange alone, a feature absent for exact Hartree–Fock exchange and described correctly by revPBE
3. Background on glycine and L-alanine Glycine and alanine are among the simplest amino acids, yet their structural evolution under pressure displays a rich phenomenology that has been studied extensively in the literature and they constitute an ideal test set for vdW-aware functionals [26–41]. Glycine, H3 N+ CH2 COO− in its zwitterionic form, has three polymorphs under ambient conditions: α (Z 0 = 4), β (Z 0 = 2), and γ (Z 0 = 3) [26–28], the most stable 3
J. Phys.: Condens. Matter 24 (2012) 424209
R Sabatini et al
Table 1. Optimized cell parameters at zero pressure for glycine phases that are stable under ambient conditions. XC
−3
˚ a (A)
˚ b (A)
˚ c (A)
β (deg)
mol cell−1
˚ ) vol mol−1 (A
5.203 5.359 5.244 5.034 5.083
12.717 14.000 12.289 11.567 11.820
5.476 5.525 5.570 5.449 5.460
109.6 109.4 111.1 112.9 111.9
4 4 4 4 4
85.32 97.76 83.74 73.06 76.0
5.167 5.414 5.216 5.017 5.094
6.768 8.231 6.410 5.990 6.286
5.411 5.546 5.503 5.399 5.383
112.0 120.0 112.5 113.9 113.2
2 2 2 2 2
87.75 106.98 85.03 74.20 79.2
7.252 7.868 7.236 6.892 6.975
7.251 7.854 7.235 6.892 6.975
5.520 5.496 5.590 5.459 5.437
90 90a 90 90 90
3 3 3 3 3
83.81 97.94 84.46 74.84 76.9
α-glycine PBE revPBE vdW-DF (revPBE) vdW-DF-c09x Exp. [48] β-glycine PBE revPBE vdW-DF (revPBE) vdW-DF-c09x Exp. [29] γ -glycine PBE revPBE vdW-DF (revPBE) vdW-DF-c09x Exp. [49] a
Upon releasing the constraint on symmetry, the angle β becomes 89.9◦ .
and γ are given in table 1. Also listed are experimental values determined for these structures. Calculated volumes are found to be sensitive to the exchange–correlation functionals used. On average, PBE is found to overestimate the volumes by +10.7%, revPBE by +30.4%, and vdW-DF (with revPBE exchange) by +9.1%, and vdW-DF-c09x is found to underestimate the volume by 4.3%. The difference between the strongly overestimated volume obtained with revPBE and the improved volume estimate obtained with vdW-DF reveals the beneficial effect of including long-range correlations. We can see that the vdW-aware c09x functional consistently underestimates the experimental volume outside of the typical GGA behavior and gives the results closest to experimental ones. Concerning the individual lattice parameters and monoclinic angle, we can see that revPBE in particular provides a poor description of the crystal shape, but inclusion of non-local correlations in vdW-DF significantly improve the situation. From our results, it is evident that crystal densities and shapes are strongly dependent on the choice of XC functional. Next, we explore the performance of these functionals in reproducing the structure of the single molecules in the crystal. In table 2, we compare the bond lengths and torsion angle for α glycine at ambient pressure. In terms of bond lengths, all functionals describe the molecule reasonably well, within a 2% range from experimental values for bonds not involving hydrogen atoms, and 6% for bonds with hydrogen atoms. We can note that in general atom–atom distances are overestimated within PBE and that this tendency is enhanced by revPBE, consistently with the observed volume behavior. Including long-range correlations in vdW-DF does not improve the description of bond lengths, at variance with the observed effect on the volume. The effect of the XC functional becomes more apparent in the torsion angle where
exchange. A further refinement of revPBE exchange, named c09x, was subsequently developed by Cooper [44] to be used in conjunction with the non-local vdW-DF correlation. The resulting vdW-DF-c09x [44] functional was found to be superior to the original vdW-DF (revPBE) functional on the benchmark S22 [45] database of weakly bound molecules. Both the original vdW-DF (revPBE) [14] functional and its refinement vdW-DF-c09x [44] are considered here. The accuracy of the q-point grid involved in the Rom´an-P´erez and Soler interpolation [24] of the vdW-DF kernel has been carefully checked and the standard 20-point grid was confirmed to be adequate. Results obtained with a denser 30-point grid did not change significantly: less than ˚ −1 0.1 mRyd for the energy per molecule, less than 10−3 eV A −3 for the forces on atoms and less than 10 GPa for the stress tensor. Ultrasoft pseudopotentials from the PSLibrary project [46] are used without further modification. A kinetic energy cutoff of 80 Ryd and a charge density cutoff of 560 Ryd are used to achieve pressure convergence within less than 0.1 GPa. Monkhorst–Pack [47] k-point grids of 4 × 2 × 4, 3 × 3 × 3, 3 × 3 × 3, 2 × 3× 3, 3 × 3 × 3 are used for the α, β, γ , δ, and phases of glycine, respectively. For alanine, an 8 × 8 × 4 grid is found sufficient. With this setting, a tight convergence of less than 0.1 mRyd in total energy is achieved. The cell parameters and atomic positions are fully optimized. Vibrational zero-point motion and finite temperature effects are not included.
5. Results at ambient pressure 5.1. Glycine Our optimized lattice parameters a, b, c, and β calculated at ambient pressure for the three monoclinic stable phases α, β, 4
J. Phys.: Condens. Matter 24 (2012) 424209
R Sabatini et al
˚ and torsion angle (deg) for α glycine at zero pressure. The numbering of the atoms is given in the Table 2. Optimized bond lengths (A) graphical representation at the bottom. Bonds
PBE
revPBE
C1–O1 C1–O2 C2–N C1–C2 C2–H1 C2–H2 N–H3 N–H4 N–H5
1.26 1.27 1.49 1.53 1.10 1.09 1.05 1.07 1.04
1.27 1.27 1.50 1.55 1.10 1.09 1.04 1.07 1.04
vdW-DF (revPBE) 1.27 1.27 1.51 1.55 1.09 1.09 1.04 1.06 1.03
vdW-DF-c09x
Exp. [48]
1.27 1.27 1.48 1.53 1.10 1.10 1.06 1.07 1.05
1.25 1.25 1.48 1.52 1.05 1.04 1.03 1.09 1.09
Torsion N–C2–C1–O1 30.2
37.6
21.5
20.1
18.6
Table 3. Optimized cell parameters at zero pressure for the L -alanine phase stable under ambient conditions.
only vdW-including functionals show a very good agreement with the experiment. We observe the same general trend also for the other glycine phases, which are stable at ambient pressure. 5.2. Alanine At ambient pressure, alanine molecules, bearing a net electric dipole, are aligned in chains along the c axis, packed in the ab plane so as to minimize electrostatic energy. In addition, a network of H bonds exists, which are rigid along chains and rather flexible in the orthogonal plane. On this basis we can expect vdW interactions to be particularly important for the description of the lateral chain packing in the ab plane. We report the optimized a, b, and c lattice parameters at ambient pressure for L-alanine in table 3, together with the experimentally determined values [40, 41]. Like for the case of glycine, we see that the calculated volume is sensitive to the choice of XC functional. While GGA functionals overestimate the volume, PBE by 8.1% and revPBE by 19%, the vdW-aware ones are able to describe the system with higher precision, vdW-DF overestimating the cell volume by 6% and vdW-DF-c09x underestimating it by 5%. The apparently better results for the crystal volume obtained with GGA functionals with respect to the glycine case are somewhat misleading, as the cell shape is very badly described by these functionals. In fact, both GGA functionals strongly overestimate the a lattice parameter while underestimating b (and for revPBE, also c), confirming the importance of the inclusion of dispersion forces for the aforementioned lateral chain packing. Indeed, the two vdW-aware functionals provide a more balanced, though not perfect, description of the crystal shape.
XC
˚ a (A)
˚ b (A)
˚ c (A)
˚ 3) Volume (A
PBE revPBE vdW-DF (revPBE) vdW-DF-c09x Exp. [40] Exp. [41]
6.831 7.241 6.119 5.886 6.102 6.049
11.615 11.939 12.637 12.055 12.339 12.345
5.871 5.525 5.910 5.781 5.780 5.780
465.91 512.64 456.99 410.17 425.13 431.6
As in the case of glycine, all functionals are able to describe with reasonable accuracy the structure of the single molecule, as we can see in table 4, where we report the characteristic bond lengths for a single alanine molecule. All bonds are correctly described within a 2% range from experimental values, including C–H and N–H ones, and also the torsion angle does not vary significantly between the different XC functionals used and agrees well with experiment.
6. Results at high pressure 6.1. Glycine Addressing the relative stability of glycine polymorphs is difficult due to the very small energy differences between the phases (less than 1 kcal mol−1 ). Calculated differences are also sensitive to the exchange–correlation functionals used. In figure 1 we report the enthalpy per molecule as a function of pressure, referenced with respect to the α phase, for all XC functionals employed. At zero pressure, vdW-DF is the only functional that predicts the γ > α > β stability order correctly, while 5
J. Phys.: Condens. Matter 24 (2012) 424209
R Sabatini et al
Figure 1. Enthalpy per molecule referenced with respect to the α phase as a function of pressure for all known phases of glycine up to 10 GPa, for all XC functionals considered. Calculations are performed with target pressure −0.5, 0, and up to 10 GPa in increments of 1 GPa. Additional calculations were performed for high-pressure phases at around ambient pressure for better convergence.
phase is well separated from the α and γ phases as expected for a metastable state; however this is not observed for the PBE functional. We see that revPBE displays a highly erratic stability order both at zero pressure and at higher pressures. The sharp increase in enthalpy for the β phase at 2.0 GPa with the revPBE functional is due to the reorganization of the hydrogen bond network taking place from the 1.0 GPa configuration to the 3.0 GPa one. At 1.0 GPa, the two molecular layers in the unit cell are stacked on top of each other, so linear hydrogen bonds between molecules of the two layers are formed, normal to the layer plane. At 3.0 GPa, the layer stacking changes such that a molecule of the top layer sits in the void between two molecules of the lower layer, forming a trigonal hydrogen bond network. The 2.0 GPa configuration is the transition point between these two, resulting in higher enthalpy. The on-top stacking at 1.0 GPa is only observed with the revPBE functional. Due to this and other anomalies in its enthalpy plot, we deduce that the revPBE functional is not reliable and it will not be discussed further. In agreement with the experimental observations, the α phase is stable up to 10 GPa for PBE, vdW-DF, vdW-DF-c09x, and after 1.1 GPa also for revPBE. The phase transition from β to δ occurs experimentally at 0.76 GPa. The δ phase is calculated to be thermodynamically more stable than the β phase after around 4.0 GPa, 1.5 GPa, and 0.7 GPa for PBE, vdW-DF, and vdW-DF-c09x, respectively. Considering the precision of our calculations, the estimates from vdW-DF and vdW-DF-c09x are equally satisfactory, while PBE overestimates the transition point. The phase transition of the γ polymorph is experimentally more complicated, as single crystals already undergo a phase transition to a polycrystalline phase at 1.9 GPa, but
˚ and torsion angle (deg) for Table 4. Optimized bond lengths (A) L -alanine at zero pressure. The numbering of the atoms is given in the graphical representation at the bottom. Bonds
PBE
revPBE
C1–O1 C1–O2 C1–C2 C2–N C2–C3 N–H1 N–H2 N–H3 C2–H4 C3–H5 C3–H6 C3–H7
1.26 1.28 1.54 1.49 1.53 1.04 1.05 1.06 1.10 1.09 1.09 1.10
1.26 1.28 1.55 1.51 1.53 1.04 1.04 1.05 1.10 1.10 1.09 1.10
vdW-DF 1.26 1.28 1.56 1.52 1.54 1.04 1.04 1.05 1.10 1.10 1.10 1.10
c09x 1.26 1.28 1.53 1.48 1.52 1.05 1.05 1.07 1.10 1.10 1.10 1.10
Exp. [38] 1.24 1.26 1.53 1.49 1.52 1.03 1.03 1.05 1.09 1.08 1.08 1.08
Torsion N–C2–C1–O1 17.0
19.5
20.0
17.2
18.7
vdW-DF-c09x also shows very similar enthalpy for the γ and α phases. Since the energy differences are extremely small, we can regard these two functionals as performing equally well at zero pressure. For both functionals, the β 6
J. Phys.: Condens. Matter 24 (2012) 424209
R Sabatini et al
Table 5. L-alanine bulk modulus B0 (GPa) and its pressure derivative B00 . Experimental values from [41] were extracted by fitting experimental data. XC
B0
B00
11.2 7.4 17.8 16.4 31.5 ± 1.4 15.1 ± 1.2
7.32 7.04 4.32 5.25 4.4 ± 0.4 5.1 ± 0.9
L -alanine
PBE revPBE vdW-DF (revPBE) vdW-DF-c09x Exp. [40] Exp. [41]
less important role, the functional approaches experimental ones. As for the non-local functionals, vdW-DF-c09x underestimates the experimental volume at all pressures, while vdW-DF overestimates it at low pressure and falls midway within experiment uncertainty at higher pressures. PBE appears to give reasonably good results for the volume, close to vdW-DF ones at all pressures. A least-squares fit of the data with the Birch–Murnaghan equation of state gives the equilibrium bulk modulus and its pressure derivative for each functional, reported in table 5. As is already evident from figure 2, compressibility from vdW-aware functionals agrees very well with the value that can be deduced from the experimental data in [41], while the agreement is less good with [40]. In any case, both GGA functionals severely overestimate the compressibility of the system. As already noted when discussing ambient pressure results, the small difference between PBE and vdW-DF volumes in L-alanine is misleading. In figure 3 we report the evolution of the a, b, and c cell parameters as a function of pressure, compared with experimental values. As could be expected, the c axis corresponds to the stiffer direction in the crystal and all functionals give a very reasonable description of its pressure dependence. It can be seen however that PBE and revPBE values for the a and b lattice parameters are completely wrong at low pressure and approach the experimental values only at higher pressure, where weak van der Walls interactions become less important. The only functionals able to correctly describe the evolution of the cell geometry in the whole pressure range are the vdW-aware ones. Experimental data from [40], reported with full symbols, show the signature of the experimentally observed transition, from orthorhombic to tetragonal, at 2.3 GPa where a and c become equal. Experimental data from [41], reported with open symbols, indicate instead a chance crossing between these two quantities, which interchange at higher pressure. Our calculations with vdW-aware functionals support this interpretation, with the original vdW-DF giving the best results.
Figure 2. Cell volume as a function of pressure for L-alanine for all XC functionals considered. Experimental data from [40, 41] are also shown.
the identification of the phase could only be obtained at 4.3 GPa. Our calculations predict that the phase becomes thermodynamically more stable than the γ one at around 4.0 GPa, 1.8 GPa, and 1.2 GPa for PBE, vdW-DF, and vdW-DF-c09x, respectively. As in the case of the β → δ transition, we can say that vdW-DF and vdW-DF-c09x perform equally well in predicting the single-crystal transition of the γ form, while PBE overestimates the transition pressure. The stability order remains the same for all XC functionals at high pressures: α > > δ > β > γ . Our calculations show that energetics and structural properties are sensitive to the choice of XC functionals in the case of glycine crystals at ambient and higher pressures. Although for a more definite discussion of the relative stability of the different polymorphs, inclusion of vibrational zero-point motion and finite temperature effects would be desirable, exchange–correlation functionals including van der Waals interactions result in better estimates for the stability order and crystal density than the functionals lacking this contribution. For the stability order at ambient pressure and transition pressures, the two vdW-aware functionals perform equally well and better than PBE and revPBE, suggesting that using vdW-aware XC functionals can greatly improve the energetics and structural properties of molecular crystals obtained by ab initio methods. 6.2. Alanine As a final case, we address the ability of the new vdW-aware functionals to correctly describe the equation of state of L -alanine and to shed light on the ongoing controversy about a debated pressure-driven phase transition [39–41]. In figure 2, we report the cell volume as a function of pressure for the four functionals considered between 0 and 10 GPa, a range where three structural phase transitions have been reported. The two available experimental volume sets [40, 41] agree at low pressure but significantly differ at higher pressure. At low pressure, revPBE fails to describe the crystal correctly, grossly overestimating the volume, while at higher pressure, where long-range correlations likely play a
7. Conclusion Our calculations show that consideration of vdW-aware XC functionals can greatly improve the description of the energetics and structural properties of molecular crystals from 7
J. Phys.: Condens. Matter 24 (2012) 424209
R Sabatini et al
Figure 3. Crystal cell parameters a (red lines and circles), b (blue lines and triangles), and c (green lines and squares) of L-alanine as a function of pressure. Lines connect the calculated points. Experimental data from [40, 41] are reported with full and open symbols, respectively.
first principles. In the case of glycine, functionals including van der Waals interactions result in better estimates for the stability order of the different polymorphs and crystal density than the functionals lacking this contribution. Both flavors of vdW-aware functionals perform equally well and significantly better than PBE and revPBE. In the case of L-alanine, only vdW-aware functionals deliver individual cell parameters that evolve correctly with pressure, while the GGA results are qualitatively wrong. The equation of state agrees reasonably well with a recent experimentally obtained one, but less so with another. The bulk modulus is much improved compared with GGA results. In conclusion, while our results confirm that vdW-DF functionals allow a significant step forward in the firstprinciples description of soft matter, they also show that significant room for improvement still remains.
[3]
[4] [5]
[6] [7]
[8] [9] [10]
Acknowledgments SdG acknowledges support by the MIUR under the PRIN 2008 program. TT acknowledges support by the Department of Energy, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, Grant No. DE-FG0208ER46491. Computational facilities were provided by CINECA and SISSA.
[11] [12] [13] [14]
References
[15]
[1] Kleis J, Lundqvist B I, Langreth D C and Schr¨oder E 2007 Phys. Rev. B 76 100201 [2] Thonhauser T, Puzder A and Langreth D C 2006 J. Chem. Phys. 124 164106
[16] [17]
8
Li S, Cooper V R, Thonhauser T, Puzder A and Langreth D C 2008 J. Phys. Chem. A 112 9031 Hooper J, Cooper V R, Thonhauser T, Romero N A, Zerilli F and Langreth D C 2008 ChemPhysChem 9 891 Sony P, Puschnig P, Nabok D and Ambrosch-Draxl C 2007 Phys. Rev. Lett. 99 176401 Mura M, Gulans A, Thonhauser T and Kantorovich L 2010 Phys. Chem. Chem. Phys. 12 4759 Kolb B and Thonhauser T 2011 Phys. Rev. B 84 045116 Cooper V R, Thonhauser T, Puzder A, Schr¨oder E, Lundqvist B I and Langreth D C 2008 J. Am. Chem. Soc. 130 1304 Cooper V R, Thonhauser T and Langreth D C 2008 J. Chem. Phys. 128 204102 Li S, Cooper V R, Thonhauser T, Lundqvist B I and Langreth D C 2009 J. Phys. Chem. B 113 11166 Li Q, Kolb B, Rom´an-P´erez G, Soler J M, Yndurain F, Kong L, Langreth D C and Thonhauser T 2011 Phys. Rev. B 84 153103 Yao Y, Nijem N, Li J, Chabal Y J, Langreth D C and Thonhauser T 2012 Phys. Rev. B 85 064302 French R H et al 2010 Rev. Mod. Phys. 82 1887 Tang K T and Toennies J P 1984 J. Chem. Phys. 80 3726 Grimme S 2004 J. Comput. Chem. 25 1463 Grimme S, Antony J, Ehrlich S and Krieg H 2010 J. Chem. Phys. 132 154104 Steinmann S N and Corminboeuf C 2010 J. Chem. Theory Comput. 6 1990 Wu Q and Yang W 2002 J. Chem. Phys. 116 515 Neumann M A and Perrin M-A 2005 J. Phys. Chem. B 109 15531 Dion M, Rydberg H, Schr¨oder E, Langreth D C and Lundqvist B I 2004 Phys. Rev. Lett. 92 246401 Kohn W, Yigal M and Makarov D E 1998 Phys. Rev. Lett. 80 4153 Mahan G D and Subbaswamy K R 1990 Local Density Theory of Polarizability (New York: Plenum) chapter 4.3 Bell R J and Zucker I J 1976 Rare Gas Solids (New York: Academic) chapter 2
J. Phys.: Condens. Matter 24 (2012) 424209
R Sabatini et al
[35] Dawson A, Allan D R, Belmonte S A, Clark S J, David W I F, McGregor P A, Parsons S, Pulham C R and Sawyer L 2005 Cryst. Growth Des. 5 1415 [36] Simpson H J Jr and Marsh R E 1966 Acta Crystallogr. 20 550–5 [37] Dunitz J D and Ryan R R 1966 Acta Crystallogr. 21 617 [38] Lehmann M S, Koetzle T F and Hamilton W C 1972 J. Am. Chem. Soc. 94 2657 [39] Teixeira A M R, Freire P T C, Moreno A J D, Sasaki J M, Ayala A P, Mendes Filho J and Melo F E A 2000 Solid State Commun. 116 405i [40] Olsen J S, Gerward L, Freire P T C, Mendes Filho J, Melo F E A and Souza Filho A G 2008 J. Phys. Chem. Solids 69 1641 [41] Tumanov N A, Boldyreva E V, Kolesov B A, Kurnosov A V and Quesada Cabrera R 2010 Acta Crystallogr. B 66 458 [42] Perdew J P, Burke K and Ernzerhof M 1996 Phys. Rev. Lett. 77 3865 [43] Zhang Y and Yang W 1998 Phys. Rev. Lett. 80 890 [44] Cooper V R 2010 Phys. Rev. B 81 161104 ˇ ˇ y J and Hozba P 2006 Phys. Chem. [45] Jureˇcka P, Sponer J, Cern´ Chem. Phys. 8 1985 [46] Adllan A A and Dal Corso A 2011 J. Phys.: Condens. Matter 23 425501 [47] Monkhorst H J and Pack J D 1976 Phys. Rev. B 13 5188 [48] Legros J P and Kvick A 1980 Acta Crystallogr. B 36 3052 ˚ Canning W M, Koetzle T F and [49] Kvick A, Williams G J B 1980 Acta Crystallogr. B 36 115
[18] Langreth D C and Perdew J P 1977 Solid State Commun. 17 1425 Langreth D C and Perdew J P 1977 Phys. Rev. B 15 2884 [19] Langreth D C et al 2009 J. Phys.: Condens. Matter 21 084203 [20] Bil A, Kolb B, Atkinson R, Pettifor D G, Thonhauser T and Kolmogorov A N 2011 Phys. Rev. B 83 224103 [21] Thonhauser T, Cooper V R, Li S, Puzder A, Hyldgaard P and Langreth D C 2007 Phys. Rev. B 76 125112 [22] White J A and Bird D M 1994 Phys. Rev. B 50 4954 [23] Nielsen O H and Martin R M 1983 Phys. Rev. Lett. 50 697 Nielsen O H and Martin R M 1985 Phys. Rev. B 32 3780 Nielsen O H and Martin R M 1988 Phys. Rev. B 37 10905 [24] Rom´an-P´erez G and Soler J M 2009 Phys. Rev. Lett. 103 096102 [25] Giannozzi P et al 2009 J. Phys.: Condens. Matter 21 395502 [26] Iitaka Y 1961 Acta Crystallogr. 14 1 [27] Iitaka Y 1960 Acta Crystallogr. 13 35 ˚ 1972 Acta Crystallogr. B 28 1827 [28] J¨onsson P-G and Kvick A [29] Pervolich G L, Hansen L K and Bauer-Brandl A 2001 J. Therm. Anal. Calorim. 66 699 [30] Chongprasert S, Knopp S A and Nail S L 2001 J. Pharm. Sci. 90 1720 [31] Park K, Evans J M B and Myerson A S 2003 Cryst. Growth Des. 3 991 [32] Sakai H, Hosogai H, Kawakita T, Onuma K and Tsukamoto K 1992 J. Cryst. Growth 116 421 [33] Ferrari E S, Davey R J, Cross W I, Gillon A L and Towler C S 2003 Cryst. Growth Des. 3 53 [34] Boldyreva E V, Ivashevskaya S N, Sowa H, Ahsbahs H and Weber H P 2004 Dokl. Phys. Chem. 396 111
9
View publication stats
