Основы МД: Анализ радиальных функций распределения

PDF version

Описание задачи

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

Цель.

Изучение работы базовых алгоритмов молекулярной динамики на примере моделирования инертного газа (аргона), построение радиальной функции распределения.
Расчеты осуществляются при помощи программного пакета Gromacs.

Подготовка к работе

Для выполнения задачи необходим компьютер, работающий под OS Linux, либо родственными ОС. Расчеты молекулярной динамики рекомендуется выполнять на высокопроизводительных компьютерных кластерах.
Необходимое ПО :
1 Gromacs - осуществляет расчеты молекулярной динамики (в данной статье использована версия 4.5.4) как в обычной, так и в MPI версиях (для распределенных вычислений на кластере)
2 Gnuplot - используется для построения графиков (в статье рассмотрена версия 5)
3 LAM-MPI - драйвер, необходимый для распределенных вычислений
Для удобства работы рекомендуем использовать файловые менеджеры, например, Midnigt Commander.

Описание системы и ее подготовка.

Система представляет собой инертный газ аргон в коробке 10х10х10 нм.
Cкачайте архив в интересующую Вас директорию

wget http://molsim.org/sites/default/files/tutor01.tar_0.gz 

Распакуйте его

tar -xvf file.tar.gz

Создаем систему с помощью файла argongen.pl, в качестве аргумента выступает необходимое количество атомов. Argongen.pl создает .gro-файл и файл топологии (*.top).
Например, мы вводим

./argongen.pl 1000

На выходе мы получим файлы argon1000.gro и topol1000.top
Рассмотрим скрипт Argongen.pl, создающий эти файлы.

#!/usr/bin/perl
$number=$ARGV[0];     #считывает аргумент, в данном случае это 1000
use POSIX;            #подключает модуль работы с числами

#считаем грань куба с округлением, если корень не целый:
$box_size=POSIX::ceil($number**(1/3));	
printf ("The size of the cub's edge = $box_size\n");

#считает исправленное количество атомов для целочисленной грани:
$num=($box_size**3); 
printf ("corrected number of atoms = $num\n");

#создаст и откроет gro-файл на запись и пишет необходимые строки
{open(TE,">argon$num.gro") || die "Can't create the main file"; 
$i=1;$h=0;$r=1;
print(TE "ARGON\n");
print(TE "$num\n");

#откроет файл топологии на запись и пишт необходимые строки
open(TT,">topol$num.top") ||  die "I can't create topology file";#создаем 
print TT ("\#include \"oplsaa.ff/forcefield.itp\" \n");
print TT (" \n");
print TT ("[ moleculetype ]\n");
print TT ("; Name            nrexcl\n");
print TT ("ARGON             5 \n");
print TT (" \n");
print TT (" \n");

#цикл для генерации координат
for ($xi=0; $xi<=($box_size-1); $xi++) 
{for ($yi=0; $yi<=($box_size-1); $yi++)
{for ($zi=0; $zi<=($box_size-1); $zi++)
{
$l="GRO";	              #имя остатка
$ar="Ar";	              #тип атома
$c=0;		              #заряд
$m=40;		              #масса
$x=($xi/$box_size*10);        #координаты в соответствии с пропорциями
$y=($yi/$box_size*10);
$z=($zi/$box_size*10);
if ($r==10){$h=$h+1;$r=1}else{$r=$r+1}

#печать данных в файл в нужном формате
print TE sprintf "%5u%5s%5s%5u%8.3f%8.3f%8.3f\n",$h,$l,$ar,$i,$x,$y,$z; 
print TT sprintf "%5u%5s%5u%6s%5s%6u%5u%5u\n",$i,$ar,$h,$l,$ar,$i,$c,$m;
$i=$i+1}}}

#печать необходимых строк в файл топологии
print TT ("[ system ]\n"); 
print TT ("; Name\n");
print TT ("ARGON\n");
print TT (" \n");
print TT ("[ molecules ]\n");
print TT ("; Compound        #mols\n");
print TT ("ARGON     1\n");
close(TT);                    #закрытие файлов
close(TE);
}

Формат gro-файла:

Номер_остатка  название_остатка  тип_атома  номер_атома координаты_x_y_z

Формат файла топологии:

Номер_атома  тип_атома  номер_остатка  название_остатка  атом  зарядовая_группа  заряд  масса

Поместите наш газ в необходимую коробку:

editconf -f argon1000.gro -o argon1000_newbox.gro -c -box 10 10 10 -bt cubic

Будет использован файл argon1000.gro, запись данных пойдет в argon1000_neewbox.gro, опция -с говорит, что система будет помещена в центр коробки, -box позволяет установить размеры коробки, -bt cubic задает форму коробки, в данном случае — куб.

Затем необходимо подготовить систему, а именно — минимизировать ее энергию. Для этого используйте следующую команду:

grompp -f en_minim.mdp -c argon1000_newbox.gro -p topol1000.top -o em1000.tpr

Опция -f указывает на файл, в котором прописаны параметры действия, совершаемого над сиcтемой, -с входной файл системы, -p файл топологии, -о выходной файл.
Для запуска МД используется:

mdrun -v -deffnm em1000

Опция -v принуждает печатать всю возможную информацию в лог, -deffnm общее имя группы файлов, которые мы используем.
Далее запустите МД в NVE-ансамбле, то есть cохраняется количество вещества, объем и энергия системы. Делается это аналогично минимизации энергии системы.

grompp -f nve_conf.mdp -c em1000.gro -p topol1000.top -o nve1000.tpr
mdrun -deffnm nve1000

Построение функции радиального распределения.

В ходе выполнения команд были получены файлы траектории. На основании них можно построить функцию радиального распределения системы.
Для этого необходимо сначала создать .tpr-файл с топологией:

grompp -f nve_conf.mdp -c argon1000_newbox.gro -p topol1000.top -o topol1000.tpr

Для построения функции радиального распределения используется следующая команда:

g_rdf -f nve1000.trr -s topol1000.tpr -o rdf1000.xvg

Опция -f указывает на входной файл траектории, -s на файл топологии, -o на выходной файл с графиком.
При выполнении команды появится запрос о необходимости выбора группы атомов, используйте группу System (0), на второй запрос ответьте так же.
На выходе будут получены файлы в формате xvg, чтобы получить из них файлы в формате tif необходимо:

gnuplot gnuplotenergy1.scr
./eps_to_img.sh

gnuplotenergy1.scr находится в папке RDF.
Для того, чтобы построить графики функции радиального распределения для системы с другим числом атомов необходимо изменить названия соответствующих файлов в скрипте.

Функция радиального распределения для системы с 1000 атомов аргона

Функция радиального распределения для системы с 8000 атомов аргона

Функция радиального распределения для системы с 64000 атомов аргона

Общий график для всех трех систем

All

На общем графике видно, что все три системы имеют разные функции радиального распределения и на разных r достигают единицы. С увеличением плотности r увеличивается, то есть увеличивается расстояние, на котором влияние частиц друг на друга отсутствует.