Project Topics
Get Complete Project Material File(s) Now! »
Static electric fields in ab initio simulations
The theoretical development and the implementation of electric fields in ab initio codes is all other than trivial. A wide literature exists in this field [63–67] and it is impossible to resume all the conceptual steps which started since the definition of polarization.
One of the key points that one has to face in dealing with this delicate aspect is the common presence of PBC. The first problem that arises is the inherent discontinuity carried by the infinite replication of the simulation box. At the edges of the boxes infinite electric field strengths clearly occur when a linear electrostatic potential is applied within. However, this problem can be more or less easily handled by confining the MD with some specific constraint even if this solution is not so rewarding.
More seriously, the periodicity in the presence of a macroscopic electric field E will lead to a change in the electron potential in each replication of the simulation box. The intrinsic problem resides in the non-periodic nature of the position operator. In particular, the electric field changes by a factor of eE · R under a traslation by a lattice vector R and even a small field varies the nature of the energy eigenstates. Moreover, because the potential is unbounded from below, the ground state is ill-defined [63, 68, 69].
Many perturbative treatments of the application of an electric field have been proposed but only within the Modern Theory of Polarization and Berry’s phases this delicate issue can be efficiently tackled [67,70]. Few of them [63] are founded on aWannier-function-based solution to the finite-field problem that was not very useful in practice [67]. By means of the use of Bloch’s functions, Nunes and Gonze [64] showed how the common perturbative treatments could be directly obtained from a variational principle based on minimizing an energy functional F of the following form F = EKS({ψkn}) − E · P({ψkn}) .
Collective variables and Metadynamics
As stated before, ab initio methods are by their nature an extremely powerful tool in predicting and analyzing chemical reactions. However, one of the most important parameters in describing activated processes is of course its free energy surface (FES). This latter is represented by a function that tipically depends on a number of variables (i.e., collective variables (CV)) which can be either few and relatively numerous depending on the specific reaction. There exist several theoretical and computational methods that are capable to reproduce the FES of a chemical process and few of them have been employed in this thesis.
Let us consider a system in the canonical ensemble; we introduce a collective variable q(R), a function of the atomic coordinates, that is able to distinguish the relevant metastable states of the system (i.e., reactants and products). The probability of finding the system in a specific configuration characterized by the reaction coordinate s is given by P(s) = 1 Q Z e−U(R) kBT δ(s(R) − s)dNR.
Path Collective Variables with a new definition of distance
By the knowledge of the dynamical evolution of a given chemical reaction – i.e., its trajectory – the problem of the a priori choise of a set of CV can, in principle, be bypassed. Branduardi et al. [35] introduced the so-called “path-CV” which are constructed from a tentative reaction path {R0} connecting the reactants and the products states. By defining a measure D of the distance between the current reaction {R} and the reference one {R0} as, e.g., the mean square difference between all pairs of the interatomic distances computed in the reference states (say A and B for the reactants and the products states, respectively).
Umbrella Sampling and the Weighted (and Dynamic) Histogram Analysis Methods
The sampling probability (1.87) can be modified by virtue of a biased potential added to the system which depends only on the CV. For the sake of simplicity now we will take only one CV but the following treatment can be straightforwardly extended to multi-dimensional FESs. The potential energy will be thus modified in U(R) → U(R) + V (s). Hence the biased distribution of the CV will be P′(s)∝ Z e−U(R)+V (s) kBT δ(s − s(R))dNR∝ e−V (s) kBT P(s) .
Dynamics of H-bond chains
With regard to the local dynamics, when the dissociation events occur we observe a shortening of the H-bond lengths that allows for a fast PT along the respective reaction coordinates, as shown in Fig. 2.6 for the two crystalline phases. In the ice XI case the average oxygen-oxygen distances are slightly shorter than in the hexagonal common phase, which suggests a more efficient transfer process of defects.
We define the H-bond wire length as the sum of the individual donor-acceptor distances be- tween H-bond connected oxygen atoms in which a molecular dissociation has occurred. Figure 2.7-a shows the behavior of an ice Ih H-bond wire, whose length has been rescaled by dividing it by the number of involved intermolecular distances and in which a combined PT between the oxygens has occurred in a few hundreds femtoseconds, followed by a ready proton recombi- nation. A pronounced relative minimum is visible just before the global one indicating a first (unsuccessful) attempt at molecular dissociation and local PT. We observe another curious aspect which characterizes all the investigated H-bond wires: after dissociation, the length of the “wire” always increases with respect to the average one and attaines a global maximum just after the global minimum, showing a sort of chain “breath”; in all cases (even in the ferroelectric one-stage mechanism of PT shown in Fig. 2.7-b) subsequent correlated molecular dissociations, produced immediately later along the same H-bond path, do not occur. This finding is consis- tent with the Jaccard theory of defects propagation [106]. Hence, we argue that the mechanism underlying molecular dissociation, which implies a certain degree of cooperation between the oxygen atoms, is the same in liquid water [22, 103] as well as in the crystalline phases that we have investigated.
Ice Ih vs. ice XI: anisotropy effects
As already noted, the field intensity value that allows for a sustained protonic current in common ice is comparable to that observed in liquid water (0.36 V/˚A) [18, 81, 97] and does not substantially depend on the field direction. The response of ice Ih to static electrical perturbations is intrinsically isotropic on account of the topological isotropy conferred to this ice phase by its proton sublattice. Incidentally, this result recalls the hypothesized isotropic nature of the G correlation parameter (also known as dipolar autocorrelation coefficient) in ice Ih; this finding is also supported by Monte Carlo simulations and is intimately related to the electrical response via the dipole moments distribution [117].
However, when the field is oriented along the b-axis of ice XI, the behavior of the material changes significantly if compared with that of the same structure under an electric field oriented along the c-axis. In these conditions the molecular dissociation events and the subsequent proton diffusion are not simultaneous, as previously noticed for ice XI under the action of a field parallel to the ferroelectric direction. In particular, as the concentration of both types of defects becomes relevant, the system starts to conduct in almost the same way as in the case of ice XI when the field lies along the ferroelectric axis, evolving, from a certain point of view, to a partial (unstable and field-induced) ferroelectricity.
In order to give a more complete description of the ionic conduction in these systems, we also calculate the ionic current-voltage diagrams corresponding to our unit cells. The results are displayed in Fig. 2.8-a for the two phases when the electric field is oriented along the b-axis; Figure 2.8-b/c shows the comparison with the previously investigated cases where the field is oriented along the c-axis for ice Ih and XI, respectively [81]. All the ionic current-voltage diagrams show a ohmic behavior after a net PT process has taken place. These results confirm those found in several experiments performed on ice Ih samples [104, 105]. It is not surprising to observe a good agreement between the two different situations (see Fig. 2.8-b); a similar conclusion follows from the calculations performed with the BLYP functional. However, as far as ice XI is concerned, we observe a different behavior below the current threshold value (Fig. 2.8-c). This result can be interpreted at relatively low voltages (up to about 5 V for our unit cell) by considering that the system under the action of a field oriented along the b-axis, generates a sufficient concentration of defects, up to a field strength of 0.29 V/˚A (ionic current threshold). Above this value, the sample is somewhat indistinguishable, as for the electrical response, from a similar sample under the effect of an electric field acting along the ferroelectric axis. The hypothetic slight differences in the derivatives (i.e., in the conductivities)
Table of contents :
Introduction
1 Theoretical Background
1.1 Density Functional Theory
1.1.1 Hohenberg-Kohn Theorem
1.1.2 The Levy approach and the Kohn-Sham scheme
1.1.3 The exchange and correlation functional
1.1.4 Plane waves and pseudopotentials
1.2 Ab initio simulations
1.2.1 Car-Parrinello Molecular Dynamics
1.2.2 Static electric fields in ab initio simulations
1.3 Collective variables and Metadynamics
1.3.1 Path Collective Variables with a new definition of distance
1.3.2 Umbrella Sampling and the Weighted (and Dynamic) Histogram Analysis Methods
2 Results on Proton Transfer in H-bonded Systems
2.1 Water ices under an electric field
2.2 Liquid methanol under an electric field
3 Results on Methanol Chemical Reactions
3.1 Methanol chemistry at high electric fields
3.1.1 Field-induced chemical reactions
3.1.2 Reactivity and (regio)selectivity
3.1.3 (Preliminary) Free Energy landscape of formaldehyde formation
Conclusions
A Theoretical background on Fukui Functions, (Regio)Selectivity, and Reactivity Descriptors
A.1 Global reactivity descriptors
A.2 Local reactivity descriptors
A.3 Another approach: the Noncovalent Interactions
B Supplemental Results on the Reactivity Descriptors
B.1 Global and (other) local indicators
References
