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

Oxygen molecule (in Angtroms) O -0.605000 0.000000 0.000000 O 0.605000 0.000000 0.000000

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

1 void writeXYZ ( const char * filename , float3 * data , int N){

2 FILE * file = fopen ( filename , "w");

3 fprintf (file , "%d\n", N);

4 fprintf (file , " Oxygen \n");

5 int i;

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

7 fprintf (file , "O\t%f\t%f\t%f\n",

8 10.0* data [i].x, 10.0* data [i].y, 10.0* data [i].z);

9 }

10 fclose ( file );

11 }

Розглядана процедура побудована таким чином, що вона може зберігати не тільки координати атомів, але і їхні швидкості. Крім того, проста зміна індикатора відкриття файлу на "a" (append - додати) дає змогу зберігати файли, що містять більше одного набору координат. Далі ми використовуємо цю процедуру в окремому класі, який будемо використовувати для читання і запису поточного стану системи.

Об'єднання процедур у працюючу програму

Для того, щоб отримати працюючу програму, необхідно об'єднати всі наведені вище процедури в один файл .cpp. Хоча ми використовували простий тип float3, у поточному варіанті програму можна компілювати звичайним компілятором C++ (наприклад g++), додавши заголовок із описом простих типів CUDA (vector_types.h). Крім перерахованих процедур, нам буде потрібен генератор псевдо-випадкових чисел Ran2 ([31] і Додаток 1). Також слід задати необхідні оператори для векторних типів float3, операцію знаходження довжини вектора (len(float3 a)) і його покомпонентного округлення (rint(float3 a)).

1 int main(int argc , char* argv[]){

2 r = (float3*)calloc(N, sizeof(float3));

3 v = (float3*)calloc(N, sizeof(float3));

4 placeAtoms();

5 assignVelocities();

6 writeXYZ("O2.coor.xyz", r, N);

7 writeXYZ ("O2.vel.xyz", v, N);

8 return 1;

9 }

У методі main(...) необхідно виділити пам'ять під координати і швидкості та викликати процедури присвоювання координат (placeAtoms()) і швидкостей (assignVelocities()) атомів, після чого зберегти ці значення у файли .xyz (writeXYZ(...)). Отриманий при запуску програми файл координат можна завантажити в програму VMD і переконатися, що атоми не виходять за межі куба LЧLЧL і розташовані досить далеко один від одного (рис. 3.1).

Атоми, що вільно рухаються

Інтегрування рівнянь руху

У наведеному прикладі ми будемо використовувати рівняння руху Ньютона:

,

,

Одним із найбільш популярних алгоритмів інтегрування рівнянь руху в молекулярній динаміці є алгоритм "стрибок жаби" (англ. "Leap-Frog"), який можна записати так:

,

,

Уважаючи, що швидкості частинок задані між кроками по часу, для однієї частинки даний алгоритм в програмному коді можна записати таким чином:

1 v += f*( dt/m);

2 r += v*dt;

Для системи в N атомів необхідно інтегрувати N рівнянь руху:

1 void integrate ( float3 * r, float3 * v, float3 * f,

2 float m, float dt , int N){

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

4 v[i] += f[i]*( dt/m);

5 r[i] += v[i]* dt;

6 f[i] = make_float3 (0.0f, 0.0f, 0.0f);

7 }

8 }

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

1__global__ void integrate_kernel ( float3 * d_r , float3 *

2 d_v , float3 * d_f , float m, float dt , int N){

3 int d_i = blockIdx .x* blockDim .x + threadIdx .x;

4 if(d_i < N){

5 float3 r = d_r[d_i ];

6 float3 v = d_v[d_i ];

7 float3 f = d_f[d_i ];

8 v += f*( dt/m);

9 r += v*dt;

10 d_r [d_i] = r;

11 d_v [d_i] = v;

12 d_f [d_i] = make_float3 (0.0f, 0.0f, 0.0f);

13 }

14 }

Дане обчислювальне ядро працює з масивами, що мають префікс "d_" (від англ. "Device"). Ми спеціально використовуємо такий префікс, щоб відокремити покажчики на змінні в глобальній пам'яті ГП від покажчиків на пам'ять ЦП. Далі ми введемо ще одне правило: якщо в масиву на ЦП є копія на ГП, ми будемо додавати до такого масиву префікс "h_" (від англ. "Host"). У нашому випадку цей префікс повинні отримати масиви r, v і f. У той час, як пам'ять під масиви з префіксом "h_" виділяють за допомогою стандартної процедури calloc(...), для виділення пам'яті на ГП необхідно використовувати процедуру cudaMalloc(...) з діалекту CUDA. Маючи два масиви "h_r" і "d_r" однакового типу (float3) і розміру (N), можна копіювати дані з ЦП на ГП і назад за допомогою процедури cudaMemcpy(...).

Індексом d_i в наведеному лістингу присвоюється значення глобального індексу потоку, за яким і визначається той елемент масиву, з яким працює даний потік. Перевірка того, що індекс не виходить за межі масиву (d_i < N), робиться тому, що кількість потоків на ГП кратна розміру блоку потоків (тобто може бути більше кількості атомів N). Далі змінні, з якими даний потік працює, зчитують у локальну пам'ять потоку, з ними проводять всі необхідні операції, і результат зберігають у глобальну пам'ять. Такий підхід зводить до мінімуму кількість звернень до глобальної пам'яті, які можуть істотно сповільнити виконання програми.

Моделювання вільного руху частинок

Наведені вище процедури роблять тільки один крок інтегрування, у той час як у молекулярній динаміці кількість цих кроків досягає мільйонів або навіть мільярдів. Наприклад, якщо інтегрувати з кроків в 1 фс (10-15 с), то для досягнення 1 нс (10-9 с) буде потрібно 106 кроків. Тому процедуру інтегрування необхідно запускати в циклі, кількість ітерацій якого визначатиме часовий інтервал, який досягається в моделюванні. Вибір кроку інтегрування визначається найбільш швидкими рухами всередині системи. У нашому випадку такими будуть вібрації ковалентних зв'язків O=O з періодом коливань близько 20 фс.

Для того, щоб стежити за динамікою системи, необхідно зберігати координати атомів у процесі інтегрування рівнянь руху. Збереження координат на кожному кроці не тільки потребує великої кількості місця на жорсткому диску, але й істотно уповільнить процес обчислення. Тому координати зазвичай зберігають кожні 1000-10000 кроків. Уведемо змінну output_freq, яка визначатиме частоту збереження координат. По закінченні output_freq, припинимо процес інтегрування і збережемо поточні координати, попередньо перенісши їх через кордони періодичних граничних умов. Після збереження координат програма повинна повертатися до інтегрування рівнянь руху. Цей процес повторюють до тих пір, поки не буде досягнута бажана кількість кроків інтегрування steps_count:

1 void computeCPU (){

2 int s, i;

3 for( step = 0; step <= steps_count ; step += output_freq ){

4 // Internal cycle - no interaptions

5 for (s = 0; s < output_freq ; s ++){

6 integrate (h_r , h_v , h_f , m, dt , N);

7 }

8 // Transfer coordinates through periodic boundary

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

10 h_r[i] = transferPBC (h_r [i], L);

11 }

12 // Save coordinates

13 appendXYZ (" O2_CPU . coor .xyz", h_r , N);

14 }

15 }

У цій обчислювальній процедурі ми використовували процедуру appendXYZ(...), яка відрізняється від writeXYZ(...) тільки модифікатором доступу до файлу (замість "w" використовується "a"), тобто замість FILE* file = fopen(filename, "w") використовується FILE* file = fopen(filename, "a"). Це дає змогу дописувати координати в файл, зберігаючи динамічну поведінку системи.

Реалізація аналогічної процедури при інтегруванні рівнянь руху на ГП має такий вигляд:

1 void computeGPU (){

2 int s, i;

3 for( step = 0; step <= steps_count ; step += output_freq ){

4 // Internal cycle - no interaptions

5 for (s = 0; s < output_freq ; s ++){

6 integrate_kernel <<<(N -1)/ BLOCK_SIZE + 1,

7 BLOCK_SIZE >>>(d_r , d_v , d_f , m, dt , N);

8 }

9 // Copy coordinates from GPU

10 cudaMemcpy (h_r , d_r , N* sizeof ( float3 ) , cudaMemcpyDeviceToHost );

11 // Transfer coordinates through periodic boundary

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

13 h_r[i] = transferPBC (h_r [i], L);

14 }

15 // Save coordinates

16 appendXYZ (" O2_GPU . coor .xyz", h_r , N);

17 }

18 }

Тут процедура integrate(...) замінена на ядро integrate_kernel<<<...>>>(...). Запис <<< N/BLOCK_SIZE + 1, BLOCK_SIZE>>> означає, що ядро буде запущено в N/BLOCK_SIZE + 1 блоках по BLOCK_SIZE потоків. Загальне число потоків (T) буде визначати розмір системи (N), і буде більшим або дорівнювати цьому розміру (T ? N). Запуск більшого числа потоків необхідний тому, що потоки на ГП виконуються в блоках, що на апаратному рівні відповідає поділу потоків між мультипроцесорами. Ядру передається той же набір аргументів, що і процедурі інтегрування на ЦП, але передаються покажчики на масиви в пам'яті ГП. Отже, всі зміни з даними будуть відбуватися виключно на ГП і, перед тим як зберегти результати у файл, необхідно скопіювати координати в пам'ять ЦП.

Об'єднання процедур у працюючу програму

Для запуску моделювання вільного руху незв'язаних атомів кисню необхідно додати в методі main(...) виклик процедури computeCPU(...). У разі моделювання на ГП пам'ять також виділяється і для змінних пристроїв (d_r, d_v і d_f), у які копіюються дані з пам'яті ЦП:

1 cudaMalloc (( void **)& d_r , N* sizeof ( float3 ));

2 cudaMalloc (( void **)& d_v , N* sizeof ( float3 ));

3 cudaMalloc (( void **)& d_f , N* sizeof ( float3 ));

4 cudaMemcpy (d_r , h_r , N* sizeof ( float3 ), cudaMemcpyHostToDevice );

5 cudaMemcpy (d_v , h_v , N* sizeof ( float3 ), cudaMemcpyHostToDevice );

6 cudaMemcpy (d_f , h_f , N* sizeof ( float3 ), cudaMemcpyHostToDevice );

Після цього повинна бути викликана процедура computeGPU(...). Ми навмисне копіюємо масив із силами, що діють на частинки, позаяк процедура cudaMalloc(...) на відміну від процедури calloc(...) не обнуляє значення масиву. У даній програмі був використаний діалект CUDA, через що скомпілювати її за допомогою стандартного компілятора не вийде і необхідно використовувати компілятор nvcc з CUDA Toolkit. Цей компілятор дає змогу також використовувати бібліотеку векторних операторів cutil_math.h з CUDA SDK. Після виконання програми в робочій директорії повинен з'явитися файл .xyz з динамікою системи. Цей файл можна завантажити в програму VMD і простежити рух частинок. Так як атоми кисню в молекулі перебувають на близькій відстані, VMD здогадається, що вони пов'язані валентно і в деяких поданнях будуть промальовані зв'язки між цими атомами. Не зв'язані атоми будуть рухатися незалежно один від одного, що візуально відіб'ється в нефізичному розтягуванні зв'язків. Це буде виправлено в наступному розділі, де ми додамо потенціал ковалентної взаємодії.

Потенціал валентних взаємодій

Формули для розрахунку сил

Реалізація будь-якого потенціалу починається з його формального запису. Зокрема, щоб отримати формулу для розрахунку сили, що діє на атом, необхідно обчислити похідну потенціалу за координатами цього атома. У разі молекули O2, кожен атом пов'язаний валентно тільки з одним іншим атомом - його сусідом по молекулі. Таким чином, від координат кожного атома буде залежати тільки один член суми валентних потенціалів Vbond з рівняння (1):

,

Введемо оператор градієнта за координатами атома i:

,

Тоді:



,

де . Використовуючи даний вираз, можна знайти градієнт потенціалу ковалентного зв'язку:

,

Таким чином, сила, що діє на атом i із-за ковалентного зв'язку з атомом j буде здаватися таким виразом:

,

Так як у випадку молекул O=O кожен атом утворює тільки один зв'язок, сила буде повною силою ковалентного зв'язку для атома i, а - для атома j, тобто і . У загальному випадку, сила буде сумою сил , порахованих для всіх ковалентних зв'язків атома i, а параметри взаємодій Ks і b0 будуть залежати від типів атомів i і j, тобто:

,



де (j1,..., jB) - набір атомів, ковалентно пов'язаних з атомом i.

Програмна реалізація розрахунку потенціалу

Як було обговорено раніше, для даної системи ковалентні зв'язки пов'язують атоми з парними індексами (N = 0, 2,..., N - 2) з атомами, у яких індекс на одиницю більше (N = 1, 3,.. ., N -1). Таким чином, розрахунок сили, що діє між двома атомами, можна записати в наступному вигляді:

1 void computeCovalent ( float3 * h_r , float3 * h_f ,

2 float Ks , float b0 , int N, float L){

3 int i;

4 float3 rij;

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

6 rij = getVector (h_r[i], h_r[i+1] , L);

7 float b = len (rij );

8 float df = Ks *(b - b0 )/b;

9 float3 fij = df*rij;

10 h_f [i] += fij;

11 h_f [i+1] -= fij;

12 }

13 }

Розрахована в даній процедурі сила залежить від координат і використовується в інтеграторі, який ці координати змінює на кожному кроці інтегрування за часом. Таким чином, на кожному часовому кроці силу необхідно розраховувати заново, що відбивається в наступній зміні процедури computeCPU():

1 void computeCPU (){

2 ...

3 for (s = 0; s < output_freq ; s ++){

4 computeCovalent (h_r , h_f , Ks , b0 , N, L);

5 integrate (h_r , h_v , h_f , m, dt , N);

6 }

7 ...

8 }

Так як кожен атом має тільки один ковалентний зв'язок, підсумовування сил на ГП можна проводити не побоюючись того, що два потоки будуть одночасно зберігати результат в одну комірку пам'яті. Кількість потоків на ГП при цьому буде дорівнює кількості ковалентних зв'язків (тобто N/2). Крім того, оскільки даний потенціал вважається першим, можливо пряме збереження розрахованої в даному ядрі сили в відповідну клітинку масиву сил, що діють на атоми системи (замість додавання відповідних значень).

1 __global__ void computeCovalent_kernel ( float3 * d_r ,

2 float3 * d_f , float Ks , float b0 , int N, float L){

3 int d_i = blockIdx .x* blockDim .x + threadIdx .x;

4 if(d_i < N /2){

5 float3 ri = d_r [2* d_i ];

6 float3 rj = d_r [2* d_i + 1];

7 float3 rij = getVector (ri , rj , L);

8 float b = len (rij );

9 float df = Ks *(b - b0 )/b;

10 float3 fij = df*rij;

11 d_f [2* d_i] = fij;

12 d_f [2* d_i + 1] = -fij ;

13 }

14 }

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

1 void computeGPU (){

2 ...

3 for (s = 0; s < output_freq ; s ++){

4 computeCovalent_kernel

5 <<<(N/2 -1)/ BLOCK_SIZE + 1, BLOCK_SIZE >>>(d_r ,

6 d_f , Ks , b0 , N, L);

7 integrate_kernel <<<(N -1)/ BLOCK_SIZE + 1, BLOCK_SIZE >>>(d_r ,

8 d_v , d_f , m, dt , N);

9 }

10 ...

11 }

В результаті виконання отриманої програми молекули кисню будуть вільно рухатися, а атоми всередині кожної молекули - коливатися навколо рівноважної відстані. Наступним кроком стане додавання взаємодій між молекулами кисню.

Розрахунок термодинамічних величин

Енергія системи

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

1 ...

2 float df = Ks *(b - b0 )/b;

3 pot_cov += 0.5f*Ks *(b - b0 )*(b - b0 );

4 float3 fij = df*rij;

5 ...

6 }

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

1 ...

2 if( step % output_freq == 0){

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

4 h_r[i] = transferPBC (h_r[i], L);

5 }

6 appendXYZ (" O2_CPU . coor . xyz", h_r , N);

7 pot_cov /= output_freq ;

8 pot_lj /= output_freq ;

9 printf ("%ld\t%f\t%f\t%f\n",

10 step , pot_cov , pot_lj , pot_cov + pot_lj );

11 pot_cov = 0.0;

12 pot_lj = 0.0;

13 }

14 ...

У першу колонку виведено крок за часом, у другу і третю - значення енергії ковалентних і нековалентних взаємодій. В останню колонку виведена повна енергія системи.

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

1 ...

2 if( step % output_freq == 0){

3

4 cudaMemcpy (h_r , d_r , N* sizeof ( float3 ), cudaMemcpyDeviceToHost );

5 cudaMemcpy ( h_pot_cov , d_pot_cov , N* sizeof ( float ),

6 cudaMemcpyDeviceToHost );

7 cudaMemcpy ( h_pot_lj , d_pot_lj , N* sizeof ( float ),

8 cudaMemcpyDeviceToHost );

9

10 pot_cov = 0.0;

11 pot_lj = 0.0;

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

13 h_r[i] = transferPBC (h_r[i], L);

14 pot_cov += h_pot_cov [i];

15 pot_lj += h_pot_lj [i];

16 }

17 pot_cov /= output_freq ;

18 pot_lj /= 2* output_freq ;

19

20 cudaMemcpy (d_r , h_r , N* sizeof ( float3 ), cudaMemcpyHostToDevice );

21 cudaMemset ( d_pot_cov , 0, N* sizeof ( float ));

22 cudaMemset ( d_pot_lj , 0, N* sizeof ( float ));

23

24 appendXYZ (" O2_GPU . coor . xyz", h_r , N);

25 printf ("%ld\t%f\t%f\t%f\n",

26 step , pot_cov , pot_lj , pot_cov + pot_lj );

27 }

28 ...

Потенційна енергія взаємодії Леннарда-Джонса ділиться навпіл, позаяк одну й та саму бінарну взаємодію між атомами i і j розраховують двічі - в i-тому і в j-тому потоках. Розрахунок потенційної енергії є тривіальним доповненням ядра, що обчислює міжатомні сили для цієї взаємодії. Наприклад, для потенціалу Леннарда-Джонса доповнене ядро набуває такого вигляду:

1 __global__ void computeLJWithPairslist_kernel ( float3 * d_r , float3 * d_f ,

2 float * d_pot_lj ,

3 int* d_pairs_count , int * d_pairslist , float lj_cutoff2 ,

4 float epsilon , float sigma , int N, float L){

5 int d_i = blockIdx .x* blockDim .x + threadIdx .x;

6 if(d_i < N){

7 int p;

8 float3 rij;

9 float3 fi = d_f [d_i ];

10 float3 ri = d_r [d_i ];

11 float pot = d_pot_lj [d_i ];

12 for (p = 0; p < d_pairs_count [d_i ]; p ++){

13 int j = d_pairslist [p*N + d_i ];

14 float3 rj = d_r [j];

15 rij = getVector (ri , rj , L);

16 float rij2 = len2 (rij );

17 if( rij2 < lj_cutoff2 ){

18 float rij2inv = 1.0 f/ rij2 ;

19 float sor2 = sigma * sigma * rij2inv ;

20 float sor6 = sor2 * sor2 * sor2 ;

21 pot += epsilon * sor6 *( sor6 - 2.0f);

22 float df = 12.0 f* epsilon * sor6 *(1.0 f - sor6 )* rij2inv ;

23 float3 fij = df*rij ;

24 fi += fij;

25 }

26 }

27 d_f [d_i] = fi;

28 d_pot_lj [d_i] = pot;

29 }

30 }

Потенційна енергія Леннарда-Джонса на ГП акумулюється окремо для кожного атома в масиві d_pot_lj. Значення, відповідне потоку, зчитують із цього масиву на початку виконання ядра, а оновлюють і зберігають по закінченні виконання ядра.

Середня температура

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

1 void integrate ( float3 * h_r , float3 * h_v , float3 * h_f ,

2 float m, float dt , int N){

3 int i;

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

5 float3 halfdv = 0.5f*h_f[i]*( dt/m);

6 h_v [i] += halfdv ;

7 temperature += len2 ( h_v[i ]);

8 h_v [i] += halfdv ;

9 h_r [i] += h_v[i]* dt;

10 h_f [i] = make_float3 (0.0f, 0.0f, 0.0 f);

11 }

12 }

Щоб із повної кінетичної енергії отримати температуру, необхідно знайти найбільш імовірну (у даному випадку - середню) кінетичну енергію атома і розділити на постійну Больцмана:

1 ...

2 if( step % output_freq == 0){

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

4 h_r[i] = transferPBC (h_r[i], L);

5 }

6 appendXYZ (" O2_CPU . coor . xyz", h_r , N);

7 pot_cov /= output_freq ;

8 pot_lj /= output_freq ;

9 temperature *= 0.5* m/N;

10 temperature /= output_freq *kB;

11 printf ("%ld\t%f\t%f\t%f\t%f\n",

12 step , pot_cov , pot_lj , pot_cov + pot_lj , temperature );

13 pot_cov = 0.0;

14 pot_lj = 0.0;

15 temperature = 0.0;

16 }

17 ...

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

На ГП розрахунок температури також здійснюють в інтеграторі, але значення зберігають для кожного потоку (атома) окремо:

1 __global__ void integrate_kernel ( float3 * d_r , float3 * d_v , float3 * d_f,

2 float * d_temperature , float m, float dt , int N){

3 int d_i = blockIdx .x* blockDim .x + threadIdx .x;

4 if(d_i < N){

5 float3 r = d_r[d_i ];

6 float3 v = d_v[d_i ];

7 float3 f = d_f[d_i ];

8 float3 halfdv = 0.5f*f*( dt/m);

9 v += halfdv ;

10 d_temperature [d_i] += len2 (v);

11 v += halfdv ;

12 r += v*dt;

13 d_r [d_i] = r;

14 d_v [d_i] = v;

15 d_f [d_i] = make_float3 (0.0f, 0.0f, 0.0f);

16 }

17 }

Де d_temperature - масив змінних типу float, виділений у пам'яті ГП. Потім, як і для енергії, отримані значення сум квадратів температур копіюють у пам'ять ЦП, де й проводять розрахунок середньої кінетичної енергії та температури:

1 ...

2 if( step % output_freq == 0){

3 cudaMemcpy (h_r , d_r , N* sizeof ( float3 ), cudaMemcpyDeviceToHost );

4 cudaMemcpy ( h_pot_cov , d_pot_cov , N* sizeof (float ,

5 cudaMemcpyDeviceToHost );

6 cudaMemcpy ( h_pot_lj , d_pot_lj , N* sizeof ( float ),

7 cudaMemcpyDeviceToHost );

8 cudaMemcpy ( h_temperature , d_temperature , N* sizeof ( float ),

9 cudaMemcpyDeviceToHost );

10 pot_cov = 0.0;

11 pot_lj = 0.0;

12 temperature = 0.0;

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

14 h_r[i] = transferPBC (h_r[i], L);

15 pot_cov += h_pot_cov [i];

16 pot_lj += h_pot_lj [i];

17 temperature += h_temperature [i];

18 }

19 temperature *= 0.5* m/N;

20 temperature /= output_freq *kB;

21 pot_cov /= output_freq ;

22 pot_lj /= 2* output_freq ;

23 cudaMemcpy (d_r , h_r , N* sizeof ( float3 ), cudaMemcpyHostToDevice );

24 cudaMemset ( d_pot_cov , 0, N* sizeof ( float ));

25 cudaMemset ( d_pot_lj , 0, N* sizeof ( float ));

26 cudaMemset ( d_temperature , 0, N* sizeof ( float ));

27 appendXYZ (" O2_GPU . coor . xyz", h_r , N);

28 printf ("%ld\t%f\t%f\t%f\t%f\n",

29 step , pot_cov , pot_lj , pot_cov + pot_lj , temperature );

30 }

31 ...

Тут також можна використовувати алгоритми редукції в пам'яті ГП і копіювати вже кінцеве значення температури.

Размещено на Allbest.ru

...