Молекулярна динаміка з використанням пакета NAMD

Размещено на http://www.allbest.ru/

Зміст

1. Молекулярна динаміка з використанням пакета NAMD

2. Програмний пакет VMD і підготовка до моделювання

3. Моделювання кисню в рівновазі Постановка завдання



1. Молекулярна динаміка з використанням пакета NAMD

Ідея молекулярної динаміки, звичайно, не нова. За останні десятиліття було розроблено кілька програмних пакетів, які дають змогу проводити молекулярне моделювання. Деякі з них безкоштовні з відкритим вихідним кодом (такі як, наприклад, пакети NAMD і GROMACS), існують також комерційні програмні продукти - CHARMM, AMBER. Ціна останніх для некомерційних організацій, таких як ВНЗ, невисока, тому багато дослідницьких груп воліють використовувати саме ці програмні пакети. У цій частині ми зупинимося на пакеті NAMD і розглянемо покроково, як даний програмний продукт може бути використаний для проведення моделювання простої білкової системи в явному розчиннику. У зв'язку з тим, що початкова конформація біомолекул виходить експериментально або з використанням численних алгоритмів, проведення моделювання в фізіологічних умовах вимагає декількох попередніх допоміжних кроків, які часто називають приготуванням системи (рис. 1.1).

Далі детальніше зупинимося на кожному з перерахованих етапів моделювання. Початкові дані

Для проведення моделювання за допомогою пакета молекулярної динаміки потрібні такі початкові дані:

Файл координат (.pdb). Цей файл містить координати атомів у декартовій системі координат.

Файл топології (.psf). У цьому файлі міститься опис всіх атомів у системі, включаючи їх типи, заряди, маси, а також інформація про ковалентні зв'язки: пари атомів, пов'язаних ковалентно, трійки атомів, що утворюють кути, четвірки атомів, що утворюють торсійні кути.

Файл параметрів силового поля (наприклад, par_all22_prot.inp). У ньому містяться параметри для потенціалів силового поля - постійні пружності ковалентних зв'язків, рівноважні відстані, кути і т.д.



Рис. 1.1 Попередні (підготовчі) етапи моделювання. Програму VMD використовують для побудови топології, додавання розчинника й іонів. У NAMD на початкових етапах проводять мінімізацію енергії, нагрівання й еквілібрацію системи. Після завершення цих етапів можна приступати до безпосереднього моделюванню з використанням пакета NAMD.

Як початкові координати атомів найчастіше постають дані експериментальних досліджень, зокрема рентгенівська кристалографія або експеримент ЯМР (ядерний магнітний резонанс). Більшість наукових проектів, пов'язаних із молекулярним моделюванням, починають із пошуку відповідних тривимірних структур на таких інтернет ресурсах, як база даних білкових структур (Protein Data Bank). На цьому ресурсі зібрана величезна кількість результатів експериментів із покликаннями на оригінальні статті дослідників, які проводили експеримент. Варто відзначити, що особливості експериментальних процедур не завжди дають змогу точно визначити місце розташування всіх атомів у системі, а положення водню в рентгенівської кристалографії взагалі не визначається. Тому, перед тим як почати процес моделювання тієї чи тієї біологічної молекули, необхідно ретельно вивчити отриману експериментаторами структуру.

Другим кроком є створення файлу топології системи (для пакета NAMD цей файл має розширення .psf). У файлі топології наявна інформація про атоми, що входять в систему, а також про зв'язким цих атомів у білку. Для того щоб отримати цей файл, треба використовувати програму VMD, створену тієї самою науковою групою, яка створила програмний пакет NAMD. У найпростіших випадках файл топології можна побудувати за допомогою автоматичної побудови через графічний інтерфейс програми VMD. При цьому буде створено два файли - .psf і оновлений .pdb. Останній буде містити недозволені в експериментах атоми водню і невеликі недозволені ділянки поліпептидного ланцюга, якщо такі є.

2. Програмний пакет VMD і підготовка до моделювання

Програмний пакет VMD пропонує зручний інтерфейс для візуалізації біомолекул, побудови біомолекулярних структур та первинної обробки результатів моделювання. Також цей пакет тісно інтегрований з програмою NAMD і дає змогу готувати всі необхідні для цієї програми дані. VMD розповсюджують безкоштовно у вигляді відкритого коду, і його можна завантажити з сайту університету Іллінойсу США. На сайті також можна знайти велику кількість документації щодо пакета для самостійного і більш глибокого його вивчення.

Файл координат .pdb

Для того щоб почати роботу з біомолекулою, необхідно отримати її тривимірну структуру. Одне з найпопулярніших джерел - база даних білкових структур. Сама база даних була заснована в 1971 році в Брукхейвенській національній лабораторії. Тоді в цій базі даних містилося всього 7 структур, однією з яких стала структура білка міоглобін, розв'язана ще в 1958 році групою Джона Кендрю (John Kendrew). Для вивчення тривимірної структури міоглобіну, Джону Кендрю і його групі довелося вручну створювати модель молекули (у статті наведено її фотографії). Для обміну результатами експерименту в ті часи використовували перфокарти. Оскільки кожен атом дозволеної структури зберігався на окремій перфокарті, передавання результатів була досить трудомістким процесом - координати міоглобіну займали близько 1000 перфокарт.

Централізоване сховище даних допомогло не тільки істотно спростити такий обмін, але і зробити розв'язані структури доступнішими. У 80_х роках, разом із можливостями експерименту, зріс і інтерес до цієї області. Кількість структур почала зростати експоненційно. Незабаром технічний прогрес дав змогу спростити і спосіб обміну розв'язаними структурами, а з появою мережі інтернет і персональних комп'ютерів - результати роботи кристалографів стали доступні кожному. На сьогодні база даних білкових структур містить близько 80000 структур різних біомолекул. Найбільш популярним форматом файлу координат атомів біомолекулярної системи донині залишається формат .pdb, у якому кожному атому системи відведений один рядок. Формат має деякі обмеження, що збереглися ще з часів перфокарт. Зокрема, довжина кожного рядка файлу не може перевищувати 80 символів, кожному з яких відведена своя роль. Оскільки основною метою формату .pdb було забезпечити можливість обмінюватися результатами експериментальних досліджень, він передбачає наявність великої кількості інформації про використану в експерименті процедуру, кількість невирішених атомів і так далі. Нас головно буде цікавити опис атомів системи, розбиття якого за стовпцями наведено в таблиця 1.1.



Таблиця 1.1 Опис формату .pdb

Стовпець	Тип	Містить
1-6	Фіксований вираз	Строка «ATOM »(з двома пробілами)
7-11	Ціле	Порядковий номер атома
13-16	Рядок (4 символа)	Назва атома
17	Символ	Індикатор неточності позиціонування атома (якщо у файлі є дублікат з іншими координатами)
18-20	Рядок (3 символа)	Ім'я амінокислоти
22	Символ	Ідентифікатор поліпептидного ланцюга
23-26	Ціле	Порядковий номер амінокислоти
27	Символ	Код додавання амінокислот
31-38	Число з плав. комою (8.3)	Координата X атома в декартовій системі (в ангстремах)
39-46	Число з плав. комою (8.3)	Координата Y атома в декартовій системі (в ангстремах)
47-54	Число з плав. комою (8.3)	Координата Z атома в декартовій системі (в ангстремах)
55-60	Число з плав. комою (6.2)	Заповненість (в кристалографічному експерименті)
61-66	Число з плав. комою (6.2)	Температурний множник
73-76	Рядок (4 символа)	Назва сегмента молекули (вирівнювання зліва)
77-78	Рядок (2 символа)	Ім'я хімічного елемента (вирівнювання справа)
79-80	Ціле	Заряд атома

Як приклад розглянемо структуру одного з доменів тітіну з індексом в базі даних білкових структур 1TIT. Рядок опису першого атома має такий вигляд:

1.......10........20........30........40........50........60........70........80

ATOM 1 N LEU A 1 -19.761 9.730 -6.893 1.00 1.53 N

Тобто це - атом азоту з амінокислоти лейцин, яка є першою амінокислотою в даному білку. Незважаючи на те, що в .pdb-файлі міститься детальна інформація про структуру білка, уявити собі цю структуру, очевидно, неможливо. Піонери кристалографії збирали тривимірні моделі білків, що, на щастя, в наші дні зовсім не обов'язково: є ціла низка програмних пакетів, які дають змогу візуалізувати структуру білка на комп'ютері.

Перегляд тривимірного зображення молекули

У програмі VMD три основних вікна - вікно терміналу, вікно зображення тривимірної структури і головне вікно (рис. 1.2). У вікні терміналу показують інформацію про роботу програми, крім того, в ньому можна безпосередньо писати команди за допомогою мови сценаріїв TCL. У вікні тривимірного зображення користувач має можливість повертати, рухати і наближати структуру молекули. Управління програмою можна здійснювати за допомогою графічного інтерфейсу головного вікна. Наприклад, щоб завантажити структуру домену тітіну, необхідно через меню "File" додати молекулу, або вказавши шлях до завантаженого файлу .pdb, або безпосередньо через інтернет, написавши у вікні чотирисимвольний PDB-код молекули (1TIT).

Рис. 1.2 Три основних вікна програми VMD: головне вікно, вікно терміналу і показу тривимірного зображення.

Після завантаження структури у вікні зображення тривимірної структури відобразиться молекула у вигляді ліній ковалентних зв'язків. При цьому ковалентні зв'язки відобразяться, ґрунтуючись як на інформації з файлу .pdb, де вони можуть бути задані явно, так і на основі радіуса обрізки (тобто якщо атоми знаходяться досить близько, VMD пов'язує їх ковалентно). Для того, щоб змінити репрезентацію молекули, у головному вікні є меню "Display" і пункт "Representations". Щоб показати елементи вторинної структури, можна використовувати метод окреслення "New cartoons" і виділити елементи вторинної структури кольором ("Secondary structure" в випадному списку "Coloring method") Кількість репрезентацій необмежена: можна намалювати частину молекули і іншим методом.

Наприклад, щоб відобразити наведений вище атом азоту, необхідно натиснути кнопку "Create rep", у рядку вибору "Selected atoms" набрати "resid 1 and name N" і встановити метод окреслення "VdW", який відобразить цей атом у вигляді сфери.

Для того, щоб перемикатися між режимами роботи миші, можна використовувати меню "Mouse" або гарячі кнопки на клавіатурі. Наприклад, можна відобразити інформацію про атом (клавіша 1) або відстань між двома атомами (клавіша 2). Детальнішу інформацію про виділені атоми або відстані можна дізнатися в пункті "Labels" меню "Graphics". Там само можна побудувати графік динаміки спостережуваних величин у випадку, коли в програму завантажено відразу багато структур (як, наприклад, при візуалізації траєкторії молекулярної динаміки).

Далі зупинимося на функціоналі цієї програми, який використовують для підготовки всіх файлів, що описують молекулярну систему і необхідні для проведення моделювання.

Підготовка до моделювання

Початкові координати (файл PDB)

Хоча зазвичай моделювання системи починається з пошуку її тривимірної структури в базах даних білкових структур (наприклад www.pdb.org), тут ми розглянемо простіший приклад, у якому часові інтервали процесів, що відбуваються, дають змогу отримати результат навіть на звичайному персональному комп'ютері. Тому ми синтезуємо структуру, що складається з восьми послідовно сполучених аланінів. Відомо, що аланін є однією з амінокислот, яка найчастіше трапляється в б-спіралях, процес утворення яких досить швидкий. Для того, щоб синтезувати структуру, використовуємо розширення Molefacture програми VMD. У цій підпрограмі виберемо утиліту побудови білків - Protein Builder, де додамо 8 аланінів у витягнутій конформації (Straight, торсійні кути: ц = -180? і ш = 180?). Коли система буде готова, координати атомів необхідно зберегти у форматі PDB через меню Файл.

Рис. 1.3 Структура поліаланіну у витягнутій конформації (торсійні кути ц = -180? і ш = 180?). Синім показані атоми вуглецю, зеленим - азоту, червоним - кисню, сірим - водню.

Створення топології системи (файл PSF)

Файл координат не несе в собі інформації про типи атомів в системі і про її зв'язаність. Ця інформація потрібна для розрахунку потенційної функції силового поля. Типи атомів, на відміну від їхніх імен, не є унікальними ідентифікаторами того чи того атома. Наприклад, атом Cб амінокислоти аланін в силовому полі CHARMM27 матиме тип CT1 - атом вуглецю з одним приєднаним воднем. Точно такі самі типи будуть і в Cб-атомів із інших аланінів у молекулі. Ім'я атома точно визначає місце розташування цього атома всередині амінокислоти, тип атома визначає його властивості при розрахунку потенційної функції, параметри якої залежать від типів атомів, утягнутих у ту чи ту взаємодію. Під зв'язаністю системи розуміють інформацію про взаємодії всередині системи - які взаємодії розраховувати для конкретних атомів у системі. Наприклад: атом номер 1 пов'язаний із атомом номер 2 ковалентно, а атоми 1, 2 і 3 утворюють ковалентний кут. У процесі моделювання цю інформацію використовують для складання списків ковалентних зв'язків, кутів і так далі. Потім типи атомів, що беруть участь у взаємодіях, зіставляють із типами у файлі параметрів силового поля. Отже, визначають потенційну функцію системи. Інформацію про типи атомів і їх зв'язаність зберігають у файлі структури білка (англ. Protein Structure File, PSF).

Створення файлу PSF може бути нетривіальним завданням. Особливо це стосується гетерогенних систем, наприклад таких, як глікопротеїни. Це відбувається тому, що в цих системах наявні нестандартні зв'язки - зв'язки між цукрами в полісахаридах і ковалентні зшивки між білками і сахаридами. У білкових системах процес створення файлу PSF значно простіший - програма має передати координати поліпептидних ланцюгів і дисульфідні зв'язки між ними, якщо такі є. Інформацію про зв'язаності всередині амінокислот, типи атомів і їхні часткові заряди беруть із файлу топологій використовуваного силового поля.

У прикладі з поліаланіном ми будемо використовувати розширення програми VMD, яке надає графічний інтерфейс для побудови файлу PSF. Для цього завантажимо отриманий раніше файл координат PDB в VMD і відкриємо вікно розширення "Automatic PSF Builder". Оскільки досліджувана система є поліпептидом, налаштування програми за замовчування повинні без проблем згенерувати файл PSF. Єдине, на що слід звернути увагу - файл топологій має відповідати використовуваному нами силовому полю CHARMM27 - "topall27_prot_na.rtf".

Формат файлу PSF

Файл PSF розпадається на кілька секцій, кожна з яких починається з кількості елементів у даній секції та її назви. Перша секція описує атоми системи. У нашому випадку їх 83, про що свідчить рядок заголовка секції:

1 83 !NATOM

2 1 P1 1 ALA N NH3 -0.300000 14.0070 0

3 ...

4 83 P1 8 ALA HB3 HA 0.090000 1.0080 0

Потім ідуть 83 рядки, кожен із яких описує один із атомів, для якого наведено його порядковий номер, ім'я сегмента молекули (поліпептидного ланцюга - "P1"), номер амінокислоти в цьому ланцюзі (від 1 до 8), ім'я амінокислоти ("ALA"), ім'я атома, тип атома, частковий заряд і масу в атомарних одиницях.

Далі у файлі PSF йде опис ковалентних зв'язків, при цьому в кожному рядку вказано до 4-х зв'язків (до 8 атомів). У наступному прикладі атом 1 пов'язаний ковалентно з атомом 5, атом 2 - з атомом 1, атом 1 - із атомом 3, атом 4 - з атомом 1 і так далі:

1 82 !NBOND: bonds

2 1 5 2 1 3 1 4 1

3 ...

4 80 82 80 83

У секціях, що описують ковалентні кути, торсійні кути і неправильні торсійні кути, атоми зібрані по 3 і 4 штуки відповідно. При цьому на кожен рядок випадає 3 кути (9 атомів) або 2 торсійних кути (8 атомів):

1 147 ! NTHETA : angles

2 1 5 6 1 5 11 2 1 5

3 ...

4 81 80 83 81 80 82 82 80 83

5

6 199 ! NPHI : dihedrals

7 1 5 7 8 1 5 7 9

8 ...

9 79 78 80 83

10

11 15 !NIMPHI : impropers

12 11 5 13 12 13 11 15 14

13 ...

14 76 71 78 77

Отже, файл PSF визначає, для яких атомів розраховувати ту чи ту взаємодію, а також вказує на набір параметрів, використовуваних для розрахунку. Останнє робиться через присвоєння типів для всіх атомів системи.

Додавання розчинника (води) та іонів

Водний розчин

Природне середовище для будь-якої біологічної системи - вода. У даному прикладі воду буде описано явно, тобто в модельованій системі будуть наявні молекули води. Оскільки ці молекули не зв'язані з пептидом і між собою ковалентно, це також вимагатиме введення граничних умов, які будуть перешкоджати розльоту молекул у просторі. Найпопулярнішим методом є введення періодичних граничних умов.

Перед тим, як додавати молекули води, необхідно оцінити лінійний розмір системи. У кінцевій системі білок повинен перебувати досить далеко від меж періодичної комірки. Вибір розмірів комірки дуже важливий. З одного боку, чим вони більші, тим більше буде потрібно молекул води, щоб її заповнити і, як наслідок, тим більше розрахунків буде потрібно на кожному кроці інтегрування. З другого боку, занадто маленька відстань від білка до межі з періодичною коміркою може призвести до взаємодії молекули зі своїм образом в наступній комірці. Зазвичай, розмір комірки вибирають, виходячи з розмірів молекули l, і беруть рівним L = l + 2d, де d - відступ від меж. На вибір d впливає безліч факторів: який саме обчислювальний експеримент буде проведено, наскільки сильно заряджена молекула. Наприклад, якщо передбачають моделювання денатурації, то розміри молекули будуть збільшуватися і d має бути достатньо великим, щоб в періодичну комірку вмістилася денатурована молекула. У разі моделювання денатурації під дією зовнішньої сили, періодичну комірку можна збільшити в одному напрямку, який збігається з напрямком дії сили, позаяк малоймовірне збільшення лінійного розміру молекули в напрямках, перпендикулярних до напрямку дії сили. У разі рівноважних симуляцій потрібно враховувати можливі обертання молекули і використовувати в оцінці максимальний лінійний розмір системи.

Виберемо для поліаланіну d ~ 6 Е і знайдемо (приблизно) максимальний лінійний розмір системи l. Для цього у вікні VMD можна показати відстань між термінальними атомами молекули - азотом першої амінокислоти та карбоксильним вуглецем останньої. Вибравши показ відстані між атомами з меню Mouse і кликнувши послідовно на ці два атоми, можна отримати поточне значення відстані в ~ 27 Е. Отже, розмір періодичної комірки буде L = 27 + 2 Ч 6 ~ 40 Е.

У розширеннях VMD є утиліта, що дає змогу додати комірку розчинника ("Extensions" > "Modelling" > "Add Solvation Box"). У вікні розширення необхідно вказати шляхи до файлів pdb і psf, створених на попередньому кроці. Також слід вказати розміри комірки: у нашому випадку мінімальні і максимальні значення будуть дорівнювати 20 і 20 Е у всіх трьох напрямках. Також варто дозволити програмі повертати молекулу для зменшення об'єму. Після виконання програми будуть створені оновлені файли координат і топології (pdb і psf), які міститимуть інформацію як про поліпептидну частину системи (молекула поліаланіну), так і про додану воду.

Іони

У фізіологічних умовах у водному розчині наявні розчинені речовини, головно іони натрію і хлору. Хоча концентрація цих речовин може відрізнятися в різних частинах організму, солі мають дуже великий вплив на фізіологічні процеси. Відхилення концентрації від норми може призвести до серйозних патологій функції білків і, в деяких випадках, навіть до смерті організму. З погляду моделювання, відсутність іонів може також призвести до некоректної роботи численних процедур. Зокрема, за використання періодичних умов необхідно, щоб повний електростатичний заряд модельованої системи залишався нульовим. Для додавання іонів в нашу систему скористаємося розширенням "Add Ions" програми VMD. Оскільки наша система нейтральна, кількість доданих негативно і позитивно заряджених іонів повинна бути однаковою, а у розширенні "Add Ions" є можливість вказати бажану концентрацію солі замість явної вказівки бажаної кількості іонів. Після вказівки концентрації і виконання програми в системі деякі молекули води будуть замінені на атоми іонів (рис. 1.4), унаслідок чого будуть згенеровані оновлені файли psf і pdb. Саме ці файли описують нашу систему (psf) і її початкову конформацію (pdb) та будуть використані далі програмою NAMD в моделюванні.

Рис. 1.4 Поліаланін, розчинений у воді без додавання іонів (ліворуч) і з ними (праворуч). Вода показана лініями ковалентних зв'язків, іони натрію і хлору - жовтими і блакитними сферами відповідно. Чорними лініями показана межа періодичної комірки.

Моделювання системи та обробка результатів

Рівноважне моделювання

Після того, як всі підготовчі етапи виконані, можна переходити до моделювання системи. У даному прикладі розглядаємо рівноважне моделювання, конфігураційний файл для якого дуже схожий на використаний при еквілібрації енергії. Як і раніше, як початкові координати будемо використовувати кінцеві координати попереднього етапу, а файл топології (psf) залишається незмінним:

1 structure pala_ionized .psf

2 coordinates pdb / pala_equil . coor

Ми будемо використовувати стандартний крок і експериментальне значення температури (1 фс і 300 К відповідно). Основною відмінністю рівноважних симуляцій від еквілібрації і нагріву є те, що швидкості атомів перепризначатися не будуть. Хоча за тривалого моделювання це може призвести до зміщення температури від заданої спочатку, у нашому випадку довжина траєкторії в 107 кроків (10 нс) дає змогу просто відключити перепризначення швидкостей. Також змінимо ініціалізуюче значення для генератора випадкових чисел.

1 timestep 1.0

2 temperature 300

3 numsteps 10000000

4 seed 322223322

У параметрах вихідних даних, крім імен файлів, можна змінити і частоту збереження енергії та координат:

1 outputenergies 1000

2 dcdfreq 1000

3 dcdfile dcd/ pala_quench .dcd

4 binaryoutput no

5 outputname pdb/ pala_quench

Утворення вторинної структури

Основною метою даної вправи було визначити час, необхідний для формування б-спіралі в короткому пептиді. Існує кілька різних алгоритмів, які дають змогу охарактеризувати вторинну структуру білка з координат його атомів. Основні параметри, що їх використовують у цих методах - значення торсійних кутів пептидного (скелетного) зв'язку та геометрія водневих зв'язків. У програмі VMD наявна вбудована можливість визначення вторинної структури, яка використовує алгоритм STRIDE.

Візуально наявність тієї чи тієї вторинної структури можна перевірити за допомогою колірного кодування, що визначається в графічних репрезентаціях молекули в програмі VMD. Варто відзначити, що в поточній версії програми, вторинна структура визначається лише для одного набору координат системи, і для того, щоб конфігурація вторинних структур перевизначалася динамічно, необхідно в TCL-консолі програми VMD виконати скрипт "sscache". Хоча вторинна структура полі аланіну постійно змінюється, по закінченні приблизно 5 нс молекула більшу частину часу перебуває в б-конформації (Рис. 40).

Рис. 1.5 б-спіраль, що утворилася в процесі рівноважного моделювання поліаланіну.

Також серед розширень VMD є утиліта "Timeline", яка дає змогу показати динаміку зміни вторинної структури для всіх амінокислот білка. У меню "Calculate" необхідно вибрати "Calc. Sec. Struc. "(Англ. Calculate Secondary Structure - розрахувати вторинну структуру). Після закінчення розрахунку у вікні утиліти з'явиться діаграма вторинних структур, колірне кодування якої наведено в нижньому лівому кутку вікна (Рис. 41).

Рис. 1.6 Динаміка вторинної структури білка, розрахована за допомогою утиліти "Timeline" програмного пакета VMD. По горизонтальній осі відкладено час моделювання, по вертикальній - вісім амінокислот молекули. Фіолетовим показані б-спіралі.

Графічні процесори (ГП)

Загальні обчислення за допомогою графічних процесорів

Завдяки зростаючим потребам ринку комп'ютерних ігор і обробки зображень програмовані графічні процесори еволюціонували у високопаралельні багатопотокові обчислювальні пристрої, які мають потужні обчислювальні можливості й велику пропускну здатність пам'яті. Основною відмінністю архітектури графічних процесорів (ГП) від центральних процесорів (ЦП) є те, що графічні процесори з самого початку проектували для високопаралельних обчислень під час обробки тривимірних зображень, тому вони відразу налаштовані на швидку обробку великих обсягів даних, а не на їх кешування і більш складний потік команд. Схематично це показано на рисунку Рис. 2.1.

Велика кількість логічних елементів кеш-пам'яті і контролера на ЦП дає змогу їм швидко виконувати складні обчислювальні процедури, маючи швидкий доступ до даних у пам'яті. На ГП більшу кількість логічних елементів відведено для обчислень, а модулі контролю процедур і кешу зменшені. Завдяки цьому щільність Арифметичних логічних пристроїв (АЛП) на ГП істотно вища, ніж на ЦП. Для того, щоб завантажити всі АЛП обчисленнями, ГП повинен виконувати багато операцій одночасно. Це здійснюється за допомогою обчислювальних потоків, кожен із яких виконує ті самі арифметичні операції, але на різних елементах масиву даних. Наприклад, на ЦП складання двох векторів розмірності M необхідно здійснювати в циклі. При цьому всі елементи складаються послідовно один за іншим. На ГП подібна процедура може бути виконана в M незалежних потоках, кожен із яких отримує один елемент результуючого вектора. Отже, всі елементи суми можуть бути отримані одночасно і час подібної операції відповідає часу виконання одного складання на відміну від М складань на ЦП.

ГП особливо добре підходять для розв'язання обчислювальних завдань, які можна розкласти на паралельні обчислення на масиві даних - одні й ті самі операції проводяться на великій кількості елементів у пам'яті. При цьому, чим більше відношення кількості арифметичних операцій до кількості операцій читання з пам'яті і запису в пам'ять, тим ефективніша буде реалізація алгоритму на ГП. Оскільки виконувані над різними елементами масиву даних операції ідентичні, складний контролер не потрібен. За великої кількості виконуваних операцій затримка в доступі до даних може бути нівельована виконанням операцій над уже доступними елементами масиву, що дає змогу істотно скоротити розмір кеш-пам'яті.

Рис. 2.1 Архітектура центрального (ліворуч) і графічного (праворуч) процесорів

На ЦП значна частина поверхні мікрочіпа зайнята кеш-пямяттю і контролером. ЦП також оснащений змінною кількістю пам'яті DRAM (Dynamic Random Access Memory). Арифметичні логічні пристрої (АЛП) на ГП, згруповані в мультипроцесори, кожен із який має свого контролера і кеш-пам'ять. Це дає змогу істотно збільшити щільність обчислювальних одиниць на мікрочіпі. ГП має окрему пам'ять DRAM, яку зазвичай називають глобальною пам'яттю пристрою. Взаємодія між ЦП і ГП здійснюється через шину PCI Express.

Платформа CUDA

Технологія CUDA, запропонована компанією NVIDIA в листопаді 2006 року, є програмною платформою для паралельних обчислень загального характеру. Вона дає змогу використовувати можливості графічного процесора в будь-яких багатопаралельних обчислювальних задачах, вирішуючи їх набагато ефективніше, ніж це можливо із застосуванням ЦП. Крім того, CUDA є розширенням мови C і поставляється з повним інструментарієм, необхідним розробнику програмного забезпечення.

Сучасні ЦП і ГП багатоядерні, тобто широко розповсюджені обчислювальні пристрої паралельні. До того ж, у той час як виробники ЦП вже давно зіткнулися з технологічними обмеженнями, які перешкоджають зростанню продуктивності одного обчислювального ядра, продуктивність паралельних пристроїв пропорційна кількості обчислювальних ядер і підлягає закону Мура. Основною складністю є розроблення програмного забезпечення, яке могло б мати пропорційний приріст продуктивності при виконанні на нових багатопаралельних процесорах із постійно зростаючим числом обчислювальних ядер. Це завдання вже було успішно розв'язане в додатках, що працюють з тривимірною графікою.

Програмна платформа CUDA побудована довкруж ієрархії обчислювальних потоків, які розпадаються на групи. На апаратному рівні це відповідає поділу ГП на мультипроцесори: що більше груп потоків, то більше мультипроцесорів може бути задіяно під час виконання програми. За наявності достатнього паралелізму в задачі, швидкість її виконання буде пропорційна обчислювальним можливостям ГП. Крім того, платформа CUDA досить проста для засвоєння людьми, знайомими із стандартними мовами програмування, такими як C і CUDA.

Програмна модель

Обчислювальні ядра

Діалект CUDA для мови C дає змогу програмісту визначати спеціальні функції, названі ядрами, які при виклику можна виконати паралельно N разів у N обчислювальних потоках. Ядро визначають за допомогою модифікатора __global__, а кількість бажаних потоків задають за допомогою синтаксису <<<...>>> - конфігурації виконання ядра. Кожен потік, що виконує ядро, забезпечений унікальним ідентифікатором потоку, який доступний усередині ядра через вбудовану змінну threadIdx. Наприклад, наступна програма виконує додавання двох векторів a і b розміру N і зберігає результат у векторі c:

1 __global__ void vecAdd ( float * a, float * b, float * c){

2 int i = threadIdx .x;

3 c[i] = a[i] + b[i];

4 }

5

6 int main (){

7 ...

8 vecAdd <<<1, N>>>(a, b, c);

9 ...

10 }

Ієрархія потоків

Індекс потоку threadIdx - тривимірний вектор, а потоки можна ідентифікувати, використовуючи три-, дво- або одновимірний індекс, формуючи три-, дво або одномірні блоки потоків. Це дає змогу розподіляти обчислення за елементами обсягу, матриці або вектора, вибираючи відповідну розмірність індексу потоку. Так як всі потоки блоку виконують на одному мультипроцесорі, доступні їм ресурси (кеш-пам'ять, регістри) обмежені. Тому кількість потоків у блоці обмежена й залежить від покоління ГП (на графічних картах з архітектурою Fermi це обмеження - 1024 потоку на блок).

Зазвичай ядро виконують у множині однакових блоків потоків. Загальну кількість потоків у такому випадку визначають як добуток кількості блоків на кількість потоків у кожному з них. Як і потоки всередині блоків, блоки об'єднуються в одно-, дво- або тривимірну сітку блоків. Блоки в ядрі можна ідентифікувати з використанням змінної blockIdx, а їхній розмір - через змінну blockDim. Отже, загальний індекс потоку на сітці в одновимірному випадку можна визначити як blockIdx * blockDim + threadIdx. Кількість потоків усередині блоку і кількість блоків може здаватися змінними типу int або dim3 всередині дужок <<<...>>> при запуску ядра. Попередній приклад із розбиттям потоків на блоки буде записаний як:

1 __global__ void vecAdd ( float * a, float * b, float * c){

2 int i = blockIdx * blockDim + threadIdx ;

3 c[i] = a[i] + b[i];

4 }

5

6 int main (){

7 ...

8 int threadsPerBlock = 128;

9 int numberOfBlocks = N /128;

10 vecAdd <<< numberOfBlocks , threadsPerBlock >>>(a, b, c);

11 ...

12 }

Блоки потоків повинні виконуватися незалежно, а порядок їх виконання заздалегідь невідомий. До того ж, у CUDA немає можливості синхронізувати виконання потоків із різних блоків. Усередині блоку потоки можуть передавати дані через пам'ять, що розділяється, а для їх бар'єрної синхронізації наявна вбудована функція __syncthreads(). Якщо при виконанні команд потік натрапляє на виклик функції __syncthreads(), його виконання зупиняється до тих пір, поки всі інші потоки всередині того самого блоку не досягнуть виклику цієї функції.

Ієрархія пам'яті

Дуже важливою частиною будь-якої програми, що використовує CUDA, є управління пам'яттю. Частину такої програми виконують на ЦП, а частину - на ГП, і ці пристрої мають свою ділянку пам'яті. Для виділення пам'яті на ЦП можна використовувати стандартні процедури мови C, такі як calloc(...). Для виділення пам'яті на пристрої використовують процедуру cudaMalloc(...) з діалекту CUDA для C. Для копіювання даних між пам'яттю ЦП і ГП використовують процедуру cudaMemcpy(...). Повернемося до прикладу зі складанням векторів. Припустимо, що вектори заповнюються на ЦП, потім дані копіюються на ГП, вектори складаються, результат повертається в пам'ять ЦП і виводиться на екран:

1 int main(){

2 // Allocate memory on CPU ( host )

3 float* h_a = (float*) calloc (N, sizeof( float ));

4 float* h_b = (float*) calloc (N, sizeof( float ));

5 float* h_c = (float*) calloc (N, sizeof( float ));

6 ...

7 // Fill arrays h_a and h_b with data

8 ...

9 // Allocate memory on GPU ( device )

10 float* d_a; float* d_b; float* d_c;

11 cudaMalloc ((void**)&d_a , N* sizeof( float ));

12 cudaMalloc ((void**)&d_b , N* sizeof( float ));

13 cudaMalloc ((void**)&d_c , N* sizeof( float ));

14 // Copy data from CPU to GPU memory

15 cudaMemcpy (d_a , h_a , N* sizeof( float ), cudaMemcpyHostToDevice );

16 cudaMemcpy (d_b , h_b , N* sizeof ( float ), cudaMemcpyHostToDevice );

17 // GPU kernel launch

18 int threadsPerBlock = 128;

19 int numberOfBlocks = N /128;

20 vecAdd <<< numberOfBlocks , threadsPerBlock >>>(d_a , d_b , d_c );

21 // Copy the result back to CPU memory

22 cudaMemcpy (h_c , d_c , N* sizeof ( float ), cudaMemcpyDeviceToHost );

23 // Print the result on the screen

24 for(i = 0; i < N; i ++){

25 printf ("%f\n", h_c[i]);

26 }

27 // Deallocate the memory

28 free (h_a ); free (h_b ); free (h_c );

29 cudaFree ( d_a ); cudaFree (d_b ); cudaFree (d_c );

30 }

У наведених процедурах указівник на дані, що знаходяться в пам'яті ЦП, позначені префіксом h_ (від англ. host), а на дані в пам'яті ГП - префіксом d_ (від англ. device). Важливо розуміти, що при прямому зверненні до даних за вказівником на пам'ять ГП з частини програми, що виконується на ЦП, отримаємо помилку. Тут можна навести аналогію з поштовими адресами: для тих, хто знаходиться в одному місті, адреса в іншому не має сенсу.

Пам'ять ГП доступна з коду ЦП за допомогою процедур діалекту CUDA. Отже, процес виконання програми управляється ЦП, а ГП відведена роль співпроцесора. Важливо відзначити, що копіювання даних із пам'яті ЦП в пам'ять ГП здійснюють через інтерфейс PCI-Express, який характеризується низькою пропускною здатністю в порівнянні з шиною доступу до пам'яті. Велика кількість операцій копіювання даних може істотно уповільнити виконання програми, тому для забезпечення високої продуктивності необхідно мінімізувати перенесення даних.

При написанні програм, що працюють на ЦП, порядок і спосіб звернення до пам'яті не відіграє ключової ролі в продуктивності програми. На ГП порядок і спосіб звернення до пам'яті є одним із головних факторів, що впливають на продуктивність. Відбувається це через малий обсяг кеш-пам'яті, а також через те, що безліч потоків може звертатися до пам'яті одночасно. Тому в діалекті CUDA програміст має більше можливостей контролювати потоки даних із однієї області пам'яті ГП в іншу.

На ГП є кілька різних областей пам'яті (таблиця 2.1). По-перше, кожен мультипроцесор забезпечений набором регістрів і кешем. Регістри розподіляються між потоками, а кеш-пам'ять може бути використана або для кешування елементів глобальної пам'яті, або як розділена пам'ять, доступна всім потокам, що їх виконують на цьому мультипроцесорі. Ці області пам'яті розташовуються на чіпі ГП. У пам'яті, розміщеній на платі графічної карти, розташовується глобальна пам'ять, доступ до якої може здійснюватися як з будь-якого потоку, що виконується на ГП, так і з ЦП за допомогою спеціальних команд CUDA (напр. CudaMemcpy(...)). На сучасних ГП доступи до глобальної пам'яті кешуються. Крім того, існує область постійної пам'яті, яка завжди кешується. Читання з глобальної пам'яті також можна здійснювати за допомогою текстурних посилань, які кешуються просторово (разом із зчитуваним елементом у кеш заносяться й елементи, які його оточують в пам'яті).

Таблиця 2.1 Типи пам'яті в CUDA.

Тип	Тип (англ.)	Розміщення	Кешування	Доступ	Видимість	Час життя
Регістр	Register	На чіпі		читання / запис	З потоку	Потік
Спільна	Shared	На чіпі		читання / запис	З блоку	Блок
Глобальна	Global	На платі	Так/Ні	читання / запис	Глобальна	Додаток
Постійна	Constant	На платі	Так	читання / запис	Глобальна	Додаток
Текстурна	Texture	На платі	Так	читання / запис	Глобальна	Додаток

Ієрархічна структура даних на ГП відкриває широкі можливості для оптимізації програмного коду. По-перше, кількість регістрів на мультипроцесорі обмежена, тому, щоб максимізувати кількість потоків, які одночасно виконуються на мультипроцесорі, програміст повинен контролювати кількість задіяних локальних величин усередині ядра. По-друге, спільна пам'ять може бути використана для явного кешування даних. По-третє, кількість звернень до глобальної пам'яті повинно бути мінімізовано або розподілено по програмному коду так, щоб усі потоки могли виконувати арифметичні операції в процесі очікування даних. Також через наявність апаратного кеша звернення до глобальної пам'яті має бути локалізоване - сусідні потоки повинні звертатися до сусідніх осередків пам'яті. Доступ до глобальної пам'яті пристрою повинен здійснюватися спряжено: наступні один за одним потоки повинні використовувати наступні один за одним елементи масиву з глобальної пам'яті пристрою. Це дає змогу пристрою склеювати декілька звернень до пам'яті в один. Вирівнювання потоків щодо масивів з даними має багато тонкощів і залежить від архітектури (покоління) ГП. Для детальнішого ознайомлення зі способами доступу до пам'яті і її вирівнювання читачеві слід вивчити інструкцію програміста, що її постачають з CUDA. Істотним фактором, який може обмежити продуктивність програми, є швидкість передачі даних від ЦП до ГП, яка обмежена шиною PCI-Express. Програміст повинен, за можливості, уникати копіювання даних між ЦП і ГП. Цього можна досягти перенесенням усього програмного коду на ГП, використовуючи ЦП лише для управління потоком обчислень і виводу результатів.

Діалект CUDA для C

Модифікатори функцій

Модифікатори функцій повідомляють компілятору, чи може задана функція виконуватися на ГП або на ЦП і чи можна її викликати з ЦП або з ГП. Модифікатор __device__ засвідчує те, що функція може виконуватися тільки на пристрої і викликається тільки з пристрою. Модифікатор __global__ ставиться перед функцією, яка є обчислювальним ядром (викликається з ЦП і виконується на ГП). Під час викликання даної функції необхідно вказувати конфігурацію її запуску (кількість блоків потоків і кількість потоків у кожному з них за допомогою синтаксису <<<...>>>). Функції з даним модифікатором повинні повертати void. Модифікатор __host__ повідомляє про те, що дана функція може викликатися тільки з ЦП, де її і виконують. Відсутність усіх трьох модифікаторів еквівалентна наявності модифікатора __host__. Комбінацію модифікаторів __host__ і __device__ можна використовувати одночасно. У цьому випадку код компілюватиметься двічі - для ЦП і для ГП.

Вбудовані векторні типи

У CUDA використовують вбудовані векторні типи для всіх типів цілочисельних змінних і змінних із плаваючою крапкою. Усі ці типи являють собою структури, перша, друга, третя і четверта компоненти яких доступні через поля х, у, z і w відповідно. Наприклад, змінна типу float3 матиме три поля - х, у і z. Для всіх векторних типів визначений конструктор, ім'я якого відповідає типу: make_<name>. Наприклад, для того щоб створити вектор з компонентами (х, у, z) можна викликати конструктор:

1 float x, y, z;

2 float3 r = make_float3 (x, y, z);

Математичні функції

У процедурах, що виконуваних на ГП, крім стандартних математичних функцій (наприклад sinf(alpha)), доступні їхні більш швидкі, але менш точні вбудовані аналоги. Останні позначаються подвійним підкресленням перед назвою функції (наприклад __sinf(alpha)). Під час використання цих функцій потрібно впевнитися, що надана ними точність достатня для виконання тих чи тих операцій. У деяких особливих випадках результати, одержувані з використанням вбудованих функцій, можуть істотно відрізнятися від бажаних, тому правильною буде заміна функцій на вбудовані в тих місцях, де це не вплине на результат. Для того, щоб компілятор використовував вбудовані функції замість стандартних, можна при компіляції програми вказати опцію -use_fast_math.

Будова ГП

Мультипроцесори

Ядро ГП компанії NVIDIA являє собою набір багатопотокових мультипроцесорів. Коли програма, написана на CUDA, викликає обчислювальне ядро, блоки потоків нумеруються і розподіляються по мультипроцесорах пристрою. Потоки з одного блоку виконуються паралельно на одному мультипроцесорі, якому можуть бути призначені відразу кілька блоків потоків. Як тільки один із блоків закінчує виконання обчислювальних процедур, мультипроцесор починає виконання наступного блоку потоків. Мультипроцесори можуть виконувати сотні потоків одночасно. Для того, щоб керувати такою великою кількістю потоків, мультипроцесори використовують унікальну архітектуру: одиночний потік команд, множинні потоки виконання (ОКМП; від англ. Single Instruction, Multiple Thread, SIMT). Інструкції виконуються паралельно з використанням власного паралелізму інструкцій, або послідовно з використанням паралелізму безлічі потоків.

Архітектура ОКМП

Мультипроцесор створює групи по 32 потоки в кожній (англ. warp). Усі потоки однієї такої групи починають виконуватися одночасно з використанням загального набору команд. Оскільки у кожної групи наявний свій набір регістрів, групи виконуються незалежно одна від одної. Якщо мультипроцесор виконує більше одного блоку потоків, то всі блоки діляться на групи по 32 потоки і далі виконуються в цих групах. Позаяк група потоків виконує той самий набір команд, для досягнення максимальної продуктивності слід уникати розгалужень усередині групи. У разі виникнення останнього, деякі потоки будуть змушені простоювати, поки інші виконують процедури, передбачені їх відгалуженням.

Основною відмінністю архітектури ОКМП від архітектури ОКМД (одиничний потік команд, множинний потік даних; від англ. Single Instruction, Multiple Data, SIMD) є те, що програмісту надана можливість контролювати не тільки загальні для всіх даних команди, а й виконувати кожен потік окремо. Таким чином перед програмістом відкриваються великі можливості для оптимізації коду: істотного приросту продуктивності можна досягти, якщо розподіляти потоки групами із 32 штук із урахуванням виникнення в них розгалужень.

3. Моделювання кисню в рівновазі Постановка завдання

Для того, щоб детально розібрати способи реалізації численних методів молекулярної динаміки на графічних процесорах, розглянемо задачу моделювання кисню в рівноважному стані. При цьому ковалентний зв'язок між атомами кисню всередині молекули O2 будемо описувати за допомогою гармонійного потенціалу, а атоми з різних молекул будуть взаємодіяти через потенціал Леннарда-Джонса (взаємодії Ван-дер-Ваальса). Оскільки атоми кисню в молекулі O2 не заряджені (частковий заряд дорівнює нулю), то електростатичні взаємодії в цій системі відсутні. Для збереження обсягу системи ми будемо використовувати періодичні граничні умови, застосування яких позначиться тільки на формулах розрахунку відстані між атомами: через відсутність далекодіючого електростатичного потенціалу, можливий прямий розрахунок невалентних взаємодій із застосуванням радіуса обрізки. Отже, потенційна енергія системи буде описана за допомогою такої функції координат атомів:

,

Де , а значення параметрів перераховані в таблиця 3.2.

Важливим моментом є вибір системи одиниць. Очевидно, що зручнішими одиницями виміру є величини, які так чи так характеризують молекулярну систему. У розгляданому прикладі ми будемо використовувати одиниці вимірювання МД, ідентичні тим, які використовують у програмі Gromacs:

Таблиця 3.1 Система одиниць МД.

Величина	Одиниці виміру
Довжина	нм = 10-9 м
Маса	г/моль = 1.6605402 Ч 1027 кг
Час	пс = 10-12 с
Температура	К

Зазначимо, що в розгляданій системі одиниць енергія вимірюється в (г/моль)нм2/пс2 = 103(кг/моль)м2/с2 = кДж/моль, а сила - в (кДж/моль)/нм. Постійна Больцмана kB = 8.314462 Ч 10-3 кДж/(моль·К). Параметри взаємодії між атомами кисню в цій системі одиниць наведені в таблиця 3.2.

Таблиця 3.2 Параметри взаємодій для атомів молекули кисню O2.

Параметр	Значення	Опис
m	15.999 г/моль	Молярна маса атома кисню
b0	0.121 нм	Рівноважна відстань між атомами в молекулі O2
Ks	518816 кДж/моль·нм2	Постійна пружності для ковалентного зв'язку O=O
у	0.302905564168 нм	Розташування мінімуму енергії потенціалу Леннарда-Джонса
е	0.50208 кДж/моль	Глибина мінімуму потенціалу Леннарда-Джонса

Далі ми детально у визначеному порядку розберемо всі обчислювальні процедури для моделювання розгляданої системи. У процесі реалізації ми буде використовувати стандартні прості типи CUDA (float2, float3 ...), а також деякі зумовлені операції з цими типами (наприклад, покомпонентна різниця двох векторів буде визначена через оператор "-" (мінус) для двох змінних типу float3).

Початкова конформація системи

Випадкове розташування атомів

Розглянемо кубічну періодичну систему, у якій межі задають як 0 < {x, y, z} < L. Нехай кількість атомів дорівнює N = 2k, де k - ціле (k - число молекул O2), і нехай кожен атом з парним індексом i (i = 0, 2, 4 ... N-2) пов'язаний валентно з наступним атомом з індексом j (j = 1, 3, 5... N-1). Для розміщення наступного парного атома i потрібно задати випадкові координати 0 < {xi, yi, zi} < L таким чином, що сфери Ван-дер-Ваальса цього атома не перетинаються зі сферами Ван-дер-Ваальса атомів, уже доданих у систему (електронні хмари атомів не перетинаються). При цьому слід ураховувати, що в системі з періодичними граничними умовами відстань між атомами i і j задається як мінімальна відстань від атома i до атома j і його періодичних відображень, тобто:

,

де квадратні дужки позначають округлення. Програмна реалізація даного алгоритму:

1 inline __host__ __device__ float3 getVector(float3 ri, float3 rj, float L){

2 float3 dr = rj -ri;

3 dr -= rint (dr/L)*L;

4 return dr;

5 }

Крім модифікатора вбудовування процедури (inline) тут використані модифікатори з діалекту CUDA - __host__ і __device__, завдяки яким процедура буде скомпільована для виконання як на ЦП так і на ГП і її getVector(...) можна буде використовувати в обох частинах коду. Відстань між двома атомами в умовах періодичних граничних умов визначається довжиною цього вектора, тобто:

1 inline __host__ __device__ float getDistance(float3 ri, float3 rj, float L){

2 float3 dr = getVector (ri , rj , L);

3 return len(dr );

4 }

Де len(float3 r) = sqrtf(r.x*r.x + r.y*r.y + r.z*r.z);. Крім того, нам буде потрібна процедура переміщення атомів в область 0 < {x, y, z} < L через періодичні граничні умови, яку можна задати як:

1 float3 transferPBC ( float3 r){

2 if(r.x > L){

3 r.x -= L;

4 } else

5 if(r.x < 0){

6 r.x += L;

7 }

8 ...

9 return r;

10 }

Розташування валентно зв'язаного з атомом i атома j при цьому можна задати випадковим кутом орієнтації зв'язку O=O, розташувавши атом j на рівноважній відстані від атома i. Як генератор випадкових чисел ми будемо використовувати ran2. молекула процесор атом валентний

1 void placeAtoms (){

2 float3 r1 , r2 , dr;

3 int i, j;

4 for(i = 0; i < N/2; i ++){

5 float mindist = 0.0;

6 while ( mindist < 2.0* sigma ){

7 r1.x = L* ran2 (& rseed ); r1.y = L* ran2 (& rseed ); r1.z = L* ran2 (& rseed );

8 dr.x = ran2 (& rseed ); dr.y = ran2 (& rseed ); dr.z = ran2 (& rseed );

9 float norm = len(dr );

10 dr /= norm ; 11 dr *= b0;

12 r2 = r1 + dr;

13 r2 = transferPBC (r2 , L);

14 mindist = FLT_MAX ;

15 for(j = 0; j < i; j ++){

16 float dist = getDistance (r1 , r[2*j], L);

17 if( dist < mindist ){

18 mindist = dist ;

19 }

20 }

21 }

22 r [2*i] = r1;

23 r [2*i+1] = r2;

24 }

25 }

У даній обчислювальної процедурі першому атому з молекули кисню присвоюються випадкові координати (r1), а другому (r2) - координати, вилучені від першого на рівноважну довжину ковалентного зв'язку O=O. При цьому для кожної пари атомів випадкове число кидається до тих пір, поки мінімальна відстань між одним із атомів пари і атомами, вже доданими в систему, не буде перевищувати двох рівноважних відстаней потенціалу Леннарда-Джонса (2.0*sigma). Початкове розташування атомів показано на рис. 3.1.

Рис. 3.1 Випадкове розташування 50 молекул кисню в кубі 4 Ч 4 Ч 4 нм. Чорними лініями показані граничні умови. Атоми показані в уявленні CPK, колір відповідає порядковому номеру атома.



Розподіл швидкостей

Початкові значення швидкостей визначаються температурою системи. Для того, щоб їх задати, ми будемо використовувати випадкові числа, розподілені по Гаусу. Зазвичай, генератори випадкових чисел повертають рівномірно розподілені випадкові числа, які потім можна перетворити у випадкові числа, розподілені нормально за допомогою перетворення Боксу-Мюллера, яке задають таким чином. Нехай u1 і u2 - незалежні випадкові величини, рівномірно розподілені на відрізку [-1, 1]. Обчислимо s = u12 + u22. Якщо виявиться, що s > 1 або s = 0, то значення u1 і u2 слід згенерувати заново. Як тільки буде виконана умова 0 < s ? 1, можна розрахувати дві незалежні випадкові величини n1 і n2, що задовольняють нормальному розподілу з математичним очікуванням 0 і дисперсією 1. Для цього необхідно використовувати такі рівняння:

,

,

Оскільки зазвичай випадкові величини потрібні одна за одною, а перетворення Боксу-Мюллера створює відразу дві нормально розподілені величини (з двох рівномірно розподілених), має сенс зберігати одну з отриманих нормально розподілених величин, щоб повернути її за повторного звернення:

1 inline double gasdev(int *idum) {

2 static int iset = 0;

3 static double gset ;

4 double fac , rsq , v1 , v2;

5 if( iset == 0){

6 do {

7 v1 = 2.0*ran2(idum ) - 1.0;

8 v2 = 2.0*ran2(idum ) - 1.0;

9 rsq = v1*v1+v2*v2;

10 } while (rsq >= 1.0 || rsq == 0.0);

11 fac = sqrt(-2.0*log(rsq)/rsq);

12 gset = v1*fac;

13 iset = 1;

14 return v2*fac;

15 } else {

16 iset = 0;

17 return gset ;

18 }

19 }

Згідно з розподілом Максвелла, компоненти швидкостей частинок розподілені нормально, з математичним очікуванням 0 і дисперсією . Отже, щоб отримати компоненту швидкості частинки масою m за температури T, необхідно отримати випадкову величину n (Рівняння 26) і помножити її на . У програмі це можна зробити за допомогою такої простої процедури:

1 void assignVelocities(){

2 for(i = 0; i < N; i++){

3 double var = sqrt(kB*T/m);

4 v[i].x = var*gasdev(&seed);

5 v[i].y = var*gasdev(&seed);

6 v[i].z = var*gasdev(&seed);

7 }

8 }

Збереження стану системи

Розглянута нами система може бути описана 6*N узагальненими координатами - трьома компонентами координат ({x, y, z}) і трьома компонентами швидкостей ({vx, vy, vz}) для кожної з N частинок. Значення цих величин необхідно зберігати для того, щоб надалі можна було проаналізувати поведінку системи. Крім того, формат даних повинен бути загальноприйнятим, щоб його можна було використовувати в стандартних програмах візуалізації та аналізу, таких як VMD. У разі збереження траєкторій, коли у файлі міститься багато станів системи, величезну роль відіграє також і розмір цього файлу. Тому зазвичай для цих цілей використовують бінарні формати (наприклад - .dcd в NAMD і CHARMM або .trr в Gromacs). У межах нашого прикладу ми зупинимося на текстовому форматі .xyz. Звичайно, можна було б використовувати й описаний раніше формат .pdb, але він охоплює багато інформації, яка в нашому випадку не обов'язкова. Формат .xyz набагато простіший: він містить один рядок із цілим числом, рівним кількості атомів в молекулярній системі (N), один рядок коментарів і N рядків з розділеними пробілами або табуляцією іменами атомів і їх координатами (в ангстремах). Наприклад файл формату .xyz для координат атомів однієї молекули кисню має такий вигляд:

...