UPTEC F 19035
Examensarbete 30 hp Juni 2019
Structural integrity of highly ionized peptides
Ibrahim Eliah Dawod
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0
Postadress:
Box 536 751 21 Uppsala
Telefon:
018 – 471 30 03
Telefax:
018 – 471 30 00
Hemsida:
http://www.teknat.uu.se/student
Abstract
Structural integrity of highly ionized peptides
Ibrahim Eliah Dawod
In order to understand the behaviour and function of proteins, their three dimensional structure needs to be known. Determination of macro-molecules’
structures is done using X-ray diffraction or electron microscopy, where the resulting diffraction pattern is used for molecular reconstruction. These methods are however limited by radiation damage.
The aim of this work is to study radiation damage of peptides in proteins using computer simulations. Increased understanding of the atomic and molecular dynamics can contribute to an improvement of the method of imaging biological molecules. To be able to describe the processes that take place as accurately as possible, the problem must treated quantum mechanically.
Thus, the simulations are performed with molecular dynamics based on first principles.
In order to capture the dynamics of the excited states of the molecule when exposed to X-rays, time-dependent density functional theory with delta self-consistent field is used. These simulations are compared to ground state simulations. The results of the thesis conclude that the excited and ground state simulations result in differences in the dynamics, which are most pronounced for lager molecules.
ISSN: 1401-5757, UPTEC F 19035 Examinator: Tomas Nyberg
Ämnesgranskare: Nicusor Timneanu
Handledare: Oscar Grånäs and Carl Caleman
Populärvetenskaplig sammanfattning
Proteiner finns i alla biologiska system och bidrar med många viktiga funktioner. Aminosy- ror är de biomoleklyer som bygger upp proteiner genom peptidbindningar, där den tre- dimensionella strukturen bestämmer den biologiska funktionen. Strukturbestämning av biomolekyler görs främst med hjälp av röntgendiffraktion eller elektronmikroskopi, som resulterar i bilder som används för molekylär rekonstruktion. Dessa metoder är dock be- gränsade av strålskadan som uppkommer vid interaktionen. Nya röntgenlasrar har medfört en revolution inom strukturbestämning de senaste 10 åren. Dessa lasrar kan producera intensiva pulser som är 10-100 femtosekunder (fs) långa med över 1012fotoner. Den kort- variga pulsen gör att en bild kan fås, innan molekylen förstörs på grund av den joniserande röntgenstrålen.
Detta arbete har som mål att studera strålskadan hos peptider i proteiner genom datorsimuleringar. Ökad föreståelse för atom- och molekyldynamiken kan bidra till en förbättring av metoden för avbildning av biomolekyler. För att kunna beskriva processerna som sker så noggrant som möjligt måste problemet beskrivas kvantmekaniskt. Därmed utförs simuleringarna med molekyldyanmik baserad på första principer. Jonisering av molekylerna introduceras dynamiskt, där informationen om när detta ska ske tagits från en databas med simuleringar av strålskadan hos proteinet lysozym.
Contents
1 Introduction 6
2 Theory 7
2.1 The biomolecules . . . . 7
2.2 Density functional theory . . . . 8
2.3 Pseudopotentials . . . . 10
2.4 Born–Oppenheimer molecular dynamics . . . . 11
2.5 Time-dependent density functional theory . . . . 11
2.6 Photon-matter interaction . . . . 11
2.7 Excited states . . . . 12
2.8 Hirshfeld charge . . . . 12
2.9 Bond-breaking . . . . 13
3 Methodology 14 4 Results 15 4.1 Glycine . . . . 15
4.2 Diglycine . . . . 16
4.3 Dialanine . . . . 20
5 Discussion 25
6 Conclusion and outlook 26
References 28
Abbreviations
BO Born-Oppenheimer
CG Conjugate gradient DFT Density functional theory
fs Femtosecond
GGA Generalized gradient approximation HOMO Higher occupied molecular orbital LDA Local density approximation
LUMO Lower unoccupied molecular orbital MD Molecular dynamics
Non-LTE Non-local thermodynamic equilibrium SCF Self-consistent field
SFX Serial femtosecond crystallography SPI Single particle imaging
TDDFT Time dependent density functional theory XC Exchange-correlation
XFEL X-ray free electron laser
1 Introduction
Proteins are essential in all biological systems, providing several different functions in cells. They are poly-peptides, built by amino-acids through peptide bonds, where the three dimensional structure decides the biological function. Solving the structure is therefore vital for increasing understanding of the behaviour and function of proteins and for guiding the design of new drugs [1]. Structure determination of macro-molecules is mostly done using X-ray diffraction or electron microscopy. These methods are however limited by radiation damage induced by the probe, due to multiple ionizations in the sample [2, 3]. The dominant cause of ionization with X-rays is from electronic transitions occurring from the absorption of photons [4].
In X-ray diffraction, the target molecule is crystallized and irradiated by X-rays, resulting in a diffraction pattern used for molecular reconstruction [1]. X-ray Free Electron Lasers (XFELs) have furthered this method of structure determination. The intense laser pulses from an XFEL are 10-100 femtoseconds (fs) long with over 1012 photons. This allows a diffraction pattern to be obtained before destroying the structure. The short duration of the pulse reduces the effect of radiation damage, enabling smaller protein crystals and experiments at room temperature. This new technique has been named Serial Femtosecond Crystallography (SFX). In SFX experiments, protein crystals of nano-meter size are provided serially using liquid jets to the area where the X-ray beam is focused. Each diffraction occurs with a single crystal in a random orientation, making it possible to acquire the entire reciprocal lattice [4, 5].
All proteins are however not easily crystallized. This motivates the need for a method of structure determination independent of periodicity. For non-crystalline samples, Single Particle Imaging (SPI) is used. In SPI, single biomolecules are exposed to X-rays [4].
This thesis aims to aid the method of imaging proteins by studying how the structural integrity of peptides in proteins are affected by ionization caused by intense femtosecond X-rays. An understanding of the dynamics of the fragmentation process and knowledge of the timescales of bond-breaking in the biomolecules will allow for an increased accuracy in protein imaging due to an improved speckle analysis.
In order to accurately describe bond-breaking of molecules, the process needs to be treated quantum mechanically. Hence, ab-initio molecular dynamics is performed. In order to capture the coupling between the excited electrons and ionic motion, a recent extension of the Siesta software is used [6, 7]. It is able to describe the ionization process induced by the probe and the excited electronic structure observed. An already developed platform used to study dissociation pathways of amino-acids will be used and adapted in order to study the ionized peptides and the impact of localized excited states [8]. Ionization of the molecules will be dynamically introduced as time evolves with data obtained from non-local thermodynamic equilibrium (non-LTE) plasma simulations of a representative molecule (Lysozyme) exposed to femtosecond X-rays [9, 10].
2 Theory
2.1 The biomolecules
Amino-acids are organic molecules, distinguished by the atoms in the side-chain R. They form larger molecules (polypeptides) through the production of a water molecule (condensation) and a peptide bond between a carbon (C) and nitrogen (N) atom, each provided from separate amino-acids [1, 11].
The biomolecules of interest in this thesis can be seen in figure (1). They are Glycine (GLY), Alanine (ALA), Diglycine (Di-GLY) and Dialanine (Di-ALA). The side-chain of
GLY is a hydrogen (H) atom, and ALA contains a methyl (CH3) group.
(a) Glycine. (b) Alanine.
(c) Diglycine. (d) Dialanine.
Figure 1: Structures of the biomolecules (Oxygen (O) - red, Carbon (C) - dark grey, Nitrogen (N) - blue and Hydrogen (H) - light grey).
2.2 Density functional theory
The dynamics of a many-body non-relativistic quantum mechanical system is described by the Schrödinger equation [12]
HΨˆ a(R1, . . . , RK; r1σ1, . . . , rNσN) = EaΨa(R1, . . . , RK; r1σ1, . . . , rNσN), (1) where the number of nuclei is K and R are the coordinates. The number of electrons is N and r, σ are the spatial coordinates and spin of the electrons. The many-body nature of most problems makes solving equation (1) computationally demanding. Molecular systems require the treatment of the coupling between electrons and nuclei. However, due to the large difference in mass between nuclei and electrons, the Born-Oppenheimer (BO) approximation can be applied. Here, the nuclei are considered to be stationary in the time it takes for the electrons to reach their ground state configuration. The many-body wave-function in equation (1) can therefore be a factor of a nuclear and electronic wave-function. The dynamics of the electrons are then determined by the electronic part of the Hamiltonian ˆHe
HˆeΨek(R1, . . . , RK; r1σ1, . . . , rNσN)
= Ek(R1, . . . , RK)Ψek(R1, . . . , RK; r1σ1, . . . , rNσN).
(2)
The basis of density functional theory (DFT) lies in the Hohenberg-Kohn theorems [13].
They conclude that the energy of a poly-atomic system can be calculated using the electronic density in the ground state n0(x, y, z),
E0 = E[n0]. (3)
The electrons’ contribution to the total energy of a molecule can then be formulated as E0 = hT [n0]i + hVN e[n0]i + hVee[n0]i. (4) with the first term being the kinetic energy of the electrons and the last two the energy contribution of the nuclei-electron and electron-electron interactions.
Kohn-Sham (KS) proposed that the correct ground state energy and density could be determined using an auxiliary reference system with electrons having no interaction between each other. This is motivated by the large computational task of determining the interaction between all electrons. The errors in the observables between the real and imaginary system are
∆hT [n0]i = hT [n0]ireal− hT [n0]iref, hVee[n0]i = hVee[n0]ireal− hVee[n0]iref.
(5)
The real kinetic and potential energy are then determined as, hT [n0]ireal= ∆hT [n0]i + hT [n0]iref,
hVee[n0]ireal= ∆hVee[n0]i + hVee[n0]iref. (6)
The Coulombic repulsion term for the reference system is formulated as hVee[n0]iref = 1
2
Z Z n1(r1)n2(r2)
r12 dr1dr2. (7)
The terms in the integral describe the electrostatic interaction between the total charge n1(r1)dr1 in the volume dr1 and n2(r2)dr2 in dr2, with r12 being their separation.
For the unknown ∆-terms an exchange-correlation (XC) functional is introduced Exc[n0] = ∆hT [n0]i + ∆hVee[n0]i, (8) which captures the dynamics of the interactions missing in the reference system. The second term in equation (4) is
hVN ei = Z
no(r)v(r)dr (9)
where v(r) is the attractive potential from the nuclei felt by the electrons. The total energy is then
E0 = Z
no(r)v(r)dr + hT [n0]iref +1 2
Z Z n1(r1)n2(r2) r12
dr1dr2+ Exc[n0]. (10) By finding the correct XC functional and ground state density, the energy would correspond to the real values. Different XC functionals have been proposed, with the simplest being Local Density Approximation (LDA) and Generalized Gradient Approximation (GGA). The basis of LDA is that the system studied should have local properties determined by a uniform electron gas. The XC term in LDA is therefore only a functional of the electronic density. For GGA, the XC functional is in addition to the density, also a functional of the gradient.
The wave function for the electronic ground state can, in the non-interacting system, be formulated as a Slater determinant [14]
Ψk = 1
√N !
ψ1(r1) ψ2(r1) . . . ψN(r1) ψ1(r2) ψ2(r2) . . . ψN(r2)
... ... . .. ... ψ1(rN) ψ2(rN) . . . ψN(rN)
, (11)
where each orbital satisfies the Kohn-Sham equation
"
−1
2∇2i − X
ion A
ZA r1A +
Z n(r2)
r12 dr2+ vxc(r)
#
ψiKS = KSi ψKSi , (12)
obtained through the minimization of the energy with respect to the KS orbitals. The exchange-correlation potential vxc is defined as
vxc[n](r) = δExc[n]
δn(r) . (13)
For solutions of the the KS equation, the kinetic energy term in (10) can be formulated with the KS orbitals [14],
T [n0] =
N
X
i=1
Z
d3rψ∗i(r)
−∇2 2
ψi(r). (14)
The density is determined by the states obtained in equation (12). For the non-interacting system it can be written as [15]
n(r) =
N
X
i=1
|ψi(r)|2. (15)
The KS molecular orbitals are expanded in terms of the basis functions φ
ψKSi =
m
X
s=1
csiφs. (16)
The functions are found by solving the KS equation for the free atoms with identical pseudopotential and XC-functional. Localization of the resulting functions is imposed, making them zero after some radius. The consequence of this confinement is that the energy is increased. In Siesta, the parameter energy shift quantifies the energy change for all basis functions. The number of functions for each angular momentum quantum number can be varied, where increasing it provides additional variational freedom in the energy optimization. Polarized orbitals are possible through the resulting solution of the KS equation when an electric field is applied to the pseudo-atom [16].
The target density in (12) is obtained by first solving the equation with an initial guess for the density. This comes from either a previous run, or a constructed density from the atoms. This results in KS orbitals which give a more accurate density. The process is repeated until the input (ni) and output (ni+1) density has a difference lower than a given tolerance [12]. In order to reduce the number of iterations, several convergence acceleration schemes exist. For this work, the basis has been linear mixing, where the input for the next iteration is a sum of both the in-and output density according to
ni+1in = αniout+ (1 − α)niin. (17) The value of α depends on the problem. For systems that are difficult to converge, values used can be in the interval [0.1, 0.001]. Larger values means that there is a higher risk for oscillation [17]. In order to further accelerate the convergence and reduce oscillations, Pulay mixing can be used. In this scheme, the values of several previous iterations are utilized to determine the new density.
2.3 Pseudopotentials
Applying a theory which treats all the electrons is computationally demanding. Furthermore, since the valence and core states need to be orthogonal, the valence orbitals are non-smooth in the core region. Instead, the approximation of treating the core electrons as stationary with respect to displacements of the valance electrons is used. A new smoother state is formulated as a linear combination of the core and valence states, where orthogonality
between the valence and core is imposed. An operator is introduced accounting for the effect of the nuclei and the core electrons on the valence states in a region defined by a radius, rc. Beyond this radius, the valence orbitals have the same form as the orbitals in the theory treating all electrons. Furthermore, the special type of pseudo-potentials used require the pseudo-states to be norm-conserving in the interval between 0 and rc[18].
2.4 Born–Oppenheimer molecular dynamics
Running DFT calculations on a molecular system and acquiring the electronic ground states enables the use of Born–Oppenheimer molecular dynamics (BO-MD). The position Ri of nuclei i with mass Mi is propagated in time with classically formulated forces
MiR¨i= Fi, (18)
determined from the Hellmann-Feynman (HF) theorem. Given normalized eigenstates of the electronic and nuclei-nuclei Hamiltonian ˆH and the nuclei positions defined by Φ, the HF theorem reads [12]
F = −dE dΦ = −
* ψΦ
d ˆHΦ dΦ
ψΦ +
. (19)
If the the basis set in (16) is dependent on atomic positions, the force in equation (19) will have extra terms called Pulay forces [19]. Furthermore, non-adiabatic coupling terms are present for systems in excited states, when the BO approximation no longer applies.
2.5 Time-dependent density functional theory
The time-dependent dynamics of the molecules in excited states is described by time- dependent density functional theory (TDDFT). In these systems, the coupling between the nuclei and electrons will be non-adiabatic, which is described by one of the terms in the force in equation (18). The time propagation of the nuclei is done using Ehrenfest dynamics. The electrons are still seen as to have their quantum mechanical properties, and their dynamics must therefore be evolved in time using the time-dependent version of the KS equation [6],
i∂ψn(t)
∂t = ˆHKS[n](t)ψn(t) (20)
n(r, t) =
N
X
i=1
|ψn(r, t)|2 (21)
2.6 Photon-matter interaction
The photon-matter interaction is governed by several possible processes. Elastic scattering is a process where the energy of the incoming photon is same as the outgoing. For inelastic scattering, energy is lost to the sample and the outgoing photon has a different wavelength.
The photon energy can be absorbed by an electron, resulting in ionization. The hole which is left tends to be filled quickly by an electron in a higher energy state. The de-excitation process can occur through the release of either an electron in a higher energy state, or a photon [4]. Since pseudo-potentials without any holes are used, the core holes are not present in the simulations. The assumption is that the lifetime of the hole is short and the
electrons fill the hole faster than the effects on the system. All the photon-matter processes are assumed to have happened at the particular value of average ionization per atom.
2.7 Excited states
In order to capture the dynamics of the excited states of the system when exposed to an intense XFEL pulse, TDDFT in combination with delta self-consistent field (∆-SCF) is used. The total charge of the system is increased using the net charge command, which removes an electron with the lowest binding energy. The molecules are then excited using
∆-SCF. The electron density is altered by taking an electron occupying a higher occupied molecular orbital (HOMO) and making it occupy a lower unoccupied molecular orbital (LUMO) according to
n(r) =
N −1
X
i=1
ψi∗(r)ψi(r) + ψ∗a(r)ψa(r). (22) The orbital ψa(r) is the LUMO state, which could be obtained by running a ground state simulation [20]. The final density is found self-consistently given the condition that the electron is forced to stay in the excited (LUMO) state. Susceptibility for bond-breaking is dependent on whether the excited electron is in a bonding or anti-bonding orbital at the particular ionization stage. Due to the increase of charge between the atoms, the former provides attraction between the atoms. The latter induces dissociation due to the reduced charge density compared to when the atoms are apart [21]
2.8 Hirshfeld charge
The atomic charge is calculated using the Hirshfeld charge [6, 22]. This approach is less dependent on the basis expansion used, compared to other charge partitioning methods.
Each atom has its fraction of the density of the molecule determined, when no charge transfer has occurred. The reference density n(r) of the molecule is defined by taking the sum over all atoms i,
n(r) =X
i
ni(r) (23)
and the weight of atom i is
wi(r) = ni(r)
n(r). (24)
As the self-consistent charge density nSCF(r) is obtained, the charge density nBA(r) of an atom i is then
nBA(r) = wi(r)nSCF(r). (25)
The Hirshfeld charge is then determined, given the atomic number Zi, as Qi= −
Z
nBA(r)dV + Zi. (26)
2.9 Bond-breaking
The timescale of bond-breaking is determined by studying the value of the chemical bond integrity BI defined as [8]
BI(A, B) = 1 NM D
NM D
X
i=1
1 + eλ(|di[A,B](t)−di[A,B](0)|−0.5)−1
. (27)
The term di[A, B](t) defines, at time t, the separation between the atoms A and B. The distance between the atoms at t = 0 can either be from the relaxed structure or taken form a MD simulation. The number of sampled structures is NM D, and λ is a constant. B ≈ 1 is the value of the bond integrity in the beginning of the simulation, when none of the atoms have moved from the initial positions. When B ≈ 0, the average distance between two atoms is large relative to the initial configuration, concluding that the bond has been broken.
3 Methodology
The results are products of computer simulations using the Siesta software and an extension to it [6, 7]. The structures of the molecules were relaxed using the Conjugate Gradient (CG) algorithm in Siesta. Several starting structures were generated for each biomolecule by running molecular dynamics (MD) simulations. Data with average ionization per atom as a function of time for the protein Lysozyme was extracted from the database FreeDam [10]. This data was taken for several different XFEL pulse parameters. The ionization was multiplied with the number of atoms for each biomolecule and a polynomial fit was done.
Python scripts were written to run the simulations. Each structure was simulated according to the chosen XFEL pulse length and the ionization was introduced at times guided by the polynomial fit. The times acquired from the fit had to be rounded to the closest integer and divided by the time-step of the MD simulation, in order to get the number of steps.
Ionization stages which lasted less than 5.0 fs (time scale of the period of the hydrogen atom) were added until the total time was approximately 5.0 fs. The average of the ionization stages during this total time was used. For each excited state, the indices of the HOMO and LUMO states were randomly chosen. The simulations were done with a time-step of 0.5 fs, mesh cutoff 150 Ry and energy shift 0.01 Ry. The XC-functional was GGA and double-zeta-polarized (DZP) orbitals were used for the expansion of the KS orbitals. The Hirshfeld charge was tracked for each MD step. Simulations were made using both ∆−SCF and ground state simulations, comparing the end results given the same starting structure.
Since difficulties with converging the system in the excited state arose, the parameter space for accelerating the convergence had to be thoroughly explored. This testing was automated by utilizing parallel processing to run several instances of the software with different mixing parameters. All simulations were made on high performance computing clusters like UPPMAX and Tetralith. Post-processing of the output files from the TDDFT simulations was made using the getcube tool and available Perl scripts were used to manipulate the resulting files. The charge densities were visualized with the software Vesta [23]. Available scripts were used in order to calculate the bond-integrity, and for extracting the atomic charges.
4 Results
The values for the XFEL pulse parameters presented in this section are taken from FreeDam [10]. They correspond to real experimental parameters that have already been achieved at XFEL facilities around the world. The following results show the bond-integrity BI defined by equation (27) as a function of time during the XFEL pulse. They correspond to various bonds in the molecule studied. The molecule and bond in question is shown in the figure.
For the dipeptides, charge density differences and Hirshfeld charges are also shown. For all simulations, the sample has been lysozyme.
4.1 Glycine
For GLY the pulse parameters used are shown in table (1).
Pulse length: 20 fs Photon energy: 9000 eV Radiant fluence: 1e6 J/cm2 Table 1: XFEL pulse parameters.
(a) Bond N0-C1. (b) Bond C1-C2.
(a) Bond C2-O3. (b) Bond C2-O8.
Figure 3: Bond-integrity for different bonds in Glycine.
4.2 Diglycine
For Diglycine two intensities have been compared during 20 fs using the parameters shown in table (2). The parameters to the left have the same total number of photons as the right, but the photons are delivered to the target in a shorter interval. This results in 5 times higher power density (W/cm2) for the left case, compared to the right case.
Pulse length: 20 fs Photon energy: 9000 eV Radiant fluence: 1e6 J/cm2
Pulse length: 100 fs Photon energy: 9000 eV Radiant fluence: 1e6 J/cm2 Table 2: XFEL pulse parameters.
(a) Peptide bond N8-C2. (b) Peptide bond N8-C2.
(c) Bond N8-C9. (d) Bond N8-C9.
Figure 4: Bond-integrity for different bonds using high (figures a and c ) and low intensity pulse (figures b and d).
(a) Bond C1-C2. (b) Bond C1-C2.
(c) Bond C9-C10. (d) Bond C9-C10.
(e) Bond O11-C10. (f) Bond O11-C10.
Figure 5: Bond-integrity for different bonds using high (figures a, c and e) and low intensity pulse (figures b, d and f).
(a) Nitrogen (N8). (b) Carbon (C2).
(c) Carbon (C9). (d) Carbon (C1).
(e) Carbon (C10). (f) Oxygen (O11).
Figure 6: Hirshfeld charges calculated according to equation (26) for different atoms in Diglycine. This value indicates how much the atom differs in charge compared to its neutral state. Positive means that the number of electrons is less than the neutral atom, while negative means that there are more electrons. The curves correspond to ∆-SCF and ground state simulations.
(a) Net charge 10. (b) Net charge 21.
(c) Net charge 25.
Figure 7: Isosurfaces for Diglycine with the charge density difference between ∆−SCF and ground state simulations, at each new ionization stage for the higher intensity pulse. The scattered spheres correspond to ionized hydrogen (H+). Yellow surface indicates increase of charge, while blue indicates a decrease. The positive and negative charge densities are shown for values above a certain threshold for the isosurface. For (c), this value is ±0.00164461.
4.3 Dialanine
Pulse length: 20 fs Photon energy: 9000 eV Radiant fluence: 5e5 J/cm2
Pulse length: 20 fs Photon energy: 9000 eV Radiant fluence: 1e6 J/cm2 Table 3: XFEL pulse parameters.
(a) Bond C1-C2. (b) Bond C1-C2.
(c) Bond C1-C6. (d) Bond C1-C6.
Figure 8: Bond integrity for different bonds in Dialanine. The fluence in the right case corresponds to twice of the left case.
(a) Bond N11-C2 (b) Bond N11-C2.
(c) Bond N11-C12. (d) Bond N11-C12.
(e) Bond C12-C17. (f) Bond C12-C17.
Figure 9: Bond-integrity for different bonds in Dialanine.
(a) Bond C12-C13. (b) Bond C12-C13.
(c) Bond C2-O3. (d) Bond C2-O3.
(e) Bond N0-C1. (f) Bond N0-C1.
Figure 10: Bond-integrity for different bonds in Dialanine.
(a) Nitrogen (N11). (b) Carbon (C2).
(c) Carbon (C1). (d) Nitrogen (N0).
(e) Carbon (C12). (f) Carbon (C13).
Figure 11: Hirshfeld charge for various atoms in Dialanine.
(a) Net charge 13. (b) Net charge 29.
(c) Net charge 34.
Figure 12: Isosurfaces for the difference between ∆-SCF and ground state simulations for Dialanine, at each new ionization stage for the pulse with the highest fluence.
5 Discussion
∆−SCF results in different dynamics in the molecular system compared to the ground state simulations. This can be seen in the charge density difference plots, which show charge transfer. The result is that the dynamics of the bonds are altered. Some bonds might be more unstable due to the decrease in charge, while others can become more stable due to excess of charge. Whether the ground or excited state simulation provides more resistance against dissociation is dependent on how the molecular orbitals of the peptide look. It is possible that the excited state is to a bonding orbital (instead of an anti-bonding), resulting in a more stable configuration. Additionally, the results show that size of the target molecule provides more stability. By the end of the pulse, GLY has all of its bonds broken. While for the same pulse parameters Di-GLY has several stable bonds.
The bond of interest is the peptide bond, which is keeping the amino-acids together, thus building up the poly-peptide. For this bond in Di-GLY (N8-C2) in figures (4a) and (4b), both the higher and lower intensity pulses predict different behaviours for the two simulation methods. For oscillatory behaviour present in Di-GLY, the difference shows primarily in the period and amplitude. For Di-ALA, the bond N0-C1 in figure (10f) is stable in ∆-SCF, while for the ground state simulation it is predicted to dissociate. This is particularly interesting since for both GLY and its dipeptide, the two methods predict similar dynamical behaviour with differences in the parameters describing the particular dynamics (e.g. period in oscillations). However, in the case of Di-ALA, the dynamics for the bond N0-C1 is significantly different. Since Di-ALA is larger with 23 atoms, compared to the 17 of Di-GLY, the difference seems to be due to the size.
For both dipeptides the charge density differences in figures (7) and (12) are not localized in most results, but instead spread over the whole molecule. This is due to the small sizes of the molecules. For figure (7c) the differences seem to be mostly translated to the right parts of the molecule. This might be indications of strong electrostatic forces due to high charge values for the atoms (relative to their initial charge values) at the right parts of the molecule. Whether the excited state provides increased charge for an atom is dependent on where the electron is excited from. Charge transfer is further confirmed by the fact that the Hirshfeld value changes as time evolves after ionization. An example is figure (6), where at ionization stage 25, the carbon (C10) has a decreasing Hirshfeld charge and carbon (C1) an increasing.
6 Conclusion and outlook
Simulations using ∆−SCF and ground state show different dynamical behaviours. For larger systems, the differences seem to be more pronounced, for instance in the the case where
∆−SCF suggests a stable bond and ground state does not. Since the infrastructure to run the simulations and automation of the convergence testing is available, it can easily be used to study larger molecular systems. This would enable the study of systems where the molecular orbitals are not as spread over the whole molecule. Furthermore, even though the results show the possible dynamics of the system, the statistics acquired in this work are not sufficient and one would ideally average over multiple simulations with different starting structures. The low statistics is a consequence of the system being difficult to converge.
With the knowledge acquired in this work, one could continue running convergence tests, which would result in sufficient converging systems with time. Finally, in order to verify the results experimentally, the fragments would need to be detected. However, the fragments which are still intact in the excited state after the duration of the pulse might have changed after the decay to the ground state.
The next step in this project in terms of exploring new dynamics, would be to introduce a core hole in the pseudo-potentials of the atoms. This is due to there being suggestions that the lifetime is not short enough to neglect. Additionally, it would be interesting to run classical MD and compare the bond-integrity with the simulation methods used here.
Acknowledgement
I’d like to thank Oscar and Calle for all the help during this project and Nic for giving great feedback on the report. Thanks to all the nice people I got to share the server room with at materials theory. To all the biophysics/molkond people, thanks for nice monday meetings and several enjoyable social events. Finally, thanks to all of my friends and family around the world.
References
[1] Bernhard Rupp. Biomolecular crystallography: principles, practice, and application to structural biology. Garland Science, 2010.
[2] Richard Henderson. The potential and limitations of neutrons, electrons and x-rays for atomic resolution microscopy of unstained biological molecules. Quarterly Reviews of Biophysics, 28(2):171–193, 1995.
[3] Henry N. Chapman, Carl Caleman, and Nicusor Timneanu. Diffraction before de- struction. Philosophical Transactions of the Royal Society B: Biological Sciences, 2014.
[4] Sébastien Boutet, Petra Fromme, and Mark S. Hunter. X-ray Free Electron Lasers, A Revolution in Structural Biology. Springer International Publishing, 2018.
[5] Henry N. Chapman et al. Femtosecond x-ray protein nanocrystallography. Nature, 2011.
[6] Grigory Kolesov, Oscar Grånäs, Robert Hoyt, Dmitry Vinichenko, and Efthimios Kaxiras. Real-time td-dft with classical ion dynamics: Methodology and applications.
Journal of Chemical Theory and Computation, 12(2):466–476, 2016. PMID: 26680129.
[7] José M. Soler et al. The siesta method for ab-initio order-n materials simulation. J.
Phys.: Condens. Matter, 2002.
[8] Oscar Grånäs et al. Femtosecond bond-breaking and charge dynamics in ultracharged aminoacids. 2019.
[9] Howard A. Scott. Cretin - a radiative transfer capability for laboratory plasmas.
Journal of Quantitative Spectroscopy and Radiative Transfer, 71(2):689 – 701, 2001.
Radiative Properties of Hot Dense Matter.
[10] H. Olof Jönsson, Christofer Östlin, Howard A. Scott, Henry N. Chapman, Steve J.
Aplin, Nicuşor Tîmneanu, and Carl Caleman. Freedam - a webtool for free-electron laser-induced damage in femtosecond x-ray crystallography. High Energy Density Physics, 26:93 – 98, 2018.
[11] Engelbert Buxbaum and SpringerLink (Online service). Fundamentals of Protein Structure and Function. Springer International Publishing, Cham, 2nd 2015.;2nd 2015;
edition, 2015.
[12] Eberhard Engel and Reiner M. Dreizler. Density Functional Theory: An Advanced Course. Springer-Verlag, 2011.
[13] Errol G Lewars. Computational Chemistry: Introduction to the Theory and Applications of Molecular and Quantum Mechanics. Springer International Publishing, 2016.
[14] Carsten Ullrich. Time-dependent density-functional theory: concepts and applications.
Oxford University Press, New York;Oxford;, 2012;2011;.
[15] Robert G. Parr and Weitao Yang. Density-functional theory of atoms and molecules.
Oxford University Press, 1989.
[16] Javier Junquera, Óscar Paz, Daniel Sánchez-Portal, and Emilio Artacho. Numerical atomic orbitals for linear-scaling calculations. Physical Review B, 64(23), 2001.
[17] Jeffrey R. Reimers. Computational methods for large systems: electronic structure approaches for biotechnology and nanotechnology. Wiley, Hoboken, N.J, 2011.
[18] David J. Singh and Lars Nordstrom. Planewaves, pseudopotentials and the LAPW method: Second edition. 2006.
[19] E. Artacho, D. S nchez Portal, P. Ordej n, A. Garc a, and J. M. Soler. Linear- scaling ab-initio calculations for large and complex systems. Physica status solidi (b), 215(1):809–817, 1999.
[20] Jeppe Gavnholt, Thomas Olsen, Mads Engelund, and Jakob Schiøtz. ∆ self-consistent field method to obtain potential energy surfaces of excited molecules on surfaces. Phys.
Rev. B, 78:075441, Aug 2008.
[21] R. P. Feynman. Forces in molecules. Phys. Rev., 56:340–343, Aug 1939.
[22] F. L. Hirshfeld. Bonded-atom fragments for describing molecular charge densities.
Theoretica Chimica Acta, 44(2):129–138, 1977.
[23] Koichi Momma and Fujio Izumi. VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data. Journal of Applied Crystallography, 44(6):1272–1276, Dec 2011.
