4.6 Models
4.6.1 Force Field
Hydrocarbon molecules were modeled by a united atom (UA) approach, as depicted in Figure 4.7, because hydrogen atoms are known to have negligible influence on adsorption126,127 and diffusion128 of guest molecules in nano-porous host materials. The important degrees of freedoms (i.e., positions and velocities of carbon atoms), and the insignificant ones (hydrogen atoms bonded to C-atoms) are comprised together to single CHi beads in the UA approach, where i= 0 – 4. A molecule consists hence of several, possibly dif-ferent CHibeads. The striking advantage of the model is that it tremendously saves computing time∗ while simultaneously maintaining a realistic physical description of the systems.129
Several UA force fields for guest molecule simulation in zeolites exist.
Most of them took as a starting point the transferable potentials for phase equilibria (TraPPE) force field which has been introduced and is still main-tained by Ilja Siepmann and co-workers.130 The present work uses however the force field developed by Dubbeldam et al. for alkanes127 and extended by Liu et al. to alkenes.131 This is, because the TraPPE force field aims
∗The calculation of energies and forces usually represents the heaviest computational burden of a simulation so that the number of floating point operations performed in the force subroutine may well serve as a measure of the total computing time. A simple system consisting ofN atoms that interact via a pair-wise additive potential such as the Lennard–Jones type will give rise toN(N−1)/2 terms to be calculated in the subroutine per MD cycle.91Reducing the amount of atoms will obviously decrease the computing time overproportionately. For example, comparing a fully atomistic simulation withN methane molecules with a corresponding UA simulation decreases the number of relevant floating point operations in first order approximation by a factor of [N(N−1)/2]/[4N(4N−1)/2]≈ 1/16. It has to be stressed that this number does yet not include the gain achieved by neglecting intramolecular interactions.
4.6 Models 71
Figure 4.8: Intramolecular potentials; here n-butane (a) interactions with parameters from Dubbeldamet al.127Bond vibration (b), bond-angle bending (c), and torsion (d). Note that subscripts “CH3,CH2” etc. have been omitted for reasons of clarity.
at predicting a wide variety of systems and, first and foremost, fluid phase equilibria. By contrast, Dubbeldam’s and Liu’s force field was specifically de-veloped to match adsorption isotherms of hydrocarbons in zeolites.126,127,131
In the remainder of this chapter, the force field will be described and the parameters relevant to the present work provided.
The force field is based on five contributions to the total potential energy, Utot:
Utot =Uvib+Ubend+Utors +ULJ+Uelec, (4.65) where superscripts vib, bend, tors, LJ, and elec refer to bond vibration, bond-angle bending, torsion, dispersive Lennard–Jones, and electrostatic
interac-72 4 Molecular Simulations tion, respectively. Each contribution is either pair-, triple-, or quadruple-wise additive. Consider a pair-wise additive interaction typeUi. All possible pairs j and k of (united) atoms participating in that interaction have their indi-vidual contribution,uij,k. The separate contributions are then simply summed up to yield the final value of Ui:
Ui =
N−1$
j=1
$N k=j+1
uij,k, (4.66)
where N is the number of atoms and/or beads relevant to the interaction type i. The functional forms of the different interaction types will be given in the following.
Vibration of two bonded beads in the same molecule (Figure 4.8a) was modeled by a harmonic potential (Figure 4.8b):
uvibi,j =1/2·kvibi,j ·(ri,j −ri,j,eqvib )2, (4.67) whereki,jvib denotes the vibrational spring constant of the bond between beads i and j,ri,j the current distance between the two beads, and ri,j,eqvib the (con-stant) equilibrium bond length of that (type of) bond.
Two bond-angle bending terms were used, depending on the type of beads:
ubend1i,j,l =1/2·kbend1i,j,l ·(ϕi,j,l−ϕbend1i,j,l,eq)2, (4.68) ubend2i,j,l =1/2·kbend2i,j,l ·(cosϕi,j,l−cosϕbend2i,j,l,eq)2. (4.69) Again, kbendmi,j,l represents a potential constant that is specific for the bond angle considered, ϕi,j,l the current bond angle between atoms or beads i, j, and l, as well as ϕbendmi,j,l,eq the equilibrium bond angle. In Figure 4.8c, the second bending potential is displayed as a function of angle. In contrast to the vibrational term, the bending potential is periodic with two equally stable states that are separated by two differently high barriers.
The torsional potential is given by a cosine power series of the dihedral angle, φ?i,j,l,m:
utorsi,j,l,m=
$5 n=0
ηn·cosnφ?i,j,l,m, (4.70)
where each ηn is a constant. Similar to the bending potential, the torsional term is periodic (Figure 4.8d), exhibiting however three stable states and thus three barriers between them. The parameters of the just described intramolecular terms are provided in Table 4.1.
4.6 Models 73 Table 4.1: Bond vibration, bond-angle bending, and torsion parameters taken from References 127 and 131. R1 and R2 represent any possible alkyl remainder [e.g., CH3(sp3) and CH3(sp3) – CH2(sp3)].
kvib/kB rvibeq potential bond pattern kbend/kB ϕbendeq ηi/kB i uvib CH3(sp3) – CH3(sp3) 96 500 K/˚A2 1.54 ˚A
CH3(sp3) – CH2(sp3) CH2(sp3) – CH2(sp3)
CH2(sp2) = CH2(sp2) 786 873 K/˚A2 1.33 ˚A CH2(sp2) = CH(sp2) 785 050 K/˚A2 1.33 ˚A CH(sp2) – CH3(sp3) 251 549 K/˚A2 1.54 ˚A ubend1 CH2(sp2) – CH(sp2) – CH3(sp3) 70 400 K/rad2 119.7◦ ubend2 CH3(sp3) – CH2(sp3) – CH3(sp3) 62 500 K 114.0◦
CH3(sp3) – CH2(sp3) – CH2(sp3) CH2(sp3) – CH2(sp3) – CH2(sp3)
utors R1 – CH2(sp3) – CH2(sp3) – R2 1204.654 K 0
1947.740 K 1
-357.845 K 2
-1944.666 K 3
715.690 K 4
-1565.572 K 5
Interactions between two beads of different molecules were described by a Lennard–Jones type potential:
uLJi,j = 4·εi,j ·
7.σi,j
ri,j
/12
− .σi,j
ri,j
/68
, (4.71)
where ri,j denotes the distance between the interacting beads i and j and εi,j and σi,j are the respective Lennard-Jones parameters (well-depth and zero-potential distance). Figure 4.9 depicts the potential as a function of inter-bead distance, using the example of two methane united atoms. The
74 4 Molecular Simulations
interaction ener. uLJ /kB / [K]
distance between beads r / [nm]
Figure 4.9: Lennard–Jones potential; here for methane-methane interaction with parameters from Dubbeldam et al.127 Note that subscripts “CH4,CH4” have been omitted for reasons of clarity.
Lennard–Jones potential has a large region of attractive character, which causes the typical strong adsorption of hydrocarbons in siliceous zeolite struc-tures and is therefore of central importance to this work.
The energies and forces due to the Lennard–Jones potential were only explicitly calculated up to a certain cutoffdistance. To avoid any singularities in the potential energy, the cut potential was shifted to zero at the cutoff distance [i.e., by a value of ∆uLJ,shifti,j =−uLJi,j(rcutoff)], thus yielding: A value of 1.2 nm for rcutoff seems appropriate, considering the vanishing interaction energy at that distance; compare zoom in Figure 4.9.
Intra- and intermolecular forces have to be calculated on the basis of the potentials chosen in order to move a given system forward in an MD simulation. The force acting on beadj due to a pair-wise interaction of type i with beadk is given by:70
4.6 Models 75 For example, the xcomponent of the force due to a Lennard–Jones interac-tion is:91
wherexj,k is the distance between the two beads along the xcoordinate and the shift part has been omitted for clarity. Functional forms of the forces corresponding to the intramolecular interaction potentials are listed in the excellent PhD thesis by Jakobtorweihen,70 which has in fact inspired a large part of the present chapter.
Interactions between molecule beads and oxygen atoms of the siliceous zeolite structure consisted also of a Lennard–Jones potential only.127 In ac-cordance with the original force field proposal,127,131 interactions between guest molecules and silicon atoms of the host structure were neglected, which is common practice in zeolite research of hydrocarbon adsorption and trans-port with molecular simulations.132 The physical reason lies in the fact that the oxygen atoms are indeed large negatively charged ions that effectively shield the small silicon ions (cf., Chapter 2). Because hydrocarbons were studied only and extra-framework cations were not present (siliceous struc-tures) a Lennard–Jones description of the guest-host interaction was suffi-cient so that the term uelec could be entirely neglected in the present work.
Furthermore, intramolecular Lennard–Jones interactions were calculated for beads separated by more than three bonds only. That is, so-called 1-4 inter-actions were excluded. Table 4.2 lists the remaining parameters used in this work.
Although electrostatic interactions and polarizability were neglected due to the nature of the guest-host systems, it should be noted that these inter-action types become important for molecules such as carbon dioxide and water. In this respect, much research effort is currently being spent on improving existing force fields by multiscale simulation approaches. As a typically example, density functional theory based electronic structure cal-culations provide partial charges of host atoms derived from a single unit cell.133 The charges are then transferred to atomistically detailed GCMC simulations with much larger super cells to determine adsorption isotherms of guest molecules (e.g., CO2).133 The need for a more accurate description of the electrostatic guest-host interaction was already indicated by the work of Clark and Snurr in 1999.134 But it was not until recently that increasing computer resources facilitated the rigorous combination of electronic struc-ture methods and classical molecular modeling.
76 4 Molecular Simulations Table 4.2: Lennard–Jones parameters taken from References 127 and 131.
bead/ CH4 CH3(sp3) CH2(sp3) CH2(sp2) CH(sp2) atom
CH4
εi,j
kB
[K] 158.50 σi,j [˚A] 3.72 CH3(sp3) εi,j
kB
[K] 108.00 77.77
σi,j [˚A] 3.76 3.86 CH2(sp3) εi,j
kB
[K] 77.77 56.00
σi,j [˚A] 3.86 3.96 CH2(sp2) εi,j
kB
[K] 93.000 70.210
σi,j [˚A] 3.685 3.710
CH(sp2) εi,j
kB
[K] 70.210 53.000
σi,j [˚A] 3.710 3.740
O(zeol) εi,j
kB [K] 115.00 93.00 60.50 82.050 55.215
σi,j [˚A] 3.47 3.48 3.58 3.530 3.502