JOHAN HATTNE
Quantum chemical
study of the UV-induced psoralen–DNA reaction
Master’s degree project
Uppsala University School of Engineering
UPTEC X 04 013 Date of issue 2004-02
Author
Johan Hattne
Title (English)
Quantum chemical study of the UV-induced psoralen–DNA reaction
Title (Swedish)
Abstract
Psoralens comprise a family of small, aromatic, and hydrophobic molecules with the ability to intercalate between base pairs of DNA. Intercalated psoralen, when irradiated with light in the UVA region, can covalently bind to both strands of DNA. Such crosslinks can form loops in DNA which eventually inhibit the replicatory machinery of the cell. This in turn is believed be of relevance in the treatment of certain skin diseases. The present work is a theoretical study of the addition of psoralen to DNA. Several experimental observations are verified in silico, and possible explanations for a few others are offered.
Keywords
Density Functional Theory, DNA, Molecular Mechanics, ONIOM, Photocycloaddition, Psoralen, Quantum Mechanics
Supervisors
Leif A. Eriksson
Dept. of Natural Sciences Orebro University ¨
Jorge Llano
Dept. of Cell and Molecular Biology Uppsala University
Scientific reviewer Johan ˚ Aqvist
Dept. of Cell and Molecular Biology Uppsala University
Project name Sponsors
Language
English
Security
Secret until 2005-03-01 ISSN 1401-2138
Classification
Supplementary bibliographical information Pages
58
Biology Education Centre Biomedical Center Husargatan 3 Uppsala
Box 592 S-75124 Uppsala Tel +46 (0)18 4710000 Fax +46 (0)18 555217
UV-induced psoralen–DNA reaction
Johan Hattne
Sammanfattning
Psoralenerna utg¨or en familj av sm˚a, plana molekyler som kan l¨agga sig mellan basparen i DNA. Denna interkalerade place- ring stabiliseras av icke-kovalenta interaktioner mellan baserna och de cykliska strukturena i psoralen. N¨ar psoralen ¨ar placerat p˚a detta vis och bestr˚alas med ultraviolett ljus kan addukter bildas, vilka korsbinder de tv˚a str¨angarna i DNA. Trots att re- aktionen ¨ar skadlig f¨or den enskilda cellen tros den ha klinisk betydelse i behandlingen av vissa hudsjukdomar. I detta projekt har teoretiska ber¨akningar baserade p˚a kvantmekanik anv¨ants f¨or att studera additionen av psoralen till DNA. Fr˚an studien kan flera experimentella observationer bekr¨aftas och i viss m˚an
¨aven f¨orklaras.
Examensarbete 20p i Molekyl¨ar bioteknikprogrammet
Uppsala universitet Februari 2004
1 Introduction . . . . 3
1.1 Introduction to the Reaction . . . 3
1.2 The Players . . . 4
1.3 Overview of the Reaction Mechanism . . . 6
1.4 Stereochemistry of the Psoralen–DNA Adducts. . . 7
1.5 Motivation . . . 8
2 Molecular Mechanics . . . . 10
2.1 Notational Conventions . . . 10
2.2 Molecular Mechanics – Force Field Methods. . . 11
2.3 Atom Types . . . 14
2.4 Force Field Parametrisation . . . 14
2.5 Geometry Optimisation . . . 15
3 Quantum Mechanics . . . . 18
3.1 Some Basic Postulates of Quantum Mechanics. . . 18
3.2 The Variational Principle . . . 21
3.3 Hartree–Fock . . . 22
3.4 Basis Functions. . . 26
3.5 SCF – Self-Consistent Field . . . 28
4 Semi-Empirical Methods, DFT and ONIOM . . 29
4.1 Semi-Empirical Methods. . . 29
4.2 Density Functional Theory . . . 31
4.3 Time-Dependent DFT . . . 33
4.4 ONIOM . . . 34
5 Methods . . . . 36
5.1 The Model Systems . . . 36
5.2 Geometry Optimisation . . . 37
5.3 TDDFT . . . 38
5.4 Basis Sets . . . 38
5.5 Some Notes on Performance . . . 39
6 Results and Discussion . . . . 40
6.1 Introduction. . . 40
6.2 Geometries . . . 40
6.3 Energies . . . 43
6.4 Excitation Energies . . . 46
6.5 General Discussion. . . 47
6.6 Acknowledgements . . . 51
Bibliography . . . . 55
Index . . . . 56
Introduction
Figure 1.1: Some of the nodes of the local Beo- wulf cluster. This clus- ter was used extensively during the course of this project.
Figure 1.2: The Lucidor cluster at PDC in Stock- holm consists of 90 Ita- nium 2 nodes. Copyright
, c° Harald Barth.
Used with permission.
Computers, such as those shown in Figs. 1.1 and 1.2, are becoming faster and faster. The natural question to ask is
“What should one do with the spare clock cycles that are becoming ever more abundant?”. For certain users, Gates’
law may offer some comfort: “The speed of software halves every 18 months”a. For users of power OS:s, Gates’ law provides little relief. A more general answer to the above posed question is presented here: Quantum Chemistry.
aFrom The Jargon Dictionary, available on the Internet at http://
info.astrian.net/jargon. Gates’ Law may have a tremendous impact on the computing community, if its balancing counterpart, Moore’s Law , is broken. With the technological principles in use today, it is indeed predicted that this will happen within a foreseeable future [1].
---
1.1
Introduction to the Reaction
This project is a study of the reaction of psoralen with DNA. The psoralens comprise an entire family of molecules and are further introduced in section 1.2.2. Their cellular reactions are not limited to DNA; psoralens are also known to react with RNA, proteins, fatty acids and a number of other molecules as well [2, 3]. In general, reactions not involving nucleic acid are minor compared to those that do [4]. The frequency at which these reactions occur is dependent on the particular psoralen involved. Angelicin, for instance, is not able to complete all the steps of the reaction studied here. This is believed to be due to its angular geometry [5].
The most important aspect of the reaction under study is that it is initiated by the absorption of energy in the form of ultraviolet light. This will transform the ground state psoralen, denoted Ps in what follows, via an excited singlet state, 1Ps?, to an activated triplet state, 3Ps? [3]. Transition from 1Ps? to 3Ps? occurs by intersystem crossing, ISC. Only photoactivated psoralen takes part in the reactions studied here. Experimental evidence suggests that
inactivated psoralen truly is inactive: in the absence of UV-light, even high concentrations of psoralen seem to have no effects on cellular function [6].
Normally, activated psoralen reacts by electron transfer to oxygen. There are two kinds of oxygen-dependent reactions, called type I and type II respec- tively, given in Eqs. 1.1a and 1.1b below [7]. In this study we instead focus on a third, oxygen-independent, mechanism, type III, which may occur when substrate and target are close together in space, see Eq. 1.1c. It will become apparent that in our case, the proximity requirement is met.
Ps hν
−−−−→1Ps? ISC
−−−−−→3Ps? Reductant
−−−−−−−−−−→Ps·− O2
−−−−→ . . .
. . . −−→O·−2 Substrate
−−−−−−−−−→Products (1.1a)
Ps hν
−−−−→1Ps? ISC
−−−−−→3Ps? O2
−−−−→1O2
Substrate
−−−−−−−−−→Products (1.1b)
Ps hν
−−−−→1Ps? Substrate
−−−−−−−−−→Products (1.1c)
In type III reactions, the psoralen interacts with the target directly – no intermediate electron transfer to oxygen is involved. This is a consequence of the [2 + 2]-cycloaddition that is believed to accurately describe the mechanism, see section 1.4. In particular this means that no singlet oxygen is produced, and the only modifications to the system brought about by the reaction are those covalent links we wish to study. This makes side-effect analysis significantly simpler.
The remaining text is a more detailed study of the reaction and the theo- retical tools used. The interested reader is referred to textbooks by Solomons and Stryer for basic biochemistry [8, 9]. The computational theory used here is covered by Atkins, Friedman and Jensen [10, 11].
---
1.2
The Players 1.2.1 DNA
O Base H H 1’
H H
3’
O P O
−O O H2C
5’
O Base H H 1’
H H
3’
O P O
−O O H2C
5’
O Base H H 1’
H H
3’
HO H H H2C
5’
O
Figure 1.3: One strand of DNA. The loose bonds at the top and the bottom would be connected to hy- drogen atoms.
DNA is a polymer of deoxyribonucleotide units. A nucleoside consists of a nitrogenous base bonded to a deoxyribose sugar; a nucleotide is a nucleoside attached to one or more phosphate groups. The nitrogenous base is either a purine or a pyrimidine derivative. The purines in DNA are adenine and guanine, the pyrimidines are thymine and cytosine [9]. The chemical structure of these largely planar bases is shown in Fig. 1.4. Only pyrimidine bases seem to be involved in the psoralen–DNA reaction. Of these, thymine is the primary base at which photoreaction occurs [12, 13].
N9 in purines and N1 in pyrimidines bond to C10 in deoxyribose to form nucleosides. These sugars, when linked by a phosphodiester bridge between the 30− hydroxyl group of one deoxyribose and the 50− hydroxyl group of another, form the backbone of DNA. Backbone and bonded bases together make up one strand of DNA. Note that a strand of DNA has polarity: one end of the strand has a 50− OH group while the other has a 30− OH group, neither of
N
3
H
2
N
1N
H26 5
4 8 H
N R
9
N
7
(a) Adenine.
N R
1
O
2
N
3N
H24 H
5
H
6
(b) Cytosine.
N
3
H2
N
2
H
N
1
O
6 5
4 8 H
N R
9
N
7
(c) Guanine.
N R
1
O
2
H
N
3
O
4
C
H25
H
6
(d) Thymine.
Figure 1.4: The bases of DNA: adenine, cytosine, guanine and thymine. In basic form, we have R = H
for all structures above. When bonded to deoxyribose in DNA, C10of deoxyribose takes the place denoted
by R.
which is linked to another nucleotide, see Fig. 1.3. Base sequences are written – and read – from the 50-end to the 30-end.
A DNA helix consists of two strands. The bases are on the inside of the structure, arranged so that the planes of the sugars are approximately orthogonal to the planes of the bases. The helical structure repeats every ten residues. The two chains are held together by hydrogen bonds between pairs of bases. The pairing is not at all arbitrary: adenine is always paired with thymine, guanine is always paired with cytosine.
1.2.2 Psoralens
Psoralens, also called furocoumarins, are a family of small, aromatic, hydropho- bic compounds naturally present in plants throughout the world. They easily penetrate both cells and virus particles [6]. All members of the psoralen family are linear fusions of a coumarin, or a dipyrone, and a furan ring – thus pso- ralens are heterocyclic. The particular molecules under study here, psoralen and 40− (hydroxymethyl) − 4, 50, 8 − trimethylpsoralen or HMT are shown in Fig. 1.6. The furan moiety of psoralens is photosensitive in the 320 to 400 nm range of radiation – UVA, see Fig. 1.5 – an interval in which cellular nucleic
acids are weakly absorbing, if at all [12]. λ/ nm
UVC UVB UVA
100 280 320 400
Figure 1.5: Ultraviolet radiation. Visible light lies to the right of UVA in the diagram.
H
8
H H 5’ 5
H
4’
O O
2
O
1
H
4 H
3
(a) Psoralen.
C
H38
H
C
5H3 5’
C
H2
HO
4’
O O
2
O
1
C
H34 H
3
(b) HMT.
Figure 1.6: The structure of psoralen and HMT. The latter is a synthetic pso-
ralen derived from the naturally occurring 4, 50, 8 − trimethylpsoralen (TMP).
---
1.3
Overview of the Reaction Mechanism
The reaction is modelled to proceed through several steps. Experimental evidence suggests that the reaction also occurs in steps conforming to this model. In the final product the psoralen will have covalently crosslinked the two strands of DNA. Isaacs proposed the mechanism in Eq. 1.2 below, where Ps · · · DNA denotes the intercalated structure of psoralen and DNA, Ps − DNA the monoadduct, further explained in what follows, and Ps = DNA the di- adduct or crosslink [4]. A photon of light is denoted hν.
Ps + DNA −−*)−− Ps · · · DNA (1.2a)
Ps · · · DNA hν1
−−−−−*
)−−−−−
hν−1
Ps − DNA (1.2b)
Ps − DNA hν2
−−−−−*
)−−−−−
hν−2
Ps = DNA (1.2c)
Ps hν3
−−−−−→Photobreakdown Products (1.2d)
The linear sequence of steps from Eq. 1.2 can be outlined as follows:
• Intercalation,Eq. 1.2a and Fig. 1.7, is a necessary step in the reaction because it will make sure that psoralen and the relevant bases of DNA are close. Closeness is required for the reaction to proceed along the steps given in Eq. 1.1c. The necessity of proximity can also be understood from the fact that the excited state of the psoralen, 1Ps?in Eq. 1.1c, has a limited lifetime. If no suitable reactant is present during this time, no reaction will occur. This is demonstrated by the finding that if psoralens are allowed to react with free nucleosides, in which case no intercalation can occur, the reaction becomes diffusion limited [4]. The precise geom- etry of intercalation strongly affects the successive photobinding to the nucleic acids [14]. In particular, intercalation will determine which one of the possible isomers is formed during the subsequent addition reactions.
The stereochemistry is discussed further in section 1.4.
Adduct Purine
Pyrimidine Purine
Pyrimidine
Furan-side Pyrone-side
Backbone Psoralen Backbone
Adduct
Figure 1.7: During intercalation, the psoralen positions itself between two base pairs of DNA. In subsequent steps, covalent bonds are primarily formed between the psoralen and pyrimidines, as is illustrated above.
• After intercalation, the first addition reaction, Eq. 1.2b, will pro- duce a monoadduct – a covalent link between one strand of DNA and the psoralen. As psoralens are heterocyclic compounds, there are two possible kinds of monoadducts: furan-side and pyrone-side. Experiments indicate that furan-side monoadduct is much more common than pyrone- side monoadduct. It is very likely that the furan-side monoadduct is also the precursor of the crosslink [6, 14]. The reason for this preference may be that the pyrone-side monoadduct will not absorb the photon energy necessary to drive crosslink formation [15]. The fact that the furan-side monoadduct may be the only monoadduct formed in nature does not make a theoretical study of the pyrone-side monoadduct any less inter- esting – through computations one may be able to find out why this is so.
• Whatever monoadduct is formed, it must be able to absorb a second photon. This will drive the crosslink or diadduct formation, Eq.
1.2c, and thus the entire reaction to completion. Two crosslinks create loops in the DNA which can be observed using electron microscopy [2].
• Free psoralens in water dimerise and photoreact with the solvent to pro- duce photobreakdown products, Eq. 1.2d. These products are unable to react with DNA any further [6]. In the context of this investigation they are also labelled “not interesting”.
In total, the full reaction, as modelled in Eq. 1.2, is a two-photon reaction.
As indicated, it is also photoreversible [4]. The cyclobutane rings formed by the addition reactions can be broken by excitation at much shorter wavelengths to reproduce the starting materials [6].
---
1.4
Stereochemistry of the Psoralen–DNA Adducts
Monoadduct and crosslink formation are addition reactions. In the case of psoralens and DNA, a [2 + 2] cycloaddition reaction as shown in Fig. 1.8 can occur between the 5, 6-double bond of a pyrimidine and either the 40, 50-double bond of the furan or the 2, 3-double bond of the pyrone ring [6, 8, 16]. In either case, cycloaddition results in the formation of a new cyclobutane ring.
C
C C
C hν
C C
C C
Figure 1.8: The [2 + 2] cycloaddition reaction. A four-membered cyclobutane ring is formed.
In general, addition reactions can occur in either a syn or anti orientation [8]. Since the particular reaction studied here occurs between doubly bonded carbons in closed rings1, all reactions will be syn with respect to the addition.
Once the reaction has reached completion by crosslink formation, the product may be further classified according to its stereochemistry.
1One could perhaps say that we are dealing with a cyclo-cyclo-cycloaddition – two rings are added to produce a total of three rings.
• The cis configuration refers to stereoisomers in which both psoralen and pyrimidine are bonded to the same side of the cyclobutane ring.
• The trans configurationrefers to rings in which psoralen and pyrimidine lie on opposite sides of the plane of the cyclobutane. Note that cis and trans are complementary to each other.
• We define syn for the furan-side adduct as the structure in which O10 of psoralen and N1 of the pyrimidine are on adjacent corners of the cyclobutane.
• We define syn for the pyrone-side adductas the structure in which C2 of psoralen and N1 of the pyrimidine are on adjacent corners of the cyclobutane. Note that syn and anti below are complementary to each other.
• The anti configurationrefers to bonding patterns in which the N1 of pyrimidine and either the 10-O or C2 of psoralen as appropriate, are on diagonal corners of the cyclobutane.
Note that syn and anti here do not refer to the stereochemistry of the reaction, but to the stereochemistry of the products. Generally the stere- ochemistry of all psoralen–thymine photoadducts has been confirmed to be cis–syn-structures [4].
Each of the four carbon atoms in the cyclobutane ring is a chiral centre, so there are 24 = 16 possible stereoisomers of each monoadduct. However, 4+ 4+ 4 = 12, of these will be trans-fused rings. These are unlikely to occur in practice due to steric constraints [6]. Theoretically, diadduct formation should increase the number of possible isomeric products. It is instead found that geometry limits the possibilities. For the thymine–thymine reaction only one isomeric pair was found to lead to sterically reasonable crosslinking. All other monoadducts will have the psoralen sticking out from the DNA helix, which makes any subsequent crosslinking improbable. Note that this may not be true for folded DNA [4].
---
1.5
Motivation 1.5.1 Theoretical Motivation
Historically, drugs were discovered by serendipity rather than by deliberate design [17]. A medical company may thus spend 75% of its research budget on projects that will never become finished products – a percentage they would be happy to cut down2. If one could use inexpensive computer simulations to screen candidate projects worth pursuing in more expensive laboratory work, much would already be gained [11]. Computer simulations may also be a way to gain understanding of processes that are difficult to observe through direct experimentation. This motivates the theoretical approach chosen for the project.
Psoralens in particular have found many scientific applications [12]. In many areas it is their ability to “freeze” DNA that makes it an interesting object of
2Attila B´erces in his talk on Computing Engines for Drug Discovery: The Next “Quantum Leap”? , given at the Fourth Workshop on Linux Clusters for Super Computing in Link¨oping, 23rd October .
study [4]. The “freezing” is a consequence of the increased rigidity introduced by crosslink formation which will be examined here. Because the reactions from Eq. 1.2 can be carried out under a wide variety of conditions, psoralens are useful for more general DNA investigations as well [4]. This in turn explains the volume of experimental evidence available.
In the field of human pathology there are still other applications of pso- ralens: autoimmune disorders, multiple sclerosis, organ transplant rejection, rheumatoid arthritis and AIDS [14]. Two others will be presented below.
1.5.2 Psoriasis
Figure 1.9: The silvery or micaceous scales char- acteristic of psoriasis [18].
Copyright , c° War- ren Piette and the Univer- sity of Iowa’s Virtual Hos- pital. Used with permis- sion.
Figure 1.10: Vitiligo typ- ically involves areas such as the ears, mouth, nose and around the eyes [18].
Copyright , c° War- ren Piette and the Univer- sity of Iowa’s Virtual Hos- pital. Used with permis- sion.
Psoriasis, Fig. 1.9, is a chronic, rarely cured, skin disease. It is characterised by hyperproliferative conditions [14]. The symptoms of the most common type of the disease are raised, thickened patches of red skin covered with silvery- white scales [19]. The molecular details of the disease are unknown; thus today, treatments in use are based on experience rather than on a thorough understanding of the disease. Therapies fall in three major categories: topical, where one applies the therapeutic agent directly to the affected areas of the skin; systemic, where the agent is instead taken internally; and finally a form of photochemotherapy [19]. The work presented here may be of relevance to the latter type of treatment.
1.5.3 Vitiligo
Vitiligo, Fig. 1.10, is a disorder which destroys pigmentation cells in the skin, the mucous membranes, and the retina. Unlike psoriasis, it is characterised by reduced proliferative conditions [14]. As a result, affected areas loose their colouring. Hair growing in vitiligo-stricken patches usually turns white [20]. As in the case of psoriasis, the cause of vitiligo is yet unknown. While treatment with psoralen and UVA, called PUVA for short, does not directly achieve the ultimate goal of therapy – permanent repigmentation of white patches – it is probably the most beneficial treatment for vitiligo available.
1.5.4 Proposed Mechanism of Treatment
Psoralens have been used in the treatment of skin diseases since approximately
BC [21]. The beneficial effect of PUVA in psoriasis is believed to be the decrease of skin cell proliferation. The psoralens when activated by ul- traviolet light will bind to DNA as was outlined in section 1.3 above. Such covalent modifications are thought to impair the action of DNA polymerases and thus inhibit the cell cycle. For vitiligo, the benefits may arise from pso- ralen cycloadducts, that together with lecithin, are able to mimic a secondary messenger which, in turn, stimulates melanogenesis. It is known that psoralens trigger the release of lysomal enzymes that degrade RNA [2]. However, the importance of lysosomes in skin disease treatment is unclear.
While the details of the cellular psoralen action are unknown, there is cer- tainly no magic involved: neither psoralen nor ultraviolet light is able to distin- guish healthy from non-healthy cells. Unsurprisingly, PUVA treatment comes with a broad set of side effects which in some cases may be severe and even counter-productive: itching, dry skin, loss of skin pigmentation [19, 20, 22]
and even skin cancer [14, 15]. PUVA therapy can also be very time-consuming.
Nevertheless, it has a documented efficacy in many patients and is one of the most promising treatments available to date [19].
Molecular Mechanics
Molecular mechanics is a generalisation of the Ball and Stick model as illustrated in Fig. 2.1, in the sense that atomic
positions are no longer fixed in space [11]. Molecular mechanics is in turn a model of the “real” quantum mechanical system, which is the subject of the next chapter.
Compared to quantum mechanics, molecular mechanics has a huge advantage: computational throughput.
Figure 2.1: Ball and stick model of psoralen. See also Fig. 1.6(a).
---
2.1
Notational Conventions
This chapter marks the beginning of the theoretical treatment of the methods used for studying the psoralen–DNA reaction. In what follows the conventions listed below are adapted:
• Vectors are written using lowercase boldface letters, such as r and x.
Vectors are distinguished from matrices, which are written using up- percase boldface letters, like A and Λ.
• Operatorsare differentiated from matrices in that they are written using uppercase italic letters, for instance P and H .
• The bar-notation is used for complex conjugates: the conjugate of a complex number α is written α. The bracket-notation is reserved for expectation values: the expectation value of a quantity β (x) is written hβi.
• The number of electronsin any given system will assumed to be n.
The number of nuclei will likewise assumed to be N .
2.1.1 The Bra–Ket Notation
Much of the mathematics of computational chemistry is really linear algebra in disguise [23]. In quantum mechanics, the standard notation for a vector in an M -dimensional complex space is:
| ϕi =
α1
α2
... αM
(2.1)
A column vector such as | ϕi above is called a ket. The vector dual to | ϕi is written hϕ | and is called a bra.
hϕ |=£
α1 α2 . . . αM
¤ (2.2)
From Eqs. 2.1 and 2.2 we see that hϕ1| ϕ2i is a compact notation for the standard inner product1. A similarly abbreviated notation is used for the tensor product:
| ϕ1ϕ2i =| ϕ1i | ϕ2i =| ϕ1i⊗ | ϕ2i (2.3) An operator can generally be thought of as a matrix – a bracket, hϕ1| ϕ2i, is often referred to as a matrix element [11]. Applying an operator A to a vector is equivalent to a matrix–vector multiplication:
hϕ1| A | ϕ2i = (hϕ1|) · (A | ϕ2i) =¡
hϕ1| A†¢
· (| ϕ2i) (2.4) In all our cases, vectors are actually one-dimensional scalars, meaning that M = 1 in Eqs. 2.1 and 2.1 above. However, the single component is a function, ϕ, of one or more spatial variables, x1, x2, . . . , xM0. The difference between a bra and its ket is thus only complex conjugation.
| ϕi = [ϕ (x1, x1, . . . , xM0)] = hϕ | (2.5) Because our bras and kets are complex functions, we depart from the stan- dard inner product and instead define the inner product in terms of integration over the domain of the spatial variables:
hϕ1| ϕ2i = Z Z
· · · Z
ϕ1(x1, x2, . . . , xM0)
ϕ2(x1, x2, . . . , xM0) dx1dx2· · · dxM0 (2.6)
---
2.2
Molecular Mechanics – Force Field Methods
To a first approximation, we may model molecules as structures composed of atoms of different sizes, held together by spring-like bonds of different strengths. This is a considerable approximation, because individual electrons are ignored. Neglecting electrons means that the method is restricted to the study of molecular ground states. Neither can we treat reactions, since the
1Which for | ϕ1i = [αi]Mi=1∈ CM and | ϕ2i = [βi]Mi=1∈ CM is defined as:
hϕ1| ϕ2i = XM i=1
αiβi
The inner product is thus written as a bracket which consists of two parts: a “bra” to the left and a “ket” to the right. The “c” in between comes at no extra cost.
bond breaking and bond making that a reaction involves is really all about electron transfer. Nevertheless, unlike quantum mechanics, the efficiency of molecular mechanics means that it has the potential to capture the dynamics of large molecules.
In order to make use of molecular mechanics, we need the energy of the molecule. This can be written as a sum of terms arising from specific molecular distortions [11]:
Etotal= Estretch+ Ebend+ Etorsion+ Eelectrostatic+ EvdW+ Ecross (2.7) The energy is actually a function of the set of all atoms that make up the molecule, {A, B, . . . , Z} and their respective locations in Cartesian space, {rA, rB, . . . , rZ}2. It can be evaluated using only classical mechanics, which makes the theory intuitively simple. The meaning of the individual terms in Eq. 2.7 is discussed below, where we chose to denote parametric constants by
k. A B
Figure 2.2: The stretch distortion.
Energy
Interatomic Distance
Figure 2.3: The bond- stretch energy as a func- tion of the interatomic dis- tance. The solid line rep- resents the harmonic po- tential from Eq. 2.8, the dashed line is the Morse potential from Eq. 2.9.
• Estretch = Estretch(A, B) is the energy function for stretching a bond between two atoms A and B, see Fig. 2.2. It is a function of the distance between the two atoms, rAB= |rA− rB|. In its simplest form it may be written as a truncated Taylor expansion around the experimental natural bond distance, r0AB:
Estretch(A, B) = kAB(rAB− r0AB)2 (2.8) Note that the first-order term in the Taylor expansion is assumed to be identically zero because the bond stretch energy should have a minimum for rAB= r0AB. For bond distances far away from the natural distance, the energy as given in Eq. 2.8 yields results that do not agree well with the fact that molecules can dissociate. As the atoms separate, one may instead apply a Morse potential:
Estretch(A, B) = k1AB
³
1 − e−k2AB(rAB−r0AB)´2
(2.9) The two potential functions are illustrated in Fig. 2.3. Even though the Morse potential yields more accurate results at long interatomic dis- tances, one may chose the harmonic potential for computational reasons;
especially if one knows that bond distances are in the vicinity of the nat- ural distance [17].
• Ebend = Ebend(A, B, C) is the energy required for distorting a bond angle, θABC, formed by three atoms A, B and C, see Fig. 2.4. It is usu- ally computed as a truncated Taylor expansion around the experimental natural bond angle, θ0AB:
Ebend(A, B, C) = kABC(θABC− θ0ABC)2 (2.10)
A
B C
Figure 2.4: The bend dis- tortion.
If one wants higher accuracy one usually includes more terms in the expansion. For atoms in special environments it may make sense to add an extra energy term for the out-of-plane bending angle [11].
2If we for a brief instant assume that the alphabet has not 26 but N unique letters from A to Z.
The notion of a “general natural bond angle” breaks down in the case of ring systems. The natural bond angle in a cyclic carbon can be very different from a the natural bond angle in a linear carbon [11]. For this reason, many force fields have separate atom types for elements when they are members of rings. Atom types are formally introduced in section 2.3.
• Etorsion= Etorsion(A, B, C, D) is the energy associated with rotation around the B-C-bond in a four atom sequence A-B-C-D, see Fig. 2.5.
The torsional energy should at least be periodic with a period of 2π; any other periods must be divisors of 2π. Because the energy introduced by torsional distortions is often low, one usually does not apply Taylor expansions, but instead chooses a Fourier series [17]:
Etorsion(A, B, C, D) =
k1
X
n=1
k2ncos (nωABCD) (2.11)
The torsion angle, ωABCD, is defined as the A-B-D angle when looking down the B-C bond. The n = 1 term in Eq. 2.11 describes rotations that are periodic 2π, n = 2 has period 2π/2 = π and so on.
A
B C
D
Figure 2.5: The torsional distortion.
• Eelectrostatic= Eelectrostatic(A, B) describes the energy that is due to the partial charge of the atoms in the molecule. The change arises from an uneven sharing of electrons. The interaction between two point charges, QA and QB, separated by a distance rAB, is given by the Coulomb potential:
Eelectrostatic(A, B) = QAQB
4π²0rAB (2.12)
Usually one does not attempt to compute partial charges in molecules, al- though this would be possible to do [24]. Instead, one assigns parameters on the basis of electronic structure calculations [11].
• EvdW= EvdW(A, B) is the van der Waals energy describing additional interaction between non-bonded atoms. Unlike the electrostatic energy in Eq. 2.12 above, the van der Waals energy is not concerned with the partial charge on the atoms. In quantum mechanical terms, it is due to the overlap of the electron clouds of two atoms A and B [11]: because of electron correlation, two electrons may attract each other at intermediate distances, while they clearly repel each other at short distances. The van der Waals energy is usually computed as a Lennard–Jones potential about the natural bond distance [11, 17]:
EvdW(A, B) = ELJ(A, B) =
= k1
(rAB− r0AB)12 − k2
(rAB− r0AB)6 (2.13) The attraction is due to the motion of the electrons which may create instantaneous polarity, even in overall non-polar molecules. It can be shown that the attraction varies as the sixth power of the distance be- tween the two atoms. The power in the repulsive part is chosen to be
12 for efficiency reasons only3. Together with the electrostatic contribu- tion, the van der Waals energy is the most time-consuming part in the evaluation of the total energy.
• Ecross= Ecross(rAB, θABC, ωABCD) attempts to cover effects arising from interaction of the three first terms. For instance, a change in bond angle may affect the natural bond length. The particular implementation is dependent on the force field used, but it is usually in the form of Taylor expansions.
A stable conformation of the molecule will be obtained by optimising Etotal
in Eq. 2.7 with respect to atomic positions. Optimisation is discussed further in section 2.5. The energy as computed using Eq. 2.7 is inherently relative – not absolute. It has no meaning by itself but only makes sense when compared to energies of molecular structures having the same chemical composition and structural units [11].
---
2.3
Atom Types
The foundation of molecular mechanics is that molecules tend to be composed of units which are structurally similar. The characteristics of, for instance, a C−H bond are approximately, but not entirely, independent of the surrounding.
As was noted in section 2.2, the bond angle of a carbon atom is different depending on whether it is cyclic or linear. The more care is taken of the particular surrounding of an atom, the better conformance to experiment one may achieve.
This motivates the conception of atom types: on a high level, different atoms are distinguished by their element types alone; on a lower level they are distinguished based on their surrounding as well. In the Amber force field for example, a nitrogen in a five-membered ring is distinguished from a nitrogen in a six-membered ring. Note, however, that there is a risk of overfitting in this process: if too many atom types are used, the generality as well as the predictive ability of the method decreases.
---
2.4
Force Field Parametrisation
For every contribution to the energy there a several parameters. These may be different among all the different atom types. The choices of atom types and parameter values are what compose the force field. The parametrisation process is built on subjective decisions. The result is that some force fields are better suited for certain tasks than others [25].
The number of parameters to determine is huge and it is not uncommon to only support a subset of all elements or possible parameters4 in a force field.
Ultimately, one would rely on experimentally, spectroscopically derived values for all parameters. In practice it is however not uncommon to make use of electronic structure computations for certain values [11].
3Once one has computed α = ą
rAB− rAB0 ć6
computation of β =ą
rAB− rAB0 ć12 is trivially β = α2= α · α.
4For instance, only ≈ 0.2% of all torsional parameters have been parametrised [11].
2.4.1 Amber
One solution to the parametrisation problem is restricting the set of molecules handled by the method. The Amber force field is parametrised for proteins, nucleic acids and other molecules with related functional groups which are of interest in organic and biological chemistry [25]. In the original paper, Amber was validated against several biomolecules such as alcohols, dimethyl phosphates, ethers, the nucleic bases, sugar peptides and sulfur compounds.
Many parameters in Amber are obtained from experimental X-ray and NMR data. Charge parameters are determined computationally. Bond stretches and angles are computed using harmonic Taylor expansions. A small set of Fourier coefficients is utilised in the evaluation of the dihedral torsion energy. Equations 2.12 and 2.13 are used in the computation of electrostatic and van der Waals energies respectively. In order to gain efficiency, non-bonded interactions are only computed if the concerned atoms are more than three bonds apart [25].
2.4.2 UFF – Universal Force Field
An alternative approach to parametrisation is taken in the Universal Force Field, UFF. The idea is to derive di-, tri- and tetra-atomic parameters from atomic constants. This leads to a reduced parameter set [11]. It also sacrifices accuracy, but the benefit is that the force field becomes applicable to a larger range of molecules. In the original paper, UFF was validated against three handful of molecules of various kinds [26].
Unlike Amber, UFF supports molecules built from all elements of the pe- riodic table. Internally, it also uses hybridisation state in computations. UFF offers a choice of Taylor expansion or Morse potential for evaluating bond stretch energy. All angular distortions are computed using Fourier expansions, including a separate term for the out-of-plane energy. Equations 2.12 and 2.13 are used for computation of the electrostatic and van der Waals energy respec- tively. Non-bonded interactions are not computed for atoms that are bonded or share a common bonded atom.
---
2.5
Geometry Optimisation
Geometry optimisation is central to the early stages of the computations per- formed here: all models are geometry-optimised prior to any further calcula- tions. Commonly geometry optimisation is done by varying the 3N Cartesian coordinates of the nuclei subject to the constraints of the force field5. This procedure requires first-, and sometimes also second-order derivatives to be computed. Fortunately the energy terms from Eq. 2.7 are easily differentiated.
The theory of calculus gives a necessary condition for an energy minimum.
In the one-dimensional case, any optimum point x? must satisfy:
E0(x?) = 0 (2.14)
But our starting point is not the energy minimum, which is yet unknown, but a point x, displaced from the optimum, x?, by a small, residual distance, δx. Mathematically, we have x? = x + δx. If the derivative in Eq. 2.14 is
5But note that geometry optimisation is not a procedure unique to force field methods.
Geometry optimisation for quantum mechanics methods is briefly discussed in chapters 5 and 6.
written as a Taylor expansion, the necessary condition for an energy minimum turns into:
E0(x + δx) = E0(x) + E00(x) δx + E000(x) δx2+ . . . = 0 (2.15) Truncating the series at second order and rearranging terms yields an ex- pression for the residual, the displacement from the true optimum.
δx = −E0(x)
E00(x) (2.16)
Truncating at second order implicitly assumes that the minimum is exactly quadratic in behaviour. Close to the minimum, this assumption is probably reasonable. Inserting Eq. 2.16 into the expression relating the energy minimum to the residual, we have:
x?= x − E0(x)
E00(x) (2.17)
This is the Newton–Raphson method for finding the optima of the function E (x) [27]. There is nothing in Eq. 2.17 that prevents the accidental finding of maxima instead of a minima. In order to find a minimum, we will have to require that:
E00(x?) > 0 (2.18)
If Eq. 2.18 is satisfied we know that we have found a minimum. In most cases this is not exactly what we want: we want to find the global minimum.
In general we would have to first find all the minima in order to know which of them is smallest. In most practical cases this is not feasible.
2.5.1 GDIIS
GDIIS is a version of the “direct inversion in the iterative subspace”, DIIS, algorithm, adapted for the special case of geometry optimisation [28]. We first note that in the more general, multidimensional case, Eq. 2.16 above can be rewritten in matrix notation:
δxi= −H−1gi (2.19)
Here H = [hij] = £
∂2E/∂xi∂xj
¤ is the Hessian matrix. The Hessian matrix and its inverse, H−1, can be derived or approximated using standard techniques known as Hessian Update Techniques [17]. The vector g = [gi] = [∂E/∂xi] is the gradient vector. Since Eqs. 2.16 and 2.19 are in principle equivalent, GDIIS also assumes that the sought minimum is quadratic. If the minimum were truly quadratic and the exact Hessian matrix were known, convergence would be obtained in one step using Eq. 2.19 [29]. Note also the introduction of indices, i, in Eq. 2.19 which mark the iteration step during which the vectors were computed.
The essence of GDIIS is that the geometry vectors, xi, and the gradient vectors, gi, at iteration i = m, are linear combinations of the corresponding vectors from previous iterations. When taking an iteration step with GDIIS, all previous steps are taken into account. Thus, at iteration m we have:
xm+1=Pm
i=1cixi
gm+1=Pm
i=1cigi
where Xm i=1
ci= 1 (2.20)
The coefficients ciare obtained by minimisation of the norm of the residual vector δxi. Doing this using Lagrangian multipliers leads to a matrix which can be inverted directly in order to solve for the coefficients ci. Because this is done in the subspace of the iterations, we are in fact “directly inverting in the iterative subspace”, which accounts for the origin of the algorithm’s name.
Iteration continues until chosen convergence criteria are met.
Quantum Mechanics
The inclusion of electrons in the molecular model is necessary if we wish to study reactions. Electrons are very light particles
which classical mechanics cannot correctly describe. We therefore make the transition to quantum mechanics.
Ehrenfest’s theorema clarifies the relationship between classical and quantum mechanics: classical mechanics deals with average values, while quantum mechanics deals with the
underlying details [10].
aIn mathematical terms, Ehrenfest’s theorem states that:
8<
:
dhpi dt = hf i
dhri dt = hpim
Where f is the force acting on a particle with mass m and momentum p, located at r.
---
3.1
Some Basic Postulates of Quantum Mechanics
The theory of quantum mechanics – at least as it is used here – is remarkable in the sense that it can be derived from a few basic assumptions:
• The state of a system of M particles at a time, t, is fully described by a mathematical function, called the wave function:
Ψ = Ψ (r1, r2, . . . , rM, t) (3.1) The spatial coordinates of the particles 1, 2, . . . , M are denoted by r1, r2, . . . , rM. Wave functions that are independent of time will be written:
ψ = ψ (r1, r2, . . . , rM) (3.2)
• The probability that a particle is found in a volume element1dτ around r is proportional to |ψ (r)|2dτ . This is the Born interpretation of the wave function.
• The wave function evolves in time according to the time-dependent Schr¨odinger equation:
H Ψ = i~∂Ψ
∂t (3.3)
In cases where the potential is independent of time, separation of vari- ables, Ψ (r1, r2, . . . , rM, t) = ψ (r1, r2, . . . , rM) θ (t), leads to the time-independent Schr¨odinger equation:
H ψ = Eψ (3.4)
• The total electron wave function, including spin, must be antisymmetric with respect to interchange of any pair of electrons. This is the Pauli principle2.
3.1.1 The Hamilton Operator
Every physical observable has an operator , O, associated with it. The value of of the observable, or its eigenvalue, ω, can in principle be obtained by solving the eigenvalue equation:
OΨ = ωΨ (3.5)
Most of the time we are interested in expectation values of the observables.
These may be extracted by integration of the wave function:
hωi =
R ΨOΨdτ
RΨΨdτ =hΨ | O | Ψi
hΨ | Ψi (3.6)
Usually the bracket notation on the expectation value of observables is omitted3. As far as we are concerned, the most important operator of quantum mechanics is the Hamilton operator , or the Hamiltonian, H :
H = T + V = −~2
2m∇2+ V (3.7)
The Hamiltonian gives the energy of the system it is applied to, see Eq.
3.4. The total energy has a kinetic, and potential contribution, T and V respectively. The kinetic energy is in turn related to the linear momentum, p, by T = p2/2m, where m is the mass of the particle. The Hamiltonian is uniquely determined by the wave function from the relationship in Eq. 3.4.
Of special interest for the upcoming derivations is the core Hamiltonian.
The core Hamiltonian is a hydrogen-like Hamiltonian for one electron i in the field of N bare nuclei located at rj with charge Zje, 1 ≤ j ≤ N .
hi= − ~2 2me∇2i −
XN j=1
Zje2
4π²0|ri− rj| (3.8)
1Conventionally the volume element is written dν = dx × dy × dz. If spin is also included in the integration we instead write dτ .
2This follows because electrons are fermions: their spin is either +1/2 or −1/2.
3But do note the similarity in notation between the expectation value of an observable and the bra–ket notation, introduced in section 2.1.1, used for the integrals to compute it.
In particular one should note that the Hamiltonian is a Hermitian operator, although this is not shown here. Furthermore, in Eqs. 3.7 and 3.8 we have implicitly ignored relativistic effects. In general, fast moving electrons close to the nuclei have speeds that are not negligible in comparison to the speed of light [11, 17].
3.1.2 The Born–Oppenheimer Approximation
It is impossible to solve the Schr¨odinger equation analytically if the system contains more than two interacting particles. The Born–Oppenheimer approxi- mation takes note of the considerable difference in mass between electrons and nuclei. Instead of attempting to solve Eq. 3.4 for all particles simultaneously, we regard the nuclei as fixed in space, and solve for the electrons in the static electric potential arising from the nuclei. This is fundamentally different from the implicit approach taken to the Born–Oppenheimer approximation in force field methods, where only nuclear motion is considered, while electrons, due to their rapid motion, give rise to a fixed electron distribution around the atoms [17].
H = − ~2
2me∇2+ V (3.9)
Note how the Hamiltonian in the Born–Oppenheimer approximation does not take the motion of the nuclei into account. The static nuclear–nuclear repulsion, however, can be absorbed into V . Motion of the nuclei can be treated at a higher level of abstraction by allowing them to move on a potential energy surface which is a solution to the electronic Schr¨odinger equation [11].
Even with the Born–Oppenheimer approximation, the only real system that can be solved analytically is the hydrogen molecule-ion, H+2. Because the interaction of the nuclei has been artificially removed, the system still only contains two interacting particles.
3.1.3 The Slater Determinant
Wave functions that satisfy the one-electron Schr¨odinger equation are called orbitals. An orbital is thus merely a synonym for a one-electron wave function [17]. This is not an approximation per se until we consider polyelectronic systems and treat the electrons separately, each with its own orbital.
In order to satisfy the Pauli principle from section 3.1, we will need a way to represent spin. We therefore introduce the spin orbital, φ, a joint spin–space state of the electron. It is a spatial orbital multiplied by a spin function. The spin function for electron i is either α (i) or β (i), corresponding to the two spin states of the electron, +1/2 and −1/2.
φα(i) = ψ (ri) α (i)
φβ(i) = ψ (ri) β (i) (3.10)
By the Pauli principle, every spin orbital can only accommodate one elec- tron. Thus, in any given system the number of electrons cannot be higher than the number of spin orbitals. For n electrons and equally many spin orbitals, φ1, φ2, . . . , φn, we may express the wave function of the ground state with a Slater determinant.
ψ (1, 2, . . . , n) = r1
n!
¯¯
¯¯
¯¯
¯¯
¯
φ1(1) φ2(1) · · · φn(1) φ1(2) φ2(2) · · · φn(2)
... ... . .. ... φ1(n) φ2(n) · · · φn(n)
¯¯
¯¯
¯¯
¯¯
¯
(3.11)
The columns of the Slater determinant are single electron wave functions:
orbitals, or molecular orbitals in case of polyatomic systems. The rows represent the electron coordinates. Writing the wave function as a Slater determinant will ensure it obeys the Pauli principle. The properties of the determinant automatically lead to the antisymmetry of the wave function: interchanging for instance electrons 1 and 2 corresponds to an interchange of the first two rows in the determinant, leading to a change of sign in the wave function. If two electrons were to occupy the same spin orbital4, the Slater determinant would have two identical columns which would make the total wave function zero. According to the Born interpretation, such a state has zero probability of existence.
We may assume that a set of spin orbitals is orthonormal, hφi| φji = δij, where δijdenotes the Kronecker delta. This is no restriction; if the spin orbitals are not orthonormal, they can be orthonormalised by a Schmidt transformation, which only changes the wave function by a multiplicative factor, lacking physical significance.
---
3.2
The Variational Principle
Variational theory is a way of assessing and improving functional guesses, f , starting with a trial function, ftrial. The variational principle is useful in this context as it can tell whether one trial wave function, ψtrial1, is better or worse than another, ψtrial2.
Suppose that the system is described by a Hamiltonian, H . Denote the unknown, true eigenfunctions of the Hamiltonian by ψi, ordered such that the eigenvalues are increasing: Ei ≤ Ej if i < j. Because the Hamiltonian is Her- mitian, the true eigenfunctions form a complete basis. Any trial wave function can thus be written as a linear combination of the these true eigenfunctions.
ψtrial= X∞ i=0
ciψi where H ψi= Eiψi (3.12) Note that Eq. 3.12 allows for the construction of infinitely many trial wave functions. We now introduce the Raleigh ratio defined as:
² =hψtrial| H | ψtriali
hψtrial| ψtriali = hψtrial| H | ψtriali (3.13) As can be seen from Eq. 3.6, ² is the expectation value of the energy of the state described by ψtrial. The last equality in Eq. 3.13 follows if we require ψtrial to be normalised. As was noted in section 3.1.3 above, we may assume the set of the ψi to be orthonormal, so that hψi| ψji = δij. Therefore:
4Meaning that they occupy the same spatial orbital and have the same spin.
hψtrial| (H − E0) | ψtriali = X∞ i, j=0
cicjhψi| (H − E0) | ψji =
= X∞ i, j=0
cicj(Ej− E0) hψi| ψji =
= X∞ i=0
cici(Ei− E0) =
= X∞ i=0
|ci|2(Ei− E0) (3.14)
By definition E0 is the smallest eigenvalue of H , so Ei ≥ E0 for all i.
Trivially |ci|2≥ 0, so the last expression in Eq. 3.14 must be greater or equal to zero. Inserting 3.14 into 3.13 we have the result:
² ≥ E0 for any ψtrial (3.15)
In order for equality to hold in Eq. 3.15, we must have ci = 0 for i 6= 0, that is, ψtrial = ψ0, the true eigenfunction with the lowest energy. We have shown that the trial function that gives the lowest energy is the optimum function of that form, because the Raleigh ratio is not less than the true ground state energy of the system: we have a lower bound on the true energy of the system.
---
3.3
Hartree–Fock 3.3.1 The Energy of a Slater Determinant
The goal of the procedure derived here is to find the set of orthonormal spin orbitals, {φ1, φ2, . . . , φn}, that yield a minimum energy, E. These spin orbitals are related to the normalised Hartree–Fock ground state wave function, ψ, by the n-electron Slater determinant.
ψ = r1
n!det |φ1(1) φ2(2) · · · φn(n)| (3.16) Equation 3.16 is just a simpler way of writing the Slater determinant as it was defined in Eq. 3.11 previously. We know from the variational principle in section 3.2 that the lowest energy will give us the wave function that lies closest to “the truth”. A suitable first step will therefore be to find an expression for the energy, E, in terms of spin orbitals5.
E = hψ | H | ψi (3.17)
In the Hartree–Fock approximation, or the mean field approximation, each electron moves independently in the average potential of the other n − 1 elec- trons and the N nuclei. This is an effect of using a single Slater determinant to describe the entire system: the electron distribution is modelled in a contin- uous, average way – electrons are not discretely localised. Electron correlation
5Other derivations are possible; in “Till¨ampad Molekylorbitalteori”, Lindner and Lunell give a version that is closer to the original derivation by Fock and Slater from .
