DIPLOMOVÁ PRÁCE

HAPTICKÝ MODEL INTERAKCE

BIOMOLEKUL

Jiří Filipovič

2008

Fakulta informatiky Masarykova univerzita

Prohlašuji, že tato práce je mým původním autorským dílem, které jsem vypracoval samo-statně. Všechny zdroje, prameny a literaturu, které jsem při vypracování používal nebo z nichčerpal, v práci řádně cituji s uvedením úplného odkazu na příslušný zdroj.

Poděkování:Děkuji vedoucímu své diplomové práce, Mgr. Aleši Křenkovi, Ph.D., za jeho cenné rady.Dále bych chtěl poděkovat své manželce Ivaně a rodičům za péči a podporu a METACentruza poskytnutí výpočetních kapacit.

Shrnutí

V rámci spolupráce Laboratoře interakcí člověka s počítačem na FI MU, Superpočítačovéhocentra ÚVT MU a Národního centra pro výzkum biomolekul PřF MU vzniká aplikace, jejímžcílem je simulovat interakci dvou molekul a umožnit uživateli pomocí haptického zařízeníjednu z těchto molekul přímo ovládat a cítit sílu danou silovými interakcemi mezi molekulami.Obsahem diplomové práce je nalezení, implementace a popis metod, které urychlí běh této

aplikace tak, aby bylo na dostupném hardware možno dostatečně precizně simulovat interakcevětších molekul a umožnit tak její nasazení na reálné chemické problémy.Metody využité v mé práci jsou mimo jiné:

• paralelizace algoritmů

• algoritmy pro rychlou detekci kolizí

• zrychlení výpočtů pomocí předpočítání

• dávkování geometrie

• dynamická úroveň detailů

Klíčová slova

molekulový doking, haptika, optimalizace, paralelizace, adaptivní mřížka, interpolace

Obsah

1 Úvod 3

2 Výpočetní chemie 42.1 Molekulové modelování . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1.1 Molekulová mechanika . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1.2 Molekulová dynamika . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2 Dokování molekul . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.3 Matematický model interakce molekul užitý v naší aplikaci . . . . . . . . . . 6

2.3.1 Van der Waalsova energie . . . . . . . . . . . . . . . . . . . . . . . . . 62.3.2 Elektrostatická energie . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.3.3 Vztah van der Waalsovy a elektrostatické energie . . . . . . . . . . . . 7

3 Aplikace pro haptickou interakci molekul 83.1 Srovnání s podobnými pracemi . . . . . . . . . . . . . . . . . . . . . . . . . . 93.2 Technologie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

3.2.1 Haptika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93.2.2 Stereovize . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103.2.3 Dynamika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113.2.4 Semiimplicitní numerická integrace . . . . . . . . . . . . . . . . . . . . 11

3.3 Architektura aplikace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.4 Možné budoucí rozšíření . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4 Optimalizace zobrazovací části 144.1 Grafický engine OGRE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144.2 Dávkování geometrie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.3 Dynamická úroveň detailů . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.4 Využití programovatelné GPU . . . . . . . . . . . . . . . . . . . . . . . . . . 164.5 Shrnutí . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

5 Optimalizace výpočetní části 185.1 Výpočet van der Waalsových interakcí . . . . . . . . . . . . . . . . . . . . . . 18

5.1.1 Nalezení blízkých atomů . . . . . . . . . . . . . . . . . . . . . . . . . . 185.1.2 Paralelní výpočet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

5.2 Výpočet elektrostatických sil . . . . . . . . . . . . . . . . . . . . . . . . . . . 245.2.1 Regulární prostorová mřížka . . . . . . . . . . . . . . . . . . . . . . . 255.2.2 Adaptivní prostorová mřížka . . . . . . . . . . . . . . . . . . . . . . . 295.2.3 Možná vylepšení adaptivní mřížky . . . . . . . . . . . . . . . . . . . . 38

1

5.3 Detaily implementace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 395.3.1 MPI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 395.3.2 Architektura paralelní aplikace . . . . . . . . . . . . . . . . . . . . . . 405.3.3 Problém využití všech procesorů . . . . . . . . . . . . . . . . . . . . . 40

6 Měření výkonu 426.1 Zobrazování . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

6.1.1 Dávkování geometrie . . . . . . . . . . . . . . . . . . . . . . . . . . . . 436.1.2 Dynamická úroveň detailů . . . . . . . . . . . . . . . . . . . . . . . . . 44

6.2 Detekce kolizí . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 446.3 Paralelní běh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

6.3.1 Prostředí . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 466.3.2 Metodika měření . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 466.3.3 Škálování algoritmu . . . . . . . . . . . . . . . . . . . . . . . . . . . . 466.3.4 Vliv latence sítě . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 476.3.5 Srovnání CPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 486.3.6 Využitelnost simulace na clusterech . . . . . . . . . . . . . . . . . . . 49

6.4 Regulární versus adaptivní mřížka . . . . . . . . . . . . . . . . . . . . . . . . 49

7 Závěr 52

2

Kapitola 1

Úvod

Výpočetní chemie nahrazuje laboratorní experiment simulací. Její výhodou je mimo jiné ušet-ření nákladů na experiment, práce chemiků i její bezpečnost. Chemické výpočty zpravidlavyžadují velmi velký výpočetní výkon a mohou pracovat s velkými objemy dat a přináší takmnoho nových, často informatických problémů. Při modelování komplexních chemických dějůlze využít pokročilé metody interakcí člověka s počítačem.Jedním z úkolů výpočetní chemie je molekulový doking. Jedná se o úlohu, která zkoumá

interakci velké molekuly receptoru (obvykle to bývá protein) s jinou molekulou, nazývanouligand (proteinem, nukleovou kyselinou, léčivem aj.). Molekulový doking má velmi závažnéaplikace, a to především v biochemii, např. při hledání léčiva upravujícího vlastnosti nějakéhoproteinu či při ekologické likvidaci toxických látek.Ve své práci se zabývám urychlením výpočtů aplikace pro hapticky řízený molekulový do-

king. Tato aplikace realizuje nekonvenční přístup k problému dokování malých molekul ligandudo molekuly receptoru. Zatímco obvyklé výpočetní řešení spočívá v automatickém prohledá-vání velkého množství možných poloh ligandu vůči receptoru [20] [13], naše aplikace umožňujeuživateli vidět jejich vzájemnou polohu a pomocí haptického zařízení vnímat sílu, kterou nasebe molekuly působí. To umožňuje jednak přinést do problému dokingu lidskou intuici, aletaké detailnější studium interakcí molekul. Hmat je totiž jediné pro člověka přirozené vnímánísíly.Obsahem této diplomové práce je nalezení a implementace metod, které by vedly ke zrych-

lení aplikace a umožnily na dostupném hardware simulovat interakce dostatečně velkých mo-lekul při zachování dostatečné přesnosti výpočtů. Fakt, že je aplikace ovládána v reálnémčase haptickým zařízením, klade netriviální nároky na rychlost chemických výpočtů. Zatímcovykreslování vzájemné polohy molekul stačí provádět cca 25× za sekundu, simulace silovéhopůsobení musí dosahovat okolo 1 000 iterací za sekundu, jinak by se uživateli nejevil pohybligandu dostatečně hladký [14].Ve druhé kapitole se věnuji stručnému popisu metod výpočetní chemie, které stojí za naší

aplikací. V kapitole třetí popisuji aplikaci bez optimalizací popsaných v této diplomové práci,věnuji se specifickým technologiím, které aplikace používá a uvádím předpokládané směrybudoucího rozšíření její funkcionality. Optimalizace zobrazovací části aplikace jsou popsányve čtvrté kapitole, v páté se pak věnuji rozsáhlejšímu popisu optimalizací výpočetní části.V šesté kapitole jsou uvedeny a diskutovány výsledky měření přínosu jednotlivých implemen-tovaných technik.

3

Kapitola 2

Výpočetní chemie

Výpočetní chemie, z anglického computational chemistry, je oborem chemie, který staví na te-oretické chemii a využívá moderní výpočetní techniku. Zatímco v chemii experimentální jeproveden experiment in vitro1, výpočetní chemie provádí experimenty in silico2 formou si-mulace. Do výpočetní chemie se řadí také problematika uchovávání chemických dat, jejichkategorizace, porovnávání atd., tato problematika však přesahuje rámec mé diplomové práce.Výhody počítačové chemie spočívají především v ceně experimentu (často je levnější poříditči pronajmout výpočetní zdroje, na kterých provedeme simulace a až jejich výsledek ověřímeexperimentem in vitro), v ušetřeném čase (viz předchozí argument) či ve vyšší bezpečnosti(vyhýbáme se práci s toxickými či jinak nebezpečnými látkami).Dále si z výpočetní chemie budeme všímat pouze molekulového modelování, pracujícího

s molekulou na úrovni jednotlivých atomů (kvantová mechanika, schopna postihnout chováníjednotlivých elektronů, je pro naši práci s velkými biomolekulami příliš výpočetně náročnýmnástrojem).

2.1 Molekulové modelování

Základní datovou strukturou je pro nás molekula, u níž definujeme strukturu (popisuje atomytvořící molekulu a jejich vazby) a geometrii (poloha atomů v prostoru). Na geometrii molekulyzávisí její potenciální energie. Ta charakterizuje stabilitu konformace molekuly (jejího prosto-rového uspořádání), umožňuje určit stabilní konformaci (lokální minimum funkce potenciálníenergie) i změnu konformace (dle gradientu potenciálové funkce – molekula s vyšší potenciálníenergií má v souladu s 2. termodynamickým zákonem tendenci změnit svou konformaci tak,aby přešla do stavu s nižší energií).

2.1.1 Molekulová mechanika

Molekulová mechanika simuluje vazby pomocí harmonického potenciálu a nevazebné interakcepomocí bodových nábojů a Lenard-Jonesova potenciálu [21]. Jedná se tedy o velmi hrubouaproximaci, kdy přiřazujeme mikrosvětu vlastnosti známé z makrosvěta. Přesto je tato metodačasto úspěšně používána. Její velkou výhodou je především nízká výpočetní náročnost, díky

1ve skle2v křemíku

4

které je možné v rozumném čase simulovat systémy o tisícovkách atomů. Potenciálová funkcezpravidla pracuje s následujícími složkami energií:

• energie pnutí vazeb

• energie ohýbání vazebných úhlů

• energie deformace torzních úhlů

• energie van der Waalsových interakcí

• energie elektrostatických interakcí

První tři zmíněné složky energie jsou energie vazebných interakcí. Jejich velikost ovlivňujeodchylka velikosti či úhlů od ideálu, kde mají vázané atomy nejnižší energii. Zbývající dvěsložky jsou nevazebné, neboli nekovalentní a budou probrány podrobě v kapitole 2.3.Do systému molekulové mechaniky jsou vneseny konstanty upravující vlastnosti bodových

nábojů a Lennard-Jonesova potenciálu reprezentujících atomy a harmonického potenciálureprezentujícího vazby tak, aby jejich chování co nejlépe odpovídalo reálnému chemickémusystému v požadovaném prostředí.

Hyperplocha potenciální energie

Hyperplocha potenciální energie (z anglického Potential Energy hyperSurface, PES) je graf po-tenciálové funkce Rn → R, kde n je počet souřadnic potřebných k popisu geometrie molekuly.Jedná se tedy o plochu v prostoru o n + 1 rozměrech, jejíž výška v bodě ~X = (x1, x2, ..., xn)představuje hodnotu potenciální energie pro molekulu, jejíž atomy mají souřadnice ~X. U PESjsou pro nás zajímavými místy lokální minima, která představují stabilní konformace mole-kuly. Často nás také zajímají tranzitní stavy, tj. lokální maxima na energeticky nejvýhodnějšícestě mezi dvěma minimy.

2.1.2 Molekulová dynamika

Zatímco molekulová mechanika se zabývá analýzou PES (hledáním minim, sedlových bodůaj.), molekulová dynamika se orientuje na cestování po PES. Prostorová derivace potenci-álové funkce vyjadřuje totiž sílu, která působí na jednotlivé atomy molekuly. Můžeme takpomocí Newtonových dynamických zákonů simulovat pohyb atomů molekuly. Celková energiemolekuly je pak dána jak potenciální energií, tak energií kinetickou, přičemž tyto dvě složkyse mohou vzájemně přeměňovat. Díky tomu může molekula překonávat bariéry potenciálníenergie a dostávat se do nových energetických minim.Výpočet molekulové dynamiky provádíme v diskrétních časových krocích. Vycházíme z po-

lohy atomů molekuly v čase t a na základě síly dané první derivací potenciální a kinetickéenergie počítáme nové souřadnice atomů v čase t + ∆t. Další výpočet děláme pro molekulus novými souřadnicemi atomů. Souřadnice atomů molekuly v určitém čase se nazývají snímek,posloupnost snímků pak trajektorie.

2.2 Dokování molekul

Molekulový doking hledá takovou vzájemnou pozici a konformaci dvou (popř. více) molekul,jejichž komplex je nejstabilnější, tj. systém se nachází v globálním minimu PES nebo v jeho

5

blízkosti. V dokovacích úlohách sledujeme buď interakci dvou velkých molekul (např. protein-protein, protein-nukleová kyselina), nebo velké molekuly s malou molekulou (např. protein-léčivo). Menší molekula se nazývá ligand, větší molekulu, do které ligand dokujeme, nazývámereceptor.Molekulový doking používáme na prozkoumání, zdali a za jakých okolností se konkrétní

ligand váže s receptorem (jakou energetickou bariéru musí překonat, jaká je nejpravděpodob-nější poloha a konformace ligandu při zadokování). Zpravidla chceme, aby došlo k navázániv aktivním místě receptoru. Pomocí automatického dokingu můžeme také z množiny kandi-dátů hledat vhodný ligand, který se s receptorem váže.Hledání ligandů padnoucích do aktivního místa receptoru můžeme využít např. při de-

signu léčiv, kdy často potřebujeme nalézt látku, která aktivuje či inhibuje funkci specifickéhoproteinu. Jedná se tedy o úlohu nalezení co nejlepšího (ideálně globálního) energetického mi-nima pro komplex ligand – protein, pro které platí, že ligand je v aktivním místě proteinu(tj. takové místo, kde námi požadovaným způsobem ovlivňuje jeho chování při chemickýchreakcích).

2.3 Matematický model interakce molekul užitý v naší aplikaci

Naše aplikace se zabývá dokingem, kdy je ligand vázán s receptorem nekovalentními silami. Toznamená, že v intramolekulární interakci uvažujeme pouze van der Waalsovy a elektrostatickésíly, struktura molekul zůstává zachována, neboť nevznikají ani nezanikají kovalentní vazby.Zároveň pracujeme s rigidním receptorem i ligandem – jejich geometrie se vlivem silovéhopůsobení nemění (to je však omezení současné implementace aplikace).

2.3.1 Van der Waalsova energie

Existence van der Waalsovy síly byla objevena při studiu zkapalněných netečných plynů.Tyto plyny jsou tvořeny samostatnými atomy, které spolu netvoří žádné vazby. Přesto všakv kapalné formě drží pohromadě, z čehož plyne, že zde musí být ještě nějaká síla, kteránení vazebná, ale působí mezi atomy. Podle svého objevitele byly tyto síly nazvány van derWaalsovy3.Zřejmě nejčastěji užívanou aproximací pro van der Waalsovu interakci je Lenard-Jonesova

12-6 funkce. Ta má pro 2 atomy tvar

EvdW = 4ǫ(

(σ

r

)12

−(σ

r

)6)

,

kde r je vzdálenost mezi atomy, σ součet jejich van der Waalsových poloměrů a ǫ koeficientLenard-Jonesovy funkce určující její minimum.Velikost síly, kterou na sebe působí dva atomy pak získáme prostorovou derivací výše

uvedeného vztahu:

FvdW = 4ǫ(

12σ11

r13− 6

σ5

r7

)

.

3V této práci pojmem van der Waalsovy interakce označujeme součet přitažlivé disperzní Londonovy síly[22] a odpudivých sil způsobených odpuzováním atomových jader.

6

2.3.2 Elektrostatická energie

V elektricky neutrálních molekulách je náboj distribuován nerovnoměrně. V okolí atomu se vy-tváří parciální kladné a záporné náboje. Vznikají tak dipóly, mezi kterými dochází k elek-trostatické interakci. V molekulové mechanice zpravidla aproximujeme distribuci náboje po-mocí bodových nábojů umístěných v pozici atomových jader. Energii elektrostatické interakcespočítáme z těchto bodových nábojů pomocí Coulombova zákona:

Ee =1

4πǫ0ǫr

q1 · q2r

,

kde q1 a q2 jsou náboje atomů, r je jejich vzdálenost, ǫ0 je permitivita vakua a ǫr je relativnípermitivita prostředí.Elektrostatickou sílu pro dva atomy vypočítáme ze vztahu:

Fe =1

4πǫ0ǫr

q1 · q2r2

.

2.3.3 Vztah van der Waalsovy a elektrostatické energie

Nyní si ještě povšimněme vztahu našich dvou energií nikoliv z chemického, ale z matematic-kého pohledu. V obou těchto funkcích je jedinou proměnnou poloha atomů (respektive jejichvzájemná vzdálenost ve skalárním vyjádření funkcí). Ostatní hodnoty jsou pro dané dvojiceatomů konstantní. Obě funkce mají singulární bod v r → 0, tedy tam, kde se překrývá polohaatomů. Rozdílné je však chování funkcí v závislosti na velikosti r. Van der Waalsova energiemá strmější průběh pro r < σ (díky 12. mocnině podílu σ

roproti 1. mocnině podílu q1·q2

r),

zatímco pro r > σ převládá u van der Waalsovy energie složka (σr)6 a konverguje k nule rych-

leji než energie elektrostatická. Toto pozorování je velmi důležité pro návrh metod urychlenívýpočtu sil a energií popsaných v této práci.

Obrázek 2.1: Průběh van der Waalsovy a elektrostatické energie.

Obrázek 2.1 zachycuje průběh van der Waalsovy a elektrostatické energie v okolí atomu.V levém grafu vidíme, že prudká energetická bariéra daná van der Waalsovou energií začínáve větší vzdálenosti od jádra, než u energie elektrostatické. V pravé části je zvětšena částgrafu od 5 A do 20 A, ve které je vidět pomalejší konvergence k nule u elektrostatické energie.

7

Kapitola 3

Aplikace pro haptickou interakcimolekul

Zatímco u dvourozměrných objektů lze silová pole vyjádřit poměrně dobře graficky, u trojroz-měrných objektů, jakým je například molekula, je grafická reprezentace silových polí obtížná.Haptické rozhraní je pro člověka intuitivním vyjádřením silových interakcí. Chceme-li zkou-mat vzájemné silové působení dvou těles v prostoru, dává nám haptické zařízení možnostpřímo cítit, jak se tělesa ovlivňují, aniž bychom museli jejich silovou interakci rekonstruovatz obrazové informace. Podrobněji se haptikou zabývám v kapitole 3.2.1.Námi vyvíjená aplikace využívá zařízení PHANToM1, které umožňuje uživateli trojroz-

měrný pohyb s molekulou ligandu. Pomocí silové zpětné vazby mu přináší informaci o vzá-jemném silovém působení molekul.

Obrázek 3.1: Screenshot aplikace.

Vzájemná poloha ligandu a receptoru je vykreslována na obrazovku spolu s informacemio elektrostatické a van der Waalsově energii aktuálního umístění ligandu – ty jsou jednakpřímo vypisovány, ale také zakreslovány do grafu, z něhož lze vidět, jaké energetické bariérypřekonával ligand při pohybu [8].

1www.sensable.com

8

3.1 Srovnání s podobnými pracemi

Naše aplikace pro simulaci haptické interakce biomolekul přistupuje k problému molekulovéhodokingu ne zcela obvyklým způsobem. Zatímco běžný dokovací software je založen na mini-malizaci potenciálové funkce pomocí automatického prohledávání stavového prostoru s heu-ristikami [20] [13], hapticky ovládaná interakce molekul je řízena přímo člověkem. Haptickéovládání pohybu molekul člověkem je účinnější než doking založený pouze na grafické infor-maci [3] a lze předpokládat, že člověk dokáže snáze vyloučit některé cesty, které by muselybýt prozkoumány při automatickém dokingu. Hapticky řízený doking by také mohl umožnitsnadné prozkoumání okolí dokovací pozice. Nezanedbatelný je také význam ve vzdělávání,kdy studenti s pomocí haptické simulace vidí a cítí vzájemné interakce molekul [1].Reálné využití haptiky při simulaci interakcí biomolekul v reálném čase je zatím stále ve

svých počátcích, ačkoliv idea haptického dokingu je poměrně stará (je popsána např. v [2]).Existuje několik aplikací realizujících hapticky řízený doking:

• [7] umožňuje pouze dokování sondy (reprezentované bodovým nábojem) a simulujepouze elektrostatickou interakci

• [4] simuluje pouze elektrostatickou interakci

• [6] dokáže pracovat s flexibilním ligandem, ale simuluje pouze elektrostatickou interakci

• [5] dokáže pracovat s flexibilním ligandem, ale simuluje pouze van der Waalsovu interakci

Naše aplikace se ve srovnání s ostatními silně orientuje na dostatečně přesné vyjádřenísilového pole, které jsme odvodili od silového pole AMBER[25]. Pracujeme jak s elektrostatic-kými interakcemi, které implementuje většina výše citovaného software, tak s van der Waal-sovými. S těmi pracuje pouze [5], a to ve zjednodušené podobě, kdy nepočítá přímo interakcemezi jednotlivými atomy, ale interakce mezi atomy ligandu a silovým polem receptoru, které jepředpočítáno na 3D prostorové mřížce. Jelikož nemůže být kvůli paměťovým nárokům mřížkadostatečně jemná, přidává [5] do modelu interakce ligandu a receptoru klasickou obsluha ko-lizí, která zamezuje průniku jednotlivých atomů do sebe.Všechny zde zmíněné aplikace mapují síly působící na ligand přímo. Naše aplikace využívá

klasické Newtonovské dynamiky pro ovládání ligandu – té se věnuje kapitola 3.2.3.

3.2 Technologie

3.2.1 Haptika

Haptika se v obecnosti zabývá zkoumáním hmatu. V kontextu této práce je podstatný po-jem haptické rozhraní, který označuje zařízení, jež pomocí silové zpětné vazby přináší dojemhmatu. Takové zařízení musí být v přímém kontaktu s uživatelem (např. rukavice, stylus,či přímo model nástroje, se kterým uživatel operuje), aby mohlo měřit svou polohu a pomocísíly, popř. i kroutícího momentu, na uživatele působit.

9

Zařízení PHANToM

Naše aplikace využívá zařízení SensAble PHANToM2, které je k dispozici v Laboratoři in-terakcí člověka s počítačem3 (HCI). Toto zařízení umožňuje snímat pohyb stylusu v šestistupních volnosti (poloha a orientace v prostoru) a zprostředkovává pomocí servo motorkůsilovou zpětnou vazbu ve třech stupních volnosti (v poloze hrotu).Pro ovládání zařízení používáme knihovnu libphantom4, vyvinutou v rámci laboratoře

HCI na základě ovladače popsaného v [15].

Obrázek 3.2: Zařízení PHANToM (zdroj: web HCI).

3.2.2 Stereovize

Člověk dokáže vnímat zrakem hloubku díky tomu, že přijímá každým okem mírně rozdílnýobraz. Jejich kompozicí lze rozlišit, které objekty jsou blíže (rozdíl pozorovacího úhlu a tedyi umístění je v obrazech vnímaných levým a pravým okem větší) a které dále.Takové vnímání dokážeme uživateli zprostředkovat, oddělíme-li zrakovou informaci pro

levé a pravé oko. K tomu existují dva základní přístupy – můžeme současně vysílat obraz proobě oči, nebo vysílat střídavě obraz pro levé a pravé oko a současně zatmívat to oko, pro kteréobraz není určen.V laboratoři HCI jsou k dispozici obě technologie, my jsme zvolili implementaci, kdy

současně vykreslujeme obraz pro levé i pravé oko. K zobrazení využíváme dvou projektorů,vysílajících rozdílně polarizovaný obraz na společné plátno. Uživatel a případní pozorovatelési nasadí jednoduché brýle obsahující polarizační filtry, které propustí obraz určený pouze prolevé oko do levého oka a naopak. Tím se tyto dva obrazy oddělí.Z hlediska aplikace je prostorové zobrazení realizováno tak, že zobrazujeme celou scénu

ze dvou kamer (každá reprezentuje jedno oko) a tyto dva obrazy vykreslujeme vedle sebe.Systém má zapnuté zobrazování celé pracovní plochy do dvou výstupů, takže se výstup zob-razovací aplikace rozdělí mezi dva projektory. Díky nutnosti vykreslovat obraz ze dvou kamersamozřejmě rostou nároky na výkon grafického rozhraní systému.

2http://www.sensable.com/haptic-phantom-premium.htm3http://decibel.fi.muni.cz/4http://www.fi.muni.cz/ xslaby/phantom.html

10

3.2.3 Dynamika

Reprezentace ligandu jako hmotného tělesa je pro člověka přirozená, jelikož odpovídá cho-vání těles makrosvěta. Ligand není napojen na haptické zařízení přímo, ale pohybuje se jakohmotné těleso přichycené k haptickému ovladači pomocí viskoelastické vazby (analogie pružinys tlumičem, popsána v [18]). Viskoelastická vazba umožňuje využít setrvačného momentu li-gandu při překonávání energetických bariér, kdy ligand při svém pohybu dokáže sám upravittrajektorii či rotaci tak, aby se vyhnul překážce.

Obrázek 3.3: Cesta ligandu k dokovací pozici.

Obrázek 3.3 zobrazuje tažení ligandu PHANToMem volným prostorem (zanedbáváme zdeúčinek elektrostatických sil a atraktivní složky van der Waalsových sil). V pravé části obrázkuligand narazil na energetickou bariéru danou repulzivní složkou van der Waalsovy síly (v gra-fické reprezentaci „kontaktem” s atomem receptoru). Pokud bychom nevyužívali dynamiku,nemohli bychom se s haptickým zařízením poskytující silovou zpětnou vazbu ve třech stup-ních volnosti tímto směrem dále pohybovat a museli bychom upravit trajektorii haptickéhozařízení.

Obrázek 3.4: Průchod ligandu do dokovací pozice.

Na obrázku 3.4 vidíme pokračování pohybu ligandu při využití dynamiky a uchyceníligandu k PHANToMu pomocí pružiny a tlumiče. I když pohybujeme haptickým zařízenímstále rovnoběžně s osou x, ligand díky tahu pružiny a získané hybnosti rotuje a zároveň ses natahováním pružiny posouvá mírně výše, čímž prochází energeticky výhodnější trajektorií(na obrázku je trajektorie PHANToMu znázorněná čárkovanou čarou, trajektorie ligandutečkovanou).

3.2.4 Semiimplicitní numerická integrace

Systém diferenciálních rovnic popisující pohyb tělesa je vyjádřen jako

y(t) = f(t, y(t)),

11

kde y je vektor popisující stav simulace (síly i rychlosti) a t je čas simulace.Pro systém rovnic je třeba nalézt numerické řešení. Můžeme využít buďto implicitní, nebo

explicitní numerickou integraci.Explicitní numerická integrace pracuje s rovnicí ve tvaru yn = yn−1 +∆tyn−1. Obecně je

pro mnoho systémů numericky nestabilní, protože chyba daná diskrétním časovým krokemvede ke kumulaci energie systému. To se může např. projevit tak, že se kolidující těleso vlivemnepřesnosti metody rozkmitá a opustí polohu, ve které by se mělo ustálit.Implicitní numerická integrace používá rovnici ve tvaru yn = yn−1 + ∆tyn. Její nevýho-

dou je nutnost vypočítat derivaci v n-tém kroku výpočtu a s tím související nutnost řešenísystému rovnic. Oproti explicitní integraci však vlivem nepřesnosti výpočtu snižuje energiiv dynamickém systému a přináší tak větší stabilitu výpočtu [16].Jako kompromis mezi rychlostí a přesností používáme semiimplicitní Eulerovu metodu,

která aproximuje f pomocí Taylorova rozvoje.Použití semiimplicitní Eulerovy metody na řešení dynamických rovnic v haptické aplikaci

je popsáno v [16]. Naše aplikace přebírá tuto myšlenku, avšak namísto přímo naprogramovanéderivace rovnic určujících pohyb ligandu využívá automatickou derivaci zprostředkovanou kni-hovnou ADOL-C [17]. Díky tomu poskytuje větší volnost v úpravách rovnic určujících pohybligandu – např. silového pole molekul. Numerický kód aplikace lze tedy snadno modifikovatbez nutnosti přepisovat vyjádření derivací užitých funkcí [24].

3.3 Architektura aplikace

Výpočetní část naší aplikace současně ovládá haptické zařízení, u něhož je zapotřebí ovládánív reálném čase, 25× za sekundu odesílá zobrazovací data klientské části, se kterou si zároveňvyměňuje řídící data a počítá interakci molekul. Z toho důvodu je rozdělena na několik vlá-ken, která mají jednotlivé úkoly na starosti. Komunikaci mezi jednotlivými částmi aplikaceznázorňuje obrázek 3.5.

Obrázek 3.5: Architektura aplikace.

Vlákno řízení PHANToMu si vytváří knihovna libphantom a slouží k ovládání fyzickéhozařízení PHANToM. Vlákno obsluha haptiky simuluje pružinu s tlumičem, kterým je virtuálníkurzor reprezentovaný PHANToMem připevněn k ligandu, čte tedy aktuální polohu PHAN-ToMu a počítá síly, které na něj působí. Je velmi důležité udržet běh těchto dvou vlákenna frekvenci 1 000Hz, aby se jevila zpětná vazba v haptice uživateli spojitá.Pozice ligandu aktualizuje vlákno nazvané výpočet simulace. To zpracovává reakci ligandu

na silovou interakci s receptorem a na tažení PHANToMem. Toto vlákno by mělo běžettaktéž na frekvenci 1 000Hz, nicméně tato frekvence zde již není tak kritická. Pokud by se

12

totiž vlákno zastavilo, bude se uživateli ligand jevit jako těleso fixované v prostoru, kterýmnedokáže pomocí haptiky hnout. V silové zpětné vazbě však nebudou nespojitosti – ty jsouzávislé na simulaci pružiny, která v tomto případě jen nebude schopna s ligandem pohnout.Takovéto chování je samozřejmě z hlediska smyslu aplikace nežádoucí, nicméně zpomaleníběhu tohoto vlákna do určité míry není pro simulaci kritické.Vlákno zasílání dat čte polohu a orientaci ligandu a pomocí protokolu UDP ji zasílá spolu

s dalšími údaji k zobrazení zobrazovací aplikaci (na obrázku je její jediné vlákno pojmeno-váno zobrazení a ovládání). Ta vykresluje vzájemnou polohu ligandu a receptoru, energieaktuální polohy ligandu i grafy jejich průběhů. Umožňuje zároveň uživateli ovládat zobra-zení (přibližování či oddalování kamery, rotace scény) i samotnou simulaci (ovládání záznamuči reprodukce [8], nastavení tuhosti pružiny aj.). Tyto údaje jsou zasílány zpět protokolemTCP do výpočetní části aplikace, kde se o jejich zpracování stará vlákno řízení.

3.4 Možné budoucí rozšíření

Současná implementace aplikace pracuje s rigidním receptorem a ligandem. To omezuje jejínasazení na chemické problémy, u nichž se nemění konformace molekul při jejich interakci.V praxi se však často stává, že ligand při průchodu do dokovací pozice změní svou kon-

formaci. Pokud máme pouze původní či novou konformaci, nemusíme být schopni ligandzadokovat (bylo-li by pro úspěšné zadokování třeba, aby ligand při průchodu do aktivníhomísta měnil svou konformaci).Samotná molekula receptoru také není nikterak prostorově rigidní. Krom toho, že může

měnit svou konformaci v reakci na přítomnost ligandu, ji mění také samovolně. Zatímcov prvním případě bychom museli zpracovávat reakce receptoru na ligand online, můžeme vedruhém zmíněném případě využít předem napočítané trajektorie receptoru, tj. jeho konfor-mačního chování v čase.Úkolem této diplomové práce není hledat metody, jak tato vylepšení aplikace realizovat,

nicméně při návrhu algoritmů urychlujících běh stávající aplikace jsem se snažil s těmitozměnami počítat tak, aby byly využitelné v momentě, kdy by byla některá z těchto rozšířeníimplementována.

13

Kapitola 4

Optimalizace zobrazovací části

Zobrazovací část aplikace se stará o vykreslování ligandu a receptoru, informací o stavu simu-lace a zpracovávání či přeposílání vstupu od uživatele. Vzhledem k faktu, že receptor můžeobsahovat tisíce atomů, se může stát jeho vykreslování omezujícím faktorem rychlosti pře-kreslování aplikace. Obrazovku je pro zdání plynulého pohybu molekul třeba překreslit cca25× za sekundu.To, díky čemu je vykreslování receptoru tak náročné, je množství atomů (ty jsou graficky

reprezentovány koulemi) a geometrická složitost každého z nich (každá koule je vytvořena z re-lativně velkého množství trojúhelníků). Způsoby, jak snížit výpočetní náročnost vykreslováníscény s velkým receptorem, se budu zabývat v této kapitole.Protože problémy, se kterými se setkává klientská část naší aplikace, jsou společné pro

mnoho aplikací, lze k jejich řešení využít již hotový software. Při hledání vhodného balíkuknihoven (v tomto kontextu se obvykle používá pojem grafický engine) padla nakonec volbana engine OGRE [10].

4.1 Grafický engine OGRE

Engine OGRE (Object-Oriented Graphics Rendering Engine) jsme zvolili především kvůlinásledujícím vlastnostem:

• implementuje všechny požadované optimalizace i mnoho jiných, které by mohly býtvyužity, pokud by se to ukázalo jako účelné

• je nezávislý jak na operačním systému, tak na používaném grafickém rozhraní (lze tedypoužít jak s OpenGL, tak s DirectX )

• je k němu k dispozici dostatečné množství nástrojů (pro nás je důležitý předevšímexportér objektů z 3ds Max1)

• je velmi dobře a čistě navržen a rozumně dokumentován

• komunita vývojářů i uživatelů je rozsáhlá a aktivní

• je vydáván pod LGPL2 licencí

1produkt firmy Autodesk – www.autodesk.com2http://www.gnu.org/licenses/lgpl.html

14

Veškeré grafické objekty v OGRE jsou ve scéně stromově uspořádány, takže mnoho akcí(např. transformace, uložení do dávky aj.) aplikovaných na jeden vrchol stromu se promítne dovšech jeho potomků. Toho je s výhodou použito při pohybování s ligandem i při vkládání re-ceptoru do statické geometrie (viz 4.2). Práci s enginem usnadňuje také možnost reprezentovatněkterá data skriptem (např. definice materiálů či zobrazování dvourozměrných informací3).

4.2 Dávkování geometrie

Geometrie4 popisující určitý 3D objekt je přesouvána na grafickou kartu v dávkách (batches).Pro tuto dávku musí platit, že je na ni aplikován stejný materiál5. Pro moderní systémy jetypické, že dokáží vykreslovat velmi rychle velké množství trojúhelníků, ale přepínání kontextumezi jednotlivými dávkami je náročné [9]. Navíc toto přepínání zatěžuje primárně CPU, takžeaplikaci, která vykresluje v příliš mnoha dávkách na jeden snímek, nelze významně zrychlitpovýšením GPU.Vzhledem k tomu, že receptor obsahuje velké množství atomů na něž je aplikován stejný

materiál, jeví se spojení jejich geometrie do jedné dávky jako vhodná optimalizace. Grafickýengine OGRE obsahuje třídu StaticGeometry, do níž lze vložit uzel grafu scény a následně me-todou void StaticGeometry::build() vytvořit dávku, která bude na grafickou kartu předávánanajednou (pokud to neuděláme, je použita extra dávka pro každý atom receptoru). Obsahuje-lipodstrom vkládaného uzlu více materiálů, vytvoří metoda void StaticGeometry::build() dávekněkolik.Tato metoda přináší několik omezení:

• dávka obsahuje pro každou instanci objektu kopii geometrie a zabírá tak více paměti

• současné vykreslování dávky znemožňuje ořezávání těch objektů, které nejsou vidět

• jakmile je geometrie jednotlivých objektů zpracována do dávky, nelze ji již transformovatnezávisle na ostatní geometrii v dávce

Z těchto omezení je pro naši aplikaci významné pouze poslední zmíněné. Zatímco pamě-ťové nároky nejsou nijak vysoké a ořezávání nutně nepotřebujeme (ořez geometrie je schopnodělat i GPU, navíc zpravidla vidíme celý receptor), nemožnost transformovat jednotlivé částigeometrie v dávce zamezuje zobrazování flexibilního receptoru. Vzhledem k tomu, že vytvo-ření nové dávky je dosti náročné, nelze dávku znovuvytvářet při každé změně konformacereceptoru. Překonáním tohoto omezení se zabývám v kapitole 4.4

4.3 Dynamická úroveň detailů

Metoda dynamické změny detailů, anglicky Level of Detail, LOD, je obecně technika snižo-vání detailnosti zobrazení objektů typicky v závislosti na jejich vzdálenosti od kamery (lzekombinovat např. s rychlostí objektů, jejich důležitostí aj.). Snižovat se může např. složitost

3to mohou být tlačítka, texty, grafy aj.4pole souřadnic vertexů, tj. prostorových bodů a jejich indexů, přiřazujících vertexy jednotlivým trojúhel-

níkům5soubor pravidel, určujících, jakým způsobem je vykreslován povrch geometrie, obsahuje např. složky od-

raženého světla, textury, aplikované shadery aj.

15

geometrie, kvalita textury, složitost materiálu aj. V našem případě je důležitá právě složitostgeometrie (atomy mají velmi jednoduchý materiál bez textury).LOD dynamicky mění detailnost atomu podle vzdálenosti od kamery, což by mělo být pro

pozorovatele prakticky neznatelné (detaily se ztrácí tam, kde je objekt už příliš malý na to,aby byly patrné). Snižuje se tím zátěž GPU, která renderuje výrazně nižší počet trojúhelníků.OGRE umožňuje buďto automatický LOD, kdy je zjednodušování geometrie v režii samot-

ného enginu, nebo LOD s předem danými objekty. V naší aplikaci je využit druhý zmíněný.

Obrázek 4.1: Receptor beta-cyklodextrin s vypnutým a zaplým LOD.

Srovnání kvality grafického výstupu s LOD a bez něj na beta-cyklodextrinu demonstrujeobrázek 4.1 (bez LOD vykreslujeme 176 322 trojúhelníků, s LOD 11 770). Napravo je receptors LOD, což lze poznat podle mírně hranatých okrajů atomů a podle méně hladkého stínování.

4.4 Využití programovatelné GPU

Jak již bylo popsáno v kapitole 4.2, je výrazně rychlejší, renderujeme-li více atomů se stejnýmmateriálem v jedné dávce. Moderní grafické akcelerátory obsahují programovatelné shaderovéjednotky, které mohou zprostředkovat pohyb geometrie přímo v GPU.Pod označením shader se skrývá program, který zpracovává přímo GPU. V současné době

rozlišujeme vertex shader (program aplikovaný na vrcholy geometrie), pixel shader (programaplikovaný na obrazové body) a geometry shader (program modifikující geometrii přidává-ním či odebíráním vrcholů). Obecně jsou shadery velmi rychlé, jelikož mohou být současněaplikovány na větší množství dat.Máme-li grafickou kartu podporující shader model 2.0 nebo vyšší, můžeme implementovat

tzv. shader instancing. CPU pouze předá grafické kartě souřadnice jednotlivých atomů v dávcea o samotnou transformaci se stará vertex shader.Tato funkcionalita zatím není vzhledem k rigidnímu receptoru v aplikaci implemento-

vána, nicméně její přidání nevyžaduje žádné větší úpravy kódu. Implementuje ji totiž přímoOGRE, a to od verze 1.4. Po vytvoření dávky je tato dávka zkopírována a přidá se k ní poletransformací o velikosti odpovídající počtu objektů v dávce. Hodnoty tohoto pole pak čtevýše zmíněný vertexový program, který již mění polohu vykreslovaných objektů, a nám stačíukládat do něj nové pozice atomů, aniž bychom museli znovuvytvářet dávku.

16

4.5 Shrnutí

S použitím enginu OGRE verze 1.2 jsem implementoval optimalizaci dávek geometrie a dyna-mickou úroveň detailů. Jejich přínos k rychlosti aplikace je popsán v kapitole 6.1. Implemen-tace flexibilního receptoru by při optimalizaci dávek vyžadovala využití programovatelnéhoGPU. Tato funkcionalita není v současné podobě aplikace implementována. K její implemen-taci je zapotřebí upravit aplikaci pro práci s OGRE 1.4.

17

Kapitola 5

Optimalizace výpočetní části

Výpočetní část aplikace má za úkol především počítat vzájemné silové působení ligandu areceptoru (van der Waalsovy a elektrostatické síly a jejich derivace) a integrovat pohybovérovnice. Zatímco sama integrace pohybových rovnic není výpočetně náročná (ligand je rigidnítěleso), výpočty sil a jejich derivací představují náročnou výpočetní úlohu. Zároveň je zapo-třebí přistupovat odlišně k výpočtu van der Waalsových a elektrostatických sil (viz 2.3.3) ajejich derivací.Kvůli připojení haptického zařízení musí být všechny algoritmy navrženy tak, aby byla

iterace výpočtu silového působení molekul dlouhá nejvýše 1 ms. Není tedy dost dobře možnévyužít některé algoritmy s dobrou amortizovanou složitostí, které někdy vyžadují složitějšíreorganizace datových struktur.

5.1 Výpočet van der Waalsových interakcí

Van der Waalsovy síly mají v blízkosti jader velmi prudký průběh, ale s přibývající vzdálenostívelmi rychle konvergují k nule (viz 2.3.3). Díky tomu si můžeme dovolit počítat pro každýatom ligandu pouze interakce s blízkými atomy receptoru [12]. Tato metoda tedy vnáší dovýpočtu nepřesnost danou ořezem vzdálených atomů.Protože i blízkých atomů receptoru může být stále příliš mnoho, je zapotřebí hledat další

možnosti urychlení výpočtu. Tou je paralelizace algoritmu počítajícího van der Waalsovy sílya jejich derivace.

5.1.1 Nalezení blízkých atomů

Přičteme-li k van der Waalsovu poloměru uživatelem zadaný ořez c, můžeme k výběru blízkýchatomů využít některý z algoritmů detekce kolizí. Klasická detekce kolizí se skládá ze dvou fází– v první vybereme velmi rychlým algoritmem potenciálně kolidující množiny objektů, vedruhé fázi počítáme precizním algoritmem detekce kolizí, zdali ke kolizi opravdu došlo. V tétokapitole probereme algoritmy první fáze detekce kolizí (druhá fáze je v našem případě triviální– pouze ověříme eukleidovskou vzdálenost atomů). Po provedení detekce kolizí počítáme prokaždý atom ligandu interakce pouze s těmi atomy receptoru, se kterými těleso tvořené atomemligandu zvětšeným o c koliduje.Počítáme-li s budoucím využitím flexibilního receptoru, je třeba zvolit takový algoritmus

detekce kolizí, který používá datovou strukturu schopnou rychlé reorganizace v případě změny

18

konformace receptoru.Naivní způsob nalezení blízkých atomů by znamenal O(r · l) výpočtů eukleidovské vzdá-

lenosti mezi l atomy ligandu a r atomy receptoru. To může být při řádově desítkách atomůligandu a tisících atomů receptoru poměrně náročná operace. Je tedy zapotřebí nalézt řešení,které nám umožní hledat blízké atomy v lepším čase a zároveň neznemožní budoucí zavedeníflexibilního receptoru.

Hierarchie obálek

Obálkou rozumíme geometrický objekt, který ohraničuje těleso, se kterým chceme počítatkolize. V první fázi pak počítáme kolize s obálkami. Kolize s tělesy, které jsou uvnitř obálek(detekce kolizí s nimi je obecně náročnější), počítáme až v případě zjištěné kolize s obálkou.Výpočet kolize s objektem tvořícím obálku musí být poměrně rychlý, proto jsou jako obálkyvoleny zpravidla jednoduchá tělesa.Při větším počtu těles je vhodné obálky hierarchicky uspořádat do stromové struktury tak,

že každý vrchol stromu obsahuje obálku, která v R3 ohraničuje obálky všech svých potomků.V listech stromu jsou obálky, které již ohraničují přímo objekty, se kterými chceme počítatkolize.V našem případě můžeme triviálně použít kulové obálky – pak je obálkou celý objekt,

jelikož atomy jsou reprezentovány jako koule (čímž úplně odpadá nutnost provézt druhou fázidetekce kolizí). Je výhodné uspořádat do hierarchie zvlášť atomy receptoru a počítat kolizebuďto mezi každým atomem ligandu s hierarchií obálek atomů receptoru, nebo vytvořit proatomy ligandu druhou hierarchii a počítat kolize těchto dvou hierarchií.

Obrázek 5.1: Hierarchie kulových obálek.

Stromovou hierarchii kulových obálek můžeme budovat dvěma způsoby – shora a zdola.Budování obálky shora je popsáno v [11], pro budováni zdola můžeme např. postupně sesku-povat nejbližší obálky do dvojic, až dokud nebude celý receptor uzavřen do jedné velké koule.Obecně zřejmě nelze říci, který přístup je pro naši aplikaci vhodnější.Výhodou tohoto prostorového uspořádání je velká rychlost, především v případě, že jsou

si ligand a receptor vzdáleny (v takovém případě je složitost hledání kolizí v O(l) pro detekcikolizí mezi každým z l atomů ligandu a obálkovou hierarchií receptoru, nebo O(1) pro detekci

19

kolizí mezi hierarchií obálek ligandu a hierarchií obálek receptoru, jelikož dokážeme vyloučitkolizi hned při kontrole kořenové obálky). Nás však zajímá především rychlost detekce kolizív případě, kdy jsou mezi receptorem a ligandem nějaké kolize detekovány (tzn. molekuly spoluinteragují). Zde je složitost v O(l · log r) pro oba zmíněné přístupy.Další výhodou hierarchií obálek je jejich prostorová složitost – ta je lineární vzhledem

k počtu objektů, které organizuje.Podstatnou nevýhodou je relativní náročnost změn ve stromové hierarchii v momentě, kdy

dojde ke změně polohy objektů, které hierarchie organizuje. Jestliže bychom měli flexibilníreceptor, museli bychom při změně jeho konformace znovu budovat hierarchii obálek, což ječasově velmi náročné.

Hierarchické dělení prostoru

Rozdílný přístup k budování hierarchií přináší hierarchie dělení prostoru. Na rozdíl od hierar-chií obálek jsou hierarchie dělení prostoru založeny na rekurzivním dělení prostoru s ohledemna objekty v něm umístěné. Určíme tedy oblast v prostoru, kterou rekurzivně dělíme a zís-káváme tak podoblasti obsahující objekty v prostoru či jejich části. Další dělení prostoruovlivňuje množství objektů, které se nachází v aktuálně dělené buňce – buňky, které neob-sahují žádné objekty, není třeba dále dělit, zatímco buňky obsahující stále velké množstvíobjektů dělíme dále.Asi nejznámějším algoritmem pro hierarchické dělení prostoru je oktalový strom (octree).

Ten pracuje s krychlemi, které dělí na poloviny v každé dimenzi – každý vrchol stromu od-povídá krychli v prostoru a může mít až 8 potomků. Každá krychle, kterou definuje vrchol,který je listem, obsahuje odkazy na objekty v ní obsažené.

Obrázek 5.2: 2D zobrazení oktalového stromu.

Domnívám se, že rychlost detekce kolizí s využitím hierarchií dělení prostoru není v našempřípadě tak vysoká, jako u hierarchií obálek (vzniklý strom zpravidla není vyvážený, nerozdělíprecizně všechny atomy, je zapotřebí spojovat seznamy atomů v jednotlivých buňkách s vylou-čením duplicit, musíme provést druhou fázi detekce kolizí). Je zde však jednodušší vyrovnáníse se změnou konformace receptoru. Pokud dojde k přesunu některých atomů, můžeme jepouze přesunout v již existující prostorové hierarchii a její reorganizaci odložit, či ji provádět

20

postupně (není-li hierarchicky dělený prostor přeuspořádán ihned, dochází pouze k tomu, ženení detekce kolizí tak efektivní, protože se může v jednom vrcholu stromu vyskytovat přílišmnoho objektů a naopak některé vrcholy nemusí obsahovat žádné objekty a zbytečně se dáledělit).

Ploché dělení prostoru

Nejjednodušším řešením detekcí kolizí je dělení prostoru bez stromové hierarchie. Tato metodadělí prostor zpravidla pravidelnou prostorovou mřížkou. Její výhodou je, že nalezení buňkyči buněk, ve které se nachází atom ligandu, lze provézt v konstantním čase. Tyto buňky pakobsahují odkazy na atomy receptoru, které se v nich nachází. Další výhodou plochého děleníprostoru je snadná reorganizace po změně polohy atomů receptoru – stačí pouze vyjmoutatomy receptoru, které změnily pozici z buněk, ve kterých byly dříve obsaženy, a vložit je dobuněk nových.Je-li příliš mnoho objektů blízko sebe, nedokáže je mřížka účinně oddělit a navíc se každý

objekt vyskytuje v několika buňkách, takže je zapotřebí sjednocovat seznamy atomů receptrouobsažených v buňkách, do kterých zasahuje atom ligandu tak, abychom nehlásili kolizi s ně-kterým z atomů receptoru vícekrát.Nevýhodou plochého dělení prostoru je také velká paměťová náročnost, která je v O(n3)

vzhledem k počtu buněk v jedné dimenzi (ty musí pokrýt celý prostor, ve kterém se můžemepohybovat). V takovémto uspořádání je navíc zpravidla většina buněk prázdných. To lze doznačné míry redukovat hashováním indexů v regulární prostorové mřížce do tabulky nižšíchrozměrů.

Detekce kolizí pomocí hashování prostoru

Při hashování prostoru dělíme prostor regulární mřížkou, ale neukládáme její buňky. Namístotoho je pomocí hashovací funkce mapuje na indexy konečné hashovací tabulky. Objekty setak vkládají na řádky hashovací tabulky odpovídající indexům buněk regulární mřížky, dokterých spadají. Obrázek 5.3 ukazuje zobrazení objektů v prostoru děleném regulární mřížkoudo tabulky. Hashovací funkce zde má (pro lepší ilustrativnost příkladu) tvar h = ⌊x

3⌋ ⊕ ⌊y

3⌋,

kde ⊕ je operace spojení řetězců.

Obrázek 5.3: Indexace prostoru.

21

Díky hashování se vyhneme problému s přílišnou prostorovou náročností, kterou trpí re-gulární mřížka, protože velikost hashovací tabulky si můžeme nastavit libovolně bez omezeníprostoru, v němž se mohou vyskytovat atomy receptoru. U hashování prostoru také můžemebalancovat rychlost detekce kolizí vůči rychlosti přesunu atomů receptoru, a to následujícímizpůsoby:

• Chceme-li co nejrychlejší úpravu hashovací tabulky při přesunu atomů receptoru, mů-žeme atomy receptoru zkolabovat do bodů (bereme tedy jejich poloměr jako nulový).Tím je zajištěno, že každý atom receptoru padne právě do jednoho řádku hashovacítabulky a pokud mění svou polohu v regulární mřížce, smažeme pouze jeden odkazz hashovací tabulky a jeden přidáme. Při detekci kolize s atomem ligandu musíme po-loměr atomu ligandu navýšit kromě hodnoty ořezu ještě o nejvyšší poloměr z atomůreceptoru. Tím zpomalíme detekci kolizí, jelikož atom ligandu zvětšený o tento poloměrzasáhne více buněk.

• Chceme-li dosáhnout co nejrychlejší detekce kolizí, zvětšíme poloměry atomů receptorujak o uživatelem zadaný ořez, tak o maximální poloměr atomů ligandu. Tím každý atomreceptoru padne do mnoha buněk, ale atom ligandu, který nyní můžeme zkolabovat dobodu, padne vždy právě do jedné buňky regulární mřížky, takže stačí vybrat pouze je-den řádek hashovací tabulky a atomy receptoru v něm obsažené vrátit jako potenciálněkolidující. Je také nutné zavést dostatečně malé buňky regulární mřížky, abychom ne-vybírali i příliš vzdálené atomy receptoru (tím bychom ztratili výpočetní čas na druhéfázi detekce kolizí, kde bychom počítali zbytečně mnoho eukleidovských vzdáleností).

Mezi těmito mezními řešeními existuje samozřejmě spousta kompromisů. Můžeme např.ukládat atomy receptoru s jejich poloměry a kolidovat s atomy ligandu zvětšenými o hodnotuořezu, ukládat atomy receptoru zvětšené o hodnotu ořezu, nebo tuto hodnotu rozdělit meziatomy ligandu a receptoru. Všechna tato řešení nekladou žádné nároky na úpravu algoritmu,vše určujeme pouze vstupními daty.Mezní řešení, kdy zkolabujeme atomy receptoru či ligandu, lze s úspěchem použít, jelikož

mají všechny atomy relativně podobné poloměry. Pokud bychom pracovali s prostorovýmiobjekty, které by se výrazněji lišily svou velikostí, jsou pro nás tato řešení nevhodná, protožebychom při přičítání maximálního poloměru atomů ligandu či receptoru často zasahovali dovíce buněk regulární mřížky, než je nezbytně nutné.Další v algoritmu volitelnou vlastností je reprezentace seznamů objektů uložených v jed-

notlivých buňkách hashovací tabulky. Pokud položky v jednotlivých řádcích nijak neorgani-zujeme, můžeme vkládat nové prvky v čase O(1), vyhledávat v čase O(n) (to potřebujeme připřesunu atomu receptoru) a sjednotit řádky o n a m prvcích v čase O(n+m) (pomocí tabulky– bude vysvětleno níže) při detekci kolizí s atomem ligandu spadajícím do těchto dvou řádků.Seřadíme-li (např. podle ID atomu) však objekty v každém řádku do vyváženého binárníhovyhledávacího stromu (např. AVL), pak přidáme či smažeme prvek v řádku délky n v časeO(log n) a sjednocení objektů z řádků délky n a m provedeme v čase O(n+m).Abychom se vyhnuli zbytečně vysokému nárůstu složitosti v závislosti na počtu řádků

hashovací tabulky, do kterých zasáhne ligand, je vhodné sjednocovat seznamy atomů vždyza použití tabulky, do které ukládáme pro každý atom informaci o tom, zda jsme jej již dosjednocení vložili. Poté, co sjednocování dokončíme, můžeme tabulku v lineárním čase vyčistit(víme, které prvky jsme vkládali do sjednocení a u kterých je tedy tato informace). Celkově jetedy časová složitost detekce kolizí mezi receptorem a ligandem o l atomech pomocí hashování

22

prostoru v O(l · c · k), kde k je nejvyšší počet atomů v jednom řádku hashovací tabulky a c

je počet řádků, do kterých zasahuje jeden atom ligandu.Přesun r atomů receptoru je v O(r · log k) při seřazených atomech v řádcích hashovací

tabulky, respektive O(r · k) při neseřazených atomech.

Srovnání

Pro detekci blízkých atomů u rigidního receptoru jsou rozumně využitelné všechny zmíněnépřístupy. Zřejmě nejvýhodnější by byla hierarchie obálek, která umožňuje velmi rychlou de-tekci kolizí a zároveň je prostorově nenáročná. Prostorová hierarchie i hashování prostoru jsouschopny rozdělit atomy receptoru obdobným způsobem (výsledek by mohl být lepší u nepra-videlných prostorových hierarchií, např. u kd-stromů), regulární prostorová mřížka přinášízbytečné plýtvání prostorem.Situace je však dramaticky odlišná v případě, kdy uvažujeme o flexibilním receptoru.

Zde se veškeré hierarchické uspořádání jeví jako nevhodné, jelikož často musíme při změněpolohy atomů receptoru vybudovat alespoň část celé hierarchie znovu. To vede k občasnémuvýraznému zpomalení výpočtu nevhodném pro haptickou aplikaci.Z tohoto důvodu jsem se rozhodl implementovat detekci kolizí pomocí výše popsaného

hashování prostoru.

5.1.2 Paralelní výpočet

I po zanedbání vzdálených atomů je výpočet van der Waalsových sil a jejich derivací stálenáročný. Naštěstí se jedná o problém, který je možné urychlit paralelizací. Máme-li k dispozicivíceprocesorový stroj nebo cluster, rozdělíme množinu atomů receptoru, se kterými počítámeinterakce, mezi jednotlivé CPU. Na každém CPU spočítáme sílu působící mezi ligandem apříslušnou podmnožinou atomů receptoru. Výsledná van der Waalsova síla, moment a jejichderivace působící na ligand je pak součtem sil, momentů a jejich derivací spočítaných najednotlivých CPU.Některé části výpočtu prakticky nelze paralelizovat (integrace dynamických rovnic), jiné

lze, ale jsou natolik nenáročné, že není paralelizace třeba (výpočty elektrostatických sil). Zvoliljsem pro paralelizaci architekturu master-slave – master proces provádí výpočty dynamiky,elektrostatiky i van der Waalsových sil, komunikuje se zobrazovací částí aplikace a s haptickýmzařízením, slave procesy pouze počítají van der Waalsovy síly.Otázkou je, jakým způsobem zajistit rozložení výpočtu mezi jednotlivé procesy. Jednou

možností je v master procesu provést detekci kolizí a množinu blízkých atomů rozdělit conejvíce rovnoměrně mezi výpočetní procesy. Nevýhodou tohoto přístupu je nutnost provéstdetekci kolizí v master procesu (detekce kolizí je i při využití algoritmů zmíněných v kapitole5.1.1 stále relativně náročné, zvláště dojde-li ke změně konformace receptoru) a také zvýšenázátěž na komunikaci mezi CPU, protože musíme přenést údaje o všech atomech, se kterýmibudeme počítat interakce.Kvůli těmto nevýhodám jsem zvolil implementaci odlišného přístupu, kdy je paralelní již

detekce kolizí. Celý receptor rozdělíme na n částí, kde jednotlivé atomy přiřazujeme výpo-četním procesům podle hodnoty jejich pořadí v receptoru modulo n. Máme-li tedy r atomůreceptoru, každý výpočetní proces počítá s ⌊ r

n⌋ až ⌈ r

n⌉ atomy, u nichž nejprve provede detekci

kolizí a pro nalezené blízké atomy spočítá van der Waalsovy síly. Tento přístup je výhodnýv tom, že je distribuována i detekce kolizí (její urychlení popisuje kapitola 6.2) a mezi výpočet-

23

ními procesy je jednodušší komunikace (master proces pouze broadcastuje informaci o polozea orientaci ligandu, slave procesy mu vrací vypočítanou sílu, moment a jejich derivace), cožje výhodné u výpočtu na clusterech (viz 6.3).Nevýhodou zvoleného řešení je, že automaticky nezaručuje rovnoměrnou distribuci vý-

počtu van der Waalsových interakcí mezi jednotlivé procesy. Stačí si představit situaci, kdymáme atomy receptoru uspořádány tak, že jednotlivým procesům přidělíme takové atomy,které jsou blízko sebe. Při přiblížení ligandu k takové skupině atomů se může stát, že budejeden výpočetní proces neúměrně zatížen (síly s téměř všemi blízkými atomy bude počítaton) a ostatní procesy nebudou mít co počítat. Ačkoliv v praxi k takovýmto jevům nedochází,protože jsou atomy receptoru indexovány tak, že blízké atomy dostávají po sobě jdoucí in-dexy a rozdělení pomocí modula je tedy přiřazuje různým procesům (to dokazují výsledkyměření v kapitole 6.3), lze řešit i situaci, kdy jsou atomy receptoru indexovány „zlomyslně”.Zde předkládám návrh k řešení takovéto situace, v současné implementaci je využito pouzerozdělení atomů mezi procesy pomocí modula.

Nevhodně indexovaná data receptoru

V případě nevhodné indexace je možno buďto indexy randomizovat (tzn. prohazovat indexynáhodně vybraných dvojic atomů), nebo formulovat algoritmus, který rozdělí atomy recep-toru do množin o velikosti lišící se maximálně o 1 tak, aby v každé množině byly co možnánejvzdálenější atomy. Odlišným řešením je rozdělit mezi výpočetní procesy namísto atomůreceptoru buňky prostorové mřížky. Tím dokážeme zajistit poměrně rovnoměrné rozděleníprostoru okolo receptoru mezi jednotlivé výpočetní procesy a pro každý proces počítáme vander Waalsovy síly s těmi atomy receptoru, které se nachází v jím vlastněných částech prostoru.V tomto případě je nevhodné, aby mohl jeden atom receptoru padnout do více buněk, jelikožby si musely jednotlivé výpočetní uzly předat informaci o tom, kdo pro něj bude počítat vander Waalsovy síly. Je proto dobré využít kolaps atomů receptoru do jednoho bodu (viz. popishashovacího algoritmu v kapitole 5.1.1).

5.2 Výpočet elektrostatických sil

Elektrostatické síly působí nezanedbatelně i na větší vzdálenosti, takže není možné počítat prokaždý atom ligandu interakce pouze s blízkými atomy receptoru (viz 2.3.3). Oproti van derWaalsovým silám však nemají tak prudký průběh v blízkosti jader (respektive byl by velmiprudký až v poloze, do které se nemohou jádra díky repulzivní složce van der Waalsových sildostat), takže si můžeme dovolit předpočítat je do mřížky, aniž by tato mřížka musela býtextrémně jemná (což je náročné jak výpočetně, tak také paměťově). Mezi předpočítanýmibody pak interpolujeme.Jak již bylo zmíněno výše, elektrostatické síly jsme schopni za cenu akceptovatelné ztráty

přesnosti zrychlit předpočítáním vektorů sil na mřížce okolo receptoru. V každém bodu mřížkyje uložen vektor síly, kterou působí všechny atomy receptoru na jednotkový náboj. Jehohodnoty získáme následujícím vztahem:

−→F norm =

14πǫ0ǫr

·r

∑

i=0

qri

|−→P −

−→R i|3

· (−→P −

−→R i),

24

kde qriznačí náboj i-tého atomu receptoru o r atomech,

−→R i polohu i-tého atomu receptoru,−→

P je poloha bodu mřížky, ǫ0 je permitivita vakua a ǫr relativní permitivita prostředí.Silovou interakci mezi atomem ligandu umístěném přímo v bodě mřížky a receptorem pak

získáme pomocí následujícího vztahu:

−→F =

−→F norm · cl.

Pro atom ligandu v obecné poloze dopočítáváme sílu lineární interpolací hodnot okolníchbodů mřížky. Tato interpolace musí být spojitá při přechodu mezi buňkami (nespojitosti bypřinášely problémy při integraci dynamických rovnic a uživatel by přes haptické zařízení cítilnepříjemné rázy). Více bude interpolace probrána v následujících kapitolách.Za předpokladu, že bychom měli v budoucnu flexibilní receptor, můžeme mít buďto před-

počítané síly v mřížce pro různé konformace receptoru, nebo paralelně počítat vzájemnépůsobení jednotlivých atomů ligandu s atomy receptoru, a to úplně bez, nebo s ořezem nahodně velkou vzdálenost (tedy analogicky výpočtu van der Waalsových sil).

5.2.1 Regulární prostorová mřížka

Zřejmě nejjednodušší je rozdělení prostoru okolo receptoru do buněk tvaru krychlí, v jejichžrozích jsou napočítány síly. Toto řešení je běžně používané v dokovacím software, jako je např.AutoDock [23].

Obrázek 5.4: 2D pohled na regulární prostorovou mřížku.

Obrázek 5.4 znázorňuje část prostoru vyplněného buňkami, ve kterém lze pohybovat ligan-dem (místa, kde jsou napočítány síly, jsou zde označena puntíky). Pro zjednodušení uvažujmetento prostor jako krychli o straně a. Je-li rozlišení mřížky (tzn. strana buňky-krychle, v je-jíž rozích počítáme síly) d, pak je potřeba (a

d)3 buněk. Časová složitost napočítání sil na

prostorové mřížce je pak O(r · (ad)3), kde r je počet atomů receptoru.

Interpolace

Jádro každého atomu ligandu je v každém časovém okamžiku právě v jedné buňce (je-liatom geometricky na rozmezí buněk, můžeme jej deterministicky zařazovat do jedné z buněk,např. do buňky s nižší hodnotou souřadnice dimenze, v níž buňky sousedí). V následujícímčasovém okamžiku může přejít do jiné buňky, a to přes stěnu, hranu a nebo roh buněk. Jetedy zapotřebí implementovat takovou interpolaci, která bude spojitá jak uvnitř buňky, tak

25

při zmíněných přechodech mezi buňkami. V případě regulární mřížky je to jednoduché, jakje vidět na následujícím vztahu.Vektor

−−−−−→cellPos obsahuje souřadnice rohu buňky, vektor

−−−−−−→cellScale obsahuje obrácenou hod-

notu velikosti buňky1. Relativní polohu atomu ligandu uvnitř buňky normovanou na velikostbuňky rovnou jedné získáme vztahem:

relPosi = (ligandPosi − cellPosi) · cellScalei; kde i ∈ {0, 1, 2}.

Pole vah jednotlivých sil v rozích buňky w je vyjádřeno:

w[x][y][z] = (x · relPos0 + (1− x) · (1− relPos0))·(y · relPos1 + (1− y) · (1− relPos1))·(z · relPos2 + (1− z) · (1− relPos2))kde x, y, z ∈ {0, 1}.

Následný výpočet interpolované síly provedeme dle vztahu:

−→F = ligandCharge ·

1∑

x=0

1∑

y=0

1∑

z=0

w[x][y][z] ·−−−−−−−−−−−−−→cellForces[x][y][z].

Lemma 5.2.1.1. Při spojitém pohybu atomu uvnitř buňky je interpolace síly−→F spojitá.

Důkaz Je-li spojitý pohyb atomu ligandu, mění se hodnoty vektoru ligandPos spojitě. Vý-počet vektoru

−−−−→relPos je také spojitý (

−−−−−→cellPos i

−−−−−−→cellScale jsou pro konkrétní buňku konstanty).

Složky vektoru vah w jsou součinem hodnot vektoru−−−−→relPos či hodnot rozdílu 1−relPosi.

Jelikož dává výpočet hodnot vektoru−−−−→relPos spojité hodnoty, je jejich součin také spojitý.

Konečně výpočet složek vektoru−→F je spojitý, jelikož se jedná o součet součinů spojitých

hodnot w a pro buňku konstantního pole vektorů−−−−−−−−−−−−−→cellForces[x][y][z] násobený konstantní

hodnotou ligandCharge.2

Lemma 5.2.1.2. V okamžiku přechodu atomu z jedné buňky do druhé přes stěnu je hodnotainterpolované síly

−→F shodná pro výpočet nad oběma buňkama.

Důkaz. Vzhledem k tomu, že jsou všechny buňky v mřížce stejně velké a zarovnané, majíjejich sousední stěny shodné pozice i hodnoty předpočítaných sil. Fixujme časový okamžik,kdy se atom nachází na rozmezí buňky C1 a buňky C2 (je tedy geometricky ve stěně C1 asoučasně v opačné stěně sousední C2). K tomu, aby byl výpočet síly

−→F shodný pro obě buňky,

potřebujeme dokázat, že váhy všech sil netvořících tyto stěny mají nulovou hodnotu a váhypřekrývajících se předpočítaných sil buněk C1 a C2 mají shodnou hodnotu pro obě buňky.Hodnoty vektoru

−−−−→relPos padnou do intervalu < 0, 1 > (pro buňky, ve kterých se na-

chází atom ligandu, nám výsledek rozdílu−−−−−−−→ligandPos−

−−−−−→cellPos udává relativní pozici ligandu

v buňce, součinem s−−−−−−→cellScale normujeme tuto hodnotu tak, aby její složky padly do< 0, 1 >).

Je-li atom ve stěně buňky, pak ∃d ∈ {0, 1, 2}; relPosd = 0 ∨ relPosd = 1. Snadno lze nahléd-nout, že má-li jedna složka vektoru

−−−−→relPos hodnotu 0 či 1, dává výpočet hodnot pole w nulu

1U regulární mřížky jsou všechny tři hodnoty tohoto vektoru shodné, a to pro všechny buňky mřížky.

26

právě pro body reprezentující protilehlou stěnu v dimenzi d (výraz d · relPosd respektive(1− d) · (1− relPosd) nabude nulové hodnoty).Zbývá dokázat, že hodnoty nenulových vah jsou shodné pro síly se shodnou pozicí v C1 i

C2. Vektor−−−−→relPos se liší pouze v dimenzi d (buňky jsou zarovnané, jejich poloha se tedy liší

pouze posunutím v dané dimenzi a poloha ligandu je pro výpočet v obou buňkách shodná).Jsou-li tedy ostatní položky

−−−−→relPos shodné, vypočtou se pro sousední stěny v C1 a C2 shodně i

hodnoty pole vah w (výpočet se liší pouze v hodnotách výrazu d·relPosd+(1−d)·(1−relPosd),který je roven 1 právě pro body na sousedících stěnách). Předpočítané síly v cellForces majíve shodných bodech shodnou velikost, takže i výpočet celkové sily

−→F je pro obě buňky shodný.

2

Lemma 5.2.1.3. V okamžiku přechodu atomu z jedné buňky do druhé přes hranu je hod-nota interpolované síly

−→F shodná pro výpočet nad oběma buňkami.

Důkaz. Lemma lze dokázat obdobně, jako 5.2.1.2, proto formuluji tento důkaz stručněji.Je třeba dokázat, že při přechodu mezi buňkami sousedícími pouze hranou (přecházíme-liprůchodem přes hranu mezi buňkami sousedícími stěnou, lze to považovat za speciální případprůchodu přes stěnu) jsou u bodů tvořících společnou hranu shodné váhy a u ostatních bodůjsou váhy nulové. Je-li atom v hraně buňky, pak ∃d1, d2 ∈ {0, 1, 2}; d1 6= d2 ∧ (relPosd1 =0 ∨ relPosd1 = 1) ∧ (relPosd2 = 0 ∨ relPosd2 = 1). Výpočet pole w je tedy nenulový právěpro indexy v dimenzích d1 a d2 rovnající se relPosd1 a relPosd2. Výpočet nenulových vahpřiřazuje silám tvořící sousední hranu váhu danou polohou ligandu na této hraně relPosd3 ,kde d3 6= d1 ∧ d3 6= d1. Ta je ovšem pro obě buňky shodná. Výpočet celkové síly je tedy proobě buňky shodný.

2

Lemma 5.2.1.4. V okamžiku přechodu atomu z jedné buňky do druhé přes roh je hodnotainterpolované síly

−→F shodná pro výpočet nad oběma buňkami.

Důkaz. Lemma lze taktéž dokázat analogicky dvěma předchozím lemmatům. Je-li atomv rohu buňky, pak ∀i ∈ {0, 1, 2} platí relPosi = 0∨ relPosi = 1. Pak ale výpočet hodnot polew je nenulový právě pro roh s indexy odpovídající

−−−−→relPos, kde má zároveň pro obě buňky

shodnou hodnotu.2

Tvrzení 5.2.1.1. Hodnota interpolované síly−→F je při spojitém pohybu atomu spojitá.

Důkaz. Atom se může pohybovat buďto v rámci jedné buňky, nebo přecházet z jedné buňkydo druhé. Jelikož jsou všechny buňky mřížky stejně velké a zarovnané, může atom přecházetmezi dvěma buňkami pouze přes jejich stěny, přes jejich hrany, nebo přes jejich rohy, majícístejnou pozici v prostoru. Je-li tedy při spojitém pohybu atomu v jedné buňce vypočtenáhodnota

−→F spojitá, jak dokazuje Lemma 5.2.1.1 a je-li v momentě, kdy atom přechází z jedné

buňky do druhé totožný výpočet−→F pro obě buňky, jak ukazují Lemma 5.2.1.2 - 5.2.1.4, je

spojitý výpočet síly−→F na celé mřížce.

2

27

Spojitost interpolace na mřížce nám tedy umožňuje využít ji v haptické aplikaci. Bezdůkazu již ponechme triviální tvrzení, že výpočet elektrostatické síly je pro každý atom liganduv O(1), jelikož interpolujeme sílu vždy z jedné buňky, která obsahuje 8 silových vektorů.Složitost výpočtu pro molekulu ligandu je pak O(l), kde l je počet atomů ligandu – zbavilijsme se tedy úplně závislosti na velikosti receptoru.

Problémy regulární mřížky

Po implementaci regulární prostorové mřížky jsme narazili na problém balancování její přes-nosti. Při příliš velkých buňkách byly výpočty v blízkosti jader atomů dosti nepřesné (elek-trostatická síla zde rychle roste) a pokud se jeden z bodů mřížky příliš přiblížil jádru atomureceptoru, blížila se hodnota síly na něm spočítané nekonečnu. Jakmile se atom ligandu dostaldo buňky s bodem o hodnotě síly blížící se nekonečnu, působila na něj díky lineární interpo-laci ihned neúměrně velká síla, což se projevilo silným rázem v haptickém zařízení a často ipádem simulace.

Obrázek 5.5: Patologický případ interpolace na mřížce.

Na obrázku 5.5 znázorňuje hustě tečkovaná čára velikost síly−→F v dimenzi d v závislosti na

poloze vůči jádru atomu receptoru (to je na ose x reprezentováno puntíkem). Řídce tečkovanéčáry znázorňují hranice buněk v této dimenzi. Interpolace na mřížce je znázorněna pomocíplné čáry. Jasně vidíme, že v dostatečné vzdálenosti od atomu receptoru je interpolace velmipřesná. S tím, jak se k němu přibližujeme, přesnost klesá, až u poslední buňky kvůli blízkostibodu mřížky s jádrem atomu receptoru jde funkce udávající elektrostatickou sílu k nekonečnu.To by nám při precizním výpočtu nevadilo, jelikož se ve skutečnosti jádra atomu nikdy taktonepřiblíží2. Pokud ovšem nejsou buňky dostatečně malé (tzn. funkci interpolujeme po přílišvelkých intervalech), můžeme se do buňky se singulárním bodem dostat. Protože je inter-akce počítána po diskrétních časových krocích, dojde vstupem do této buňky k neúměrnémunárůstu hodnoty interpolované síly, která je vysoko nad precizně počítanou hodnotou.V případě regulární prostorové mřížky lze řešit situaci těmito způsoby (které je možno

spolu kombinovat):

2zabrání jim v tom van der Waalsova síla

28

• Umělým ořezem sil blížících se singularitě (buďto nastavíme pro každý bod mřížkynejvyšší sílu, která na něm může být, nebo nastavíme práh minimální vzdálenosti jádraatomu receptoru od bodu mřížky a u všech bližších atomů počítáme jako vzdálenosttento práh). Toto řešení odstraní problémy s rázy v haptickém zařízení a pády simulace,přetrvává ale nepřesnost výpočtu, pokud se přiblíží atom ligandu atomu receptoru, ataké je sporné, jaké prahy má mít aplikace nastaveny. Přesto je toto řešení použito např.v programu AutoDock 3.0 [23].

• Zmenšením buněk mřížky. Toto řešení odstraňuje problémy s nepřesností i nestabilitou,avšak je omezeno možnostmi hardware (mřížka se nemusí vlézt do paměti, nebo můževýpočet sil v ní uložených trvat neúměrně dlouho).

5.2.2 Adaptivní prostorová mřížka

Při analýze problémů regulární mřížky jsem si uvědomil, že ukládá data velmi nehospodárně.Zatímco v blízkosti jader atomů receptoru má předpočítávaná elektrostatická síla velký gra-dient a pro dostatečně přesný výpočet potřebujeme velmi jemnou mřížku, ve větší vzdálenostije naopak gradient velmi nízký a pro obdobnou přesnost nemusí být mřížka tak jemná.

Vygenerování adaptivní mřížky

Inspirován octree algoritmem pro dělení prostoru (využíván zpravidla pro detekci kolizí, viz5.1.1) jsem implementoval algoritmus adaptivní prostorové mřížky pro ukládání předpočíta-ných hodnot elektrostatických sil (takovéto algoritmy se používají např. pro numerické řešenírovnic, jejich nasazení na předvýpočet silového pole vyžadující spojitou interpolaci jsem ne-našel). Dělení prostoru je zde poměrně jednoduché – analogicky octree algoritmu rekurzivnědělíme prostor do krychlí tak, že receptor a jeho okolí uzavřeme do jedné krychle, pro jejížrohy napočítáme silové vektory. Je-li třeba mřížku zjemnit (kritérium bude popsáno níže),rozdělíme krychli rovnoměrně na 8 krychlí o poloviční délce strany a rekurzivně pokraču-jeme v algoritmu pro každou z těchto menších krychlí (obecně by u octree algoritmu nebylozapotřebí buďto dělit buňku na všech 8 částí nebo dělení zastavit, ale mohli bychom dělitjen některé části buňky – zde nám ale takovéto řešení významně zjednodušuje situaci přiimplementaci interpolace).Algoritmus musí vytvářet jemnější mřížku (tzn. sestupovat v rekurzi hlouběji) úměrně ve-

likosti gradientu předpočítávané elektrostatické síly. Abych se vyhnul analýze funkce udávajícíelektrostatickou sílu (bylo by zapotřebí hledat maximum první derivace funkce na omezenémprostoru), zavedl jsem jednodušší kritérium pracující se vzdáleností nejbližšího jádra atomureceptoru. Vzhledem k tomu, že funkce udávající elektrostatickou sílu klesá s druhou moc-ninou vzdálenosti od jádra, dokážeme přesně odhadnout na základě vzdálenosti její průběhpro jednotlivé atomy receptoru. Kritérium nedokáže postihnout případné interference silovéhopůsobení jednotlivých atomů, které však mohou významným způsobem měnit hodnotu sílypouze směrem k nule (velikost síly klesá s druhou mocninou vzdálenosti, proto vzdálenějšíatomy nejsou schopny sílu významně zvýšit), což nevede k neočekávané nepřesnosti mřížky,pouze k zbytečnému dělení v místech, kde je síla vlivem interference velmi malá.Na vstupu algoritmu je zadána maximální hloubka rekurze maxDepth, násobek délky

hrany buňky radMul, který určuje, jak musí být daleko nejbližší atom receptoru, aby se buňkadále dělila a okolí jader atomů receptoru inner, ve kterém již buňky nedělíme. Buňka se tedydělí na 8 potomků tehdy a jen tehdy, když není ve větší hloubce než maxDepth a zároveň

29

existuje alespoň jeden atom receptoru, jehož jádro je blíže, než radMul krát délka stranybuňky a zároveň neexistuje atom, který je blíže než inner krát polovina délky strany buňky.Díky tomuto kritériu lze balancovat přesnost adaptivní mřížky v závislosti na vzdálenosti odatomů receptoru – s tím, jak se ligand přibližuje k receptoru může přesnost výpočtu klesat,zůstávat stejná, ale i růst (což může být zajímavé pro dokovací úlohy, jelikož zadokovanýligand se zpravidla přímo „dotýká”3 svými atomy atomů receptoru). Hodnota inner, je-linenulová, může zastavovat rekurzi tam, kde je celá buňka příliš blízko jádra některého atomureceptoru, tedy je v pozici, kam se nemůžeme atomem ligandu dostat díky repulzivní van derWaalsově síle.Při prvním volání algoritmu zadáváme také pozici kořenové buňky −→pos a její velikost l.

Algorithm 1 Vytvoř buňku adaptivní mřížky

1: if maxDepth > 0 ∧ ∃i ∈ N ; |−→Ri −−→pos| < l · lmul ∧ ∀j ∈ N ; |

−→Ri −−→pos| > inner − l

2then

2: for i := 0 to 8 do3: vlož do −→p polohu i-tého rohu buňky4: children[i] :=Vytvoř buňku adaptivní mřížky(−→p , l

2,maxDepth − 1)

5: end for6: charges = ∅7: else8: for i := 0 to 8 do9:

−−−−−→force[i] := −→0

10: vlož do −→p polohu i-tého rohu buňky11: for r := 0 to r − 1 do12:

−−−−−→force[i] =

−−−−−→force[i] +

qi

|−→p −−→Ri|3

· (−→p −−→Ri)

13: end for14: end for15: children := ∅16: end if

V případě adaptivní mřížky již není možné tak přesně omezit složitost, jako u mřížkyregulární. Rychlost růstu časové a prostorové složitosti je možné určit v závislosti na maxi-mální hloubce rekurze maxDepth, násobku poloměru radMul, hodnoty inner a počtu atomůreceptoru r.Prostorovou složitost určíme jako O(r · (maxDepth · radMul3 − inner3)), jelikož pro r

atomů receptoru dělíme mřížku do hloubky maxDepth, přičemž na každé úrovni rekurzedělíme buňky vzdálené nejvýše radMul-násobek velikosti buňky, odtud radMul3 (těch je nakaždé úrovni stejně, jelikož radMul násobíme velikostí buňky, která s každou úrovní zanořeníklesne na polovinu). Složitost snižuje pro každý atom hodnota inner3, která vyjadřuje početbuněk, které již dále nevytvoříme kvůli přílišné blízkosti jádru některého atomu receptoru.Pro každý bod mřížky je nutno spočítat sílu vůči r atomům receptoru, časová složitost je

tedy O(r2 · (maxDepth · radMul3 − inner3)).V případě, že jsou atomy receptoru blíže, vychází reálná složitost nižší (o konstantní

násobek, nikoliv asymptoticky), protože se překrývají okolí atomů, ve kterých se buňky dělí(pro každý atom tedy nevytvoříme celých maxDepth · radMul3 − inner3 buněk).

3vzdálenost jader atomů odpovídá součtu van der Waalsových poloměrů

30

Interpolace

Zatímco generování adaptivní mřížky je poměrně jednoduché, interpolace na této mřížce nenípři požadavku na její spojitost triviální. Na rozdíl od regulární mřížky, kde jsou všechny buňkystejně velké, můžeme u adaptivní mřížky přecházet mezi buňkami různé velikosti. Tím neníautomaticky zajištěno, že interpolovaná síla na hraniční pozici dvou buněk je shodná, jelikožtyto dvě buňky nemusí sousedit shodnými stěnami či hranami (tzn. že interpolujeme prokaždou buňku z různě umístěných interpolačních bodů). Obr. 5.6 demonstruje výše zmíněnýpřípad (kříž reprezentuje polohu atomu ligandu, zvětšené puntíky jsou pak ty předpočítanésíly, které by byly zapojeny do interpolace).

Obrázek 5.6: 2D zobrazení přechodu atomu mezi různě velkými buňkami.

Tento problém lze řešit tak, že pro každou prostorovou dimenzi kontrolujeme, zda jealespoň některá sousední buňka menší. V takovém případě interpolujeme pro kvádr, kterývytvoří průmět menší sousední buňky přes buňku, ve které hledáme interpolovanou hodnotu.Hodnoty na opačné straně dopočítáme interpolací ve stěně menší z dvojice aktuální buňkaa soused. Případ, kdy bereme síly pro levou stěnu z menší sousední buňky a pro pravoudopočítáváme síly ze stěny buňky, ve které interpolujeme, je zobrazen na obr. 5.7.

Obrázek 5.7: 2D zobrazení přechodu atomu mezi různě velkými buňkami s opravenou inter-polací.

To nám ale ještě pro zajištění spojité interpolace nestačí. Zúžíme-li buňku využitím před-počítaných sil ze sousední buňky v dimenzi d1 a přesuneme-li se do jiné buňky v dimenzid2; d2 6= d1, můžeme ztratit souseda s přesnější silou. Proto je zapotřebí kontrolovat i buňkysdílící hranu s buňkou, ve které interpolujeme sílu, jak ukazuje obrázek 5.8.Náš algoritmus tedy musí pro danou polohu atomu ligandu

−−−−−−−→ligandPos nalézt buňku, ve

které se tento atom nachází a zmenšit ji, pokud má přesnější sousedy. Výstupem algoritmujsou údaje o buňce, tedy její pozice

−−−−−→cellPos, velikost

−−−−−−→cellScale a pole vektorů sil v rozích

cellForces. Z těchto údajů již můžeme počítat sílu působící na ligand způsobem popsanýmv 5.2.1.

31

Obrázek 5.8: Zobrazení přechodu atomu mezi buňkou zúženou dle svého souseda přes stěnudo buňky zúžené dle souseda přes hranu.

Algoritmy jsou zapsány poněkud méně exaktně, než bývá zvykem – autor zvolil tentopopis z důvodu jejich zpřehlednění.Nejprve najdeme v oktalovém stromě buňku, ve které se atom ligandu nachází (od ko-

řenové buňky sestupujeme do nelistových potomků obsahujících atom ligandu až dokud ne-narazíme na list). Následně pustíme algoritmy 2 a 3, abychom buňku zmenšili a zpřesnili nazákladě buněk z jejího okolí.Pro pochopení popisu algoritmů si ještě zavedeme několik pojmů (ty jsou, kromě formál-

ních definic, zachyceny na obrázku 5.9).Bod a náleží buňce C o počátku C0, C1, C2 a velikosti Cl0, Cl1, Cl2 právě tehdy, když leží

geometricky uvnitř buňky (tedy ∀d ∈ {0, 1, 2} platí ad ≥ Cd ∧ ad ≤ Cd + Cld.Nechť buňka C má pro bod −→a v dimenzi d předchůdce C ′ právě tehdy, když a leží v buňce

C a pro C ′ platí, že C ′

d + C ′

ld= Cd ∧ ∀i ∈ {0, 1, 2}; i 6= d platí C ′

i ≤ ai ∧ C ′

i + C ′

li≥ ai.

Nechť buňka C má pro bod −→a v dimenzi d následníka C ′ právě tehdy, když a leží v buňceC a pro C ′ platí, že C ′

d = Cd + Cld ∧ ∀i ∈ {0, 1, 2}; i 6= d platí C ′

i ≤ ai ∧ C ′

i + C ′

li≥ ai.

Nechť buňka C má pro bod −→a v dimenzi d souseda přes hranu C ′ právě tehdy, kdyžC ′

d ≤ ad ∧ C ′

d + C ′

ld≥ ad a ∀i ∈ {0, 1, 2}; i 6= d platí C ′

i +C ′

li= Ci ∨ C ′

i = Ci + Cli .

Obrázek 5.9: Ilustrace zavedených pojmů.

Interpolaci sil ze stěny provedeme analogicky interpolaci v regulární prostorové mřížce.Stěna je v R2, tudíž si vystačíme s osami x a y. Polohu bodu, pro který interpolujeme, udávávektor

−−−−→relPos. Jeho souřadnice jsou již normovány na velikost stěny a mají počátek v rohu

stěny, ze které interpolujeme.Hodnoty pole vah w získáme vztahem:

w[x][y] = (x · relPos0 + (1− x) · (1− relPos0))

32

·(y · relPos1 + (1− y) · (1− relPos1))kde x, y ∈ {0, 1}

Následný výpočet interpolované síly provedeme dle vztahu:

−→F =

1∑

x=0

1∑

y=0

w[x][y] ·−−−−−−−−−−−→cellForces[x][y]

Interpolace sil z hrany je již zcela triviální – sílu vypočítáme jako vážený součet dvou silna hraně, ze které interpolujeme. Váhou je zde relativní poloha bodu pro který interpolujeme.

Algorithm 2 Upřesni buňku pomocí sousedů přes stěnu

1: cell := Najdi buňku(−−−−−−−→ligandPos)

2: for d = 0 to 2 do3: oldCell := cell

4: do lowCell vlož předchůdce buňky cell v dimenzi d pro−−−−−−−→ligandPos

5: do hiCell vlož následníka buňky cell v dimenzi d pro−−−−−−−→ligandPos

6: if ∃lowCell ∧ (lowCell.scaled > cell.scaled) ∧ (∄hiCell ∨ (lowCell.scaled ≥hiCell.scaled)) then

7: for ∀ dimenze i; i 6= d do8: cell.scalei := lowCell.scalei

9: cell.posi := lowCell.posi

10: end for11: překopíruj síly z pravé stěny lowCell do levé stěny cell

12: if hiCell.scaled > cell.scaled then13: interpoluj síly z levé stěny buňky hiCell do pravé stěny buňky cell v dimenzi d14: else15: interpoluj síly z pravé stěny buňky oldCell do pravé stěny cell v dimenzi d16: end if17: end if18: if ∃hiCell∧(hiCell.scaled > cell.scaled)∧(∄lowCell∨(hiCell.scaled > lowCell.scaled))

then19: for ∀ dimenze i; i 6= d do20: cell.scalei := lowCell.scalei

21: cell.posi := lowCell.posi

22: end for23: překopíruj síly z levé stěny hiCell do pravé stěny cell

24: if lowCell.scaled > cell.scaled then25: interpoluj síly z pravé stěny buňky lowCell do levé stěny buňky cell v dimenzi d26: else27: interpoluj síly z levé stěny buňky oldCell do levé stěny cell v dimenzi d28: end if29: end if30: end for

Funkci algoritmu 2 si můžeme představit tak, že vytváří pro každou dimenzi jakési prů-měty buňky cell a jejích sousedů. Existuje-li předchůdce či následník buňky4, která pro polohu

4to nemusí nastat, pokud jsme již buňku zmenšili průmětem v některé předchozí dimenzi

33

−−−−−−−→ligandPos poskytuje v tomto průmětu přesnější síly, tzn. obaluje ji v průmětu menším čtver-cem (to testují řádky 6 a 18), zúžíme i buňku cell. Mají-li předchůdce či následník přesnějšísíly, interpolujeme nové síly na stěnách buňky cell ze sil předchůdce či následníka, v opačnémpřípadě k interpolaci využijeme síly buňky cell. Pokud zpřesňujeme síly na stěně buňky cell

pomocí jejích vlastních sil (řádky 15 a 27), využíváme kopii buňky vytvořenou na řádku 3,abychom si pro interpolaci udrželi původní velikost a polohu buňky.

Algorithm 3 Upřesni buňku pomocí sousedů přes hranu1: for d := 0 to 2 do2: oldCell := cell

3: do pole neighbourCell vlož sousedy přes hranu buňky cell pro−−−−−−−→ligandPos v dimenzi d

4: nechť i ∈ 0, 1, 2, 3;neighbourCell[i].scaled = max(neighbourCell[0].scaled, ...,

neighbourCell[3].scaled)5: if neighbourCell[i].scaled ≤ cell.scaled then6: pokračuj začátkem nového cyklu7: end if8: cell.scaled := neighbourCell[i].scaled

9: cell.posd := neighbourCell[i].posd

10: for i := 0 to 3 do11: if ∃neighbourCell[i] ∧ neighbourCell[i].scaled > oldCell.scaled then12: interpoluj síly z hrany buňky neighbourCell[i] sousedící s cell do hrany buňky cell

sousedící s neighbourCell[i]13: else14: interpoluj síly z hrany buňky oldCell sousedící s cell do hrany buňky cell sousedící

s oldCell

15: end if16: end for17: end for

Algoritmus 3 slouží ke zmenšení a zpřesnění buňky cell pomocí sousedů přes hranu. Prokaždou dimenzi d hledá souseda přes hranu, který vede k největšímu zpřesnění (tj. pro dimenzid má nejvyšší hodnotu scale[d]) a zúží podle něj buňku cell v dané dimenzi. Pro každou hranurovnoběžnou s osou dimenze d použije k zpřesnění sil na této hraně v cell buďto síly souseda(řádek 14), nebo síly uložené v cell (řádek 16), podle toho, kdo poskytuje větší přesnost.Každá buňka má konstantní počet sousedů (nejvýše 6 sousedů přes stěnu a 12 sousedů

přes hranu). Celkově je tedy složitost interpolace pro jeden atom v O(maxDepth), jelikožmusíme nalézt nejmenší buňku, která obsahuje jádro atomu a její sousedy.

Vlastnosti interpolace

Výše zmíněná interpolace nedává shodné výsledky vzhledem k pořadí probíraných dimenzí(cykly na řádku 2 algoritmu 2 a na řádku 1 algoritmu 3). Při zmenšení buňky průmětemv jedné dimenzi totiž buňka ztrácí sousedy v jiné dimenzi, jak demonstruje obrázek 5.10.V levé části obrázku zúžíme buňku v průmětu v 1. dimenzi tak, že v průmětu ve 2. dimenzinemůžeme buňku dále zužovat (svým zmenšením ztratila souseda, který by ji v 2. dimenzizužoval). V pravé části je pořadí dimenzí obrácené a buňku zde nejprve zúžíme ve smysluvertikální osy, čímž ale neztrácí souseda ve 2., horizontální dimenzi a zúžíme ji podruhé.

34

V implementovaném algoritmu je ovšem pořadí dimenzí fixováno, takže nevznikají ve výpočtužádné nespojitosti a výpočet je pro stejné parametry mřížky vždy konzistentní.

Obrázek 5.10: Rozdíl v interpolaci podle pořadí probíraných dimenzí

Tento přístup může také přinášet zbytečně nepřesnou interpolaci (zpřesnění dané sousedyv později probíraných dimenzích je uplatněno až při blízkosti interpolovaného bodu těmtosousedům, viz obrázek 5.11), avšak v naší aplikaci lze předpokládat sousedství buněk, jejichžvelikost se neliší více než o 1

2(to je dáno přímo kritériem pro dělení buněk mřížky při její

konstrukci – kolem každého jádra se snižuje přesnost buněk v jakýchsi kruzích, které vymezujepodmínka ∃i ∈ N ; |

−→Ri −−→pos| < l · lmul na řádku 1 algoritmu 1). V takovém případě nám dává

zde formulovaná interpolační metoda dostatečně přesné výsledky.

Obrázek 5.11: Nevhodná interpolace v případě velkých změn v měřítku sousedních buněk

Spojitost interpolace

V této kapitole vyslovím několik vět o vlastnostech interpolace na mřížce, ke kterým uvedutaké zdůvodnění jejich pravdivosti. Je zapotřebí terminologicky odlišit pojem buňka a buňkazmenšená sousedem. Zatímco v prvním případě se jedná o buňku tak, jak je reprezentovánav mřížce, buňka zmenšená sousedem je jakousi dočasnou datovou strukturou vzniklou přiinterpolaci, která přenáší přesnost danou sousedními buňkami tak, aby nedošlo k porušeníspojitosti interpolace. Nebude-li explicitně zmíněn pojem buňka zmenšená sousedem, je nutnointerpretovat slovo buňka jako buňka nativně reprezentovaná v mřížce.

Definice. Korektní zmenšení buňky C vzhledem k sousedům C1, ..., Cn je takové zmenšení,kdy nastavíme buňce C nové rozlišení a novou polohu a její síly interpolujeme z původníbuňky C a nebo sil sousedních stěn či hran buněk C1, ..., Cn, podle toho, která buňka poskytnejemnější rozlišení.

35

Lemma 5.2.2.1. Nechť stěna S′;S′ ∈ S má uloženy síly vzniklé interpolací ze stěny S.Hodnota síly v bodě a náležícím do S i S′ je shodná pro interpolaci z S i S′.

Důkaz. Pravdivost lemmatu přímo plyne z linearity použité interpolace. Pro bod, na kterémhledáme interpolovanou hodnotu, sčítáme síly uložené v rozích stěny násobené váhou, která jeúměrná vzdálenosti bodu, na kterém hledáme interpolovanou hodnotu, od jednotlivých rohů.Pokud vytvoříme ve stěně S novou stěnu S′, jejíž rohy mají hodnoty vzniklé interpolací na S,dává díky linearitě interpolace interpolace v S′ stejné hodnoty, jako interpolace téhož boduv S.

2

Pozn.: Argumentace je vedena ne příliš formální cestou – přiměl mě k tomu fakt, že se jednáo poměrně triviální tvrzení, u něhož by však precizní argumentace založená na algebraickémodvození interpolace na S′ pomocí interpolace na S byla neúměrně dlouhá. Ilustrace výšezmíněného argumentu je zachycena pro dva body na obrázku 5.12, rozšíření pro čtyři bodyje přímočaré.

Obrázek 5.12: Lineární interpolace z bodů stěny S a S′.

Lemma 5.2.2.2. Má-li některý z předchůdců či následníků buňky C v některé dimenzijemnější rozlišení, algoritmus 2 korektně zmenší buňku C alespoň na úroveň tohoto rozlišení.

Důkaz. Existuje-li předchůdce s jemnějším rozlišením, než má buňka C, a zároveň máalespoň tak jemné rozlišení jako následník C, je podmínka na řádku 6 kladná. Pak je v cykluna řádcích 7-10 buňka C zmenšena a zarovnána v dimenzích různých od d s předchůdcem.Následně jsou přepočítány její síly (z předchůdce jsou pouze překopírovány, na opačné stěněbuňky jsou interpolovány z buňky C nebo následníka, podle toho, kdo poskytuje jemnějšírozlišení).Analogický kód je proveden v bloku podmínky na řádku 18, má-li následník jemnější

rozlišení než C i předchůdce.Algoritmus tedy provede pro každou dimenzi zmenšení buňky na jemnější z rozlišení před-

chůdce a následníka, pokud je jemnější než rozlišení C a korektně upraví síly uložené v buňce.2

Lemma 5.2.2.3. Má-li některý ze sousedů přes hranu buňky C jemnější rozlišení, algorit-mus 3 korektně zmenší buňku C alespoň na úroveň tohoto rozlišení.

36

Důkaz. Algoritmus najde pro každou dimenzi všechny 4 sousedy přes hranu. Je třeba uká-zat, že po provedení algoritmu nezůstane žádný soused přes hranu, který by měl jemnějšírozlišení. Taková situace nemůže nastat, jelikož na řádku 4 vybereme souseda s nejnižšímrozlišením, na řádku 8 a 9 nastavíme rozlišení a zarovnání buňky C podle tohoto souseda.Pro korektní zmenšení buňky je zapotřebí ještě přepočítat síly v ní uložené – to provedeme

v cyklu na řádcích 10-16, kde pro všechny 4 hrany rovnoběžné s osou dimenze d přepočítámesíly buďto ze sousedící hrany souseda přes hranu, nebo z hrany buňky C, podle toho, kdoposkytuje jemnější rozlišení.

2

Tvrzení 5.2.2.1. Hodnota síly F interpolovaná na adaptivní prostorové mřížce je vzhledemk pohybu atomu spojitá.

Důkaz. Na adaptivní mřížce můžeme zaznamenat tyto případy, pro které je třeba ukázat,že nedojde k porušení spojitosti:

• pohyb uprostřed buňky

• přechod mezi buňkami shodné velikosti

• přechod mezi buňkou s jemnějším a buňkou s hrubším rozlišením

• přechod mezi buňkou zmenšenou sousedy a buňkou s hrubším rozlišením

• přechod mezi dvěma sousedy zmenšenými buňkami

Pohybujeme-li se uvnitř buňky, nebo přecházíme-li mezi buňkami stejné velikosti, dokazujespojitost interpolace tvrzení 5.2.1.1.Existuje-li soused C ′, který má jemnější rozlišení, než buňka C, pro kterou interpolujeme

sílu F , je buňka C korektně zmenšena dle svého souseda, jak ukazují lemma 5.2.2.2 a lemma5.2.2.3. Proto při přechodu mezi C a C ′ procházíme přes shodné stěny, hrany či body, cožzachovává spojitost interpolace.Je-li buňka C zmenšena sousedem C ′ a jiný soused C ′′ má rozlišení shodné či jemnější než

původní buňka C, ale hrubší než zmenšená buňka C, není při přechodu mezi C a C ′′ buňkaC ′′ zmenšována a přecházíme tedy mezi dvěma různě velkými buňkami. Síly na stěně či hraněsousedem zmenšené buňky C sousedící s C ′′ však vznikly interpolací ze sil z C ′′ nebo z C,pokud mají stejné rozlišení, takže při přechodu z C do C ′′ nedojde k porušení spojitosti, cožukazuje lemma 5.2.2.1.Analogicky má-li soused C ′′ hrubší rozlišení než C před zmenšením,vznikly síly na sousední stěně či hraně mezi C a C ′′ interpolací z C. Při přechodu do C ′′ jeC ′′ zmenšena dle původní C, ovšem interpolace na stěně zmenšené C i původní C je shodná.Sousedí-li spolu dvě buňky vzniklé zmenšením pomocí sousedů hranou či stěnou rozdílné

velikosti, jsou nezarovnané body hrany či stěny interpolovány z hrany či stěny o shodnévelikosti pro obě tyto buňky (pokud by nebyly, jedna z buněk by byla sousedem zmenšujícím tudruhou). Spojitost interpolace takovéhoto přechodu ukazuje lemma 5.2.2.1. U buněk sousedícístejně velkou hranou či stěnou plyne spojitost interpolace z lemmatu 5.2.2.2 a 5.2.2.3, jelikožty provádí zmenšení buněk korektně, musí mít sousedící hrany i stěny shodné síly.

2

37

5.2.3 Možná vylepšení adaptivní mřížky

V této kapitole letmo rozeberu další možné pokračování práce na adaptivní mřížce. Jedná seo pouhé myšlenky možného vývoje, proto zde neuvádím konkrétní algoritmy ani nepředkládámdůkazy o správnosti mých úvah, které se mohou v praxi ukázat jako kontraproduktivní. Jelikožvšak tuto problematiku považuji za přínosnou a v implementaci prezentované v této diplomovépráci je stále hodně ke zlepšování, považuji za účelné popsat, jak vidím další možnosti jejíhozlepšení.

Zlepšení prostorové a časové složitosti vygenerování mřížky

Současná implementace adaptivní mřížky ukazuje její potenciál (viz kapitola 6.4), má všakještě značné rezervy, spočívající především ve zbytečně vysoké multiplikativní konstantě pro-storové a časové složitosti. Tato konstanta je skryta v předpočítávání a ukládání sil pro všechnyrohy každé buňky – ty jsou v případě shodně velkých sousedních buněk zbytečně propočítá-vány vícekrát pro stejnou pozici. V případě, že jsou všechny buňky stejně velké (např. přivelmi hustém receptoru a velmi vysoké hodnotě radMul), je mimo hranice mřížky každá sílavypočítána a uložena 8×.Tomuto omezení lze předejít (za cenu zhoršení multiplikativní konstanty složitosti inter-

polace) ukládáním pouze těch sil, které nejsou obsaženy v sousedních buňkách. Buňka pakmusí nést informaci, které síly jsou uloženy v ní a které se mají hledat u sousedů, což všaklze uložit do 8-bitvého pole.

Paralelizace výpočtu mřížky

Generování octree mřížky je možno paralelizovat. Paralelizace však není zcela přímočará,protože je mřížka silně nepravidelná – pokud bychom např. provedli na jednom výpočetnímuzlu dělení mřížky do určité úrovně a následně pouze nechali jednotlivé výpočetní uzly dáledělit vzniklé buňky, hrozí riziko, že některé uzly skončí s dělením příliš brzy a jiné budou ještědlouho pracovat. Je tedy zapotřebí, aby se uzly vzájemně informovaly o dokončení své práce,aby jim mohly ty uzly, které pracují, přidělit buňky k dalšímu dělení.

Zrychlení interpolace

V naší aplikaci nepůsobí rychlost interpolace žádné potíže (výpočet je stále výrazně rych-lejší nežli výpočet van der Waalsových interakcí). Pokud by však měla být mřížka nasazenav software, kde je předpočítána i van der Waalsova síla, byla by rychlost interpolace proaplikaci velmi důležitá. V takovém případě je možno urychlit hledání sousedních buněk ulože-ním seznamu ukazatelů na ně pro každou buňku. Tím se nám poněkud zhorší multiplikativníkonstanta časové i prostorové složitosti jejího generování, nicméně získáme možnost provéztrekurzivní sestup v mřížce pouze jednou pro jeden výpočet interpolace, zatímco v současnéimplementaci musíme hledat buňku obsahující polohu interpolované síly, jejích 8 sousedů přesstěnu a 12 sousedů přes hranu.Domnívám se, že pro zrychlení interpolace je také možno předpřipravit buňky, které nor-

málně zmenšujeme pomocí jejich sousedů. Buňky, které v současné implementaci zmenšujemepři interpolaci, by tedy byly ve zmenšeném stavu vytvořeny už při vytváření mřížky. Sílyv jejich rozích by byly vytvořeny stejnou metodou, jako při interpolaci (tedy kopírováním silsousedů, interpolací ze sil sousedů nebo sil vlastních). Takovéto buňky by nepřidávaly mřížce

38

na přesnosti a zvyšovaly by její prostorovou složitost. Při interpolaci bychom však již nemuselivůbec přihlížet k okolním buňkám, protože měli-li by mít na buňku nějaký zpřesňující efekt,byl tento efekt do mřížky zanesen už při její konstrukci.

Zlepšení interpolace pro mřížku s velkými rozdíly hustoty sousedních buněk

Jak již bylo vysvětleno v části Vlastnosti interpolace kapitoly 5.2.2, implementovaná inter-polace nevykazuje příliš uspokojivé chování v případě, kdy spolu sousedí buňky s velkýmirozdíly v rozlišení. Pokud bychom používali adaptivní mřížku v aplikaci, v níž by k takovétosituaci docházelo, bylo by zapotřebí modifikovat algoritmus zužování buňky tak, aby při kaž-dém zúžení vytvořil novou buňku. Výsledná buňka by pak obsahovala síly vzniklé interpolacígradientů ve všech buňkách vzniklých zúžením.

Odlišná reprezentace adaptivní mřížky

S reprezentací adaptivní mřížky pomocí rekurzivně se dělících krychlí se poměrně snadnopracuje, ale při detailnějším rozboru problematiky narazíme na nutnost ošetření speciálníchsituací daných sousedstvím různě velkých buněk. To vede jednak k nutnosti implemento-vat složitější algoritmus, ale také k větší výpočetní náročnosti interpolace. Domnívám se, žeje možné se tomu vyhnout, pokud opustíme myšlenku uzavírání atomů do krychlí. Pokudbychom např. prostor rekurzivně dělili stejným způsobem jako u současné implementace, aleukládali síly pouze pro středy buněk, pak uzavřeme každý atom ligandu např. do čtyřstěnu.Interpolaci v něm by pak bylo třeba vážit podle vzdálenosti atomu od stěn tak, aby stěna,pokud se v ní nachází atom, měla plnou váhu (čímž zajistíme spojitou hodnotu interpolovanésíly při přechodu mezi čtyřstěny). Je však třeba zvážit, jak obtížné by bylo nalézt mezi před-počítanými silami čtyřstěn, který obsahuje atom, pro který interpolujeme; zdali by takovýtočtyřstěn byl nalezen jednoznačně (tzn. algoritmus by garantoval, že nebude atom liganduv další iteraci výpočtu uzavřen do jiného čtyřstěnu, aniž by opustil původní); zdali by novýčtyřstěn, do kterého by atom vstoupil z původního, sdílel s původním stěnu, přes kterou byatom přecházel. Konečně také není jisté, jestli by hledání čtyřstěnů šlo provézt rychleji, nežliinterpolace s prozkoumáváním sousedů v mé implementaci.Odlišná reprezentace adaptivní mřížky (např. pomocí čtyřstěnů) by se tedy mohla ukázat

jako vhodnější, než ta, kterou popisuji v rámci své diplomové práce (přinesla by o konstantnínásobek lepší prostorovou složitost i časovou složitost interpolace). Je však zapotřebí naléztuspokojivé odpovědi na výše zmíněné otázky.

5.3 Detaily implementace

5.3.1 MPI

Rozlišujeme tři základní typy paralelního programování:

• Datový paralelismus, kdy shodné instrukce pracují nad různými daty. Toto řešení jetypické pro vektorové počítače.

• Sdílená paměť představuje model, kdy každý proces vidí všechna data. Toto řešenílze nejsnáze využívat na strojích realizujících symetrický multiprocesoring (zde majívšechny CPU rovnoprávný přístup do celé paměti).

39

• Při předávání zpráv má každý proces svou lokální paměť a s ostatními procesy komuni-kuje explicitním zasláním či přečtením zprávy.

MPI, z anglického Message Passing Interface, je standardizované rozhraní pro komunikaciprocesů paralelního programu pomocí předávání zpráv. To má výhodu oproti modelu datovéhoparalelismu a sdílené paměti ve snadném nasazení na skupině počítačů propojenou sítí (sou-časně lze ale použít i na multiprocesorových strojích, kde se zprávy předávají přes sdílenoupaměť namísto sítě). Velkou výhodou je také snažší laditelnost systémů využívající předávánízpráv (komunikace přes zprávy nemůže nepozorovaně, jako v případě sdílené paměti, přepsatdata patřící jednomu procesu druhým procesem).Rozhraní MPI poskytuje základní funkce pro komunikaci mezi procesy, od zjištění, kolik

procesů běží ve výpočtu a jaký má který proces identifikátor (kolikátý z kolika procesů je),přes jednoduché blokující i neblokující operace odeslání a přijmutí zprávy až po složitějšíkolektivní operace, realizující broadcast či skládání zpráv od různých procesů.Při zasílání zpráv explicitně definujeme jejich datový typ, takže je možné do výpočtu

zapojit stroje různé architektury, aniž bychom museli ručně zajišťovat konverzi dat.Podrobný popis MPI lze nalézt např. v [19].MPI definuje pouze standard, pro potřeby vývoje software je tedy zapotřebí zvolit kon-

krétní implementaci. Těch existuje několik, např. LAM5, MPICH6 či OpenMPI7.

5.3.2 Architektura paralelní aplikace

V kapitole 3.3 je popsána architektura neparalelizované verze aplikace. V paralelní verzi zů-stává zachována původní struktura výpočetní části aplikace v master procesu, ve slave pro-cesech je speciální verze vlákna výpočet simulace, jak ukazuje obrázek 5.13.Protože slave procesy vůbec nekomunikují se zobrazovací částí (rozdělení práce viz 5.1.2),

zasílá vlákno řízení master procesu některé údaje do slave procesů pomocí MPI. Velká částoperací je však prováděna pouze v master procesu (záznam a reprodukce, nastavování pružinyaj.). Výpočetní vlákna slave procesů dostávají polohu ligandu, hledají blízké atomy, pro kterépočítají van der Waalsovy síly a ty zasílají zpět výpočetnímu vláknu master procesu.

5.3.3 Problém využití všech procesorů

Při simulaci na víceprocesorovém stroji jsme narazili na výkonové problémy aplikace při alo-kaci všech procesorů. Jednak docházelo k zadrhávání PHANToMu, kdy zařízení nějakou dobukladlo neúměrně velký odpor, a navíc se nedeterministicky stávalo, že celý výpočet běžel ex-trémně pomalu (aplikaci stačilo spustit několikrát po sobě a rychlost výpočtu byla i o 2 řádyrůzná).Nakonec jsme identifikovali dva problémy:

• Knihovna libphantom8 si pro ovládání zařízení PHANToM vytváří vlákno, které trávívětšinu času neaktivním čekáním na IRQ zařízení. Dojde-li k saturaci všech dostupnýchCPU náročnými procesy, probouzí někdy plánovač toto vlákno na příliš nízké frekvenci.To se projeví tak, že se PHANToM zasekne v jedné pozici a uživatel pocítí nepříjemné

5http://www.lam-mpi.org6http://www.mcs.anl.gov/research/projects/mpich27http://www.open-mpi.org8http://www.fi.muni.cz/ xslaby/phantom.html

40

Obrázek 5.13: Architektura paralelní aplikace.

cuknutí. V případě, že je výpočetními procesy využito méně procesorů, je thread po-chopitelně spouštěn dostatečně často. Rozhodli jsme se tedy tomuto vláknu a vláknuobsluha haptiky přidělit realtime prioritu, která zmíněný problém odstraňuje.

• Slave procesy trávily mimo výpočet van der Waalsových sil čas v blokujících MPI funk-cích, kdy čekaly, až budou moci odeslat výsledky výpočtu či přijmout data udávajícípolohu ligandu. Jelikož MPI funkce při aktivním čekání předávají neustále řízení já-dru, očekávali jsme, že plánovač dokáže při zátěži nějakého CPU procesem mimo našiaplikaci přidělit polovinu CPU času našemu výpočetnímu procesu a druhou polovinutomuto cizímu procesu. V praxi však docházelo (opět nedeterministicky) k tomu, že sejednomu nebo více slave procesů téměř nedostával CPU čas. Tím se celý výpočet zablo-koval, protože všechny ostatní procesy aktivně čekaly na jeden proces, který se nemohldostat k CPU času. Jako řešení se nakonec ukázalo přepsat MPI operace na neblokujícía dokud se nedokončí, tak cyklicky volat funkci sched yield()9, aby se všechny procesys dokončeným výpočtem vzdaly CPU času a dostal se na řadu proces, na který čekají10.

9ta zařadí proces, ze kterého je volána, na konec plánovací fronty10volání funkce usleep() není možné, protože by bylo probouzení procesů závislé na časovači, který je pronaše účely příliš pomalý

41

Kapitola 6

Měření výkonu

K měření výkonu jsem používal dvě dvojice receptor-ligand:

• lektin PA-IIL (3 258 atomů) a alfa-L-fukosa (23 atomů)

• beta-cyklodextrin (147 atomů) a katechin (35 atomů)

Receptory se liší významně jak svou velikostí, tak tvarem (PA-IIL je poměrně prosto-rově kompaktní, beta-cyklodextrin tvoří prstenec) a jsou proto vhodnými reprezentanty protestování aplikace.

Obrázek 6.1: Receptory beta-cyklodextrin a PA-IIL.

42

6.1 Zobrazování

K testu jsem použil tyto stroje:

• phantom.fi.muni.cz, 2x Quad Core Intel Xeon E5335 (2,0GHz), 8GB RAM, nVidia Ge-Force 8600GT

• lynkos.fi.muni.cz, 2x AMD Athlon MP 2600+ (2,133GHz), 2GB RAM, nVidia GeForce7600GT

• osobní notebook, AMD Turion 64 x2 TL-52 (1,8GHz), 1GB RAM, nVidia GeForce Go7300

Notebook jsem zvolil jako reprezentanta stroje s relativně slabou grafickou kartou, neboťse domnívám, že aplikace může být často využita (alespoň pro přehrávání zaznamenanýchvýsledků simulace) na kancelářských strojích, které obvykle nemívají příliš výkonné grafickéakcelerátory. Ze srovnání jeho výsledků s výsledky naměřenými na ostatních strojích je taképatrné, kdy zobrazování limituje výkon CPU a kdy GPU. Stroj lynkos je připojen k dvěmaprojektorům s polarizačními filtry a testoval jsem na něm výkon při zapnuté stereovizi.Výkon zobrazování byl měřen při iniciálním pohledu na receptor a ligand (pohled nastaví

aplikace po svém startu, přiblížení je nastaveno tak, aby byl vidět celý receptor) v rozlišení1280 × 1024 obrazových bodů.

6.1.1 Dávkování geometrie

Rychlost zobrazování jsem měřil dvěma způsoby:

• u výpočetní části aplikace jsem nastavil, aby zasílala data k zobrazení 1000× za sekundua následně měřil počet snímků, které dokázala za jednu sekundu zobrazit zobrazovacíčást aplikace

• při normální rychlosti odesílaných dat jsem měřil, kolik procent CPU spotřebuje zobra-zovací část aplikace

Zatímco první zmíněný způsob demonstruje celkovou propustnost systému a umožňujenám tak odhadnout, jak velkou molekulu dokážeme zobrazovat, druhý se zaměřuje pouze nazátěž CPU a je pro nás důležitý, pokud spouštíme výpočetní i zobrazovací část aplikace najednom systému.Z tabulky 6.1 vidíme, že optimalizace dávek má na výkon značný vliv, a to především

při zobrazování lektinu PA-IIL, který je tvořen řádově více atomy. U notebooku není nárůstvýkonu tak dramatický, zřejmě tedy narážíme na omezení výkonu daném grafickou kartou.Při vypnuté optimalizaci dávek nedokázal u lektinu PA-IIL žádný ze strojů zobrazovat při

vypnuté optimalizaci dávek dostatečnou rychlostí (viz tabulka 6.1) a CPU tedy bylo vytíženovždy na 100 %. U výsledků naměřených na stroji phantom je vidět, jak při dostatečně rychlégrafické kartě významně ulehčí optimalizace dávek procesoru. Na tomto stroji není problémspustit výpočetní část aplikace paralelně pro všech jeho osm jader a současně na spustitzobrazování.

43

Tabulka 6.1: Dávkování geometrie – snímků za sekunduStroj optimalizace dávek PA-IIL beta-cyklodextrin

ano 117,18 317,05notebook

ne 19,65 209,37

ano 252,50 723,28phantom

ne 22,95 391,22

ano 115,19 222,78lynkos – mono

ne 13,53 134,87

ano 62,69 145,27lynkos – stereo

ne 6,97 78,06

Tabulka 6.2: Dávkování geometrie – zátěž CPUStroj optimalizace dávek PA-IIL beta-cyklodextrin

ano 32 % 12 %notebook

ne 100 % 20 %

ano 3 % 0 %phantom

ne 100 % 2 %

ano 19 % 9 %lynkos – mono

ne 100 % 15 %

ano 35 % 16 %lynkos – stereo

ne 100 % 30 %

6.1.2 Dynamická úroveň detailů

Při přiblížení receptoru roste počet trojúhelníků na atom a výkon aplikace se mění rozdílně připoužití dynamické úrovně detailů (LOD) a bez něj. Zároveň jsou ale některé atomy již mimopohled a nejsou renderovány, což kompenzuje přibývání detailů u jednotlivých atomu. Protojsem se rozhodl měřit výkon pouze při iniciálním pohledu na receptor, stejně jako v případědávkování geometrie.Při zobrazování PA-IIL s alfa-L-fukosou bylo vykreslováno 186 378 trojúhelníků na jeden

snímek při LOD a 3 151 362 trojúhelníků bez LOD. Zde je přínos LOD velmi významný, cožodráží i výsledky naměřené v 6.3. Při zobrazování beta-cyklodextrinu s katechinem bylo zapoužití LOD vykreslováno 11 770 trojúhelníků a bez LOD 176 322 trojúhelníků, což je takévýznamný rozdíl, nicméně na výkonu se příliš neprojevuje (grafická karta při tomto počtutrojúhelníků není úzkým hrdlem).

6.2 Detekce kolizí

Rychlost detekce kolizí jsem měřil pouze na stroji phantom, jelikož zde není příliš zajímavéporovnání mezi různými stroji. Údaje jsou získávány při přehrávání pohybu ligandu v okolízadokované polohy nepřekračující 0,5 A, která je na detekci kolizí nejnáročnější – je-li ligand

44

Tabulka 6.3: Dynamická úroveň detailů – snímků za sekunduStroj LOD PA-IIL beta-cyklodextrin

ano 252,50 723,28phantom

ne 88,73 589,41

ano 115,19 222,78lynkos – mono

ne 33,30 213,15

ano 62,69 145,27lynkos – stereo

ne 17,89 129,74

daleko od receptoru, je hashovací algoritmus velmi rychlý, nicméně nás zajímá rychlost při conejvyšším počtu nalezených kolizí. Počítal jsem průměrný a nejhorší čas detekce kolizí za 100iterací dynamiky.

Tabulka 6.4: Detekce kolizí mezi PA-ILL a katechinemOřez detekce kolizí průměrný čas nejhorší čas

hashování 69µs 141µs2 A

naivní 478µs 524µs

hashování 201µs 236µs4 A

naivní 494µs 565µs

hashování 317µs 457µs6 A

naivní 512µs 533µs

Z tabulky 6.4 jasně vidíme, že čas běhu algoritmu využívajícího hashování prostoru jezávislý na počtu nalezených kolizí (ten roste s rostoucí hodnotou ořezu). Zatímco pro nízkéhodnoty ořezu je algoritmus velmi rychlý, pro vyšší hodnoty (okolo 6 A) se již blíží jehorychlost rychlosti naivní implementace testování všech dvojic. Máme-li tedy ořez, který připřiblížení ligandu zasáhne velkou část celého receptoru, není již algoritmus hashování prostorupříliš účinný (je s určitými režiemi vybrán velký počet atomů receptoru, u nichž musímetestovat exaktní kolize). Pokud bychom však pracovali s větším receptorem, než je lektin PA-ILL, pak by pracoval hashovací algoritmus obdobnou rychlostí (je závislý pouze na hustotěatomů v blízkosti ligandu), zatímco kontrola všech párů by se zpomalila. To ostatně dokládáměření na beta-cyklodextrinu, kde je při ořezu 2 A hashovací algoritmus obdobně rychlý jakou PA-ILL (u větších ořezů je rychlejší, jelikož detekuje kvůli nedostatečné velikosti beta-cyklodextrinu méně kolizí).Jelikož větší ořez počítáme obvykle s použitím více CPU a detekce kolizí je taktéž pa-

ralelizována, provedl jsem ještě měření při alokaci CPU odpovídající náročnosti úlohy1. Tozachycuje tabulka 6.5.Čas detekce kolizí pomocí hashování prostoru se zde drží na zhruba konstantní hodnotě.

Přizpůsobujeme-li tedy počet alokovaných procesorů náročnosti celého výpočtu, nezačne se

1Tj. počet CPU, který je třeba alokovat, abychom drželi rychlost iterace výpočetní části aplikace pod 1ms

– tento počet vychází z empirického měření rychlosti, neodpovídá tedy přesně asymptotické složitosti výpočtu

van der Waalsových interakcí O( c3

k), kde c je hodnota ořezu a k počet výpočetních uzlů.

45

Tabulka 6.5: Detekce kolizí mezi PA-ILL a katechinemOřez #CPU detekce kolizí průměrný čas nejhorší čas

hashování 69µs 141µs2 A 1 CPU

naivní 478µs 524µs

hashování 52µs 114µs4 A 4 CPU

naivní 117µs 140µs

hashování 43µs 65µs6 A 8 CPU

naivní 62µs 90µs

vlivem zvyšování hodnoty ořezu detekce kolizí zpomalovat.

6.3 Paralelní běh

6.3.1 Prostředí

Aplikace byla kompilována na strojích manwe3 a clusteru konos pomocí gcc a g++ verze4.0.2, na stroji phantom ještě pomocí verze 4.2.1. Na všech strojích byla zapnuta volba -O5.Pro testy jsem použil MPI implementaci LAM 7.1.Běh aplikace jsem měřil jak na víceprocesorových strojích, tak také na clusteru – konkrétně

na těchto:

• phantom.fi.muni.cz, 2x Quad Core Intel Xeon E5335 (2,0GHz), 8GB RAM

• manwe3.ics.muni.cz, 8x Dual Core AMD Opteron 885 (2,6GHz), 64GB RAM

• cluster konos.civ.zcu.cz, složen ze strojů 2x Dual Core AMD Opteron 270 (2,0GHz),4GB RAM

Aplikace ve své současné implementaci není optimalizovaná pro běh v heterogenním pro-středí – pokud bychom ji spouštěli na clusteru s různě rychlými stroji, bude její výkon omezenrychlostí nejpomalejšího stroje, jelikož rozděluje výpočet mezi stroje co nejvíce rovnoměrně.

6.3.2 Metodika měření

Čas běhu algoritmu jsem zkoumal u simulace, kde byl katechin zaveden do dokovací pozicelektinu PA-IIL. Toto zavedení bylo provedeno pomocí uložené trajektorie ligandu, aby ne-došlo k ovlivnění výsledku uživatelem. V dokovací pozici je ligand jednak dostatečně blízkoreceptoru, ale také je v místě, které je pro nás nejvíce zajímavé.Časy výpočtu byly měřeny jednou za 100 iterací. Z nich jsem vybíral nejlepší z posledních

pěti od momentu, kdy byl ligand zadokován. Nejlepší čas bereme kvůli odstínění případnýchkrátkých vytížení sítě či stroje.

6.3.3 Škálování algoritmu

Měření škálování algoritmu jsem prováděl na stroji manwe pro různé hodnoty ořezu van derWaalsových sil.

46

Obrázek 6.2: Výpočet na stroji manwe při různých hodnotách ořezu.

Obrázek 6.2 znázorňuje průběh rychlosti výpočtu dle alokovaných procesorů pro různéhodnoty ořezu van der Waalsových sil. Jak vidíme, pro zachování frekvence výpočtu 1 000Hzje pro ořez 4 A zapotřebí alespoň 6 CPU, pro ořez 6 A 12 CPU a u ořezu 8 A je zapotřebízapojit prakticky všechny procesory. Zároveň je z grafu patrné lepší škálování algoritmu přivyšší hodnotě ořezu. Toto pozorování platí i pro jiné stroje a je dáno uniformnějším rozdělenímjednotlivých atomů receptoru mezi výpočetní uzly.

6.3.4 Vliv latence sítě

Obrázek 6.3: Výpočet na strojích manwe a konos, ořez 10A.

Graf 6.3 srovnává multiprocesorový stroj manwe s clusterem konos. Ten má pomalejšíprocesory než manwe, pro srovnání škálování nám to však nevadí – podstatná je zde efektivita.U výpočtu na clusteru je vidět omezení latencí sítě, které u většího počtu alokovaných CPUprakticky vyváží přínos z rozložení výpočtu (do čtyř procesorů může aplikace běžet pouzena jednom stroji z clusteru konos – spomalení vlivem latence se tedy projevuje až u vyššího

47

počtu alokovaných CPU). Cluster začne v absolutním čase výpočtu předbíhat manwe až přivýrazně větší výpočetní náročnosti dané zvětšením hodnoty ořezu, jak je vidět v grafu 6.4.

Obrázek 6.4: Výpočet na strojích manwe a konos, ořez 20A.

U takto velké hodnoty ořezu algoritmus velmi dobře škáluje – na stroji manwe i clusterukonos dokáže do čtyř alokovaných CPU zkracovat dobu výpočtu s malou odchylkou úměrněpočtu CPU. U vyššího počtu CPU se již negativně projeví síťové spojení clusteru konos, takžese u 8 a 16 CPU dostane efektivita na konosu pod efektivitu na manwe. Úbytek efektivity jepřitom v případě většího ořezu nižší.Měření ukázalo pokles efektivity při čtyřech alokovaných CPU u stroje manwe, který

nenastal u konosu. Důvod tohoto poklesu jsme dále nezkoumali, je možné, že je způsobenoneuniformním přístupem do paměti či složitější synchronizací cache procesorů u manwe.

6.3.5 Srovnání CPU

Obrázek 6.5: Výpočet na strojích phantom a manwe při hodnotě ořezu 8A.

Na grafu 6.5 vidíme srovnání výpočtu s ořezem 8A na strojích manwe a phantom.

48

Při kompilaci pomocí gcc verze 4.0.2 spotřebuje stroj phantom při jednom alokovanémCPU cca 1,57-násobek času výpočtu na stroji manwe, přičemž procesory manwe běží na 1,3-násobku frekvence. Pokud však použijeme novější verzi kompilátoru gcc 4.2.1, výkon strojůna architektuře AMD se prakticky nezmění2, zatímco stroje postavené na architektuře Corefirmy Intel získají v naší aplikaci 1,82-násobný výkon. Tím se dostává Intel Xeon E5335 dovedení před AMD Opteronem 885 (manwe potřebuje při jednom alokovaném procesoru najednu iteraci 1,16-násobek času, za který zvládne výpočet stroj phantom).

6.3.6 Využitelnost simulace na clusterech

U výpočtu na clusteru je limitující latence sítě spojující jednotlivé výpočetní uzly. Díky tomuse u běžně používané sítě (1Gbit ethernet) nedá dosáhnout 1 000 iterací výpočtu silovéhopůsobení za sekundu, což nám ale při zachování frekvence 1 000Hz u obsluhy PHANToManemusí nutně vadit. Pokud se zpomalí výpočet sil mezi ligandem a receptorem, jeví se liganduživateli těžší, protože pomaleji reaguje na silové působení haptiky. Díky tomu je manipulaces ligandem obtížnější a méně intuitivní, subjektivně však lze konstatovat, že frekvence okolo300Hz ještě přináší použitelné chování. Výpočty na clusterech lze tedy využít, potřebujeme-livyšší výpočetní výkon (např. při požadavku na velkou hodnotu ořezu van der Waalsových sil)a nemáme-li k dispozici dostatečně výkonný multiprocesorový stroj.Je také možné, že využitím rychlejšího síťového rozhraní, jako je např. Infiniband3, se

sníží latence komunikace mezi výpočetními uzly natolik, že bude možné počítat molekulovéinterakce rychlostí 1 kHz a zcela se tím vyhneme nutnosti používat velké multiprocesorovéstroje.

6.4 Regulární versus adaptivní mřížka

Podle předpokladu by měla být adaptivní mřížka při nejmenších buňkách shodných s buňkamiregulární mřížky stejně přesná v blízkosti receptoru a s rostoucí vzdáleností by měla růst jejínepřesnost. Zároveň by měla být menší než mřížka regulární (což znamená i rychlejší propředpočítání, jelikož většinu prostoru zabírají uložené vektory sil, které se počítají u oboumřížek shodně).K měření přesnosti jsem vytvořil záznam, který nejprve zadokoval ligand a následně jej táhl

pomalou rychlostí zpět z dokovací pozice. Samotné měření jsem prováděl od momentu, kdy bylligand zadokován a naměřenou nepřesnost ukládal odděleně pro intervaly vzdálenosti 0–1 A,1–2 A až 24–25 A. Přesnost interpolace na mřížce byla ověřována vůči přesnému výpočtu síly(tedy interakce každého atomu ligandu s každým atomem receptoru). Napočítané nepřesnostipro každý interval jsou průměrem nepřesností ze všech výpočtů sil v tomto intervalu.U hodnoty nepřesnosti interpolace je třeba vzít v úvahu, že zcela korektní porovnání

mřížek by bylo dosti náročné. Současná implementace přehrávání pohybu trajektorie neza-ručuje determistický pohyb ligandu (např. vlivem nepřesnosti výpočtu u konkrétní mřížkyse v jednom kroku simulace dostane ligand do jiné polohy, než by se dostal u jiné mřížky av následujícím kroku již simulace vychází z jiných parametrů). Díky tomu měříme nepřesnostinterpolace pro ligand v jiném umístění. To může mít vliv na zjištěnou přesnost mřížky, proto

2Ověřeno na stroji vanda.ics.muni.cz s procesorem Athlon 64 X2 5000+, na strojích METACentra jsemneměl možnost novým gcc kompilovat.

3http://www.infinibandta.org/home

49

nelze brát naměřené hodnoty absolutně, ale pouze z nich interpretovat trendy a vzájemněmřížky porovnávat.

Obrázek 6.6: Adaptivní vs. regulární mřížka.

Na obrázku 6.6 vidíme srovnání adaptivní a regulární mřížky na lektinu PA-ILL.Adaptivní mřížka pokrývá krychli o velikosti 1 (v měřítku simulace, přepočteno na An-

gströmy cca 142 A). Mřížka se rekurzivně dělí do maximální hloubky 10, přičemž koeficientvzdálenosti lmul, který určuje, v jakém okruhu se musí nacházet atom receptoru, abychombuňky nadále dělili, má hodnotu 8. Hodnota inner zastavující dělení buněk v blízkosti jádraatomu je 1,5. Při takovýchto parametrech má adaptivní mřížka 95MB.Pokud bychom chtěli pokrýt adekvátní prostor regulární mřížkou obdobné jemnosti, které

může dosahovat adaptivní mřížka, museli bychom mít pro každou dimenzi 1 000 buněk. Ta-ková mřížka zabere 16GB. Jelikož má použitý receptor poměrně kompaktní tvar, který jemožno poměrně přesně uzavřít do krychle, lze kolem něj vytvořit mřížku o velikosti strany0,25 – máme tedy na každou dimenzi 250 buněk o velikosti 0,001. Taková mřížka má již ak-ceptovatelných 250MB, ovšem musíme zcela zanedbávat síly působící ve větší vzdálenosti odreceptoru. Je třeba si ale uvědomit, že při méně šťastném tvaru receptoru není možné mřížkutakto zmenšit a museli bychom použít výrazně větší mřížku. Obdobné zmenšení prostoru po-krytého mřížkou nepřinese u adaptivní mřížky příliš významné úspory místa, jelikož ve většívzdálenosti od receptoru jsou pouze buňky s nízkým rozlišením, které zabírají malé množstvípaměťového prostoru.Je vidět, že s rostoucí vzdáleností ligandu od dokovací pozice klesá chyba regulární mřížky,

zatímco u adaptivní mřížky dojde nejprve k poklesu, ale následně (se snížením rozlišenímřížky) chyba vzroste. Pokud se však podíváme na hodnotu absolutní chyby, vidíme, žerůst relativní chyby nastane až v momentě, kdy je hodnota síly velmi malá. V takových mís-tech pro nás není silové působení natolik důležité, protože se ztrácí v síle, kterou liganduudělujeme haptickým zařízením.Rozdílná hodnota chyby v malé vzdálenosti ligandu od dokovací pozice je zřejmě dána

šťastnějším rozmístěním buněk adaptivní mřížky (ty nejsou s regulární mřížkou zarovnány),jinak by měla být hodnota chyby v této vzdálenosti velmi podobná.Obrázek 6.7 srovnává přesnost mřížek s nižšími prostorovými nároky. Zatímco u regulární

mřížky můžeme jedině zvětšit buňky a tím zmenšit jejich počet, u adaptivní mřížky můžemenechat běžet rekurzi do stejné hloubky a pouze snížit hodnotu lmul, čímž snížíme přesnostinterpolace v místech vzdálenějších od receptoru, ale v dokovací pozici jsme schopni přesnostzachovat.Při těchto parametrech zabírá adaptivní mřížka 41MB, regulární v plné velikosti 2GB a

při zúžení okolo receptoru 31MB.

50

Obrázek 6.7: Adaptivní vs. regulární mřížka.

51

Kapitola 7

Závěr

Úkolem této diplomové práce bylo nalézt, implementovat a popsat metody, které by zrych-lily běh aplikace pro haptickou interakci biomolekul a tím přiblížily její využitelnost reálnýmchemickým problémům. Domnívám se, že v rámci funkcionality aplikace a jejích omezení(silové pole odvozené od pole AMBER, rigidní konformace ligandu i receptoru) bylo tohotocíle dosaženo. Aplikace nyní umožňuje práci s receptorem o tisících atomů, přičemž se lzedomnívat, že by tento počet mohl být i řádově vyšší. Na běžně dostupném hardware si mů-žeme dovolit ořez van der Waalsových sil na 6-8 A, což pro naše účely poskytuje dostatečnoupřesnost. Elektrostatické síly předpočítáváme do mřížky, což je běžné i v ostatních dokova-cích programech. Taktéž zobrazovací část aplikace dokáže vykreslovat všechna potřebná datadostatečnou rychlostí i na slabším hardware.Původní implementace aplikace dokázala pracovat s receptorem o cca 150 atomech a

ořezem van der Waalsových sil na 2 A, zcela bez elektrostatických sil. Větší receptor byl i přinízké hodnotě ořezu náročný na detekci kolizí i zobrazování.Kromě průzkumu a aplikace známých algoritmů, využitých v zobrazovací části a detekci

kolizí, jsem navrhl a implementoval také několik vlastních řešení. Paralelizace výpočtu van derWaalsových interakcí a algoritmu detekce kolizí nebyla technicky náročná, bylo však zapotřebídůkladně analyzovat možné aspekty chování paralelního algoritmu, aby bylo možné vyhnoutse případnému nechtěnému chování danému nevhodně zadanými daty. Při implementaci jsemnarazil na nečekané problémy s rychlostí paralelního výpočtu a ovládáním haptiky při saturacivšech procesorů v systému, které se však podařilo překonat.Velká část práce je věnována popisu a argumentaci o korektnosti mého algoritmu pro efek-

tivní uložení mřížky s předpočítanými elektrostatickými silami a spojitou interpolaci v rámcitéto mřížky. Domnívám se, že mnou prezentovaný algoritmus může být (při definici odlišnéhokritéria pro dělení buněk mřížky) využitelný u velké škály aplikací pracujících s vektorovýmipoli a vyžadující spojitou interpolaci na celém prostoru reprezentovaném polem.V práci jsou také nastíněna některá témata k dalšímu vývoji či výzkumu. Pro případ

využití flexibilního receptoru jsem popsal způsob, jak měnit jeho zobrazovanou geometriipomocí vertex shaderu, ostatní algoritmy jsou přímo navrhovány tak, aby byly využitelné ipo přidání možnosti pracovat s flexibilním receptorem. Další vývoj vyžaduje také adaptivníprostorová mřížka, která v současné implementaci nevyužívá veškerý svůj potenciál a v prácijsou nastíněny možnosti odlišné implementace uložení mřížky i interpolace.

52

Literatura

[1] Petter Bivall Persson, Matthew D. Cooper, Anders Ynnerman, Lena A. E. Tibell, Bengt-Harald Jonsson, Evaluating the Effectiveness of Haptic Visualization in BiomolecularEducation - Feeling Molecular Specificity in a Docking Task. 12th IOSTE Symposium.2006.

[2] Ouh-young M, Pique M, Hughes J, Srinivasan N, Brooks Jr FP, Using A ManipulatorFor Force Display In Molecular Docking. Proceedings of IEEE Robotics and AutomationConference, Vol. 3. 1988.

[3] Ouh-young M, Beard DV, Brooks Jr FP, Force Display Performs Better than VisualDisplay in a Simple 6-D Docking Task. Proceedings of IEEE International Conferenceon Robotics and Automation. 1989.

[4] Weghorst S, Augmenting Tangible Molecular Models. The Fourth International Workshopon Spatial Media, Japan. 2003.

[5] Susana K. Lai-Yuen, Yuan-Shin Lee, Interactive Computer-Aided Design for MolecularDocking and Assembly. Computer-Aided Design and Applications. 2006.

[6] Petter Bivall Persson, Matthew D. Cooper, Lena A.E. Tibell, Shaaron Ainsworth, AndersYnnerman, and Bengt-Harald Jonsson, Designing and Evaluating a Haptic System forBiomolecular Education. IEEE Virtual Reality. 2007.

[7] Nagata H, Mizushima H, Tanaka H., Concept and prototype of protein-ligand dockingsimulator with force feedback technology. Bioinformatics. 2002.

[8] Tomáš Golembiovský, Záznam a reprodukce průběhu řízené simulace v haptické aplikaci[bakalářská práce]. Fakulta informatiky Masarykovy univerzity v Brně. 2006.

[9] Matthias Wloka, ”Batch, batch, batch:” What Does It Really Mean? Game DevelopersConference. 2003.

[10] www.ogre3d.org

[11] I. J. Palmer and R.L. Grimsdale. Collision Detection for Animation using Sphere-Trees.Computer graphics forum. 1995.

[12] Radka Svobodová Vařeková, Jaroslav Koča and Chang-Guo Zhan. Complexity and Con-vergence of Electrostativ and van der Waals Energies within PME and Cutoff Methods.International Journal of Molecular Sciences. 2004.

53

[13] Gisbert Schneider and Hans-Joachim Böhm. Virtual screening and fast automated doc-king methods. Drug Discovery Today, Vol. 7. 2002.

[14] J. Kenneth Salisbury and Mandayam A. Srinivasan. Phantom-Based Haptic Interactionwith Virtual Objects. IEEE Computer Graphics and Applications, Vol. 17. 1997.

[15] Zdeněk Kabeláč, Aleš Křenek. Studying Conformational Behaviour with Phantom Device.The First PHANToM Users Research Symposium (PURS’99). 1999.

[16] Otaduy, M.A. and Lin, M.C. Stable and responsive six-degree-of-freedom haptic manipu-lation using implicit integration. Proceedings of the First Joint Eurohaptics Conferenceand Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems.2005.

[17] Andreas Griewank and David Juedes and Jean Utke. ADOL-C: A Package for the Auto-matic Differentiation of Algorithms Written in C/C++. ACM Transactions on Mathe-matical Software (TOMS), Vol. 22. 1996.

[18] J.E. Colgate, M.C. Stanley and J.M. Brown. Issues in the haptic display of tool use.International Conference on Intelligent Robots and Systems. 1995.

[19] Gropp, William. MPI - the complete reference: the MPI-2 extensions. Vol. 2. Cambridge:MIT Press. 1998.

[20] Sergio Filipe Sousa, Pedro Alexandrino Fernandes, and Maria Joao Ramos. Protein-Ligand Docking: Current Status and Future Challenges. PROTEINS: Structure,Function, and Bioinformatics, Vol. 65. 2006.

[21] Lennard-Jones, J. E. Cohesion. Proceedings of the Physical Society. Vol. 43. 1931.

[22] London F. Zur Theori und Systematik der Molekularkrafte. Zeitscgrift fur Physik, Vol.63. 1930.

[23] Garrett M. Morris, David S. Goodsell, Robert S. Halliday, Ruth Huey, William E. Hart,Richard K. Belew, Arthur J. Olson. Automated Docking Using a Lamarckian GeneticAlgorithm and an Empirical Binding Free Energy Function. Journal of ComputationalChemistry, Vol. 19. 1998.

[24] Aleš Křenek and Jiří Filipovič. Haptická simulace interakce molekul se semiimplicitníintegrací Newtonovské dynamiky. Nepublikované výsledky. 2006.

[25] Yong Duan, Chun Wu, Shibasish Chowdhury, Mathew C. Lee, Guoming Xiong, WeiZhang, Rong Yang, Piotr Cieplak, Ray Luo, Taisung Lee, James Caldwell, Junmei Wang,Peter Kollman. A point-charge force field for molecular mechanics simulations of proteinsbased on condensed-phase quantum mechanical calculations. Journal of ComputationalChemistry Vol. 24. 2003.

54