H2 Adsorption in a Porous Crystal: Accurate First-Principles Quantum

Subscriber access provided by UNIV OF CAMBRIDGE

Article 2

H Adsorption in a Porous Crystal: Accurate First Principles Quantum Simulation Jordan H D'Arcy, Meredith J.T. Jordan, Terry James Frankcombe, and Michael Anthony Collins J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.5b06074 • Publication Date (Web): 31 Aug 2015 Downloaded from http://pubs.acs.org on September 2, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

H2 Adsorption in a Porous Crystal: Accurate First Principles Quantum Simulation

Jordan H. D'Arcy1, Meredith J. T. Jordan1, Terry J. Frankcombe2 and Michael A. Collins2

1

2

School of Chemistry, The University of Sydney, Sydney NSW 2006, Australia

Research School of Chemistry, Australian National University, Canberra ACT 0200, Australia.

Abstract A general method is presented for constructing, from ab initio quantum chemistry calculations, the potential energy surface (PES) for H2 absorbed in a porous crystalline material. The method is illustrated for the metal-organic-framework material MOF-5. Rigid body quantum diffusion Monte Carlo simulations are used in the construction of the PES and to evaluate the quantum ground state of H2 in MOF-5, the zero-point energy, and the enthalpy of adsorption at 0 K.

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 55

1. INTRODUCTION Molecular hydrogen is a possible alternative fuel for use in transport applications.1 In this context, it must be safely and efficiently stored. One possible storage mechanism is via adsorption of the gas in porous materials. The US Department of Energy has set benchmarks that indicate the required H2 storage capacity, at moderate temperatures and pressures, for such porous materials.2 Metal-organic-framework materials (MOFs),3,4 graphene,5 carbon nanotubes6 and clathrates7-10 are being investigated as suitable storage materials. In general terms, the aim is to find a stable, easily manufactured material that maximises the hydrogen storage capacity. Theoretical calculations provide a useful complement to experimental investigations of storage capacity, as not yet realised potential candidate materials can be investigated in silico. The theoretical study of adsorption of H2 in such materials requires an accurate description of the energetics of H2 interactions with the material as well as the interaction between H2 molecules adsorbed within the material. To understand, quantitatively, the nature of the adsorbed H2, its quantum mechanical behaviour must also be accurately calculated. Thus, quantitative studies of hydrogen adsorption require a potential energy surface (PES) and a quantum method appropriate to the dimensionality of the system. In this paper, we present a general approach to the construction of a PES which describes the interaction energy of H2 with a porous crystal lattice. The method is illustrated here for MOF5.4,11 This material is the prototypical metal-organic-framework: it is composed of (Zn4O)6+ cores coordinated with six 1,4-benzenedicarboxylate (BDC) linkers. The MOF-5 unit cell is illustrated in Figure 1. While this material does not have a very high capacity for H2 adsorption at room temperature, 12,13 it is a relatively simple structure and exhibits a number of key features required by potentially useful storage materials: (i)

a moderately large, stable, unit cell with connected void spaces,

(ii)

highly charged metal cations with large coordinated anions, and

(iii) highly polarizable moieties that provide large dispersion interactions with adsorbates.


2

Page 3 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


These properties make MOF-5 a particularly useful case study for the development of methods to determine the energetics of H2 adsorption in porous materials.14,15 Several experimental studies have probed the nature of the H2 adsorption sites in MOF-5,4,16-19 and the vibration/rotation energy levels of adsorbed H2.20-23 A force-field model has been developed, for which translation-rotation eigenstates have been calculated.15 In these calculations, the H2 nuclear wavefunction was observed to be delocalized over several adjacent MOF-5 binding sites. However, it is suggested that this was an artefact of the potential, which underestimated the barrier height between the binding sites. These results provide additional motivation to develop chemically accurate potentials for H2 in MOF-5. The PES construction method we present here is general for any crystalline, porous, nonconducting material. Indeed, with minor modifications, the method could be employed to describe the energetics of any small molecule in such a porous crystal lattice. Our application to H2 adsorption in MOF-5 shows proof-of-principle, that is, that our PES construction strategy is valid. It also allows us to study the quantum nature of the H2 adsorbate and accurately predict enthalpies of adsorption. Given the light mass of H2, classical dynamics is unlikely to provide an adequate description of the motion, and a quantum method is required. Significant quantum effects include nuclear spin isomerism,24 tunnelling through energy barriers25 and the delocalisation of the nuclear wavefunction.26,27 In particular, delocalisation may result in different binding enthalpies and H2 loadings than would be predicted by classical mechanics. Comparisons between classical and quantum-mechanical simulations of H2 uptake in microporous structures suggest that the classical treatment of H2 will over-estimate H2 adsorption, even at room temperature, where quantum effects might be thought to be negligible.28-30 Infrared and neutron spectroscopy have been used to probe the behaviour of H2 in MOFs and carbon-based structures. Using diffuse reflectance infrared techniques, FitzGerald and coworkers23 predicted a zero-point energy (ZPE) of 1.5 kJ mol−1 associated with the centre-of-mass motion of H2 adsorbed to MOF-5. The ZPE for H2 adsorbed to MOF-74 was also reported to be of


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 55

this magnitude.31 These ZPEs are a substantial fraction of the total binding energies, so quantum effects are clearly significant. In this work, we employ the quantum diffusion Monte Carlo (QDMC) technique to estimate the lowest energy state of H2 in MOF-5. This provides a useful case study of how an ab initio PES can be used to accurately predict the lowest energy state of an adsorbate in a general porous material. This quantum technique is also used to explore the most significant configurations of H2 in the crystal lattice; an exploration that is an integral part of the process of constructing the PES (as explained in detail below). Previous quantum studies of H2 physisorption28-30,32,33 have treated H2 as a point particle. Here we account for H2 as a diatomic molecule, to take account of orientation effects on binding with the lattice. The paper is set out as follows. In the following section, we describe the several techniques used to construct the PES and the QDMC method used to evaluate properties. Section 3 presents results from the PES, including binding energies, ZPE and probability density histograms. Concluding remarks are contained in the final section.

2. METHODS The construction of a potential energy surface (PES) for a small molecule like H2 in a porous crystalline material involves the combined use of several disparate techniques: ab initio quantum chemistry, perturbation approaches to long range molecular interactions, an energy-based fragmentation method, Shepard interpolation, the Grow methodology and QDMC simulations. Moreover, some modification of the Grow methodology is necessary to account for the periodicity of the crystal lattice. Summaries of the various methods, and more details of the necessary modifications, are given here.

2.1 Systematic fragmentation


4

Page 5 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Systematic molecular fragmentation by annihilation (denoted as SMF) is a procedure for breaking molecules, crystals and crystal surfaces into relatively small molecular fragments, from which the energy and properties of the whole system can be accurately estimated. The method has been described in detail previously,34-40 so only a brief sketch of the application of SMF to crystals is included here. Molecules and crystals are treated as collections of functional groups connected by single bonds. A functional group might contain a single atom, or some larger object, such as a benzene ring. The entire molecule or crystal is broken into relatively small sets of connected groups by following a systematic set of rules that determine which single bonds are broken. Where bonds are broken, the valence is restored by “capping” the molecular fragment with hydrogen atoms.36 The systematic set of rules contains a hierarchy, whereby at each successive level in the hierarchy, the bonds that are broken are farther and farther apart in terms of bonded connectivity. The resultant sets of fragments are labelled as a “Level 1” or “Level 2” or “Level 3” fragmentation, and so on. In the case of crystals, the periodicity of the lattice is retained. As an illustration, consider a hypothetical one-dimensional lattice, composed of groups A, B and C in each unit cell: (1D) Ccrys ≡ ....A n-1Bn-1Cn-1An BnCn An+1Bn+1Cn+1.... .

(1)

In a Level 1 fragmentation, the lattice is broken into fragments containing at most two groups: (1D ) (1D ) Ccrys → CLevel 1

= ...+An-1Bn-1 + Bn-1Cn-1 + Cn-1An + An Bn .... ... − A n-1 − Bn-1 − Cn-1 − ...... . ∞

=

∑ n=-∞

(2)

 Bn-1Cn-1 + C n-1An + A n Bn   −C − A − B   n-1 n n 

At Level 2 fragmentation,


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 55

(1D ) (1D ) Ccrys → CLevel 2

= ...+A n-1Bn-1C n-1 + Bn-1Cn-1A n + Cn-1A n Bn + A n Bn Cn .... ... − A n-1Bn-1 − Bn-1Cn-1 − Cn-1A n − ......

.

(3)

∞  Bn-1C n-1A n + C n-1A n Bn + A n Bn Cn  =∑   n=-∞  −C n-1 A n − A n Bn − Bn C n 

Each fragment is capped with hydrogen atoms. These fragmentations preserve the number and type of each atom in the crystal, and the hydrogen atoms in the fragments with negative signs “cancel” those in the fragments with positive coefficients. As the value of Level increases, the interactions between each group and its neighbouring groups are better and better preserved. The total energy of this one-dimensional lattice is therefore approximated at Level 1 by:

 E ( Bn-1Cn-1 ) + E ( C n-1An ) + E ( An Bn )  =∑  . −E C − E A − E B ( ) ( ) ( ) n=-∞   n-1 n n  ∞

Level 1 1D

E

(4)

Since the lattice is periodic, the energy per unit cell (which is all that is needed) is: Level 1 Eunit cell = E ( Bn-1C n-1 ) + E ( C n-1 A n ) + E ( A n Bn )

− E ( Cn-1 ) − E ( An ) − E ( Bn )

.

(5)

This method can be generalised to a crystal lattice in three dimensions and we can apply SMF to obtain a set of fragments at the specified value of Level.36,37 For later convenience, we define some quantities that are needed to describe the crystal and an adsorbate. Each atom in a crystal has Cartesian coordinates:

x(i,l1 ,l2 ,l3 ) = x(i,0,0, 0) + l1a(1) + l2 a(2) + l3a(3),

(6)

where i numbers the atoms in a unit cell, the {a(j)} represent the lattice vectors, and ( l1 ,l2 ,l3 ) are integers that enumerate the unit cells in the crystal. When the crystal is fragmented, some fragments (sets of connected groups with hydrogen caps) can be associated with the “central unit cell”. For example, at Level 1, these are the terms in the square brackets of eq 3 with n = 0, and the corresponding terms in eq 4 at Level 2. A fragment associated with the central unit cell is denoted

Fn (0, 0, 0) . A fragment that is related to this fragment by a simple lattice translation is denoted


6

Page 7 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fn (k1 , k2 , k3 ) . By definition, if the atom with coordinates x(i,l1 ,l2 ,l3 ) is contained in Fn (0, 0, 0) , then the atom at x(i, k1 + l1 , k2 + l2 , k 3 + l3 ) is contained in Fn (k1 , k2 , k3 ) . Corresponding to eq 3, the fragmentation of a crystal lattice can be denoted by

C→

∞

∞

∞

N frag

k1 = −∞

k2 = −∞

k3 = −∞

n =1

∑ ∑ ∑ ∑

cn Fn (k1 , k2 , k3 ) ,

(7)

where the cn are simple integers, and there are a total of Nfrag fragments that describe a unit cell. The fragment, Fn (k1 , k2 , k3 ) , is a set of connected atoms, that have a corresponding set of coordinates

{x [ m(i),l(i) + k ],i = 1,...., N

F

(n)} , where fragment Fn contains NF(n) atoms, and where m(i) denotes

that the ith atom in the fragment is the m(i)th atom in the unit cell. 2.1.1 Examples It is useful to consider some illustrative examples: silicon and MOF-5. At standard temperature and pressure, silicon has the diamond structure, with 8 Si atoms per unit cell, each tetrahedrally bonded to 4 other Si atoms. At Level 1 fragmentation the unit cell is represented by 16 Si2 fragments (each with cn = 1 in eq 7) and 8 fragments with a single Si atom (each with cn = −3 in eq 7). Since hydrogen atom caps are used to restore the valence, the Level 1 fragments are actually disilane, Si2H6, and silane, SiH4. It is important to note that due to the high symmetry of the silicon lattice, all 16 Si2H6 fragments have the same structure, and all 8 SiH4 fragments have a common structure. Hence, eq 7 represents the silicon crystal lattice at Level 1 fragmentation in terms of just two unique fragment structures. Similarly, at Level 2 fragmentation, the unit cell is represented by 8 Si5 fragments (each with cn = 1 in eq 7) and 16 Si2 fragments (each with cn = −2 in eq 7). Again, due to the high symmetry, there are only two unique structures, Si5H12 and Si2H6. Figure 1 depicts the 424 atoms in the unit cell of MOF-5. The crystal structure can be viewed as an array of functional groups connected by single bonds. In SMF, bonds between oppositely charged atoms are considered to be stronger than normal single bonds, so that such connected atoms are taken to be in a common functional group. Hence, (Zn4O)6+(CO2-)6 comprises ACS Paragon Plus Environment

7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 55

a single (neutral) functional group. Phenylene, C6H4, is the only other functional group in MOF-5. The crystal structure is then viewed as a cubic array of (Zn4O)6+(CO2-)6 groups, connected along each edge by C6H4 groups. There are 8 (Zn4O)6+(CO2-)6 groups and 24 C6H4 groups in the unit cell, giving a total of 32 groups in the unit cell (denoted as Ngroup). For Level 1 fragmentation, the unit cell is represented by a set of 112 fragments, with hydrogen caps: 48 (Zn4O)6+(HCO2-)5 CO2-C6H5 molecules (each with cn = 1 in eq 7), 40 (Zn4O)6+(HCO2-)6 molecules (each with cn = −1 in eq 7) and 24 C6H6 molecules (each with cn = −1 in eq 7). Because of the high symmetry of the crystal structure, there is only one unique structure for each of these three molecules, as illustrated in Figure 2.

2.2 Interaction energy From eq 7, the energy of the entire crystal lattice is given by a sum of these fragment energies: ∞

Ecrys =

∞

∞

N frag

∑ ∑ ∑ ∑ c E {x [ m(i),l(i) + k ],i = 1,...., N n

k1 =−∞

k2 =−∞ k3 =−∞

n

n=1

F

(n)}  .

(8)

As described in previous work,36,37 eq 8 neglects the contributions to the total lattice energy that arise from the interaction between atoms that are sufficiently far apart that they never occupy a common fragment. We neglect such interactions herein, because we are not concerned with the total lattice energy, but only with the interaction of an adsorbate molecule, H2, with the lattice. This interaction energy, denoted ∆E(H 2 ) , is given by

∆E ( H 2 ) = Ecrys ( H 2 ) − Ecrys ,

(9)

where the first term on the right hand side of eq 9 indicates that H2 is included in the crystal. Using eq 8, this interaction energy is given by a sum of interaction energies between H2 and the fragments that comprise the crystal:

∆E ( H 2 ) =

∞

∞

∞

N frag

k1 =−∞

k2 =−∞

k3 =−∞

n=1

 En { x(i),i = 1,.., N F (n)} , { x ( H(i),i = 1,2 )}     . n −En {x(i),i = 1,.., N F (n)}  

∑ ∑ ∑ ∑c


(10)

8

Page 9 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Here, x(H(i)) represents the coordinates of the adsorbed hydrogen atoms, and the inclusion of

{x ( H(i),i = 1,2 )} denotes the addition of the H2 molecule in the ab initio calculation of the fragment energy. According to eq 10, we need to evaluate the interaction of H2 with each fragment in every unit cell. 2.2.1 Example As noted above, each unit cell of MOF-5 is represented at Level 1 SMF by 112 fragment molecules (Nfrag = 112 in eq 10). The unique fragment structures are illustrated in Fig. 2. Since there are only three unique fragment geometries, the electronic energies of these three structures provide all the data needed to evaluate the negative term in eq 10. However, 112 electronic structure calculations would be necessary to evaluate the energy of H2 in the presence of each of the fragments in one unit cell, since the relative position of the H2 is different for each fragment. These 112 calculations should be repeated an infinite number of times, since the relative position of the H2 is different for each fragment in each unit cell.

2.3 Perturbation approximation Fortunately, it is not necessary to evaluate all the electronic energies in eq 10 using ab initio quantum chemistry. The unit cell of a porous crystalline material may be very large (at least thousands of Å3 in volume), and largely devoid of atoms. This is the case we consider herein. The H2 molecule is small. Hence, if H2 is in one unit cell, it must be well separated from the atoms in all other unit cells. Thus the interaction of H2 with fragments in other unit cells should be accurately described by perturbation theory estimates of the interaction between molecules at large distances. The principal types of interaction that must be accounted for are electrostatic interactions (a first order term in the perturbation expansion), the associated induction effect (second order) and the dispersion interaction (also second order). Before discussing the evaluation of this perturbation energy in detail, we invoke a further approximation. Since the unit cell volume is large, if the position of the H2 molecule is close to one


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 55

fragment in a unit cell, it must be far from most other fragments in that unit cell. For example, in MOF-5 the cubic unit cell has a length of 25.669 Å. Hence if the H2 molecule is within a few Angstrom of some atom in the unit cell, it is of the order of 20 Å distant from most other atoms in the unit cell. Figure 3(a) illustrates that, if the H2 is closest to an atom in a (Zn4O)6+(CO2-)6 group, then it is normally only close to three SMF Level 1 fragments. Figure 3(b) illustrates that, if the H2 is closest to an atom in a C6H4 group, then it is only close to three (different) SMF Level 1 fragments. The hydrogen caps shown in Figure 3 indicate the fragment to be subtracted in the first term of eq 10 to avoid overcounting of the H2 interaction. This fragment has cn = −1 in eq 10. In general, we can define a subset S of fragments which are close to the H2 adsorbate. For MOF-5, we define S to be the set of three closest fragments. Hence, we arrange eq 10 as

 En { x(i),i = 1,.., N F (n)} , { x ( H(i),i = 1, 2 )}     ∆E ( H 2 ) = ∑ cn   n −E  x(i),i = 1,.., N (n)  { }   n F    Fn ∈S ∞

+∑

k1 =−∞

∞

N frag

∞

∑ ∑ ∑c

n

k2 =−∞

k3 =−∞

n=1 Fn ∉S

 En {x(i),i = 1,.., N F (n)} , { x ( H(i),i = 1,2 )}       −En { x(i),i = 1,.., N F (n)}   ,

(11)

 En { x(i),i = 1,.., N F (n)} , { x ( H(i),i = 1, 2 )}     ≡ ∑ cn   n −En {x(i),i = 1,.., N F (n)}     F ∈S n

+∆E pert ( H 2 ) where all the interaction energies in the second (infinite) sum are evaluated perturbatively. Although ∆E pert ( H 2 ) in eq 11 contains an infinite sum, the interaction of a H2 molecule with lattice fragments decays rapidly with distance to the fragment, so that the sum can be truncated to a reasonable number of terms, which are easily evaluated. Thus, the interaction energy of a H2 molecule with a MOF-5 lattice mainly requires the ab initio calculation of the interaction of H2 with just a few molecular fragments. This approach is quite general, although the number of ab initio calculations may differ for different materials. For example, the same approach for the covalentorganic-framework material COF-10241 requires ab initio calculations for five fragments with H2, rather than three fragments for MOF-5.


10

Page 11 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The perturbation component of the interaction energy is evaluated in two parts, an electrostatic component and a part due to dispersion:

∆E pert ( H 2 ) = ∆Eele ( H 2 ) + ∆Edisp ( H 2 ) .

(12)

The interaction energy due to the induction effect is incorporated into the electronic structure calculations, as described below. 2.3.1 Electrostatics The electrostatic interaction between the lattice and the H2 molecule is evaluated as a sum of the interaction of an electric quadrupole moment located at the midpoint of the H2 bond with charges, dipoles and quadrupoles located on each atom in the lattice. That is, a sum of two-body interactions. Ab initio calculations, for the three unique structures obtained from Level 1 fragmenation of MOF-5, were carried out at the required level of ab initio theory (see Section 2.7). The electric charge density for each fragment is then represented as a set of distributed charges and multipoles on each atom, via Stone's method 42,43. The net charge on each atom is then given by q(i): N frag

q(i) =

∑

cn q(i;n) ,

(13)

n=1 i∈Fn

where q(i;n) is the charge on atom i in fragment Fn. The distributed charges and multipoles on any fragment can be evaluated from the data for the corresponding unique structure. Similar formulae apply to the dipole and quadrupole on each atom. The ab initio calculation for each of the unique fragment structures is now repeated, except that each fragment is embedded in a set of point charges. These point charges are located at all the atoms in a 3 × 3 × 3 block of unit cells surrounding the fragment of interest. Charges at the locations of the fragment atoms are excluded from the array of point charges. The atomic charges evaluated by eq 13 are used. The resultant new charge densities are used to evaluate new charges and multipoles, via eq 13 and Stone's method. The point charges are updated with the new charges, and the ab initio calculations repeated. This self-consistent iterative procedure gives well-converged charges after approximately four iterations. ACS Paragon Plus Environment

11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 55

The iterative procedure accounts for the mutual polarisation of the electron density on each group by all the other groups in the crystal. The final set of atomic charges, dipoles and quadrupoles are used in the evaluation of the electrostatic interaction of the lattice with H2. The quadrupole moment of H2 was evaluated at the required level of ab initio theory for discrete H2 bond lengths. This quadrupole moment was fitted to a quadratic polynomial function of the bond length, and this function was used to evaluate the quadrupole for arbitrary values of the bond length and orientation of H2 (see the Supporting Information for details). The electrostatic energy is then evaluated as a sum over the lattice atom charge (q), dipole (µ) and quadrupole ( Θ ) interacting with the quadrupole on H2: Eele ( H 2 ) =

∞

∞

∞

N uc

k1 =−∞

k2 =−∞

k3 =−∞

i=1 i∉S

∑ ∑ ∑ ∑

EqΘ  x(i, k), x ( H 2 )  + EµΘ  x(i,k), x ( H 2 )  + EΘΘ  x(i, k),x ( H 2 )  .

(14)

The form of these electrostatic interactions is well known.44 2.3.2 Dispersion The dispersion interaction is evaluated as a sum of interactions between H2 and each functional group, G, in the lattice, excepting the groups that are contained in the few fragments for which ab initio interactions are evaluated. Letting Ngroup denote the number of groups in the unit cell: ∞

Edisp =

∞

∞

N group

∑ ∑ ∑ ∑

k1 =−∞

k2 =−∞ k 3 =−∞

U disp G ( n, k1, k2 , k3 ) ,H 2  .

(15)

n=1 n∉S

The dispersion interaction, Udisp, between H2 and each group can be approximated as follows. The detailed derivation has been presented previously,45 so only a summary is given here. Let Z(H2) represent the midpoint of the H2 bond, and X(n,k1,k2,k3) represent the geometric centre of the group G(n,k1,k2,k3). Let −1

D = Z ( H 2 ) − X ( n, k1, k2 , k3 ) ,

(16)

Dαβ = ∇α ∇ β D .

(17)

and


12

Page 13 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Then, 44

U disp G ( n, k1 , k2 , k3 ) ,H 2  = −

∞ 1 G (n ) Dαβ Dγδ ∫ α αγH2 (iν )α βδ ( iν ) dν , 0 2π

(18)

H2 where α αγ ( iν ) is the Cartesian polarizability tensor for H2 as a function of imaginary frequency, G (n) and α βδ ( iν ) is the corresponding quantity for group G(n). Repeated indices in the tensors in eq

18 are summed. As previously derived, eq 18 can be accurately approximated by45 H2 U disp G ( n, k1 , k2 , k3 ) ,H 2  = −2Dαβ Dγδ α αγ ( 0 )α βδG (n) ( 0 )

P H 2 P G (n) . P H2 + P G (n )

(19)

Here, α (0) denotes the usual static polarizability tensor, and 2

1 ∞  α (iν )  P= dν , 2π ∫0  α (0)  where

(20)

α is the isotropic component of the polarizability. G(n) The P factor was evaluated, using the DALTON program package,46,47 for each

functional group (just two groups in the case of MOF-5) at the required level of ab initio theory, and H H stored. P 2 and α 2 ( 0 ) have been calculated, using GAUSSIAN09,48 at a series of H2 bond

lengths and fitted to polynomial functions for later evaluation at arbitrary bond lengths (see the Supporting Information for details). Given some position of H2 in the lattice, the Dαβ tensor components are readily evaluated, and eq 19 is calculated for each group in eq 15. The dispersion interaction energy falls rapidly as the distance between H2 and a group increases, so the summation in eq 15 is accurately approximated by setting −1 ≤ k1 , k2 , k3 ≤ 1 . 2.3.3 Induction Induction is a fundamentally many-body effect. The electron distribution in the H2 molecule is polarised by the net electric field at the location of the H2 molecule, due to all the charges and electric multipoles in the crystal. However, polarisation of the entire crystal by H2 is only treated approximately, herein. Short range interactions, including polarisation of the nearby lattice molecules by H2, are explicitly included in the electronic energy calculations for the three “closest” ACS Paragon Plus Environment

13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 55

MOF-5 fragments. Because of the small polarizability of the H2 molecule and the relatively large distances between it and the “distant” MOF-5 fragments, polarization of these distant lattice molecules is likely to be small in this case, and such effects are neglected. This approximation significantly simplifies the calculation of the potential energy. Incorporation of the induction energy has been achieved by a number of authors for calculations of molecular and cluster energies, by employing embedded charges in the ab initio calculations.49-57 That is, all ab initio calculations are carried out in the presence of the point charges, which were evaluated in the iterative procedure above. In practise, point charges are placed at the positions of all lattice atoms in a K × K × K block of unit cells surrounding the fragment (except at locations of atoms in the fragment). Simple trials have shown that, for MOF-5, we can take K = 3 without loss of accuracy; higher values of K only change the energy by a value of order 10−7 au. In principle, a similar method might used to include the induction effect of H2 on the MOF lattice. Note: Since the ab initio calculations implicitly account for the interaction of the quadrupole of H2 with the crystal atom charges, only the H2 quadrupole–lattice-atom dipole and –lattice-atom quadrupole interactions are included in the electrostatic calculation above; EqQ  x(i, k),x ( H 2 )  is set to zero.

2.4 Total interaction energy It is useful to summarize the total interaction energy of the H2 molecule with the lattice:

 En {x(i),i = 1,.., N F (n)} , { x ( H(i),i = 1,2 )}; {q}embedded   ∆E ( H 2 ) = ∑ c j   n −En { x(i),i = 1,.., N F (n)}; {q}embedded     Fn ∈S ∞

∞

∞

N uc

k1 =−∞

k2 =−∞

k 3 =−∞

i=1 i∉S

∞

∞

∞

N group

k2 =−∞

k 3 =−∞

n=1 n∉S

+∑

+∑

k1 =−∞

∑ ∑ ∑ ∑ ∑ ∑

EµΘ  x(i, k),x ( H 2 )  + EΘΘ  x(i, k),x ( H 2 ) 

(21)

U disp  G ( n, k1, k2 , k3 ) ,H 2  ,


14

Page 15 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where the presence of embedded charges, {q}embedded, in the ab initio evaluation of fragment energies is explicitly noted. The sum over unit cells in each direction in eq 21 is restricted to a 3 × 3 × 3 block of cells without significant loss of accuracy. It is important to note that ∆E ( H 2 ) is a periodic function of the Cartesian coordinates of H2. That is, translation of the H2 molecule by an integer number of lattice vectors does not change the interaction energy. Moreover, the symmetry of the crystal lattice is reflected in the symmetry of

∆E ( H 2 ) . The space group of MOF-5 is Fm-3m, which has 192 elements.58,59 This means that any configuration of H2 within a unit cell is one of 192 energetically equivalent positions within that cell. Note also, that the “internal” energy of the H2 molecule has not been removed from ∆E ( H 2 ) , so that the internal energy associated with changes to the H2 bond length is included in ∆E ( H 2 ) .

2.5 Construction of the PES Equation 21 gives the interaction energy for H2 in some given location within the crystal. To describe the H2 wavefunction in the lattice, we need to evaluate this interaction energy at any relevant configuration. That is, we need a global PES. This global PES for the energy of H2 in a crystal is a function of six degrees of freedom, three for each of the hydrogen atoms. Since the crystal lattice is held rigid, the positions of the crystal atoms do not appear explicitly in the PES. The H2 wavefunction calculated on this PES may occupy a significant volume in this six dimensional space, so that evaluation of eq 21 (including electronic structure calculations) at every required configuration would be prohibitively time consuming. Hence, we construct an explicit function for the PES, which can subsequently be evaluated at all necessary configurations at manageable computational cost. We express the PES in the form of a modified Shepard interpolation, as in numerous previous applications to molecules and crystal surfaces.60-69 Significant modifications of this


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 55

approach are necessary to properly account for the lattice periodicity and symmetry. The general form of the interpolated PES is: .

(22)

Here, the index n enumerates “data points”, which are configurations of the H2 molecule in the lattice. At each data point, we evaluate ∆E ( H 2 ) and the first and second derivatives of ∆E ( H 2 ) with respect to convenient coordinates (see below). A second order Taylor expansion, Tn, for the value of ∆E ( H 2 ) at nearby configurations can then be evaluated. A set of Ndata data points is created (see Section 2.5.5). At any arbitrary configuration of H2, each Taylor series will result in a different value of ∆E ( H 2 ) . The final value of ∆E ( H 2 ) in eq 22 is given by a weighted average of the Taylor series values, where wn denotes the relative weight assigned to each Taylor series. Because of the symmetry of the lattice, any configuration is one of an equivalent set of configurations (for example, consider a configuration of H2 near any one of the eight equivalent “corners” in the unit cell of Figure 1). These symmetry equivalent positions are transformed from one to another by the operations of the symmetry group, which for crystals is denoted as the space group. In addition, the hydrogen atoms in H2 are indistinguishable, so permutation of the two hydrogen atoms leaves ∆E ( H 2 ) unchanged. Hence, in eq 22, GS denotes the symmetry group, which is the product of the space group and S2 (the permutation group for two objects). The symbol g in eq 22 denotes an element of GS. The data set is expanded (at no cost in electronic structure calculations) to include all symmetry equivalent data points. Finally, we note that a suitable choice of the weight function, w, ensures that eq 22 is a strict interpolation of the value of ∆E ( H 2 ) . That is, if eq 22 is evaluated at a configuration that is exactly at a data point, then the value of ∆E ( H 2 ) at that data point is obtained. 2.5.1 Coordinates


16

Page 17 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The coordinates describing the position of the two hydrogen atoms in the lattice are chosen to be periodic functions, so translational symmetry of the PES is ensured. This approach follows from a similar treatment for molecular collisions with a crystal surface.64 We recall that a(1), a(2), and a(3) denote the lattice vectors which define the unit cell. Let z(1), z(2), and z(3) denote the corresponding reciprocal lattice vectors. Then, the position of each hydrogen in the unit cell, X(i), is given by the fractional coordinates (x,y,z) as:

X(i) = x(i)a(1) + y(i)a(2) + z(i)a(3); i = 1,2,

(23)

where

X(i)T z(1) = x(i); X(i)T z(2) = y(i);

(24)

X(i)T z(3) = z(i). Translation of X(i) by one lattice vector changes the corresponding fractional coordinate by one. We define 12 periodic coordinates, ρ i :

ρ i = sin ( 2π X(1)T z(i)) ,i = 1,.., 3; ρ 3+i = cos ( 2π X(1)T z(i)) ,i = 1,.., 3; ρ 6+i = sin ( 2π X(2)T z(i)) ,i = 1,.., 3;

(25a)

ρ 9+i = cos ( 2π X(2)T z(i)) ,i = 1,.., 3. Note that we employ both cosine and sine functions to produce a redundant set of coordinates (the space is only six dimensional). In addition, we define a 13th coordinate:

ρ13 = X(1) − X(2)

−1

.

(25b)

The last coordinate, a reciprocal atom-atom distance, has been found to be very useful in Taylor expansions of molecular PES. 63 ρ13 is also (trivially) periodic with respect to lattice translation. The reason for defining a redundant set of coordinates is made apparent by the construction of the Taylor series and the operation of the space group. 2.5.2 Taylor series.


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 55

Construction of the Taylor series requires that the first and second derivatives of ∆E ( H 2 ) , with respect to the coordinates, have the correct periodicity. The derivatives of ∆E ( H 2 ) are initially obtained from the electronic structure calculations as derivatives with respect to the Cartesian coordinates of the two hydrogen atoms. The corresponding Cartesian derivatives of the electrostatic and dispersion energies are obtained, using the fact that the quadrupole moment, polarizability and P factor for H2 are known functions of the bond length. The derivatives of the total interaction energy, ∆E ( H 2 ) , are transformed into local periodic quantities for each Taylor series, using the established methodology63 with some modifications. To understand this, we write the matrix of derivatives of the { ρ (i)} coordinates with respect to the Cartesians (Y), at the data point, [Y(n)] as B:

Bij =

∂ ρi ∂Y j

i = 1,...,13; Y1 = X1 (1),......,Y6 = X 3 (2).

;

(26)

Y=Y(n)

This “Wilson B matrix” 70 can be written in a singular value decomposition as

B = UΛVT ,

(27)

where U is a 13 × 13 unitary matrix, V is a 6 × 6 unitary matrix, and Λ is a diagonal matrix of socalled singular values. The variation of the { ρ i } corresponding to variation of the Cartesian coordinates about the data point is given from eqs 26 and 27 as

δρ = UΛVT δ Y .

(28)

Since there are only six independent degrees of freedom, there can only be six linear combinations of the { ρ i } coordinates that vary independently as the Cartesian coordinates vary, and so there are only six non-zero singular values in Λ . We assume these are the first six diagonal elements of Λ . Thus, we define six coordinates in the locality of each data point: 13

ξi = Λ ii−1 ∑ ( UT ) ρ k , i = 1,..., 6; k=1

ik

(29)

13

≡ ∑ Rik ρ k . k=1


18

Page 19 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


From eq 28, we see that variations of these coordinates about some configuration correspond to variations of a unitary set of Cartesian coordinates:

δξ = VT δ Y .

(30)

The Cartesian energy derivatives can be transformed into derivatives with respect to the {ξ } coordinates:

∂∆E ( H 2 ) ∂∆E ( H 2 ) ∂∆E ( H 2 ) ∂Y . = =V ∂ζ ∂Y ∂Y ∂ζ ξ =ξ (n) ξ =ξ (n) Y=Y(n) Y=Y(n )

(31)

Note that the {ξ } are “curvilinear coordinates”, so that in transforming the Cartesian second derivatives to derivatives with respect to the {ξ } coordinates, we have to correct for the second derivatives of these coordinates with respect to the Cartesians:

∂ 2 ∆E ( H 2 ) 6 =∑ ∂Yi ∂Y j k=1

6

∑ l=1

2 ∂ξ k ∂ ∆E ( H 2 ) ∂ξl 6 ∂∆E ( H 2 ) ∂2 ξ k . +∑ ∂Yi ∂ξ k ∂ξl ∂Y j k=1 ∂ξ k ∂Yi ∂Y j

(32)

We can now define a Taylor expansion for ∆E ( H 2 ) about each data point in terms of periodic coordinates that are defined uniquely for each data point (see eq 29): 6

Tn = ∆E ( H 2 )ξ =ξ (n) + ∑ ξ i − ξi (n)  i=1

6 ∂∆E ( H 2 ) ∂ 2 ∆E ( H 2 ) 1 6 . (33) + ∑ ξ i − ξi (n) ∑ ξ j − ξ j (n)  ∂ζ i 2 i=1 ∂ξi ∂ξ j ξ =ξ (n) j=1 ξ =ξ (n)

There are important reasons for defining the coordinates and associated Taylor series in this complicated way (involving redundant coordinates). The { ρi } coordinates are periodic functions, so that ∆E ( H 2 ) has the periodicity of the lattice. The use of both cosine and sine functions in eq 25a ensures that each position in the unit cell is uniquely identified by the { ρi } coordinates, and that the B matrix of eq 26 has six non-zero singular values. Use of (say) cosine functions only would result in zero derivatives in eq 26 at some positions in the unit cell, and would also not uniquely define the hydrogen atom positions. Moreover, the { ρi } coordinates form a representation of the space group for the crystal. This would not be the case if only sine or only cosine functions were employed in eq 25a, as some symmetry operations in some space groups transform one ACS Paragon Plus Environment

19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 55

configuration to another which is one quarter of a unit cell away [for example, transforming

cos(2π x) to cos(2π x + π / 2) = − sin(2π x) ]. Thus, we require a redundant set of coordinates. 2.5.3 Symmetry The advantage of using coordinates that form a representation of the symmetry group, is that .64 The

it is a simple matter to form Taylor series at all symmetry-related configurations, space group operations are available as transformations of the fractional coordinates.58,59 It is

relatively straightforward to evaluate the corresponding transformations of the { ρ i } coordinates as matrix multiplications: .

(34)

The Appendix shows that only the definition of the local coordinates changes from a data point to a symmetry-related data point, so that . (35)

2.5.4 Weight function In this first application of Shepard interpolation to PES in crystals, we have employed the simple form of weight function, which was used initially in applications to gas phase molecular PES.60 The weight of data point n at some configuration, ρ , is defined in terms of a primitive weight factor, vn :

vn ( ρ ) = ρ − ρ (n)

−2 p

,

(36)

where p is a positive integer. The relative weight in eq 22 is then ,

(37)

so that .

(38)


20

Page 21 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A positive value of p ensures that wn ( ρ ) → 1, as ρ → ρ (n) , so that eq 22 is an interpolation. For sufficiently large p, it can be shown that eq 22 converges to the exact PES as the number of data points approaches infinity. 63 We have set p = 8 in the calculations on MOF-5 presented below. 2.5.5 Growing the PES The data set for the PES of eq 22 is “grown” iteratively. A small set of configurations for H2 in the crystal is chosen as the initial data set. Data points approximating the four MOF-5 binding sites, α, β, γ and δ,15,16 were generated by optimizing the position and orientation of the H2 molecule relative to a rigid (Zn4O)6+(HCO2-)5 CO2-C6H5 fragment. An additional data point was generated by placing the H2 at a ‘bridging’ site between the α and δ binding sites. Thus, the initial H2 + MOF-5 PES data set comprised five symmetry-unique data points. This initial PES is then “explored” to discover optimum locations for additional data points. In previous applications of the Grow scheme 63,66,67 classical trajectories and/or QDMC simulations71-74 have been used to explore configuration space. After each new data point is added, further simulations on the now modified PES, explore the configuration space in a slightly different way. This process of exploration, followed by choice and addition of a new data point, is repeated over and over to grow the size of the data set in locations determined by the classical dynamics or QDMC sampling. The PES is taken to be “converged” when the value of calculated observables of interest do not change significantly with increasing number of data points. As in previous studies,6669,75

QDMC simulations alone were used to explore the PES and to calculate the ground state

wavefunction and ZPE. A new data point, at each iteration, was determined by the variance in the PES, σ 2 : .

(39)

This variance is a measure of the uncertainty in the estimate of ∆E ( H 2 ) . During each iteration, over 4 million configurations are sampled and σ 2 is evaluated for each. The configuration with the largest value of σ 2 is chosen as a new data point. In addition, previous experience with growing ACS Paragon Plus Environment

21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 55

molecular PES, suggests that the number of data points required to converge the PES is reduced if prompt action is taken to place data points in locations where the PES of eq 22 is spuriously low in value. The lowest energy configuration is estimated from the lowest energy data point in the existing data set. If the QDMC calculation reveals one or more locations where eq 22 gives a value lower than the known lowest energy, then the location of the lowest energy QDMC configuration is also added to the data set.

2.6 Quantum diffusion Monte Carlo The QDMC simulations employed here have been discussed in detail previously.69,75 Rigid body quantum diffusion Monte Carlo (RBQDMC) is used to explore the configuration space of the adsorbed H2 molecule within MOF-5. Similarly to our previous work, we have fixed the H2 bond length. Here it is set at its optimized value when adsorbed to the α binding site of the metal oxide fragment, approximately 0.75 Å. This site is illustrated in Figure 3(a) and the coordinates of all of the atoms in the optimized fragment are given in the Supporting Information. Whilst we can, in principle, perform six-dimensional simulations that include the H2 bond length coordinate, the intrinsic time-scale for H2 vibration is considerably faster than that of the intermolecular modes; QDMC parameters suitable for H2 vibration result in slow convergence of the intermolecular components of the ground state wavefunction. Conversely, if we optimise the QDMC parameters to the intermolecular degrees of freedom, the high frequency H2 vibration introduces significant statistical noise into the final results. Because we are primarily interested in the intermolecular degrees of freedom associated with H2 adsorption, we have frozen the H2 bond length. The RBQDMC simulations are thus only five-, rather than six-dimensional. Walker configurations are generated as random displacements of the H2 centre-of-mass, followed by random rotations around its principal axes. The magnitude of the displacements and rotations are taken from appropriate Gaussian distributions, in terms of either mass or moment of inertia, as described in Reference 73.


22

Page 23 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In constructing the PES, we use two types of RBQDMC simulation: “growing” and “convergence”, the parameters of which are given in Table I. “Growing” RBQDMC simulations are performed for each Grow iteration. They are designed to explore the configuration space relevant to the ground state H2 wavefunction within MOF-5. Sampled configurations are used to choose additional points to add to the PES data set, using the criteria given in Section 2.5.5. “Convergence” RBQDMC simulations are used to evaluate the ZPE for the intermolecular degrees of freedom of the adsorbed H2 molecule as a function of the number of data points in the PES data set. TABLE I. Parameters used in the five-dimensional RBQDMC simulations. Parameter

“Growing”

“Convergence”

Time step (a.u.)

5.0

5.0

Feedback parameter (α)a

0.2

0.2

Maximum initial displacement (a.u.)

9.0

9.0

Number of equilibration steps

3000

4000

Number of sampling steps

3000

30000

Number of walkers per simulation

1350

13500

Number of walkers per sub-cell

50

500

Number of independent simulations

1

20

Number of generations for descendant weighting simulations

N/A

60

Delay between successive generations (in time steps)

N/A

500

Number of steps per generation

N/A

200

Descendant weighting parameters:

a

Feedback parameter, α, used as defined in Ref. 73

The time step and feedback parameters used in the RBQDMC simulations were chosen as those previously used for simulations for the H2-Li+-benzene complex. 69,75 Preliminary optimisation of these parameters for this system showed that the ZPE was insensitive to these values. Within the MOF-5 motif, we can consider two ‘sub-cells’: one where the benzene plane of the BDC linkers faces inward, the other with benzene facing outward. In both “growing” and


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 55

“convergence” RBQDMC simulations, the initial walker configurations were distributed throughout 3 × 3 × 3 of the sub-cells in order to test for equilibration of the walker populations and demonstrate the periodicity of the system. Growing and convergence simulations were initiated with 50 and 500 walkers per sub-cell, respectively. Walkers were initially placed at the centre of each sub-cell, again as a test for equilibration. The large initial displacements given in Table I ensured that walkers were distributed throughout the sub-cell without an initial bias to any given binding site. Descendant weighting calculations were performed as part of the convergence simulations.73 These provide an approximation to the square of the ground state nuclear wavefunction, that is, they provide approximate ground state probability distributions. The descendant weighting parameters are given in Table I. Because the H2 molecule can move within MOF-5, the probability density histograms were defined in terms of the nearest (Zn4O)6+ core, and the origin redefined accordingly. That is, if descendant weighting walkers were beyond the MOF-5 central unit cell, they were translated to the equivalent position within the unit cell, and the nearest O atom of the eight (Zn4O)6+ clusters determined. Figures 4(a) and 4(b) depict the coordinate system used. The central O atom of the nearest (Zn4O)6+ is defined as the origin; the z-axis is the vector joining the central O atom to the centre of the sub-cell containing the H2 molecule. The x- and y-axes make up a righthand set, where the x-axis is defined as the projection of one of the O-Zn vectors onto the plane containing three of the Zn atoms (the fourth Zn atom is on the z-axis; see Figure 4). The position of the H2 centre-of-mass was binned into probability density histograms in three degrees of freedom: r, the radial distance from the origin to the H2 centre-of-mass; θ, the polar angle between r and the z-axis, and ϕ, the azimuthal angle between the x-axis and the H2 centre-ofmass projected in the xy plane. The {r, θ, ϕ} coordinates are shown in Figure 4. We define two additional coordinates to describe the orientation of the H2 molecule: the ‘ferris wheel’ angle between the H-H bond and r, and the ‘helicopter’ angle, the dihedral angle between the z-axis, r and the H-H bond. Thus, the ferris wheel angle describes the orientation of the H2 molecule in a plane defined by the origin and the two H atoms, while the helicopter angle describes the orientation of


24

Page 25 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the H2 molecule in a plane perpendicular to r. The histograms presented below project the probability density onto the various internal coordinates. They were monitored as a function of Ndata, the number of symmetry unique data points in the PES data set.

2.7 Electronic structure calculations For each data point, the energy, gradient and Hessian must be calculated for the necessary Level 1 fragments, H2 + (Zn4O)6+(HCO2-)5CO2-C6H5, H2 + (Zn4O)6+(HCO2-)6 and H2 + C6H6, for MOF-5. The largest of these molecules contains 41 atoms and 308 electrons, which puts a constraint on the size of the basis set and level of theory that can be readily applied. The quantum chemistry approaches that have been applied to MOFs have been recently reviewed.76 Here we briefly discuss the previous studies of MOF-5. Second order Møller-Plesset perturbation theory (MP2) calculations have been used to estimate the binding energy of H2 with some moieties in MOF-5.77 In particular, calculations were performed for H2 + (Zn4O)6+(HCO2-)6. A triple zeta valence basis (TZVP) was used for the Zn atoms, and the aug-cc-pVTZ basis78 for all other atoms. The binding energy was estimated to be −6.9 kJ mol−1. A more recent paper,79 gives the basis set limit of the MP2 binding energy at the α site of the H2 + (Zn4O)6+(HCO2-)5 CO2-C6H5 fragment as −7.6 kJ mol−1. MP2 calculations with very large basis sets are not feasible for the construction of a global PES, as presented herein. However, since dispersion is likely to make a large contribution to the binding energy of H2 in the crystal, standard DFT calculations are also not appropriate. Hence, we have adopted the use of dispersion-corrected DFT methods. We have used the def2-TZVP basis set80,81 with the GAUSSIAN09 program package,48 and the generalized gradient approximation functional, PBE0.82 We include an empirical treatment of dispersion, as suggested by Grimme83 (“EmpiricalDispersion=GD3BJ” in GAUSSIAN09). In all electronic structure calculations, the fragments have the same structure as in the crystal11 (CCDC number 256965), with capping hydrogens attached as previously described.36 The PBE0-D3/def2-TZVP calculations yield slightly


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 55

stronger binding for H2 + (Zn4O)6+(HCO2-)5 CO2-C6H5 than the MP2/CBS calculations: −7.8 vs −7.6 kJ mol−1. As each data point in the PES requires the second derivatives of the energy with respect to the six coordinates that specify the position of the H2, it is more tractable to consider these six degrees of freedom than to constrain the H2 bond length. We have found that the most computationally efficient method to compute these second derivatives is by central difference of the Cartesian gradients. Hence, 13 gradient calculations (one for the original position of H2 and 12 displaced geometries) are performed for each of the required fragments. For MOF-5, there are three “close” fragments, so that a total of 39 gradient calculations are carried out in parallel. Each calculation is performed in the presence of embedded charges, obtained using the iterative procedure described in Section 2.3.

3. RESULTS Figure 5 shows the ZPE calculated by “convergence” RBQDMC simulations, as a function of the number of symmetry-unique data points, Ndata, in the PES. The total ZPE for H2 adsorbed to MOF-5, estimated by Matanović and coworkers,15 is shown with a dotted line. The ZPE shown in Figure 5 increases from an initial value of approximately 2.77 kJ mol–1 at Ndata = 5 to a maximum of 3.26 kJ mol–1 at Ndata = 35, before converging to approximately 2.84 kJ mol–1 for Ndata ≥ 105. The statistical errors in these calculations are considerably smaller than the variation arising from the different individual PES. Error bars, representing twice the standard error of the mean of 20 independent RBQDMC simulations, are within the size of the symbols used in Figure 5. It is worth nothing that the permutation group for one H2 molecule adsorbed to MOF-5 contains 384 elements: 192 symmetry operations of the Fm-3m space group, and the permutation of the H atoms. Thus, our PES with 135 symmetry-unique data points comprises a total of 51,840 data points in the MOF-5 unit cell. Matanović and colleagues15 estimate a total ZPE of 2.9 ± 0.3 kJ mol−1 for H2 adsorbed to ACS Paragon Plus Environment

26

Page 27 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MOF-5. Our converged value is well within error of this estimate. Indeed, the RBQDMC ZPE calculated for Ndata ≥ 45 is within 0.1 kJ mol–1 of this number. Furthermore, our RBQDMC ZPE is fully anharmonic and includes coupling between the translational and rotational degrees of freedom of the H2 molecule. The ultimate goal of our simulations is the accurate prediction of H2 adsorption enthalpies, ∆Hads, to a porous material such as MOF-5. We must understand the factors governing ∆Hads if we are to design microporous materials suitable for hydrogen storage at near to room temperature. Our RBQDMC simulations provide anharmonic, coupled values of the ZPE for the intermolecular degrees of freedom, and therefore allow accurate calculation of the 0 K H2 adsorption enthalpy for MOF-5, ∆Hads (0 K). Estimates of the electronic adsorption energy, ∆Eads, and ∆Hads (0 K) are given in Table II. TABLE II. PBE0-D3/def2-TZVP electronic H2 adsorption energies, ∆Eads (kJ mol–1), for H2 adsorbed to the four MOF-5 binding sites in a single (Zn4O)6+(HCO2-)5 CO2-C6H5 fragment are shown in column 2. The corresponding values for the converged H2 + MOF-5 PES are shown in column 3. The 0 K H2 adsorption enthalpy, ∆Hads (0 K) (kJ mol–1), calculated on the converged H2 + MOF-5 PES is also given. Fragment Species

∆Eads

MOF-5

MOF-5 ∆Eads

∆Hads (0 K)

−9.4

−6.5

MOF-5 α binding site

−7.8

−9.4

MOF-5 β binding site

−4.5

−5.5

MOF-5 γ binding site

−4.3

−4.0

MOF-5 δ binding site

−4.9

−3.6

Table II shows electronic adsorption energies, calculated at the PBE0-D3/def2-TZVP level of theory, for H2 adsorbed to the α, β, γ and δ binding sites of MOF-5. These were calculated both for a molecular fragment, H2 + (Zn4O)6+(HCO2-)5 CO2-C6H5, and on the converged H2 + MOF-5 PES. The total RBQDMC ZPE we have calculated for H2 adsorbed to MOF-5, 2.84 kJ mol−1, is with respect to the minimum energy on the converged H2 + MOF-5 PES, that is, with respect to ∆Eads = –9.4 kJ mol−1, as shown in Table II. This minimum corresponds to the electronic energy of ACS Paragon Plus Environment

27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 55

H2 bound to the α binding site. The 0 K H2 adsorption enthalpy, ∆Hads (0 K), calculated on the converged H2 + MOF-5 PES, is also shown in the table. It should be noted that this value of ∆Hads (0 K) refers to the lowest energy state of the H2 molecule, that is, the state with total angular momentum, J = 0, corresponding to para-H2. In the simplest approximation, ∆Hads (0 K) for orthoH2 can be calculated using the lowest energy H2 nuclear wavefunction for J = 1. A fixed-node formulation of the RBQDMC method would allow the inclusion of coupling and anharmoniicty between the translational and rotational degrees of freedom of the adsorbed H2 molecule84 but this is beyond the scope of this work. As shown in Table II, and consistent with previous studies,77,79 H2 binds most strongly to the α site. In our RBQDMC simulations we have a single H2 molecule in the periodic MOF-5 potential, and we anticipate it will adsorb only to the α site. This is based on the fact that the ZPE is smaller than the difference in electronic binding energies between the α site and the β, γ and δ binding sites. The ∆Eads for the β, γ and δ binding sites, calculated for the H2 + (Zn4O)6+(HCO2-)5 CO2C6H5 fragment, have magnitudes of 4.3–4.9 kJ mol–1 and they can be considered approximately isoenergetic for the fragment. The magnitudes and relative ordering of ∆Eads for the four binding sites differ from the fragment to the converged H2 + MOF-5 PES. The magnitude of ∆Eads increases for the α site, by 1.6 kJ mol−1, to 9.4 kJ mol−1. A 1.0 kJ mol−1 increase is also observed for the β site. The magnitude of ∆Eads for the γ site remains approximately the same, while the magnitude of ∆Eads for the δ site decreases by 1.3 kJ mol−1. Whilst the most favorable binding site is the α site for both the fragment and MOF-5 potentials, the δ site is least favored in the periodic potential. This is consistent with experiment.16 Our calculated ∆Hads (0 K), for one H2 molecule adsorbed to MOF-5, was −6.5 kJ mol−1. In this case, the rigid-body ZPE of the H2 molecule accounts for 31 % of the total electronic adsorption energy, indicating that quantum effects cannot be ignored. If the H-H bond length were allowed to relax from its value in the gas phase to its value for the adsorbed hydrogen, the change in ZPE for this degree of freedom (0.34 kJ mol–1 in the harmonic limit) would increase the magnitude


28

Page 29 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of ∆Hads (0 K). This gives an estimate for ∆Hads (0 K) of −6.8 kJ mol−1. Heats of adsorption for gas sorption in porous materials at finite temperature can be estimated from experimental uptake data. For H2 adsorbed to MOF-5 at 77 K, estimates of ∆Hads range from −3.8 kJ mol−1 85 to −4.8 kJ mol−1.86 Our ∆Hads (0 K) value is broadly consistent with these results. As temperature increases, thermal contributions will reduce the magnitude of ∆Hads. Ground state nuclear probability density histograms were obtained from convergence RBQDMC simulations, performed on the PES with Ndata = 135. Figure 6 shows RBQDMC probability density histograms for the {r, θ, ϕ} spherical polar coordinates defining the location of the H2 centre-of-mass: r, panel (a); θ, panel (b) and ϕ, panel (c). These histograms implicitly include the spherical polar volume elements and this is apparent in Figure 6(b). The r probability density histogram shows a single distribution in r, indicating that, at low H2 loading (here, one adsorbed H2) adsorption occurs at a single site. The H2 is indeed adsorbed at the MOF-5 α site, and this is described further below. The r histogram has maximum probability density at r ≈ 3.97 Å, which is 0.30 Å longer than the global minimum energy H2-oxygen distance. The full-width-half-maximum (FWHM) of the r histogram is 0.84 Å, and the vibrationally averaged expectation value, 〈〉, is 4.06 Å. This is in good agreement with results obtained from neutron diffraction experiments, which indicate a vibrationally averaged r distance of 3.95 Å.17 The relatively large FWHM indicates that there is significant delocalisation of the H2 nuclear wavefunction in this coordinate. Figure 6(b) displays the probability density histogram for the polar angle, θ, between the H2 centre-of-mass and the z-axis. The histogram peaks at θ ≈ 6°, with a FWHM of approximately 9°. The expectation value, 〈〉, is approximately 7°. At the global minimum energy geometry, the H2 centre-of-mass is on the z-axis; that is, θ = 0°. However, the maximum of the θ histogram and the value of 〈〉 reflect the fact that there is no probability of finding the wavefunction at θ = 0°, as the volume element at this point is zero. The relatively small FWHM suggests that the H2 wavefunction is comparatively localised in the θ coordinate. The probability density histogram for the azimuthal angle, ϕ, is shown in Figure 6(c). The ϕ


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 55

coordinate describes the location of the H2 centre-of-mass with respect to the x- and y- axes, shown in Figure 4(b). The ϕ probability density histogram exhibits three-fold symmetry with maxima at 0°, 120° and 240°. These angles correspond to the three closest Zn atoms to H2, which form the “cup” of the MOF-5 α binding site, with the O atom at the coordinate origin. The maxima in the ϕ histogram indicates that H2 preferentially adsorbs over a Zn atom. It is important to note, however, that the probability density in ϕ is never zero; there is finite probability of finding the H2 molecule at any azimuthal angle ϕ. The H2 molecule is not localised in the ϕ coordinate. The three-fold structure of the ϕ probability density histogram is apparent even from the initial PES, defined with Ndata = 5. The amplitude of the periodicity, however, increases markedly as Ndata increases from 35 to 40, and is effectively constant for Ndata > 40. This is an indication that at least 40 symmetry-unique data points are needed to accurately describe the interaction of H2 at the α binding site. In contrast, there is relatively little change in the θ histogram as Ndata increases, reflecting the comparative localisation of the H2 wavefunction in θ. The width of the r probability density histogram is approximately constant with Ndata, although the maximum varies from r ≈ 3.5 Å, for Ndata = 5 to r ≈ 4.0 Å for Ndata > 40. Figure 7 shows the probability density histograms for the ‘ferris wheel’ angle, panel (a), and the ‘helicopter’ angle, panel (b), that describe the relative orientation of the H-H bond. The probability density for the ferris wheel angle peaks at 90°, and has a FWHM of approximately 84°. The maximum at 90° shows that the H2 molecule is preferentially oriented such that the H-H bond is perpendicular to the r vector. This is consistent with the minimum energy geometry of H2 bound to the frozen fragment, shown in Figure 3(a), where the optimised angle is approximately 90°. Notably, the large FWHM of this histogram indicates that there is considerable delocalisation of the H2 in this coordinate. There is also very little probability (~10−5) of finding the H2 bond parallel to the r vector. The potential energy in this coordinate is relatively flat in the region about 90°, but increases markedly as the ferris wheel angle approaches 0° (or equivalently, 180°). This is similar to the interactions we have previously observed between H2 and Li+-doped benzene.69


30

Page 31 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7(b) shows the probability density histogram for the helicopter angle, describing the orientation of the H2 molecule perpendicular to the ferris wheel angle. The probability density is almost constant for all helicopter angles, with a slight preference shown for angles of approximately 90º and 270º. The H2 molecule is, essentially, a free rotor in this coordinate. This is also consistent with previous work on H2-Li+-benzene.69 The ferris wheel probability density histogram shows very little change with the value of Ndata. All PES considered accurately represented this degree of freedom. The two-fold symmetry seen in the helicopter angle probability density histogram does change with Ndata, becoming stable and robust for Ndata > 45. The anisotropy in the PES in this degree of freedom is small and is not accurately described until sufficient data has been added to the PES data set. The probability density histograms for the ferris wheel and helicopter rotations demonstrate anisotropy in the PES for the H2 rotational degrees-of-freedom. This is consistent with the two-dimensional PES scans performed by Matanović and colleagues,15 for the two H2 rotational coordinates. Figure 8 compares the location of D2 molecules in MOF-5, as determined by neutron diffraction,16 with the results of our RBQDMC simulations. The experiments used D2 rather than H2, as it has a large neutron scattering cross section, and neutron powder diffraction patterns were recorded at a loading of 32 D2 molecules per MOF-5 unit cell. Our initial RBQDMC simulations used H2 because quantum effects are more significant in the lighter isotopologue. The contour diagrams shown in Figure 8 clearly illustrate the locations of the adsorbed hydrogen molecules, shaded with red-yellow-green. The primary adsorption site is the α site, with the highest H2 adsorption energy (Table II). These are located in ‘inward’ MOF-5 sub-cells described in Section 2.6. Within each inward sub-cell there are eight α sites, with four shown in Figure 8. The secondary H2 adsorption sites, the β sites, are located in the ‘outward’ sub-cells, in the centre of Figure 8(a) and (b). Because of the higher D2 loading in the neutron diffraction experiments, primary and secondary adsorption sites are apparent in Figure 8(a), whereas only the primary adsorption sites are seen in Figure 8(b). At higher H2/D2 loading the γ and δ sites are


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 55

occupied.16 The comparison with experiment serves to validate our methodology. For a single H2 molecule we see adsorption only at the α sites in the inward sub-cells. This is consistent with experimental results.4,16-19 No adsorption is seen in the outward sub-cells, even when the RBQDMC walkers are initially distributed evenly throughout both types of sub-cell. We anticipate H2 will adsorb at the other binding sites as we increase the H2 loading in our RBQDMC simulations. Our accurate, anharmonic, fully coupled H2 + MOF-5 PES makes these simulations possible.

4. CONCLUDING REMARKS A general method has been presented for evaluating the energy of interaction of a H2 molecule with a porous crystalline material. The crystal lattice is held rigid. The interaction energy is calculated using a combination of ab initio quantum chemistry and perturbation theory. Systematic molecular fragmentation (SMF) was used to greatly simplify the quantum chemistry calculations to that of a few molecular fragments interacting with H2. Weaker, longer range, interactions of the H2 molecule with other parts of the lattice are evaluated using perturbation theory. All the molecular properties needed for these perturbation calculations are obtained from ab initio quantum chemistry calculations. There are no empirical or adjustable parameters employed. The use of SMF and perturbation theory greatly facilitates the energy calculation. A global PES was constructed from this ab initio data using the Grow methodology with Shepard interpolation. The established methodology has been modified to employ periodic coordinates to ensure that the resultant PES has the correct lattice translational symmetry. The space group of the lattice, together with the permutation group for H2, has been employed to ensure that the PES has the correct symmetry under all symmetry operations. The high symmetry of such crystals greatly increases the number of effective ab initio “data points” which define the PES, and thereby greatly improves the accuracy of the PES interpolation. ACS Paragon Plus Environment

32

Page 33 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Comparison with experimental neutron diffraction results shows that our RBQDMC simulations on the H2 + MOF-5 PES reproduce the key features of hydrogen adsorption. Whilst this work is proof-of-principle, our methods can be applied, within the same methodology, to additional H2 loading and to other isotopologues. The method, as a whole, can also be applied to other guest molecules within the MOF-5 lattice, or any other non-conducting crystal.


33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 55

APPENDIX (i) Calculation of the symmetry matrices For each element of the space group, g, the fractional coordinates of the two hydrogen atoms, represented by the six-dimensional column vector R, are transformed as ,

(A.1)

The elements of the unitary matrix M( R) (g) are typically 0 or ±1. The vector b(g) is given by a sum of the lattice vectors multiplied by fractions like 0, 1/2 or 1/4. M( R) (g) and b(g) are easily obtained from the crystallographic tables.58,59 To evaluate the corresponding unitary matrix representation for the { ρi } coordinates, we take a set (Nrandom = 6) of random values for R, {R(i)}, and minimise the function (A.2) with respect to the elements of the unitary matrix M( ρ ) (g) . Setting the derivatives of L with respect to the elements of M( ρ ) (g) to zero, gives a simple set of linear equations for M( ρ ) (g) . These equations are easily solved for M( ρ ) (g) , and the process is repeated for each symmetry operation, g. The permutation of the two identical hydrogen atoms permutes the elements of { ρi } in an obvious fashion, allowing straightforward evaluation of the symmetry matrices involving this permutation. (ii) Coordinates and the Taylor series The Cartesian coordinates for each hydrogen atom are related to the fractional coordinates by eq 23. In eq 26 we wrote the six Cartesians for the two hydrogen atoms as a six-dimensional column vector, Y. If we combine the lattice vectors a(1), a(2), and a(3) in a nine-dimensional column vector, A, then we can rewrite eq 23 for the two hydrogen atoms as

Y = CA ,

(A.3)


34

Page 35 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where C is a 6 × 9 matrix, with elements given by the fractional coordinates. In this format, it is easy to see that eqs (A.1) and (A.3) give the symmetry operations on the Cartesian coordinates as ,

(A.4)

where b '(g) is a combination of the elements of the vector A. Importantly, the same symmetry matrix, M( R) (g) , applies to both the Cartesian and fractional coordinates. Then, from eqs 26 to 28 and eq A.4, the effect of the symmetry operations on the B matrix can be determined:

(A.5)

. Then,

(A.6)

implies ,

(A.7)

.

(A.8)

and

The singular values are unchanged by the symmetry operations as both Cartesian and local coordinates are correspondingly transformed by the operations. From eq 29, the definition of the local coordinates about a symmetry transformed data point is changed by the symmetry operation: 13

−1 (ρ ) ξ (g) (g)U )ik ρ k , i = 1,..., 6. i = Λ ii ∑ ( M T

(A.9)

k=1


35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 55

However, the value of the local coordinates at the symmetry transformed data point is unchanged:

(A.10)

Similar derivations show that the value of the first and second derivatives of the energy, evaluated at the transformed data point, are the same as the corresponding values at the original data point, because the ρ coordinates and the definition of the local coordinates are transformed in complementary fashion.

Acknowledgements J.H.D. acknowledges the financial support of the Australian Postgraduate Award. The authors wish to thank Dr. Lars Goerigk (University of Melbourne) for his valuable assistance regarding electronic structure methods for MOF-5. This research was undertaken with the assistance of resources from the National Computational Infrastructure (NCI), which is supported by the Australian Commonwealth Government. References (1)

Balat, M. Potential Importance Of Hydrogen as a Future Solution to Environmental and Transportation Problems. Int. J. Hydrogen Energy 2008, 33, 4013–4029.

(2)

"Hydrogen and Fuel Cells Program, 2014 Annual Progress Report" U.S. Department of Energy, Washington D. C, 2014; http://www.hydrogen.energy.gov/annual_progress14.html, accessed 8/11/2015.

(3)

Suh, M. P.; Park, H. J.; Prasad, T. K.; Lim, D.-W. Hydrogen Storage in Metal-Organic Frameworks. Chem. Rev. 2012, 112, 782–835.


36

Page 37 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(4)

Rosi, N. L.; Eckert, J.; Eddaoudi, M.; Vodak, D. T.; Kim, J.; O’Keeffe, M.; Yaghi, O. M. Hydrogen Storage in Microporous Metal-Organic Frameworks. Science 2003, 300, 1127– 1129.

(5)

Pumera, M. Graphene-Based Nanomaterials for Energy Storage. Energy Environ. Sci. 2011, 4, 668–674.

(6)

Orikov, R.; Orik, A. Recent Applications of Carbon Nanotubes in Hydrogen Production and Storage. Fuel 2011, 90, 3123–3140.

(7)

Struzhkin, V. V.; Militzer, B.; Mao, W. L.; Mao, H.-k.; Hemley R. J. Hydrogen Storage in Molecular Clathrates. Chem. Rev. 2007, 107, 4133–4151.

(8)

Florusse, L. J.; Peters, C. J.; Schoonman, J.; K. C. Hester; Koh, C. A.; Dec, S. F.; Marsh, K. N.; Sloan, E. D. Stable Low-Pressure Hydrogen Clusters Stored in a Binary Clathrate Hydrate. Science 2004, 306, 469–471.

(9)

Frankcombe, T. J. The Importance of Vibrations in Modelling Complex Metal Hydrides. J. Alloys Compds. 2007, 446-447, 455–458.

(10) Matsumoto, Y.; Grim, R. G.; Khan, N. M.; Sugahara, T.; Ohgaki, K.; Sloan, E. D.; Koh, C. A.; Sum, A. K. Investigating the Thermodynamic Stabilities of Hydrogen and Methane Binary Gas Hydrates. J. Phys. Chem. C 2014, 118, 3783–3788. (11) Li, H.; Eddaoudi, M.; O’Keeffe, M.; Yaghi, O. M. Design and Synthesis of an Exceptionally Stable and Highly Porous Metal-Organic Framework. Nature 1999, 402, 276–279. (12) Kaye, S. S.; Dailly, A.; Yagh, O. M.; Long, J. R. Impact of Preparation and Handling on the Hydrogen Storage Properties of Zn4O(1,4-benzenedicarboxylate)3 (MOF-5). J. Am. Chem. Soc. 2007, 129, 14176–14177. (13) Panella, B. Hydrogen Physisorption in Metal–Organic Porous Crystals. Adv. Mat. 2005, 17, 538–541. (14) Belof, J. L.; Stern, A. C.; Space, B. A Predictive Model of Hydrogen Sorption for MetalOrganic Materials. J. Phys. Chem. C 2009, 113, 9316–9320.


37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 55

(15) Matanović, I.; Belof, J. L.; Space, B.; Sillar, K.; Sauer, J.; Eckert, J.; Bačić, Z. Hydrogen Adsorbed in a Metal Organic Framework-5: Coupled Translation-Rotation Eigenstates from Quantum Five-Dimensional Calculations. J. Chem. Phys. 2012, 137, 014701. (16) Yildirim, T.; Hartman, M. R. Direct Observation of Hydrogen Adsorption Sites and Nanocage Formation in Metal-Organic Frameworks. Phys. Rev. Lett. 2005, 95, 215504. (17) Spencer, E. C.; Howard, J. A. K.; McIntyre, G. J.; Rowsell, J. L. C.; Yaghi, O. M. Determination of the Hydrogen Absorption Sites in Zn4O(1,4-benzenedicarboxylate) by Single Crystal Neutron Diffraction Chem. Commun. 2006, 3, 278–280. (18) Rowsell, J. L. C.; Eckert, J.; Yaghi, O. M. Characterization of H2 Binding Sites in Prototypical Metal-Organic Frameworks by Inelastic Neutron Scattering. J. Am. Chem. Soc. 2005, 127, 14904–14910. (19) Mulder, F. M.; Dingemans, T. J.; Schimmel, H. G.; Ramirez-Cuesta, A. J.; Kearley, G. J. Hydrogen Adsorption Strength and Sites in the Metal Organic Framework MOF5: Comparing Experiment and Model Calculations. Chem. Phys. 2008, 351, 72–76. (20) Centrone, A.; Siberio-Perez, D. Y.; Millward, A. R.; Matzger, O. M. Y. A. J.; Zerbi, G. Raman Spectra of Hydrogen and Deuterium Adsorbed on a Metal–Organic Framework. Chem. Phys. Lett. 2005, 411, 516–519. (21) Bordiga, S.; Vitillo, J. G.; Ricchiardi, G.; Regli, L.; Cocina, D.; Zecchina, A.; Arstad, B.; Bjorgen, M.; Hafizovic, J.; Lillerud, K. P. Interaction of Hydrogen with MOF-5. J. Phys. Chem. B 2005, 109, 18237–18242. (22) Vitillo, J. G.; Regli, L.; Chavan, S.; Ricchiardi, G.; Spoto, G.; Dietzel, P. D. C.; Bordiga, S.; Zecchina, A. Role of Exposed Metal Sites in Hydrogen Storage in MOFs. J. Am. Chem. Soc. 2008, 130, 8386–8396. (23) FitzGerald, S. A.; Allen, K.; Landerman, P.; Hopkins, J.; Matters, J.; Myers, R.; Rowsell, J. L. C. Quantum Dynamics of Adsorbed H2 in the Microporous Framework MOF-5 Analyzed Using Diffuse Reflectance Infrared Spectroscopy. Phys. Rev. B 2008, 77, 224301.


38

Page 39 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(24) Santos, S. F. d.; Balakrishnan, N.; Lepp, S.; Quéméner, G.; Forrey, R. C.; Hinde, R. J.; Stancil, P. C. Quantum Dynamics of Rovibrational Transitions in H2-H2 Collisions: Internal Energy and Rotational Angular Momentum Conservation Effects. J. Chem. Phys. 2011, 134, 214303. (25) Kisvarsanyi, E. G.; Sullivan, N. S. Quantum Tunneling in Solid D2 and Solid H2. Phys. Rev. B 1995, 51, 3462–3468. (26) Fernandez-Alonso, F.; Cabrillo, C.; Fernández-Perea, R. F.; Bermejo, F. J.; González, M. A.; Mondelli, C.; Farhi, E. Solid Para-Hydrogen as the Paradigmatic Quantum Crystal: Three Observables Probed by Ultrahigh-Resolution Neutron Spectroscopy. Phys. Rev. B 2012, 86, 144524. (27) Warnecke, S.; Sevryuk, M. B.; Ceperley, D. M.; Toennies, J. P.; Guardiola, R.; Navarro, J. The Structure of Para-Hydrogen Clusters. EPJ D 2010, 56, 353–358. (28) Wang, Q.; Johnson, J. K. Molecular Simulation of Hydrogen Adsorption in Single-Walled Carbon Nanotubes and Idealized Carbon Slit Pores. J. Chem. Phys. 1999, 110, 577–586. (29) Garberoglio, G.; Skoulidas, A. I.; Johnson, J. K. J. Adsorption of Gases in Metal Organic Materials: Comparison of Simulations and Experiments. Phys. Chem. B 2005, 109, 13094– 13103. (30) Martínez-Mesa, A.; Yurchenko, S. N.; Patchkovskii, S.; Heine, T.; Seifert, G. Influence of Quantum Effects on the Physisorption of Molecular Hydrogen in Model Carbon Foams. J. Chem. Phys. 2011, 135, 214701. (31) FitzGerald, S. A.; Hopkins, J.; Burkholder, B.; Friedman, M.; Rowsell, J. L. C. Quantum Dynamics of Adsorbed Normal- and Para-H2, HD, and D2 in the Microporous Framework MOF-74 Analyzed Using Infrared Spectroscopy. Phys. Rev. B 2010, 81, 104305. (32) Teufel, J.; Oh, H.; Hirscher, M.; Wahiduzzaman, M.; Zhechkov, L.; Kuc, A.; Heine, T.; Denysenko, D.; Volkmer, D. MFU-4 - a Metal-Organic Framework for Highly Effective H2/D2 Separation. Adv. Mater. 2013, 25, 635–639.


39


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 55

(33) Martínez-Mesa, A.; Zhechkov, L.; Yurchenko, S. N.; Heine, T.; Seifert, G.; Rubayo-Soneira, J. Hydrogen Physisorption on Carbon Foams upon Inclusion of Many-Body and Quantum Delocalization Eﬀects. J. Phys. Chem. C 2012, 116, 19543–19553. (34) Deev, V.; Collins, M. A. Approximate Ab Initio Energies by Systematic Molecular Fragmentation. J. Chem. Phys. 2005, 122, 154102. (35) Collins, M. A.; Deev, V. A. Accuracy and Efficiency of Electronic Energies from Systematic Molecular Fragmentation. J. Chem. Phys. 2006, 125, 104104. (36) Netzloff, H. M.; Collins, M. A. Ab initio Energies of Nonconducting Crystals by Systematic Fragmentation. J. Chem. Phys. 2007, 127, 134113. (37) Collins, M. A. Ab initio Lattice Dynamics of Nonconducting Crystals by Systematic Fragmentation J. Chem. Phys. 2011, 134, 164110. (38) Collins, M. A.; Cvitkovic, M. W.; Bettens, R. P. A. The Combined Fragmentation and Systematic Molecular Fragmentation Methods. Acc. Chem. Res. 2014, 47, 2776–2785. (39) Collins, M. A. Molecular Forces, Geometries, and Frequencies by Systematic Molecular Fragmentation Including Embedded Charges. J. Chem. Phys. 2014, 141, 094108. (40) Collins, M. A.; Bettens, R. P. A. Energy-Based Molecular Fragmentation Methods. Chem. Rev. 2015, 115, 5607–5642. (41) El-Kaderi, H. M.; Hunt, J. R.; Mendoza-Cortés, J. L.; Côté, A. P.; Taylor, R. E.; O’Keeffe, M.; Yaghi, O. M. Designed Synthesis of 3D Covalent Organic Frameworks. Science 2007, 316, 268–272. (42) Stone, A. J. Distributed Multipole Analysis; or How to Describe a Molecular Charge Distribution. Chem. Phys. Lett. 1981, 83, 233–239. (43) Stone, A. J. Distributed Multipole Analysis: Stability for Large Basis Sets. J. Chem. Theory Comput. 2005, 1, 1128–1132. (44) Stone, A. J. The Theory of Intermolecular Forces; Clarendon: Oxford, 1996.


40

Page 41 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(45) Addicoat, M. A.; Collins, M. A. Accurate Treatment of Nonbonded Interactions Within Systematic Molecular Fragmentation. J. Chem. Phys. 2009, 131, 104103. (46) Aidas, K., et al. WIREs Comput. Mol. Sci. 2014, 4, 269. (47) Dalton. Dalton, a molecular electronic structure program, Release Dalton2015.X 2015. (48) Frisch, M. J., et al. Gaussian 09; Revision A.1 ed.; Gaussian, Inc.: Wallingford CT, 2009. (49) Jiang, N.; Ma, J.; Jiang, Y. Electrostatic Field-Adapted Molecular Fractionation with Conjugated Caps for Energy Calculations of Charged Biomolecules. J. Chem. Phys. 2006, 124, 114112. (50) Li, W.; Li, S. H.; Jiang, Y. S. Generalized Energy-Based Fragmentation Approach for Computing the Ground-State Energies and Properties of Large Molecules. J. Phys. Chem. A 2007, 111, 2193–2199. (51) Dahlke, E. E.; Truhlar, D. G. Electrostatically Embedded Many-Body Expansion for Large Systems, with Applications to Water Clusters. J. Chem. Theory Comput. 2007, 3, 46–53. (52) Dahlke, E. E.; Truhlar, D. G. Electrostatically Embedded Many-Body Correlation Energy, with Applications to the Calculation of Accurate Second-Order Møller−Plesset Perturbation Theory Energies for Large Water Clusters J. Chem. Theory Comput. 2007, 3, 1342–1384. (53) Xie, W.; Gao, J. The Design of a Next Generation Force Field: The X-POL Potential. J. Chem. Theory Comput. 2007, 3, 1890–1900. (54) Ji, C. G.; Mei, Y.; Zhang, J. Z. H. Developing Polarized Protein-Specific Charges for Protein Dynamics: MD Free Energy Calculation of pKa Shifts for Asp26/Asp20 in Thioredoxin. Biophys. J. 2008, 95, 1080–1088. (55) Le, H. A.; Lee, A. M.; Bettens, R. P. A. Accurately Reproducing. Ab Initio Electrostatic Potentials with Multipoles and Fragmentation. J. Phys. Chem. A 2009, 113, 10527–10533. (56) Le, H. A.; Bettens, R. P. A. Distributed Multipoles and Energies of Flexible Molecules. J. Chem. Theory Comput. 2011, 7, 921–930.


41


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 55

(57) Collins, M. A. Molecular Forces, Geometries, and Frequencies by Systematic Molecular Fragmentation Including Embedded Charges. J. Chem. Phys. 2014, 141, 094108. (58) Hahn, T. International Tables for Crystallography. Vol A., 5th ed.; John Wiley & Sons: Hoboken, 2005. (59) Bilbao Crystallographic Server, http://www.cryst.ehu.es/ (accessed 11 August, 2015). (60) Ischtwan, J.; Collins, M. A. Molecular Potential Energy Surfaces by Interpolation. J. Chem. Phys. 1994, 100, 8080–8088. (61) Jordan, M. J. T.; Thompson, K. C.; Collins, M. A. Convergence of Molecular Potential Energy Surfaces by Interpolation: Application to the OH+H2→H2O+H Reaction. J. Chem. Phys. 1995, 102, 5647–5657. (62) Thompson, K. C.; Jordan, M. J. T.; Collins, M. A. Polyatomic Molecular Potential Energy Surfaces by Interpolation in Local Internal Coordinates. J. Chem. Phys. 1998, 108, 8302– 8316. (63) Collins, M. A. Molecular Potential-Energy Surfaces for Chemical Reaction Dynamics. Theor. Chem. Acc. 2002, 108, 313–324. (64) Frankcombe, T. J.; Collins, M. A.; Zhang, D. H. Modified Shepard Interpolation of GasSurface Potential Energy Surfaces with Strict Plane Group Symmetry and Translational Periodicity. J. Chem. Phys. 2012, 137, 144701. (65) Frankcombe, T. J.; Collins, M. A. Growing Fragmented Potentials for Gas–Surface Reactions: The Reaction between Hydrogen Atoms and Hydrogen-Terminated Silicon (111). J. Phys. Chem. C 2012, 116, 7793–7802. (66) Crittenden, D. L.; Thompson, K. C.; Chebib, M.; Jordan, M. J. T. Efficiency Considerations in the Construction of Interpolated Potential Energy Surfaces for the Calculation of Quantum Observables by Diffusion Monte Carlo. J. Chem. Phys. 2004, 121, 9844–9854. (67) Thompson, K. C.; Crittenden, D. L.; Jordan, M. J. T. CH5+: Chemistry's Chameleon Unmasked. J. Am. Chem. Soc. 2005, 127, 4954–4958.


42

Page 43 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(68) Kolmann, S. J.; Jordan, M. J. T. Method and Basis Set Dependence of Anharmonic Ground State Nuclear Wave Functions and Zero-Point Energies: Application to SSSH. J. Chem. Phys. 2010, 132, 054105. (69) Kolmann, S. J.; D’Arcy, J. H.; Jordan, M. J. T. Quantum Effects and Anharmonicity in the H2-Li+-Benzene Complex: A Model for Hydrogen Storage Materials. J. Chem. Phys. 2013, 139, 234305. (70) Wilson, E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations; Dover: New York, 1955. (71) Anderson, J. B. A Randomwalk Simulation of the Schrödinger Equation: H3+. J. Chem. Phys. 1975, 63, 1499–1503. (72) Coker, D. F.; Watts, R. O. Quantum Simulation of Systems With Nodal Surfaces. Mol. Phys. 1986, 58, 1113–1123. (73) Suhm, A.; Watts, R. O. Quantum Monte Carlo Studies of Vibrational States in Molecules and Clusters. Phys. Rep. 1991, 204, 293–329. (74) Buch, V. Treatment Of Rigid Bodies By Diffusion Monte Carlo: Application To The Para-H2--H2O And Ortho-H2---H2O Clusters. J. Chem. Phys. 1992, 97, 726–729. (75) D'Arcy, J. H.; Kolmann, S. J.; Jordan, M. J. T. "Plug-and-Play" Potentials: Investigating Quantum Effects in (H2)2-Li+-benzene. J. Chem. Phys., 2015, in press, accepted 8/7/2015. (76) Odoh, S. O.; Cramer, C. J.; Truhlar, D. G.; Gagliardi, L. Quantum-Chemical Characterization of the Properties and Reactivities of Metal–Organic Frameworks. Chem. Rev. 2015, 115, 6051–6111. (77) Sagara, T.; Klassen, J.; Ganz, E. Computational Study Of Hydrogen Binding By MetalOrganic Framework-5. J. Chem. Phys. 2004, 121, 12543–12547. (78) Dunning, T. H.; Peterson, K. A. Basis Sets: Correlation Consistent Sets. In Encyclopedia of Computational Chemistry; Schleyer, P. v. R., Ed.; Wiley: Chichester, 1998; Vol. 1; pp 88. (79) Sillar, K.; Hofmann, A.; Sauer, J. Ab Initio Study of Hydrogen Adsorption in MOF-5. J. Am. Chem. Soc. 2009, 131, 4143–4150. ACS Paragon Plus Environment

43


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 55

(80) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. (81) Weigend, F. Accurate Coulomb-Fitting Basis Sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057–1065. (82) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158–6170. (83) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comp. Chem. 2011, 32, 1456–1465. (84) Buch, V. Treatment of Rigid Bodies by Diffusion Monte Carlo: Application to the ParaH2…H2O and Ortho- H2…H2O clusters. J. Chem. Phys. 1992, 97, 726–729. (85) Panella, B.; Hirscher, M.; Pütter, H.; Müller, U. Hydrogen Adsorption in Metal-Organic Frameworks: Cu-MOFs and Zn-MOFs Compared. Adv. Funct. Mater. 2006, 16, 520–524. (86) Rowsell, J. L. C.; Yaghi, O. M. Effects of Functionalization, Catenation, and Variation of the Metal Oxide and Organic Linking Units on the Low-Pressure Hydrogen Adsorption Properties of Metal-Organic Frameworks. J. Am. Chem. Soc. 2006, 128, 1304–1315.


44

Page 45 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure Captions: Figure 1. The MOF-5 unit cell. Zinc atoms are coloured blue, oxygen atoms are red, carbon atoms are black and hydrogen atoms are white.

Figure 2. The three unique molecular fragments obtained by Level 1 SMF of the MOF-5 unit cell. Zinc atoms are coloured blue, oxygen atoms are red, carbon atoms are black and hydrogen atoms are white.

Figure 3. The three nearest MOF-5 fragments for (a) H2 adsorbed at the α binding site of MOF-5, and (b) H2 adsorbed at the δ binding site. Hydrogen caps on the fragments in panels (a) and (b) indicate the fragment to be subtracted in the first term of eq 10; that is, where cn = −1. Figure 4. (a) The r and θ spherical polar coordinate system for H2 adsorbed to the α binding site in MOF-5, used to define probability density histograms in the descendant weighting simulations; (b) The coordinate system, viewed along the z-axis, showing the r and ϕ coordinates. Some framework atoms have been omitted for clarity.

Figure 5. The calculated RBQDMC ZPE as a function of the number of symmetry-unique data points in the PES data set. The dotted line is the total ZPE estimated by Reference 15. Error bars, representing twice the standard error of the mean of 20 independent RBQDMC simulations, are within the size of the symbols used.

Figure 6. Probability density histograms for the coordinates (a) r (Å), (b) θ (°), and (c) ϕ (°), for H2 adsorbed at the α binding site in MOF-5. Error bars, calculated as twice the standard error of the mean of 20 independent RBQDMC simulations, are within the thickness of the lines for all histograms.

Figure 7. Probability density histograms for the coordinates describing the orientation of the H2 molecule adsorbed at the α binding site in MOF-5: (a) the ‘ferris wheel’ angle (°), and (b) the ‘helicopter’ angle (°), as defined in Section 2.6. Errors, calculated as twice the standard error of the mean of 20 independent RBQDMC simulations, are within the thickness of the lines.

Figure 8. (a) The real space Fourier difference D2 neutron scattering-length density MOF-5-4D2, superimposed with the MOF-5 structure, reproduced with permission from Reference 16; (b) our RBQDMC probability density histograms projected onto the plane defined by the a- and b- lattice vectors, superimposed with the MOF-5 structure.


45


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 55

Table of Contents Graphic:


46

Page 47 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1 71x72mm (600 x 600 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Page 48 of 55

Page 49 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60





1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4 107x182mm (600 x 600 DPI)


Page 50 of 55

Page 51 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60





1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6 131x308mm (600 x 600 DPI)


Page 52 of 55

Page 53 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7 84x127mm (600 x 600 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 8 137x276mm (600 x 600 DPI)


Page 54 of 55

Page 55 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents Graphic 50x50mm (300 x 300 DPI)


H2 Adsorption in a Porous Crystal: Accurate First-Principles Quantum

Recommend Documents