Applying an Electric Field via The Modern Theory of Polarisation
To perform AIMD calculations in an electric field it would seem that we only need to add an electric potential into the calculation of the force. However we have: ~E = rVE(r) (2.23).
which means that to apply a constant electric field in a given direction we need a linear potential, which is not periodic. The problem is that in our framework the Hamiltonian needs to have the same periodicity as the system. We can therefore only apply sawlike potentials, which have the right periodicity but which represent unrealistic, discontinuous electric fields. The solution is to include an electric field term in the energy functional. We obtain the new energy functional : EE[f kg] = E[f kg] EP[f kg] (2.24) where is the unit cell volume, and P[f kg] refers to the polarisation in the direction of the field, as given by the modern theory of polarisation explained briefly below. The sole dependency on the electronic density is lost because the polarisation s calculated using a quantum phase of the wavefunction, requiring information which cannot be obtained from the density alone. The macroscopic polarisation, also referred to as polarisation density P, is a vector quantity which expresses the electric dipole moment per unit volume of a material. This notion is useful to describe the response of materials to an electric field, such as inducing a switch of polarisation in ferroelectric materials. The microscopic description of polarisation is not as obvious. The historical ClausiusMossotti model defines the macroscopic polarisation PCM in a periodic solid as the sum of the dipole moments in a given cell divided by the volume of the cell.
This definition implies that we can use the electronic density n(r) to find areas of more positive or more negative charge densities to define dipole moments within a unit cell. In reality this is problematic because even if there are well defined areas of higher and lower charge density in a material, the choice of the unit cell is arbitrary and can lead to different values of P. This point is illustrated in figure 2.1.
Enhanced Sampling Techniques
The AIMD framework requires a large amount of computational resources as a tradeoff for accuracy in the description of polarisation and making and breaking of bonds in complex environments, such as the oxide/water interface which is studied in this work. However this means that within the available computational resources today, systems studied in this framework may contain at most several hundreds of atoms, with trajectory lengths of several tens of picoseconds.
These small sizes and timescales make it impossible to correctly sample the thermodynamics of a system, without resorting to some clever tricks known globally as enhanced sampling techniques. Many of these techniques are available in the literature, we focus here only on the ones used in this work – metadynamics and umbrella sampling. The choice of the reaction coordinate is also important, as is outlined below.
Parametrisation of the Free Energy Landscape
A free energy landscape can be described in onedimension if the reaction going from state A to state B has a clear reaction coordinate. The reaction coordinate can be defined in terms of the positions of the atoms in the system; it could be a physical distance between atoms, an angle, torsional angle, or any combination of these things. However such a reaction coordinate is not always obvious to find, and in complex cases can lead to restricting the system to follow a strict path in configuration space, which may not necessarily reflect the path of lowest free energy.
Path collective variables (pathCV) are a way to parametrise the free energy landscape according to configurations along a reaction path, and using these variables we can also define a second reaction coordinate in which the system may stray from the path. If we call the set of coordinates describing a state in the transition R, the vector describing the path is a set of Rk where the first ‘frame’ R1 = RA is the initial state and RN = RB is the final state. From this definition we can introduce variables s and z .
One monolayer of H2O on MgO(001)
A monolayer of water (1ML) on MgO(001) is known to exist in several configurations depending on temperature and pressure . For this study we chose to use the (2 3) reconstruction, which is stable between 185 K and 235 K in UHV conditions. In this reconstruction two out of six water molecules are dissociated as shown in figure 3.10. Initially the six water molecules were placed on top of MgO(001) and the positions allowed to relax. This resulted in spontaneous dissociation of two water molecules, with no energy barrier, in a configuration which was slightly different from the one in figure 3.10, but which exhibits the same number of OH groups on the surface and has the same symmetry (2 3 with a glideplane).
Nevertheless the configuration chosen is the one in figure 3.10 because it is accepted as the most stable configuration in the literature and because its energy was lowest. It is likely that several metastable analogues of this configuration exist, especially at nonzero temperatures, but due to the small size of the simulation box we chose for the following to remain with the most stable configu44
ration seen experimentally, even later on when several layers of water were added. occurs within the adsorbed layer at the height of the adsorbed H+ ions, and the second peak seems analogous to that in the bare surface case, i.e. intense, highly localised, but permeating a few Ångstrom beyond the surface. This peak is not as intense as in the bare surface case, which is a sign of screening by the H2O molecules. Figure 3.13 shows the profile in the xy plane at a height of 4 Å from the surface, so about 2 Å from the adsorbed monolayer. At this height the underlying surface structure is not longer ‘visible’ in the field profile, and is instead screened by the less ordered adsorbed monolayer. We notice zones of positive and negative z fields corresponding to local water dipoles, as well as small scale oscillations which would probably become irrelevant when taking into account the dynamics of the system. Some of these oscillations are probably also noise in the electric field, which being a gradient is more numerically noisy than the potential or charge density. Such complexity calls for a coarsegrained description of the electric field, especially when we consider its effect on the structural and dynamical properties of the water above the first adlayer. For the case of the water/MgO interface with more than a water adlayer, we will thus retain the profile along z: ¯Ez(z). Figure 3.11 shows selected vibrations taken from a simulation of the adsorbed monolayer at 300 K. The vibration frequency spectra were obtained from taking a Fourier transform of the selected OH bond lengths as a function of time. A velocity autocorrelation function of all atoms could also be used to be able to directly compare such a spectrum to experimental data, however in this case the statistics were not sufficient to obtain discernible peaks, and it was preferable to isolate specific bonds and apply the autocorrelation function method to get specific vibration modes. Preselecting the characteristic OH bonds enables us to check which frequency contributions come from the surface OH groups. The blue peak corresponds to the OH in vacuum frequency at about 3900cm1 and the green peak corresponds to the OHs which are directly adsorbed. We see that these are shifted towards smaller frequencies due to the higher coordination of the oxygen atoms, and the peak isn’t as sharp because of the increased interaction with surrounding oxygen atoms.
Surface assisted proton transfer
The green histogram in figure 3.15 was obtained by counting the number of times an oxygen atom had a switch in one of its two nearest neighbour hydrogen atoms, thus counting the number of proton transfer, or ‘proton hopping’ events. In the (2 3) reconstruction shown in figure 3.10, proton transfer is susceptible to occur, simply because the high density of the adsorbed layer means that the oxygen atoms are near enough to each other, and also because the layer is partly dissociated which means there are OH groups free to receive the hopping proton from a neighbouring water molecule. Thus it is possible to imagine further metastable configurations from figure 3.10, and these were observed frequently in our 300K simulation. Such frequent proton transfer at an oxide surface has also been observed in simulations of ZnO . There are two peaks on each side of the histogram. The peak closest to the surface counts the proton transfers within the adsorbed layer, parallel to the surface. The second peak is less intense, and corresponds to hydrogen bonds in the z direction. As shown in the third snapshot in figure 3.3 this corresponds to proton transfer events between the first adsorbed layer and the second layer. It has been recently shown that strong external electric fields promote molecular dissociation and proton hopping in bulk water , ice [89, 90], liquid methanol , as well as a strong chemical reactivity in molecular mixtures. In this case we notice that proton transfer occurs in the first and second layers above the surface corresponding to the zones of highest electric field.
Proton transfer in water has also been shown to be correlated to short oxygenoxygen distance. For instance in the simplest case, the Zundel ion, there is a threshold of oxygenoxygen distance under which the excess proton prefers to lie in the centre of the two oxygen atoms. In the DFT framework this distance is found to be 2.5 Å although more precise treatments of the problem finds this distance to be closer to 2.4 Å . Classically, proton transfer is expected to occur when the oxygen atoms are closer than this minimum distance, apart from processes of tunnelling, where a quantum treatment of the proton is needed . Figure 3.17 shows distributions of the oxygenoxygen distances for hydrogen bonded water molecules, in the first and second layers, and in the other water molecules. With regard to proton transfer, it is important to note that the OO distances concerned are those which are below the threshold mentioned above. It is clear that there are more very short OO distances near the surfaces than in the middle of the simulation box. The layers nearest the surfaces have more than 10% of their OO distances below the DFT protonhopping threshold, while those in the next slice have between 5% and 10%. Although some bulk water properties seem to be recovered just after the first layer, the heightened density seems to permeate through to the second layer and promote proton transfer. Although proton transfer is more likely when there are more OO distances well below the threshold, this is not the only driving force for the transfer. Thus the percentage of OO distances lying below the threshold is not linearly representative of the number of hops.
Table of contents :
1.1 Prebiotic Context
1.2 Mineral Surfaces and Confined Water
1.3 The Electric Field
2 Theoretical Background
2.1 Ab Initio Molecular Dynamics
2.1.1 Molecular Dynamics
2.1.2 Density Functional Theory
2.1.4 Applying an Electric Field via The Modern Theory of Polarisation
2.2 Accessing the Electric field
2.2.1 Application to a model surface
2.3 Enhanced Sampling Techniques
2.3.1 Parametrisation of the Free Energy Landscape
2.3.3 Umbrella Sampling
3 The Electric Field at the Surface of MgO
3.1 Introduction and Context
3.2 Static Results on Bare MgO
3.2.2 MgO(111) and MgO(015)
3.3 One monolayer of H2O on MgO(001)
3.4 Dynamic Results on the MgO/water Interface at 300K
3.5 Mimicking the Surface Electric Field
3.5.1 Description and implementation
3.5.2 Tests and Remarks
4 Surface Effects on a Prebiotic Reaction
4.1 Introduction and Context
4.2 Preliminary Study of the Configuration Space
4.3 Simulation Procedure
4.4 Temperature and Electric Field Effects
4.5 Effect of the Mineral Surface
4.6 Further Analysis
4.6.1 Entropic vs Enthalpic Contributions
4.6.2 Transition States and Isomers
5 Conclusions and Perspectives
A Simulation Details
A.1 Computational Details
A.1.1 AIMD with Quantum Espresso
A.1.2 Metadynamics and Umbrella Sampling with the Plumed Plugin
A.2 System Descriptions