Possibility of Coherent Delocalized Nuclear ... - ACS Publications

pubs.acs.org/JPCL

Possibility of Coherent Delocalized Nuclear Quantum States of Protons in Li2NH ~ a† and Daniel Sebastiani*,‡ Guillermo A. Luduen †

Max-Planck-Institut f€ ur Polymerforschung, Ackermannweg 10, 55128 Mainz, Germany, and Physics Department, Freie Universit€ at Berlin, Arnimallee 14, 141965 Berlin, Germany

‡

ABSTRACT We analyze the possibility of quantum delocalization in lithium imide (Li2NH) in the condensed phase using ab initio path-integral molecular dynamics. Our results provide evidence that the effective potential felt by the protons in the material has a toroidal shape. The virtually flat potential may lead to quantum delocalization of the NH protons over the torus. SECTION Molecular Structure, Quantum Chemistry, General Theory

U

nder ambient conditions, most atoms can normally be approximated as point particles, and this approximation yields a good level of accuracy for many properties. However, in the case of hydrogen, nuclear quantum effects can be observable due to its low atomic mass. In some cases, hydrogen has been reported to exhibit quantum properties, such as tunneling and quantum delocalization.1-5 These effects might also play a role in industrially relevant situations, such as hydrogen storage.6 Besides this special aspect, understanding the nature of quantum effects of hydrogen in bulk materials is a matter of high fundamental and general relevance. As of today, no production-ready material has been found to meet the required per-weight capacity of hydrogen storage.6-8 In this work, we focus on the particular case of lithium imide (Li2NH), which is a promising hydrogen storage system due to its low weight. The lithium imide/amide system can undergo a reversible hydrogenation with high H storage densities (6.5%) via the reaction Li2NH þ H2 S LiNH2 þ LiH. A lot of research, both computational and experimental, has been done on this system, particularly to understand its structure and phase transition.9,10 In this context, the elementary steps of the reaction have recently attracted considerable interest.11 In particular, the possibility of quantum delocalization of the protons over several possible sites has been reported, as a consequence of the relatively flat potential energy landscape in the highly symmetric structure.12 In our previous publication,10 we observed that the effect of nuclear quantum delocalization for the hydrogen atoms is considerably reduced compared to a full delocalization of each proton on the six octahedral sites around its nitrogen atom. In the present work, we analyze the potential energy landscape of the protons in the condensed phase by means of ab initio path-integral molecular dynamics simulations (PIMD).13,14 We show that this surface is surprisingly flat within a torusshaped volume, which might lead to coherent delocalization at low temperatures. We model the initial configuration of the Li2NH protons as delocalized particles, which are distributed over six

r 2010 American Chemical Society

Figure 1. Left: Li2NH converged structure. The Li lattice is disordered, and the N lattice is face-centered cubic (FCC). The N-H bonds precess around a preferred N-N axis while keeping 30 respect to that axis. Right: The quantum particle density of the Li2NH system, as obtained from PIMD simulations at 300 K. Li atoms are colored green, N atoms blue, and H atoms gray/ white.

crystallographically symmetric sites, using the setup presented in our previous work.10 Such a proton configuration is plausible because of the low potential barrier between these sites.12 As reported previously,10 the equilibrated structure at 300 K features a somewhat disordered Li lattice. The NH bonds form angles of 30 from a preferred N-N axis15 (Figure 1, left), and the protons remain in a continuous precession around the preferred axis. This structure is observed both in PIMD (P > 1) and Car-Parrinello molecular dynamics (CPMD) (P=1) simulations at ambient temperature.10,15 This continuous motion suggests a low potential energy surface along the proton trajectory. In the case of protons, such a flat potential can lead to important quantum effects. In the following, we study the Received Date: September 1, 2010 Accepted Date: October 15, 2010 Published on Web Date: October 25, 2010

3214

DOI: 10.1021/jz1012388 |J. Phys. Chem. Lett. 2010, 1, 3214–3218

pubs.acs.org/JPCL

Figure 3. Two-dimensional correlation functions of pairs of azimuth angles φ around the precession axis between two NH bonds in the Li2NH bulk. Two of the six possible NH bond correlation plots are shown for illustration.

Figure 2. Dependence of the averaged Kohn-Sham energy EhKS(φ) on the azimuth angle φ (sketched in Figure 1 and Figure 5). The error bars show the statistical standard deviation for a given angular range.

possibility of quantum delocalization by looking in the PIMD particle densities. The quantum particle density resulting from our PIMD simulations is presented in Figure 1, right. The figure shows the superposition of the beads from 5 ps of a PIMD trajectory (32 beads, 300 K). Each individual quantum proton density forms a torus of radius r ≈ 0.5 Å with an homogeneous proton distribution. However, from this distribution alone, it cannot be derived whether the proton is being delocalized or moving as a localized particle in a torus-shaped potential well. In order to discriminate the two scenarios, we have computed the potential landscape along the precession angle of the protons (φ, sketched in Figure 1). The resulting potential EhKS(φ) was obtained by averaging the energies EKS of about 50 000 individual configurations over all degrees of freedom, for each value of φ. Formally, this can be expressed as P EKS ðφÞ ¼

t , p, i

Figure 4. Quantum particle density obtained from PIMD calculations for (a) CH3OH molecule, (b,c) CH3OH view from the CO axis, and (d) Li2NH molecule. Colors: H atoms in black, Li in green, N in blue, O in red, and C in cyan.

NH bonds. In order to check this assumption, we calculated the two-dimensional correlation function of both azimuthal NH bond angles (φi, φj). In Figure 3, we present two plots of the correlation function for the φ angles of two NH bonds. The plots show that this joint probability is scattered over all possible angle pairs. This indicates that the orientations of different NH bonds are essentially uncorrelated. Hence, we conclude that the individual NH protons precess in an uncorrelated way on a mainly flat potential energy landscape. As a mean for comparison with simpler systems, additional PIMD simulations of the isolated Li2NH and CH3OH molecules at 300 K were performed. In Figure 4, the obtained quantum particle densities are shown for both molecules. The NH bond in the isolated Li2NH molecule shows a different behavior to the bulk Li2NH, as the molecule converges to an essentially flat and rigid configuration. The final proton density is localized for all particles in the molecule (see Figure 4d). The case of CH3OH is particularly interesting since the OH proton can freely rotate around the C-O axis, and therefore it experiences a situation similar to the NH protons in bulk Li2NH. As shown in Figure 4a, the particle density of the OH proton resembles a torus of radius ∼0.9 Å. In this respect, the proton density distribution looks very similar to the one observed in bulk Li2NH. Conversely, the protons forming CH3 are localized due to the repulsion exerted between them. There is, however, a small but noticeable correlation between the density of the OH proton and the orientation of the CH3 protons (Figure 4b,c), due to the proton-proton repulsion.

t, p EKS δðφt, p, i - φÞ

P

t , p, i

δðφt, p, i - φÞ

ð1Þ

where φt,p,i is the azimuth angle of the ith individual NH bond in the pth replica at time t, and Et,p KS is the energy of the configuration in replica p at time t. Using this procedure, we effectively average the orientation of the other NH bonds in the system. In Figure 2, the resulting potential EhKS(φ) is shown. It varies by less than 0.8 mHa as a function of φ, which is on the order of the thermal energy kBTat T = 300 K. We have also computed the statistical standard deviation of the energies Et,p KS at fixed angle φ, shown as error bars in Figure 2. We observe that this standard deviation is virtually constant and independent of the value of the angle. A completely flat potential would be perfectly contained within these statistical error bars. Together with the very small variations of the averaged potential EhKS(φ), we thus conclude that the effective potential for hydrogen motion along the torus is most likely very flat. In the procedure followed to obtain the potential energy (eq 1), we have assumed that the energy of an individual NH bond azimuthal angle is uncorrelated from the orientation of the rest of the bonds, i.e., there is no concerted rotation of the


3215


pubs.acs.org/JPCL

OH proton, which is not subjected to a confining potential in that angle. When comparing the arc half-widths in the tori for both materials, the angular distribution of the protons in bulk Li2NH turned out to be 26% wider than that of the OH proton in CH3OH. This narrower distribution in CH3OH may be a consequence of the previously mentioned correlation of the OH with the CH3 protons. More importantly, the much wider distribution in Li2NH indicates that the bulk Li2NH protons are exposed to a considerably more flat potential landscape than those in CH3OH. The flatness of the potential in which Li2NH protons are found is important evidence supporting proton quantum delocalization in that material. In this work, we have looked at nuclear quantum effects in the condensed phase by means of first-principles path integral molecular dynamics simulations. Specifically, we have modeled the quantum delocalization of protons in Li2NH, an inorganic compound that is highly relevant as a potential hydrogen storage material. Our simulations provide evidence for a virtually flat effective potential for the motion of the protons within a toroidal region corresponding to the precession of the NH bond around the crystal c axis. The flatness of the effective potential with respect to the azimuthal angle might lead to a coherent quantum delocalization at low temperatures. While there is most probably no such coherent delocalization at ambient temperatures, we are currently performing calculations of the nuclear momentum distribution n(k) at low temperature. This momentum distribution will show whether the protons exhibit only zero point broadening or a coherent delocalization.

Figure 5. Distribution of arc distances to the centroid in the azimuth angle φ (see inset) for the protons in bulk Li2NH, in the CH bonds of CH3OH, and in the OH bond of CH3OH.

The proton distributions in CH3OH (Figure 4) have the expected shape. They illustrate qualitatively that even a small interaction (such as between CH3 and OH protons) can lead to a visible signature in the proton density distribution, as seen in the anticorrelation of the OH proton density with the CH3 protons in Figure 4c. Turning back to the condensed-phase Li2NH system, the remaining question is whether the motion of the protons in imaginary time is a fully classical precession (with an additional zero-point broadening) or a coherent delocalized quantum state. Such a coherent state would form a wavepacket at high temperatures, but it would turn into a plane wave state close to zero temperature. While this question is nontrivial to answer, we will shed some light on the matter by studying the delocalizaton of the ring polymer in PIMD imaginary time, with the aim of finding evidence supporting one of these two possibilities. In order to quantify the proton ring polymer (de)localization with respect to the azimuthal angle, we compute the arc distance of bead p to its centroid as

COMPUTATIONAL DETAILS The quantum treatment of the nuclei is provided by the path integral technique (PIMD)13,14,16 running on top of CPMD simulations.17-19 The PIMD technique provides a way to compute the quantum properties of a material from a statistical ensemble of quasi-classical trajectories. These trajectories are generated by quasi-classical nuclei, called beads, governed by the classical Hamiltonian of the nuclei. The ensemble of coordinates for a given (imaginary) time is called a “replica”. Each nucleus in a replica (i.e., each bead) is connected by an additional harmonic force to the same nucleus in another replica, conforming a “ring polymer”. All replicas are simultaneously propagated in imaginary time, eventually forming a set of path integral trajectories. Then, the ensemble averaged properties of the bead distribution for a given nucleus correspond to that of the stationary quantum problem, when the number of replicas P f ¥. We have used Trotter decompositions into P = 4,6,16,32, and 64 replicas. It turns out that the configurational energies are converged for P = 16, which led us to base all further calculations on this number of beads. The calculation of the forces exerted over the nuclei is based on density functional theory (DFT),20 in the generalized gradient approximation using the BLYP exchange-correlation functional.21,22 The orbitals were expanded in plane waves, using Goedecker pseudopotentials for N and Li.23 The planewave cutoff was set to 120 Ry, which turned out to be

Δarcφ ðpÞ ðtÞ ¼ rðpÞ ðtÞ sinðθðpÞ ðtÞÞφðpÞ ðtÞ - rC ðtÞ sinðθC ðtÞÞφC ðtÞ

ð2Þ

)

)

where C indicates the centroid of the ring polymer, and r(p), θ(p), and φ(p) are the coordinates of the NH bond in a given replica p. In Figure 5, we show the distribution of arc distances during the whole PIMD trajectory for the protons in bulk Li2NH and CH3OH. We extract from these distributions the arc half-width of the proton ring polymers, which we define as the standard deviation of the arc distances distribution, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ÆΔarc2 æ - ÆΔarcæ2 . We want to remark that the width of the ring polymers at high temperature is dominated by the strong harmonic forces between beads, which are proportional to mPT2 R(p) - R(pþ1) . The effect of the external potential in the ring polymer width is therefore expected to be subtle. The arc half-widths obtained are 0.117 Å for the protons in bulk Li2NH (16 beads, 300 K), 0.093 Å for the OH of CH3OH (32 beads, 300 K), and 0.087 Å for the (localized) organic protons in CH3OH (32 beads, 300 K). It is remarkable that the ring polymer angular distribution of the organic protons in CH3OH is only 6.5% narrower than that corresponding to the mobile


3216


pubs.acs.org/JPCL

Table 1. Proton Bead Distribution Converged Width W (t . 0) in Lithium Imide at 300 K for Different Numbers of Beads (in Å)a beads W (t . 0) a

6 0.11

16 0.10

32 0.11

64 0.12

The statistical error amounts to (0.1 Å.

occurs within 0.3 ps. For t > 0.3 ps, the width W (t) remains essentially constant, and its final average value is independent of the number of beads (see Table 1). This instantaneous width of the path integral ring polymer does not have a direct physical meaning; only the ensemble average of all beads over a converged trajectory has. Nevertheless, we believe that the average width in imaginary time still provides a qualitiative measure for discriminating between different degrees of delocalization, especially because of the fast convergence of the bead distribution and the stable width during the MD run.

Figure 6. Bead distribution width W (t) for the initially delocalized configuration of Li2NH, during the first picoseconds of a simulation (P = 32, T = 300 K). The inset is a zoom to the first steps, showing how the initially highly delocalized proton beads collapse rapidly to a more localized distribution.

AUTHOR INFORMATION Corresponding Author: *To whom correspondence should be addressed. E-mail: daniel. [email protected].

necessary to achieve convergence in the frequencies of the vibrational normal modes of molecular Li2NH. The energy hypersurface generated by DFT was used for the PIMD simulations via the CPMD scheme,17,18 thus applying realistic atomic motion at ambient temperatures. The time step was set to 0.085 fs to ensure an overall stability. Our unit cell for Li2NH has dimensions of 5.007 5.007 5.007 Å3 and contains 16 independent atoms (i.e., four Li2NH molecules). The system was equilibrated during 5 ps of PIMD simulation, by stepwise increasing temperatures every 1 ps, until the target temperature of 300 K was reached. The subsequent production runs were between 4 and 9 ps long. For the study of isolated molecules, a single molecule was placed in a unit cell of 15 15 15 Å3, and the remaining setup was left unchanged with respect to the bulk Li2NH simulations. For the calculation of the particle densities, the translational and the rotational degrees of freedom of the molecules as a whole were eliminated a posteriori from the generated trajectories. For the calculation of EhKS(φ), a full wave function optimization was performed on each snapshot from the PIMD trajectories, because the Car-Parrinello orbitals are energetically slightly above the Born-Oppenheimer surface. In order to test the convergence in the dynamics of the bead distribution width, we calculated the width W (t) of a particle's ring polymer at simulation time t, vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uP u P u jjRðpÞ ðtÞ - RC ðtÞjj2 tp ¼ 1 ð3Þ W ðtÞ ¼ P

ACKNOWLEDGMENT This work was financially supported by the

Deutsche Forschungsgemeinschaft (DFG) under Grants Se 1008/5-1 and Se 1008/6-1. The authors are grateful to A. Alavi, A. B. Poma, and H. W. Spiess for valuable discussions and continuous support.

REFERENCES (1) (2)

(3)

(4)

(5)

(6) (7)

(8)

where P is the number of beads, and R(p)(t) and RC(t) are the position of the pth bead and the centroid at time t, respectively. The time evolution of the width W (t) during the initial stage of the simulation is presented in Figure 6. The convergence of the (initially delocalized) distribution of proton beads


(9)

(10)

3217

Bakker, H. J.; Nienhuys, H.-K. Delocalization of Protons in Liquid Water. Science 2002, 297, 587–590. Uras-Aytemiz, N.; Devlin, J.; Sadlej, J.; Buch, V. HCl Solvation in Methanol Clusters and Nanoparticles: Evidence for ProtonWires. Chem. Phys. Lett. 2006, 422, 179–183. Li, X.-Z.; Probert, M. I. J.; Alavi, A.; Michaelides, A. Quantum Nature of the Proton in Water-Hydroxyl Overlayers on Metal Surfaces. Phys. Rev. Lett. 2010, 104, 066102. Ranea, V. A.; Michaelides, A.; Ramírez, R.; de Andres, P. L.; Verg es, J. A.; King, D. A. Water Dimer Diffusion on Pd{111} Assisted by an H-Bond Donor-Acceptor Tunneling Exchange. Phys. Rev. Lett. 2004, 92, 136104. Kumagai, T.; Kaizu, M.; Hatta, S.; Okuyama, H.; Aruga, T.; Hamada, I.; Morikawa, Y. Direct Observation of HydrogenBond Exchange within a Single Water Dimer. Phys. Rev. Lett. 2008, 100, 166101. Schlapbach, L.; Zuttel, A. Hydrogen-Storage Materials for Mobile Applications. Nature 2001, 414, 353–358. Gregory, D. H. Lithium Nitrides, Imides and Amides As Lightweight, Reversible Hydrogen Stores. J. Mater. Chem. 2008, 18, 2321–2330. Klerke, A.; Christensen, C. H.; Norskov, J. K.; Vegge, T. Ammonia for Hydrogen Storage: Challenges and Opportunities. J. Mater. Chem. 2008, 18, 2304–2310. Ichikawa, T.; Isobe, S. The Structural Properties of Amides and Imides As Hydrogen Storage Materials. Z. Kristallogr. 2008, 223, 660–665. Ludue~ na, G. A.; Wegner, M.; Bjålie, L.; Sebastiani, D. Local Disorder in Hydrogen Storage Compounds: The Case of Lithium Amide/Imide. ChemPhysChem 2010, 11, 2353–2360.


pubs.acs.org/JPCL

(11)

(12)

(13)

(14) (15)

(16)

(17) (18) (19)

(20)

(21)

(22)

(23)

Miceli, G.; Cucinotta, C. S.; Bernasconi, M.; Parrinello, M. First Principles Study of the LiNH2/Li2NH Transformation. J. Phys. Chem. C 2010, 114, 15174–15183. Zhang, C.; Dyer, M.; Alavi, A. Quantum Delocalization of Hydrogen in the Li2NH Crystal. J. Phys. Chem. B 2005, 109, 22089–22091. Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. Efficient Molecular Dynamics and Hybrid Monte Carlo Algorithms for Path Integrals. J. Chem. Phys. 1993, 99, 2796– 2808. Ceperley, D. Path Integrals in the Theory of Condensed Helium. Rev. Mod. Phys. 1995, 67, 279–355. Ara ujo, C. M.; Blomqvist, A.; Scheicher, R. H.; Chen, P.; Ahuja, R. Superionicity in the Hydrogen Storage Material Li2NH: Molecular Dynamics Simulations. Phys. Rev. B 2009, 79, 172101. Tuckerman, M.; Marx, D; Klein, M.; Parrinello, M. Efficient and General Algorithms for Path Integral Car-Parrinello Molecular Dynamics. J. Chem. Phys. 1996, 104, 5579–5588. Car, R.; Parrinello, M. A Combined Approach to DFT and Molecular Dynamics. Phys. Rev. Lett. 1985, 55, 2471. Hutter, J. Computer code CPMD, version 3.12, Copyright IBM Corp. and MPI-FKF Stuttgart 1990-2008, www.cpmd.org. Iftimie, R.; Minary, P.; Tuckerman, M. Ab Initio Molecular Dynamics: Concepts, Recent Developments, and Future Trends. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6654–6659. Jones, R. O.; Gunnarsson, O. The Density Functional Formalism, Its Applications and Prospects. Rev. Mod. Phys. 1989, 61, 689–746. Becke, A. D. Density-Functional Exchange-Energy Approximation With Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098. Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the Electron-Density. Phys. Rev. B 1988, 37, 785–789. Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B 1996, 54, 1703.


3218


Possibility of Coherent Delocalized Nuclear ... - ACS Publications

Recommend Documents