Application of simplified models of electron density for estimation of the electrostatic properties of molecules

Application of simplified models of electron density for estimation of the

electrostatic properties of molecules

PhD Thesis

Author: Sławomir A. Bojarowski Supervisor: Paulina M. Dominiak

Faculty of Chemistry

University of Warsaw

To my wife and mom.

Acknowledgements

Special thanks to my supervisor Paulina M. Dominiak for a great support through all my studies.

Also, I would like to thank my colleague and friend Pra- shant Kumar for a perfect collaboration.

Many greetings and thanks to all members of our lab super-group.

And last but not least, to my family members and spe- cially to my wife Beata, for all patience, support and indul- gence.

Also, the National Science Centre for PRELUDIUM-UMO

/2015/19/N/ST4/00010 grant and Foundation for Polish Science

for POMOST/2010-2/3/styp4 stipend are gratefully acknow-

ledged. Financial support from the Ministry of Science and

Higher Education of Poland is also acknowledged

5

Summary (in Polish)

Przybliżenie Borna-Oppenheimera, a także jego dokładniejsza wersja - przybliżenie adiaba- tyczne, jest podstawą dzisiejszej chemii. Wynika ono z obserwacji, że elektrony są dużo lżejsze, a zatem szybsze od jąder atomowych. Co za tym idzie, jądra poruszają się w uśrednionym poten- cjale od elektronowym (PES). Bez takiego podejścia ciężko byłoby w prosty sposób zdefiniować wiele podstawowych własności chemicznych i fizycznych. Jedną z takich własności jest energia oddziaływania międzycząsteczkowego - E

int

. Z definicji E

int

= E

_AB

− (E

A

+ E

_B

), czyli energia odziaływnia układu AB odjąć energię monomerów A i B, przy tych samych geometriach co w AB.

Są dwa podstawowe podejścia do liczenia energii odziaływań. Metoda supramolekularna bez- pośrednio odnosi się do definicji. Druga korzysta z rachunku zaburzeń Reileigha-Schroedingera, w którym, jako niewielkie zaburzenie, zakłada się oddziaływanie jednej molekuły na drugą. Drugie podejście daje głębszy wgląd w charakter oddziaływań. Kolejne poprawki to: energia oddziały- wań elektrostatycznych E

_es

, energia indukcyjna E

_ind

, energia dyspersyjna E

_dysp

, a także energia oddziaływań wymiennych E

_exch−rep

obecna nie tylko w pierwszym, ale także w wyższych rzędach poprawek.

Dla polarnych molekuł, np. dla makromolekuł o znaczeniu biologicznym: białek, nici DNA/

RNA, na odległościach równowagowych, największe znaczenie będzie miała energia oddziaływań elektrostatycznych. W takich przypadkach, E

es

zazwyczaj wykazuje największy wkład oraz największą korelację z E

_int

.

E

_es

, często nazywana też oddziaływaniem kulombowskim, swoją nazwę bierze stąd, iż fak- tycznie dotyczy oddziaływania zamrożonych układów, analogicznie do klasycznego oddziaływa- nia dwóch ładunków. W teorii obliczenie tego typu energii nie powinno nastręczać problemów.

Znając funkcję falową obu układów lub ciągły rozkład gęstości ładunku, należy przeprowadzić procedurę całkowania po całej przestrzeni dla układów A i B. W rzeczywistości zazwyczaj tego typu całkowanie w sześciowymiarowej przestrzeni lepiej jest uprościć. Oddziaływanie ciągłego rozkładu ładunku, dla się podzielić na szereg wkładów oddziaływań: jądra A - jądra B, elektrony A - elektrony B, elektrony A - jadra B i vice versa. Niemniej jednak, nadal stanowi to złożony problem i dalsze uproszczenia są pożądane.

Pewnym uproszczeniem jest rozwinięcie operatora oddziaływań w szereg potęg odległości między układami. Rozwinięcie to, zwane multipolowym, zakłada dobrą separację układów, co jest możliwe jedynie na większych odległościach. Jest to na tyle dobre uproszczenie, że energię sumuje się z oddziaływań pomiędzy współczynnikami przy kolejnych potęgach odległości (E

mtp

).

Współczynniki te nazywa się multipolami, jako że są one odchyleniami od symetrii sferycznej.

Kolejne multipole przyjmują nazwy zgodnie z 2

^k

, przy czym k=0,1,2,3, itd: monopole, dipole, kwadrupole, oktupole, itd.

Gdy, przybliżenie/rozwinięcie multipolowe zastosujemy niezgodnie z założeniem, czyli przy niewielkiej odległości układów, E

_mtp

nie będzie równało się E

_es

. Błąd przy tym powstały nazywa się energią penetracji, gdyż dotyczy penetrujących się rozkładów gęstości ładunku:

E

_pen

= E

_es

− E

mtp

. E

_pen

jest zawsze ujemne i szybko maleje, wraz ze zmniejszającymi się

odległościami miedzy układami.

dów o znaczeniu biologicznym. Ich prostota, parametry empiryczne zastosowane do funkcji zna- nych z fizyki klasycznej, pozwoliła zredukować koszt obliczeń do minimum, w stosunku do metod kwantowych. Prawo Kulomba zastosowane do ładunków punktowych prowadzi jednak do zbyt uproszczonego opisu właściwości elektrostatycznych molekuł. Można powiedzieć, że w tym przy- padku pola siłowe korzystają z rozwinięcia multipolowego zredukowanego tylko do pierwszego członu, do odziaływań monopoli. A zatem, w klasycznych polach siłowych nie ma uwzględnienia anizotropowości ładunku czy efektów penetracyjnych, rzutujących na dokładność estymacji E

_es

. Ostatnio gwałtownie wzrasta liczba modeli, które starają się szacować powyższe właściwości.

Część z nich opracowywana jest tak, by przy jednoczesnej dużej dokładności względem obliczeń ab initio, nie generowały dużych kosztów obliczeniowych. Modele te często polegają na adaptacji już istniejącego modelu, szeroko wykorzystywanego w innych dziedzinach chemii. Przykładowo popularnym źródłem punktowych momentów multipolowych są: distributed multipole analy- sis (DMA) Stone’a, cummulative atomic multipole moments (CAMM) Sokalskiego, stockholder partitioning (HD/HDI/BS-ISA) w pierwotnej wersji zaproponowane przez Hirshfelda oraz mo- del multipolowy Hansen’a-Coppens’a (HCM) (publikacje są wymienione w części References). Z drugiej strony modele pól siłowych nowej generacji nawiązują pośrednio, lub bezpośrednio, do technik chemii kwantowej, a przede wszystkim do DFT-SAPT.

O ile takie zjawiska dotyczące anizotropowego rozkładu ładunku są uwzględniane wraz z dodawaniem kolejnych oddziaływań multipoli. O tyle problem efektów penetracyjnych nie ma jednolitego rozwiązania. Część modeli stara się nawiązywać bezpośrednio do mechaniki kwanto- wej szacując E

_pen

przy pomocy tzw. całek nakrywania, bądź szacując bezpośrednio z definicji albo po prostu całkując oddziaływania z ciągłego rozkładu ładunku. Tego typu podejście cha- rakteryzuje się dokładnymi wynikami porównując do metod ab initio jednakże, często wiąże się z dużym czasem obliczeniowym, co utrudnia adaptację modelu bezpośrednio do modelowania molekularnego.

Z drugiej strony, częstymi składnikami modeli są różnego rodzaju funkcje wygaszające (ang.

damping functions). W większości przypadków można je przedstawić jako funkcje typu: (1 − exp( −αr)), gdzie r to odległość między układami, a α to współczynnik wyliczany przy użyciu różnego rodzaju obliczeń, bądź w najprostszym przypadku dopasowywany do jakiejś wartości z mechaniki kwantowej, bądź z innego źródła. To podejście nie powoduje znacznego zwiększenia kosztów obliczeniowych, natomiast często wprowadza elementy empiryczne, które z kolei mogą być specyficzne tylko dla jednego typu układów.

Należy przy tym dodać, że znaczna większość modeli, z obu podejść, jest tzw. testowana na danym zestawie molekuł. To znaczy, parametryzacja i weryfikacja przebiegają na tych sa- mych układach, co z definicji daje najmniejszy możliwy błąd, przez co fałszuje prawdziwy obraz zastosowań modelu.

Model, który byłby dokładny, przenaszalny, a do tego prosty jest wciąż w obszarze poszukiwań w obrębie modelowania molekularnego.

University at Buffalo DataBank jest bankiem parametrów modelu HCM, korzystającym z rozwinięcia multipolowego aż do heksadekapoli i powstał z myślą o udokładnieniach danych po- chodzących z rentgenografii. Jednak zakres jego zastosowań okazał się dużo szerszy. Wraz progra- mem LSDB, UBDB jest w stanie odtworzyć gęstość elektronową zadanej molekuły. Taka gęstość elektronowa może być wykorzystana do szacowania np. energii oddziaływań elektrostatycznych, potencjałów elektrostatycznych czy całkowitych molekularnych momentów dipolowych.

7

E

_es

przy użyciu UBDB można szacować korzystając z metody EPMM. Metoda ta polega na połączeniu bezpośredniej całki kulombowskiej przy małych odległościach, oraz wyliczania oddzia- ływań punktowych multipoli typu Buckinghama na dalszych odległościach. W celu wyznaczania błędów i limitów, a także korzyści płynących z metody UBDB+EPMM, przeprowadziłem szereg badań. Badania te w większość przeprowadziłem na małych dimerach cząsteczek organicznych, skonstruowanych tak, by w zbalansowanych sposób mogły odtwarzać oddziaływania obecne w kompleksach biomakromolekularnych. Użyłem zestawów, w których dimery te obecne były przy geometriach równowagowych, a także przy ośmiu różnych odległościach wzdłuż krzywej dysocja- cyjnej oraz przy ośmiu różnych kątach wokól geometrii równowagowej.

Metoda UBDB+EPMM daje błąd średniokwadratowy (RMSE) na poziomie 1,0 kcal mol

⁻¹

, porównując do SPDFG/aug-cc-pVTZ. Niski błąd jest możliwy, ponieważ UBDB+EPMM szacuje zarówno energię penetracyjną, jak i efekty związane z anizotropią ładunku. Dla cząsteczek polarnych E

_pen

może stanowić 50% E

_es

, dla niepolarnych nawet 100%. Energia penetracyjna do niedawna była uznana za zjawisko występujące na odległościach znacznie mniejszych niż geometria równowagowa. Jednakże moje badania pokazują, że w wielu przypadkach E

_pen

znika całkowicie dopiero przy odległościach dwukrotnie większych niż geometria równowagowa.

Wracając do anizotropii ładunku, przy użyciu UBDB, do dobrego opisu E

_mtp

, w większo- ści przypadków nie potrzeba heksadekapoli, jednak rozwinięcie przynajmniej do kwadrupoli jest niezbędne. Błąd związany z uśrednianiem parametrów HCM w UBDB jest nieznaczny i mniej- szy niż błąd metody. Największą przeszkodą w bezpośrednim zaadaptowaniu UBDB+EPMM do modelowania molekularnego jest potrzeba korzystania z procedury całkowania przy bliskich odległościach, co znacznie wydłuża czas obliczeniowy.

W celu obejścia tego ostatniego, ale bardzo istotnego problemu, zaproponowałem model aug-PROmol. Aug-PROmol polega na uproszczeniu modelu HCM do tak zwanej promolekuły, zmodyfikowanej o ładunki punktowe z dopasowywania do potencjału elektrostatycznego. Model promolekuły polega na modelu gęstości elektronowej wyrażonym poprzez sumę sferycznych gę- stości od pojedynczych atomów wchodzących w skład tej molekuły. Tego typu podejście gwaran- tuje płynne przejście przy ewentualnym zbliżaniu się układów. Szacowanie energii oddziaływań elektrostatycznych przy dalekich odległościach, to po prostu szacowanie E

_mtp

używając samych ładunków punktowych. Przy bliższych odległościach efekty penetracyjne uwzględniane są przy całkowaniu oddziaływań sferycznych gęstości elektronowych układów.

Wyniki E

_es

dla aug-PROmol, korzystając z EPMM, sa bardzo zbliżone do tych z UBDB+

EPMM. RMSE, na podobnym poziomie ≈1,0 kcal mol

⁻¹

, dla odległości równowagowych, a także przy różnych kątach przy odległości równowagowej. Przy tym, tak samo jak przy UBDB+EPMM, zauważalna jest niezależność kątowa. Należy przy tym nadmienić, że oddziaływanie ładunków punktowych (tych samych użytych do budowy aug-PROmol), bez uwzglednienia penetracji, daje wyniki obarczone bardzo dużym błędem przy odległościach równych i bliższych niż równowagowa.

Dodatkowo widać także zależność kątową, szczególnie w przypadku cząsteczek polarnych.

Dodatkowo, uproszczenie modelu HCM do gęstości sferycznych dało możliwość pominięcia

procedury całkowania. Zaproponowałem analityczny odpowiednik aug-PROmol, dopasowany

do energii oddziaływań elektrostatycznych z aug-PROmol. Badania zostały przeprowadzone

na homodimerze benzenu w ustawieniu tzw. pi-pi stacking. Błąd tej funkcji przy odległości

równowagowej był niewielki w stosunku zarówno do oryginalnego aug-PROmol+EPMM, jak i

mol

⁻¹

, 0.07 kcal mol

⁻¹

). Oczywiście należy przy tym zauważyć, że jest to tylko zarysowanie możliwości aug-PROmol, gdyż pełna analiza nie może być przeprowadzona w oparciu o tylko jeden przypadek.

Dalsze badania nad możliwościami aug-PROmol obejmowały próbę rozszerzenia aug-PROmol o ładunki punktowe z innych, mniej kosztownych obliczeniowo metod. Dopasowywanie do poten- cjału ab initio daje dobrej jakości ładunki punktowe, ale było wąskim gardłem metody. Przetesto- wałem aug-PROmol z ładunkami punktowymi z popularnej metody półempirycznej AM1-BCC, oraz z tabelaryzowanych ładunków punktowych związanych z bazą General Invariom Database (GID) tzw. zbankowanych RESP’ów. RMSE dla metody półempirycznej jest na tym samym po- ziomie co w przypadku wyjściowego aug-PROmol (< ≈ 1,0 kcal mol

⁻¹

) i nieznacznie większy dla zbankowanych ładunków punktowych (> ≈ 1,0 kcal mol

⁻¹

). Co więcej, uzyskane wyniki wska- zują, że dokładność szacowania E

es

zależy tylko od jakości użytych ładunków punktowych. E

pen

szacowana z aug-PROmol praktycznie nie zależy od źródła badanych ładunków punktowych.

Dodatkowo przy użyciu aug-PROmol, z ładunkami punktowymi z badanych metod, oszacowa- łem też potencjały elektrostatyczne molekuł. Wyniki są dalekie od tych jakie są generowane przez same ładunki punktowe. Aug-PROmol szacuje potencjały elektrostatyczne w sposób zbli- żony do metod chemii kwantowej, uwzględniając anizotropię od wolnych par elektronowych czy elektronów pi.

Aug-PROmol może być niezwykle pomocne w dziedzinie modelowanie molekularnego. Model ten pochodzi z uproszczonego formalizmu Hansena-Coppensa poszerzonego o ładunki punktowe dopasowywane do potencjału elektrostatycznego. Użycie tych dwóch, dobrze znanych i szeroko stosowanych metod zapewnia możliwość łatwej aplikacji w już istniejących technikach. Pro- stota modelu zapewnia także wysoką przenaszalność w obrębie atomów będących w podobnych środowisku chemicznym. Średni błąd kwadratowy dla dimerów molekularnych, przy odległości równowagowej wynosi ok. 1,0 kcal mol

⁻¹

w odniesieniu do DFT-SAPT/B3LYP/aug-ccpVTZ.

Odtwarzane potencjały elektrostatyczne są zbliżone do tych otrzymywanych przy użyciu me- tod chemii kwantowej. Ograniczenie się do sferycznych gęstości elektronowych pozwoliło na opracowanie analitycznej wersji aug-PROmol, dzięki której możliwe jest pominięcie kosztownych obliczeniowo procedur całkowania. O ile mógłby pojawić się zarzut, że ładunki punktowe dopa- sowywane do potencjału ab initio wpisują się we wspomnianą wyżej procedurę parametryzacji i weryfikacji na tym samym zestawie molekuł. O tyle aug-PROmol ze zbankowanymi ładun- kami punktowymi jest już całkowicie niezależne od zadanych molekuł. To wszystko sprawia, że aug-PROmol może być ważnym narzędziem w modelowaniu molekularnym.

Przykładem może być skonstruowanie funkcji szacującej (en. scoring function) w procedu- rze dockingu. Według Langnera dokładna elektrostatyka, tzn. taka, która zawiera też efekty penetracyjne, jest najlepszym wskaźnikiem stabilności kompleksów, gdy współrzędne obarczone są dużym błędem. W jego pracy przetestował wiele członów wchodzących w skład E

_int

oraz sam człon elektrostatyczny. Wyniki wskazywały, że w okół geometrii równowagowej największy succes rate miała sama elektrostatyka.

Drugim przykładem zastosowania aug-PROmol jest bezpośrednie zastosowanie modelu do pól siłowych jako gotowy człon E

es

. Do tego potrzebne byłoby dokończenie parametryzacji analitycznego odpowiednika aug-PROmol. Adaptacja do istniejącego pola siłowego nie powinna nastręczyć większych problemów. Na pewno w takich przypadkach istnieje potrzeba ponownej parametryzacji innych członów odpowiedzialnych za oddziaływania międzycząsteczkowe.

9

Abstract

Classical force fields are still in mainstream of biomacromolecular studies. Empirical parameters applied to functions known from classical physics provides less computational effort, but also are encumbered with large errors. The Coulomb law applied to point charges is not enough to estimate the electrostatic properties with a good accuracy. At shorter distances, charge anisotropy and charge penetration start playing signif- icant role. An accurate, transferable and simple model, which can be applied in large biosystems is still a demanding task.

Results obtained with UBDB model showed a dominance of charge penetration, over other effects, at the equilibrium geometry and around.

UBDB produces errors at reasonable level, however, its complexity em- phasized a demand for a more robust model.

The aug-PROmol model of an electron density became the solution.

It originates from the combination of simplified Hansen-Coppens model of electron density with point charges from restrained fitting to elec- trostatic potential. These well known and widely used theories let the aug-PROmol model be ready for applications in already existing tech- niques. The simplicity of the model provide a good transferability among atoms in similar chemical environment.

Aug-PROmol estimates the electrostatic interaction energies includ- ing the penetration effects. The root mean square error for molecular dimers, at equilibrium distance, equals ≈ 1.0 kcal mol

⁻¹

referring to DFT-SAPT at B3LYP/aug-cc-pVTZ level of theory. Estimated molecu- lar electrostatic potentials are in good agreement with these from quan- tum theory. Limiting the molecular electron density just to spherical parts created a possibility to omit the computationally costly integra- tion procedure. Corresponding error of analytical version of aug-PROmol equals 0.07 kcal mol

⁻¹

for benzene-benzene homodimer.

The aug-PROmol analytical representation may be used directly in

molecular modeling.

List of Abbreviations

AMOEBA atomic multipole optimized energies for biomolecular applications force field BS-ISA basis space iterated Stockholder analysis

CAMM cumulative atomic multiple moments DFT density functional theory

DMA distributed multipole analysis E

_es

electrostatic interaction energy E

_int

total interaction energy

E

_mtp

energy from multipole moments interactions E

_pen

penetration energy

EFP effective fragment potential

EPMM exacts potential multipole moments FF force fields

GEM Gaussian electrostatic model

HCM/HCMM Hansen Coppens multipole model HD Hirschfeld density

HD Hirschfeld iterative density MAE mean absolute error

MEP molecular electrostatic potential MSE mean signed error

ODS outer density screening

QM/MM quantum mechnics/molecular modeling RMSE/RMSD root mean square error

SAPT symmetry adapted perturbation theory

SIBFA sum of interactions between fragments ab initio computed force field SPDFG a program for fast ab initio electrostatic interaction energy calculation TAAM transferable aspherical atom model

UBDB University at Buffalo databank

aug-PROmol a promolecule model of electron density augmented with point charges

11

Contents

Summary (in Polish) . . . . 6

Abstract . . . . 10

List of Abbreviations . . . . 11

1 Introduction 14 1.0 Introduction to the intermolecular interactions . . . . 14

1.1 Born-Oppenheimer approximation . . . . 14

1.2 Intermolecular interactions . . . . 14

1.3 Multipole expansion . . . . 15

1.4 Quantum theory to calculate the intermolecular interactions . . . . 17

2.0 Electrostatic interaction energy . . . . 19

2.1 DFT-SAPT a quantum mechanic reference method . . . . 21

2.2 Distributed multipoles . . . . 21

3.0 Force fields . . . . 23

3.1 New models for an accurate electrostatic term in force fields . . . . 24

3.2 UBDB . . . . 31

2 Objectives 34 References . . . . 40 3 Interplay of point multipole moments and charge penetration for intermolec-

ular electrostatic interaction energies from the University at Buffalo pseu-

doatom databank model of electron density 42

4 A universal and straightforward approach to include penetration effects in

electrostatic interaction energy estimation 56

5 Universal method for electrostatic interaction energies estimation with charge

penetration and easily attainable point charges 66

6 Conclusions 78

7 Supporting information 80

8 Attachments 158

1. Introduction

1.0 Introduction to the intermolecular interactions

From the ancient times the discrete nature of matter was presumed. However, the centuries had passed till atoms and molecules were considered the basic building block of matter. Even without the modern techniques it should be easy to predict that there must be some kind of forces that keeps the matter together. Liquids and solids are not squeezable same as they do not vapor instantly. There must be some kind of repulsive forces at the shorter distances and on the other hand, there must be attraction at long range. This simple observation let to assume that interactions can be described by a potential curve over the inter- atomic/molecular distance.

1.1 Born-Oppenheimer approximation

Observation that electrons are much lighter than nuclei became a basis of modern chemistry. The Born-Oppenheimer approximation, and/or its more accurate version adiabatic approximation, relies on difference in the motion between electrons and nuclei. Since electrons are much faster than nuclei, the latter can be treated as fixed. With this assumption, electronic wave function depends parametrically on a nuclear coordinates - (R). The total wave function of a molecule can be split into electronic Ψ

^el_k

(r; R) and nuclei f

_kⁿ

(R) parts:

Ψ

^tot

(r, R) ∼ = X

k

Ψ

^el_k

(r; R)f

_kⁿ

(R) (1.1.1)

Solving the electronic Schrodinger equation of the k-state for fixed i-th molecular geometry R

⁽ⁱ⁾

gives the electronic energy:

H

^el

(r, R

⁽ⁱ⁾

)Ψ

^el_k

(r, R

⁽ⁱ⁾

) = E

_k

(R

⁽ⁱ⁾

)Ψ

^el_k

(r, R

⁽ⁱ⁾

) (1.1.2) Obtained E

_k

(R) for infinite set of R

⁽ⁱ⁾

become the representation of the potential energy surface (PES). Nuclei dynamics is determined by the average electronic field described by PES.

1.2 Intermolecular interactions

Without above considerations, it would not be possible to define rigorously most of the crucial chemical and physical properties. The intermolecular interaction energy can be easily defined only for the fixed geometries:

E

_int

= E

_AB

− E

A

− E

B

(1.2.3)

where E

AB

is energy of AB molecular dimers and E

A

, E

B

are energies of the separate molecules

additivity of the intermolecular interactions. There are always some contributions from three-, four-, etc. body interactions. In many cases the many body effects can be treated as a small correction to E

int

. In present work there will be no need for further consideration of this problem.

Stone classifies the forces into two categories: long and short range and enumerates them as follows: long (U ∼ R

⁻ⁿ

) electrostatic, induction, dispersion, resonance, magnetic) and short (U ∼ e

^−αR

) exchange-repulsion, exchange-induction, exchange-dispersion, charge transfer. [1]

This variety of forces are present only as a consequence of applications of particular theories.

According to Hellmann-Feynman theorem, [2] all the forces in the system can be calculated using classical electrostatics. This should accompany to everyone to prevent from drawing far-reaching conclusions.

1.3 Multipole expansion

Interaction between two atoms - A and B can be illustrated as below:

a

z

a

x

a

y

_a

i r

_ai

R b

z

_b

x

_b

y

_b

j

r

_bj

r

_ij

j

r

_bi

i

r

_aj

Figure 1.1: Coordinates centered at the nuclei of interacting atoms and distances occurring in the interaction operator.

where a and b are nuclei of atoms A and B respectively. The i and j are i-th and j-th electrons of atoms A and B respectively. To achieve the simplicity, the origins are centered at the nuclei, and the atomic units are applied. The ˆ V operator consists of operators of interactions between individuals from atoms of different molecules: nucleus-nucleus, electron-electron and nucleus-electron. For this simple system the operator ˆ V is:

V = ˆ Z

a

Z

_b

R − X

j

Z

a

r

_aj

− X

i

Z

_b

r

_bi

+ X

ij

1

r

_ij

(1.3.4)

The R is distance between the a and b nuclei and the Z is the charge of particular nucleus.

The r

ij

represents the distance between electrons of different molecules, the r

_bi

and r

aj

are distances between electron from one molecule and nucleus of another one. When the r

_ij

distance is expressed using only one origin, for example centered at nucleus a it can be written as:

r

ij

= [(x

ai

− x

ai

)

²

+ (y

ai

− y

aj

)

²

+ (z

ai

− z

aj

)

²

]

^1/2

. According to Fig. 1.1:

x

_a

= x

_b

y

_a

= y

_b

z

_a

= z

_b

+ R

(1.3.5)

15

and the r

_ij

can be expressed using both origins:

r

_ij

= [(x

_ai

− x

bi

)

²

+ (y

_ai

− y

bj

)

²

+ (z

_ai

− (z

bj

+ R))

²

]

^1/2

(1.3.6) so the operator

_r¹

ij

, from the eq. 1.3.4, can be expressed as follows:

1 r

_ij

= 1

R

n 1 − 1

R (z

_ai

− z

bj

) + 1

R

²

[(x

_ai

− x

bj

)

²

+ (y

_ai

− y

bj

)

²

+ (z

_ai

− z

bj

)

²

] o

(1.3.7) It is to be noticed, that coordinates: x

_ai

, y

_ai

, z

_ai

define the size of an A atom, similarly for x

_bj

, y

_bj

, z

_bj

and atom B. When R >> r

_a

and R >> r

_b

the ˆ V operator can be expanded into infinite series of R

⁻¹

according to the formula (1 + x)

^1/2

= 1 − 1/2x + 3/8x

²

+ ... . Expanding the

_r¹

ij

up to R

⁻³

we get:

1 r

_ij

= 1

R

n 1 + 1

R (z

_ai

− z

bj

) − 1

2R

²

[(x

_ai

− x

bj

)

²

+ (y

_ai

− y

bj

)

²

− 2(z

ai

− z

bj

)

²

] + ... o

(1.3.8)

When x

ai

= y

ai

= z

ai

= 0 then

_r¹

ij

become

_r¹

aj

: 1

r

_aj

= 1 R

n 1 + 1

R z

_bj

− 1

2R

²

(x

²_bj

+ y

_bj²

− 2z

bj²

) + ... o

(1.3.9)

and similar for

_r¹

bi

: 1 r

_bi

= 1

R

n 1 + 1

R z

_ai

− 1

2R

²

(x

²_ai

+ y

²_ai

− 2z

ai²

) + ... o

(1.3.10) When the equations 1.3.8, 1.3.9 and 1.3.10 are substituted to the eq. 1.3.4, and also assuming that x

a

= P

na

i=1

x

ai

, x

b

= P

nb

i=1

x

bi

, etc. and (Z

c

− n

c

) = q

c

where q

c

is and effective charge of an atom, we get:

V = q ˆ

_a

q

_b

R

⁻¹

−(q

b

z

_a

−q

a

z

_b

)R

⁻²

+ 

q

_b

(x

²_a

+y

_a²

−2z

a²

)+q

_a

(x

²_b

+y

_b²

−2z

²b

) 

R

⁻³

+(x

_a

x

_b

+y

_a

y

_b

−2z

a

z

_b

)R

⁻³

+...

(1.3.11) finaly, it can be rewritten into the common form as follows:

V = q ˆ

a

q

_b

R

⁻¹

−(q

b

µ ˆ

az

−q

a

µ ˆ

_bz

)R

⁻²

+ (q

_b

Q ˆ

_bz²

+ q

a

Q ˆ

_bz²

)R

⁻³

+ (ˆ µ

ax

µ ˆ

_bx

+ ˆ µ

ay

µ ˆ

_by

−2ˆµ

az

µ ˆ

_bz

)R

⁻³

+ ...

(1.3.12) where µ, ˆ ˆ Q are the multipole moments operators, dipole and quadrupole, respectively.

Here, because of the simplicity of the system (Fig. 1.1), the multipole moments operators are expressed in a simplified form:

ˆ

µ

_Az

= −z

a

Q ˆ

_Az^z

= −(2z

a²

− x

²a

− y

a²

) (1.3.13)

For many-atoms molecules the origins are usually centered at the center of the mas. Then the

multipole operators must have the non-zeroed coordinates of the nuclei.

ˆ

µ

_Az

= − X

i

z

_i

+ X

a

Z

_a

z

_a

Q ˆ

_Az2

= − X

i

1

2 (2z

²_ai

− x

²ai

− y

²ai

) + X

a

Z

_a

1

2 (2z

_a²

− x

²a

− y

²a

)

(1.3.14)

where lower case coordinates refers to particles, and capital A refers to system, where the origin is centred. Capital Z refers to the charge of the particular nucleus. The subtraction of 1/r operators in the eq. 1.3.4 should be executed only when the R >> r, where r is any of r

_ai

or r

_bj

. Whole equation 1.3.12 is called a multipolar expansion. Here, it is truncated at the R

⁻³

level. In literature it is often described as a truncation at quadrupole level or expansion up to quadrupoles (when interaction µQ and QQ are accounted). In theory the multipolar expansion should be expanded to infinity. In practise, the higher order multipoles are taken into account, the better accuracy is achieved. The greater power of

_R¹

the range of interactions are shorter, and at larger distances will be giving only small contributions to the E

_int

. Generally, multipolar expansion describe the interaction of separated systems by the series of interacting multipoles moments, a deviations from the spherical symmetry.

1.4 Quantum theory to calculate the intermolecular interactions

There are two main approaches for the calculation of the intermolecular potentials. First one, called supramolecular method, involves the equation 1.2.3 directly. The E

int

is calculated by a difference between energies of dimer and separate molecules. Though its simplicity, it does not give deeper insight into the nature of intermolecular forces. To obtain it, the decomposition scheme is necessary. For good accuracy the energies should be calculated using a theory which includes the correlation of the electrons, such as Moller-Plesset, Coupled Clusters or Configura- tion Interaction. The level of complexity and as well as the need to use a good basis-set make this approach hardly available for larger molecules.

The second one uses the perturbation approach, with the assumption that interactions are relatively weak. Another far-reaching assumption is good separation of interacting molecules, that the whole system can be described by a product of two wave functions of isolated molecules.

It is called a polarization approximation. For ground state it is:

ψ

₀⁽⁰⁾

= ψ

_A,0

ψ

_B,0

(1.4.15)

it is nothing but statement that we can assign a set of electrons to molecule A and B separately.

Then function in this form is eigenfunction of Hamiltonian H

⁽⁰⁾

= H

^A

+ H

^B

. Then:

H

⁽⁰⁾

ψ

₀⁽⁰⁾

= (E

^A

+ E

^B

)ψ

₀⁽⁰⁾

(1.4.16) Naturally, this can work only if the molecules are located at large distances. The cases at shorter distances will be discussed later. According to the Rayleigh-Schroedinger perturbation theory, the interaction energy for first and second order is given by:

E

_int

= E

_pol⁽¹⁾

+ E

_pol⁽²⁾

=< ψ

₀⁽⁰⁾

| ˆ V ψ

₀⁽⁰⁾

> + X

nA

X

nB

| < ψ

A,n_A

ψ

B,n_b

| ˆ V ψ

A,0

ψ

B,0

>

(E

_A,0

− E

A,nA

) + (E

_B,0

− E

B,nb

) (1.4.17)

17

naturally (n

_A

n

_B

) 6= (0, 0), and operator ˆ V is the operator from eq. 1.3.4, generalized to many nuclei molecules.

The meaning of the first part will be introduced further in the text. The second part is usually divided into another three:

E

₀⁽²⁾

= X

n_A=0

X

nB6=0

· · · + X

nA6=0

X

n_B=0

· · · + X

nA6=0

X

nB6=0

· · · (1.4.18)

where the particular parts are usually presented as follows:

E

₀⁽²⁾

= E

_ind

(A → B) + E

ind

(A ← B) + E

disp

(1.4.19) 1.4.1 Induction energy

The first two parts of eq. 1.4.19 describe the induction energy (E

_ind

). It occurs when permanent multipole from one molecule induce the multipoles at the second one. Hereby, both permanent and induced multipoles interact with each other. The strength depends on the polarizabilities of the second molecule. The eq. 1.4.17 shows that the polarizabilities are high when the energies of excited states are close to energies of ground states. When induction energy is described in multipole expansion it is proportional to R

^−2(k+2)

, where k is the smaller from the lowest multipoles of molecule A and B.

1.4.2 Dispersion energy

Third part of eq. 1.4.19 (E

_disp

) can be treated as energy of a interaction between dynamic multi- pole polarizabilities of monomers, since there are involved only excited states of both monomers (eq. 1.4.17). In case of dispersion forces, the polarizabilities are much harder to be accurately calculated, the larger basis set are required. However, the leading contribution comes from the first non vanishing - dipole dipole term. Thus, dispersion term is proportional to the square of R

⁻³

i.e. E

_disp

≈ −C

6

R

⁻⁶

. The C

₆

is the dispersion coefficient. The E

_disp

≈ −C

6

R

⁻⁶

is approxi- mation of orientation-averaged dispersion energy and it is the bigger, the more polarizable atoms are. Here, the strength depends on the polarizabilities of both molecules. Dispersion interaction is always present unless one of the molecules has no electrons at all.

Many approximations are present for the dispersion energy. One worth to mention, is Drude model, [3] since it is often used in the area of molecular modeling. It relies on a one-dimensional harmonic oscillator. The interaction between electron and nucleus of a particular atom is treated as classical potential of harmonic spring. Movement of an electron, around the nucleus, causes the presence of a dipole moment. So obtained atomic dipole moment correlates with the others, causing lowering the energy of the system.

1.4.3 Exchange-repulsion energy

It must be noticed, that the second order correction, describes always attractive forces. Before

the first part of the equation 1.4.17 will be discussed, it is necessary to introduce another impor-

tant term. Above considerations were following the assumption of good separation between the

function of a system must be antisymmetric due to the Pauli condition of antisymmetry of the total wave function:

Z

Ψ

^A

(1, 2, ..., n

_A

)

^∗

Ψ

^B

(1

⁰

, 2

⁰

, ..., n

⁰_B

)

^∗

... × Ψ

^A

(1

⁰

, 2, ..., n

_A

)Ψ

^B

(1, 2

⁰

, ..., n

⁰_B

)dτ (1.4.20) Applying antisymmetry into the eq. 1.4.17 creates several new correction terms to the E

_int

. In the first order, when the single determinant wavefunction is applied, the correction is described as exchange-repulsion and consists of two parts:

E

_exch−rep

= E

_exch

+ E

_rep

(1.4.21)

The first part of eq. 1.4.21 is attractive. It presents the energy associated to the motion of the electrons which are free to move over whole interacting system. The second term describes the repulsion between electrons fo the same spin. According to the Pauli principle they cannot occupy the same region. Obviously, the E

_exch−rep

arises between interacting closed-shell systems.

In such a case, chemical bonding cannot emerge and, as a consequence of further approaching of the interacting systems, the E

_int

must raise fast. The exchange-repulsion energy can be approximated with the exponential function:

E

_exch−rep

≈ K e

^{[−α(R−ρ)]}

(1.4.22)

where K, α and ρ are arbitrary parameters, responsible for the shape of the function. In many theories the E

_exch−rep

is presented as just E

_exch

. It must be noted that, exchange-repulsion force occurs already in the first order correction and further in higher orders.

2.0 Electrostatic interaction energy

Finally, leaving behind the exchange-repulsion forces, we can define the first term of the eq.

1.4.17:

E

es

=< Ψ

^A

Ψ

^B

|V |Ψ

^A

Ψ

^B

> (2.0.1) as a Coulomb interaction between charge distributions of separated molecules. The name electro- static interaction refers to the classical Coulomb interaction of static point charges. It is highly visible when the equation 2.0.1 is rewritten as the interaction between unperturbed charge dis- tributions:

E

es

= Z Z

ρ

^tot_A

(r

_A

) 1

|r

A

− r

B

| ρ

^tot_B

(r

_B

)dr

_A

dr

_B

(2.0.2) and total charge distribution for monomer A may be presented as below:

ρ

^tot_A

(r) = X

a∈A

Z

_a

δ(r − R

a

) − ρ

A

(r) (2.0.3)

The first part of eq. 2.0.3 represents the positive potential from the nuclei at R

_a

positions, whereas the second one represents the electron charge distribution. According to eq. 2.0.1 the E

_es

is an average value of the interaction operator, when the approximation of the Ψ

⁽⁰⁾₀

as a Ψ

_A,0

Ψ

_B,0

is valid. The eq. 2.0.2 can be explicitly divided into the nuclei-nuclei, electrons-nuclei

19

and electrons-electrons contributions:

E

_es

= X

a∈A

X

b∈B

Z

_a

Z

_b

R

_ab

+

Z

A

ρ

_A

(r

_A

)V

_B^nuc

(r

_A

)d(r

_A

)

+ Z

B

ρ

B

(r

B

)V

_A^nuc

(r

B

)d(r

B

) + Z

A

Z

B

ρ

A

(r

A

)ρ

B

(r

B

)

|r

A

− r

B

| dr

A

dr

B

(2.0.4)

where the Z is a charge of an particular nucleus, ρ(r) are electron densities and the V

^nuc

are potentials of the nuclei. The A or B indexes refers to a particular molecule.

On the other hand, at large distances R, the ˆ V operator can be applied as in the form shown in the eq. 1.3.12. The first order correction uses the unperturbed wavefunction of the isolated systems. Moreover the interaction operator consists of operators depending on electrons A or B separately. Then, it is easy to divide eq. 2.0.1 into two: A and B parts and perform the integration of both systems separately. Extending the ˆ V operator to larger systems, to many nuclei molecules A and B, the equation 2.0.1 can be transformed into the form:

E

_mtp

= q

_A

q

_B

R

⁻¹

− (q

B

µ

_A,z

− q

A

µ

_B,z

)R

⁻²

+(q

_B

Q

_A,z2

+ q

_A

Q

_B,z2

)R

⁻³

+(µ

_A,x

µ

_B,x

+ µ

_A,y

µ

_B,y

− 2µ

A,z

µ

_B,z

)R

⁻³

+ ...

(2.0.5) Now the µ, Q, etc. are averaged values of the operators in the ground state and they are called point multipole moments of a particular molecule. The electrostatic interaction energy can be calculated by the simple summation over interacting point multipoles from both molecules A and B. Since it relies only on the simple interactions between multipoles it will be abbreviated as E

mtp

.

In most of the tested approaches the multipole moments expansion is described using the interaction tensor notification:

E

mtp

= X

a∈A

X

b∈B

T q

a

q

_b

+ T

^α

(q

a

µ

^α_b

− q

b

µ

^α_a

) + T

^αβ

( 1

3 q

a

θ

^αβ_b

+ 1

3 q

_b

θ

^αβ_a

− µ

^αa

µ

^α_b

) + ... (2.0.6) where T

^αβ

is interaction tensor: ∇

α

∇

β

R

⁻¹

, etc.. R and remaining symbols correspond to these introduced in eq. 2.0.5, with the modifications described in the 2.2 section.

It is to be noted, that the eq. 2.0.5 was constructed with the assumption of the good separation of interacting systems. In theory, this assumption is impossible to maintain since the wavefunctions extend to infinity. However, at larger distances this error is not significant.

When molecules are approaching each other the error is the greater, the closer molecules are.

At some distances the electron clouds cannot be described by dimensionless elements and the error becomes significant. In such a case E

_es

is not identical with E

_mtp

. This error it is called a penetration energy and is defined as a difference between exact electrostatic interaction energy and energy from multipole moments interactions:

E

_pen

= E

_es

− E

mtp

(2.0.7)

When the separation of the systems is not fulfilled the exact electrostatic interaction energy

should be calculated directly from the eq. 2.0.2 or estimated in another way.

2.1 DFT-SAPT a quantum mechanic reference method

Along with the progress of computational power of the super fast computers the symmetry adapted perturbation theory (SAPT) becomes more and more popular. There are couple of the- ories which uses the perturbation scheme with function adapting the antisymmetry. However, SAPT refers to particular one proposed by Jeziorski et al.. [4] Originally it is very complicated method which uses high-order perturbation theory to account the correlation corrections. To simplify its use, Williams and Chabalowski proposed a combination with Density Functional Theory (DFT), where problem of the correlation effects is set in the correlation-exchange func- tional. [5, 6] The E

_int

calculated up to second order correction with the DFT-SAPT looks as below:

E

int

[XC] = E

_es⁽¹⁾

+ E

_exch⁽¹⁾

+ E

_ind⁽²⁾

+ E

_exch⁽²⁾ _−ind

+ E

_disp⁽²⁾

+ E

_exch⁽²⁾ _−disp

(2.1.8) where [XC] denotes that total interaction energy relies on a exchange-correlation functional.

Hesselmann and Jansen improved the DFT-SAPT performance by the usage of asymptotically corrected functional PBE0AC. [7]

DFT-SAPT by its functionality is very good approach to become a reference method. It provides less complex theory than other SAPT variants. The interaction energies are in excellent agreement with the CCSD(T). [8] However, to calculate intermolecular interactions with a good accuracy there is a need to use larger basis set, at least triple zeta. Finally, the most valuable fea- ture of any SAPT theory is the decomposition scheme (eq. 2.1.8) of interaction energy accounted by the definition.

For the purpose of this work the first part of the eq.2.1.8, the E

es⁽¹⁾

, will be treated as a reference for electrostatic interaction energy estimation.

2.2 Distributed multipoles

The most rough version of multipolar expansion (molecular multipole moments - MMM) relies on a multipole moments representation of an electron density for whole molecule. In most of such cases, the origin is centered at the center of the mas of molecule. Then, whole MMM is expanded about this point. If there is a possibility to close both molecules in non overlapping spheres the MMM may converge. However, in most of the important situations (from the chemical or biological point of view), the interacting molecules are too close and the molecular multipole moments approximation must diverge.

Moreover, such an approach strongly relies on the choice of the origin, since only the first non vanishing term is origin independent. Finite numbers of multipoles can be infinite at other origin. Centering an origin at the atomic ion makes that the charge distribution can be described be a single monopole, whereas other center force to use infinite numbers of multipoles.

One of the solution may be the partitioning the molecule into smaller contributions. Usually the most obvious choice is the atom partitioning. Then, the multipole moments can be expanded separately about the nuclei. Also, instead of single atom, it can be group of atoms or a bond.

This approach should improve the convergence.

In distributed multipole analysis (DMA) the molecular charge density is divided into several pieces. [9] Each piece is characterized by the multipole moment expansion. The molecular charge

21

density can be presented in term of primitive Gaussians [10]:

ρ(r) =

AOs

X

µ,ν

P

_µν

χ

_µ

(r − R

a

)χ

_ν

(r − R

b

) =

AOs

X

µ,ν P G∈ν

X

u

P G∈µ

X

t

[P

_µν

P

_ut⁰

φ

_ut

(r − R

⁰z

(α

_u

+ α

_t

))] (2.2.9)

where φ

_u

(r − R

a

, α

_u

) are primitive Gaussians (PGs) centered on atom a with contraction ex- ponents α

_u

. The P

_ut⁰

are products of contraction coefficients of primitive Gaussians. The PGs are used to build the basis functions χ(r − R

a

), where P

µν

is the RHF density matrix element for atomic orbitals (AOs): µ and ν. R

_z

is Gaussian overlap point. Expressly, it is centered at the line between R

_a

and R

_b

(α

_u

R

_a

+ α

_t

R

_b

/α

_u

+ α

_t

). Apart from the radial distributions the PGs contains also the angular part represented by spherical harmonics. At the R

z

center the spherical harmonics are rank from 0 to l

_a

+ l

_b

, where l

_a

is rank at center a. Each products of primitive Gaussians, and further each piece of charge density, can be represented by a set of multipoles. The greater the number of z centers the more accurate charge density is imitated by a multipole moments expansion. The lower number of z the expansion will diverge easier.

Although, it is to be stressed that above analysis is much dependent on a chosen basis set.

Another similar approach is present for example in the cumulative atomic multipole moment (CAMM) model. [11] The molecules are divided into the atomic regions. The multiple moments are obtained directly from a wavefunction using the Muliken population analysis. [12] The au- thor stress the facts that the model can be easily referred to experiment, for example X-ray crystallography.

The other multiple models can be adapted from, for example, crystallography. A popular method is stockholder partitioning introduced by Hirshfeld. [13] It relies on a partitioning the electron density of a molecule to its atomic contributions according to densities of the isolated atoms. The partitioned electron densities of an a atom are weighted:

ρ

_a

(r) = ρ

⁰_a

(r)

ρ

_promol

(r) ρ(r) (2.2.10)

the ρ

⁰_a

(r) is spherically averaged electron density of an atom and ρ

_promol

(r) is a promolecule model. [14] Promolecule model of a molecular electron density relies on a sum of spherical, neutral, free atomic electron densities. Point multipole moments on atom a can be obtained from ρ

_a

(r).

Another very popular approach from crystallography is the Hansen-Coppens multipole (HCM) model of an electron density: [15]

ρ

a

(r) = P

c

ρ

c

(r) + P

v

κρ

v

(κr) + κ

⁰

X

4

l=0

R

_l

(κ

⁰

r) X

l m=−l

P

_lm

Y

_lm

(θ, φ) (2.2.11)

where the molecular electron density is divided into pseudoatomic contributions (ρ

_a

). The pseu- doatomic densities are partitioned into three: core (c), valence (v) and aspherical (lm) parts.

First two uses the free-atom, spherically averaged Hartree-Fock densities normalized to one elec-

tron.

N

_s

ρ

_s

(r) = 1 4π

µs

X

j=1

N

_j

νj

X

i=1 νj

X

k=1

C

_j,i

C

_j,k

× h(2ζ

j,i

)

^(2nj,i+3)

(2ζ

j,k

)

^(2nj,k+3)

(2n

_j,i

+ 2)!(2n

_j,k

+ 2)!

i

1/2

×r

^(nj,i+nj,k)

exp[ −(ζ

j,i

+ ζ

_j,k

)r]

(2.2.12)

where s stands for the core or valence shells, µ

_s

is the number of atomic orbitals, ν

_j

is the number of basis-set functions for orbital j, N

_j

is the occupancy of the jth orbital, N

_s

= P

µs

j=1

N

_j

. n

i

= n

_k

= n − 1 are the exponents of r, where n is the principal quantum number. The values of expansion coefficients C

_j,i

, C

_j,k

and the exponents ζ

_j,i

, ζ

_j,k

are taken from the self- consistent-field (Hartree-Fock) energy-optimization of Clementi and Roetti for the ground-state of chemical elements from He up to Xe. [16] For the hydrogen atom, single Slater function ρ

v

= (ζ

³

π)exp( −2ζr) with ζ = 1.0 a

⁻¹0

(a

₀

= 0.529177 Å) is used to describe the electron density.

Third part of eq. 2.2.11 refers to aspherical deformation represented by a multipolar expansion.

The Y

lm

are spherical harmonics, R

l

is the radial function R

l

(r) =

_(n^ζnl+3

l+2)!

r

^nl

exp( −ζ

l

r), the P

c

, P

_v

and P

_lm

are populations of core, valence and multipoles, respectively. The κ and κ

⁰

are expansion/contraction variable radial parameters.

The HCM became a basis for a new transferable aspherical atom model (TAAM). [17] Pre- viously, to refine crystal structure only the independent atom model (IAM) could have been applied. IAM relies on a assumption that most of electron density is located around atomic cores and centered at nuclei. The plainness of this model reduce its usage just for refinement of the average atomic position and atomic displacement parameters. Well known artifacts are, for example, underestimation of X-H bond distances, and overestimation of atomic displacement parameters.

TAAM can replace IAM in procedure of X-ray data refinement at standard resolution. It is to be noted, that HCM model uses multipolar expansion up to hexadecapoles to account the aspherical deformations of electron density. Such a representation allows to account subtleness of electron density which are highly visible for example at the deformations map. The usage of TAAM at standard resolutions, improves the statistical parameters indicating the consistence of modeled and experimental data. Also, many other parameters, as average geometries, X-H bonds, atomic displacement parameters are ameliorated.

The TAAM was created due to the observations that HCM parameters are almost identical within similar chemical environment. Thus, electron density of a molecule can be divided into contributions of its atoms and represented by the set of HCM parameters. The parameters from a particular pseudoatom can be transfer to reconstruct the electron density around a corresponding pseudoatom. Such a reconstructed electron density can be used not only for its main purpose, a X-Ray data refinement. It may be used for direct calculation of the molecular properties, for example the E

es

applying the eqs. 2.0.2 and 2.0.3, or it may serve as a source of point multipole moments.

3.0 Force fields

Far different approach is present in a class of methods called force fields. Force fields rely on empirical parameters applied to simple functions known from the classical physics. The most

23

simplified versions define potential energy function as below: [18]

V (r) = X

bonds

k

_b

(b −b

0

)

²

+ X

angles

k

_θ

(θ −θ

0

)

²

+ X

torsions

k

_φ

[cos

²

(nφ + δ) + 1] + X

non− bonded

hA

ij

r

_ij¹²

− C

_ij

r

_ij⁶

i + q

_i

q

_j

r

_ij

(3.0.1) The first three parts represent, respectively, the 1-2 (bonds), 1-3 (angles) and 1-4 (torsions) interactions, where b, b

₀

, k

_θ

, θ

₀

, etc. are empirical parameters. The non-bonded term apply to the interaction between 1-5 neighbors and further. It consists of two components: van der Waals (E

_vdw

) and electrostatic (E

_es

). The former component relies on a model of exchange-repulsion mixed with dispersion forces, a so-called 6-12 potential. The latter one relies on a simple Coulomb law with point charges.

Naturally, for intermolecular interactions the non-bonded terms matter the most. The main advantage of such a simplification is extreme acceleration in computational time, comparing to quantum methods. Since there is no need for wavefunction calculation, integration procedure, etc. it can be applied to large scale systems as biomacromolecules. However, such a robust approach, must be encumbered with a large error. Additionally, the eq. 3.0.1 neglects the induction energy, charge anisotropy or penetration effects. At the very beginning there have been many attempts to improve the potential energy.

It must be noticed, that main purpose of force fields is to predict the total interaction energy rather than its individual contributions. However, it is not by accident that non-bonded terms are compatible with the perturbation approach. It is specially visible with the modern force fields, which try to account effects listed above. These force fields are developed in a way to loose the empiricism and to mimic the SAPT scheme. Thus, it is justified to compare such components directly to theirs equivalents from SAPT theory.

3.1 New models for an accurate electrostatic term in force fields

As a definition of new generation force fields McDaniel and Schmidt give the relation to quantum theory in terms of physical interactions. [19] In other words modern force fields should leave the Lennard-Jones and Buckingham approximation. There, the empirical parameters are responsible for description of multiple physical interactions at once. Such a superposition makes the descrip- tion of physical properties rather ineffective. Also, in terms of transferability, authors, noticed that only parameters which are based on a physical ground, as for example on SAPT, can be transferable to other systems.

Such an approach is developed, for example, by Szalewicz and Podeszwa. [20, 21] The basic potential relies on electrostatic, exchange and dispersion forces:

V

_ab

(r) = q

a

q

_b

r + A

_ab

exp( − r

B

_ab

) − C

₆^ab

r

⁶

(3.1.2)

where r is interatomic distance between atoms a and b. A, B, and C are exchange and dispersion coefficients. Whole equation 3.1.2 is fitted to the DFT-SAPT data.

In recent years the number of models, as above, has increased rapidly. Only few of them are

presented in this work. It is so, not only due to its brief character. This studies are dedicated to

3.1.1 Hirshfeld partitioned atomic multipoles

A method proposed by Elking et al. applies the partitioning method (HD) defined by the eq. 2.2.10 to study the atomic multipoles for particular molecules. [22] The densities can be obtained using any of the theory: HF, DFT, MP2 or CCSD. Slight modification of this approach is accounted in iterative partitioning (HDI). It relies on a replacement of isolated atomic charge densities by fractional charge densities. The latter one, in a iterative way, combines the atomic charge with the ab initio atomic charge densities. The authors adapted so-obtained charge densities for calculation of atomic properties, with the special attention paid on higher multipole moments. The results were compared to values taken from ab initio method. The molecular dipole moment and molecular electrostatic potentials obtained with HD and HDI were quite positively verified on a set of small molecules. It is proposed to use point charges obtained with the HD or HDI method since the most popular force fields uses only these parameters to calculate the electrostatics. However, authors noticed that molecular dipole moments calculated using only the HD/HDI point charges produce relatively large errors, so the possible usage in the force fields may be limited. Apart from that, all the properties calculations using the proposed two methods require time consuming calculations. For this reason the HD and HDI can be used only as a support in the area of force filed development.

3.1.2 Effective fragment potential

Efective fragment potential with multipole moments is another strictly theoretical approach introduced by Bertoni et al.. [23] Authors declare no single empirical parameter. The main purpose of the method is to calculate total interaction energy. Here, E

_int

consists of four pairwise additive contributions: Coulomb, exchange-repulsion, dispersion, charge-transfer and one many- body term: a polarization. For calculations of Coulomb, polarization and charge transfer energy the multipole moments representation of a Coulomb field is used. Gordon et al. did a wide studies of two multiple moments sources: the distributed multipole analysis (DMA referring to eq. 2.2.9) and basis-space iterated stockholder atoms (BS-ISA; referring to eq. 2.2.10). Both are roughly described in the 2.2 section, however, authors applied several variations. The energy calculations with EFP using both DMA and BS-ISA require two steps. First one is initial quantum mechanic calculations to generate the parameters of the EFP method. Second one is proper calculations of requested molecular properties. The Coulomb interaction energy (E

_es

) are taken from the multipole moments interactions. However, for systems at shorter distances, a charge penetration correction is applied. In the EFP, it is based on a overlap of localized molecular orbitals (LMOs).

The overlap integral was estimated to approximate the E

_pen

. The E

_es

energies were estimated for group of small organic molecules and compared to corresponding SAPT2+/aug-cc-pVTZ energies components. The mean unsigned error (MUE) for particular EFP variations ranges from 1.595 kcal mol

⁻¹

to 2.545 kcal mol

⁻¹

. The penetration energy correction is at the same level for all tested methods. The complexity of the calculations and the quite high level of errors make the EFP method hard to implement in molecular modeling. However, it must be noticed, the EFP method is targeted on the whole E

_int

estimation. The lowest MUEs for E

_int

equals 1.210 kcal mol

⁻¹

(2.842 kcal mol

⁻¹

the highest).

25

3.1.3 Multicenter multipolar expansion

The main aim of multicenter multipolar expansion (MME) method is to account for both charge- anisotrophy and penetration correction. [24] As the name of the method suggests, it relies on a multipoles moments centered not only on the nuclei, but also in the middle of the chemical bonds.

The multipole moments are expanded up to octupoles from ab initio Hartree-Fock at 6-31G**

and 6-31G**+ wave functions for atomic and bond centers respectively. [25] The character of MME is similar to EFP as it follows the QM/MM approach. The penetration energy correction is evaluated with a damping function (1 −exp(−α

i

r)) proposed by Freitag et al.. [26] The proposed elaborate expression for total electrostatic interaction energy between two centers is:

E

_es^{M M E}

∼ =

N1

X

i=1 N2

X

j=1

( q

i

q

j

r

_ij

"

1 − (e

^−αirij

+ e

^−αjrij

)

2 − (α

²_i

+ α

²_j

)

2(α

_i²

− α

²_j

) (e

^−αjrij

− e

^−αirij

)

#

+ q

_i

Z

_j

r

_ij

[1 − e

^−αirij

] + q

_j

Z

_i

r

_ij

[1 − e

^−αjrij

] )

+

N1

X

i=1 N2

X

j=1

Z

_i

Z

_j

r

_ij

+

X

3 n=1

E

_0,i^n,j

+ X

3 n=1

E

_n,i^0,j

+ X

3 n=1

E

_n,i^1,j

+ X

3 n=1

E

_0,j^n,i

+ E

_2,i^2,j

!

(3.1.3)

where i and j are expansion centers, r

_ij

is distance between them, q

_i

and q

_j

is the charge as- sociated to i and j. Z is the nuclear charge (Z=0 for center at bond) and N

₁

and N

₂

are the number of expansion centers. The second summation concerns interactions between higher mul- tipole moments. Only the first summation of the charge-charge and charge-nuclei interactions, accounts the charge penetration error. Finally, the α are the damping parameters, which are crucial to the MME method. The clue of MME method is to achieve whole set of the optimal α parameters. They are obtained in the procedure of minimization of the obtained molecular electrostatic potential to corresponding one from ab initio method. The procedure runs over given number of grid points (N

p

). The larger N

p

the more complex optimization is. Authors dedicated a vast amount of work to find the best optimization procedure, which made the issue even more complex. The optimizations and tests were performed over 20 aminoacids and 8 small organic molecules at the equilibrium distance. The yielded error, comparing to ab initio method is low and equals 0.26 kcal mol

⁻¹

(1.0 kcal mol

⁻¹

for undamped MME).

3.1.4 Screened charges

Thrular et al. [27] proposed charge model called the outer density screening (ODS). ODS relies on

a partitioning of the electron density into two components. The point charges in the inner region

which mimics the potential resulting from nuclei and closely surrounding electrons. The outer

region is represented by smeared charge of n

_screen

electrons. The smeared charge is described

by a single spherical Slater-type orbital (STO), normalized to n

_screen

. The n

_screen

equals 1 -

q for hydrogens (q is a partial charge on hydrogen), n

_screen

= 1 or 2 for the rest of remaining

nonmetals (depends on studies) and n

_screen

= 0 for metals. The ζ exponents used in STO were

parametrized with QM/MM studies for hydrogen atom and by fitting to interaction energies

scheme for other elements (C, N, O F, Si, P, S and Cl). Authors emphasized that without higher