Работа 3

Материал из md
Перейти к: навигация, поиск
Моделирование простой жидкости

Скачать файлы для выполнения работы: 03_argon.tar.gz

Скачать описание работы в pdf-формате: 03_argon.pdf

Цель

Провести МД-моделирование одноатомной жидкости с потенциалом Леннарда-Джонса. Научиться работать с основными программами анализа моделей. Изучить основные структурные и динамические характеристики простой жидкости.

Задание

  1. Создать начальную конфигурацию леннард-джонсовских атомов.
  2. Провести релаксацию системы
  3. Получить равновесную модель.
  4. Рассчитать функцию радиального распределения (ФРР).
  5. Рассчитать автокорреляционную функцию скорости (АФС).
  6. Рассчитать среднеквадратичное смещение r2(t) и коэффициент самодиффузии D для атомов модельной жидкости.

Вопросы

  1. Зачем нужны два этапа релаксации модели?
  2. Чем отличаются кривые, полученные в ваших расчётах, от соответствующих кривых, приведенных на Рис. 3.1–3.3?
  3. Найти соотношение между положениями первого и второго максимумов ФРР.
  4. Сформулировать основные структурные черты простой жидкости.

Общие замечания

Проще всего метод молекулярной динамики реализуется для сферически-симметричного потенциала взаимодействия между атомами. В этом случае, с точки зрения механики, мы имеем дело с системой взаимодействующих точек без каких-либо дополнительных наложенных ограничений и связей. Обычно подразумевается, что взаимодействие парное, т.е. сила, действующая на данный атом, рассчитывается как простая сумма сил от других атомов системы. Атомы благородных газов (кроме гелия) являются хорошим примером таких систем. Однако это приближение работает и в других случаях, в частности, для металлов типа меди, свинца, натрия (т.е. для достаточно тяжелых атомов с простыми внешними электронными оболочками), для органических молекул, имеющих компактную форму, как например, метан, молекулы которого можно считать сферическими. Жидкости из таких молекул обычно называют простыми, все они имеют общие структурные черты. Это отличает их от систем с направленными (ориентационно-зависимыми) взаимодействиями, как например между атомами кремния, или молекулами воды.

Конкретный вид сферически-симметричного потенциала межатомного взаимодействия может быть различным. Общим свойством является наличие сильного отталкивания при сближении атомов, т.е. наличие "жесткой сердцевины" атомов. Известные ван-дер-ваальсовские радиусы атомов характеризуют эти расстояния. Поэтому в качестве базовой структурной модели простых жидкостей часто используют неупорядоченную систему твердых шаров.

Компьютерное моделирование простых жидкостей позволило понять сущность жидкого состояния, его отличия от кристаллической и газообразной фазы. Прежде всего было осознано, что в жидкости, как и в твердой фазе, существует структура. Однако она понимается иначе, чем в кристалле. (Напомним, что кристаллографы под структурой понимают "конкретное расположение атомов в пространстве", которое задается кристаллической решеткой). В жидкости, а также в стекле, т.е в некристаллических системах, структура проявляется в виде статистических корреляций в расположении атомов. При этом корреляции существуют только на расстояниях, которые не превышают нескольких атомных диаметров.

Основной структурной характеристикой жидкостей является функция радиального распределения (ФРР), она задаёт изменение числа атомов жидкости (или плотности) в зависимости от расстояния от центра заданного атома. Заметим, что в жидкости для любого конкретного атома в заданный момент времени такое распределение соседей будет каким-то своим. Однако после усреднения по всем "центральным" атомам (или вокруг одного атома, но в течение долгого времени), получается функция, несущая информацию о характерном взаимном расположении атомов в пространстве для данной системы. Обычно ФРР нормируют на среднюю плотность системы и делят на квадрат расстояния, чтобы она стремилась с ростом r к единице. В теории жидкостей такую функцию называют парной корреляционной функцией и обозначаются как g(r). Можно сказать, что она характеризует "вероятность" обнаружить атом на расстоянии r от данного атома. На Рис. 3.1 показана g(r) жидкого аргона. Имеется четкий пик, соответствующий ближайшим соседям, и последующие затухающие осцилляции, что типично для всех простых жидкостей.

Жидкий аргон был одним из первых объектов применения метода молекулярной динамики в статистической физике. Взаимодействие атомов аргона хорошо описывается леннард-джонсовским потенциалом . Значения параметров σ=0,34 нм и ε=9,7 КДж/моль (или 119.8 K) были оценены в свое время из экспериментов с газообразным аргоном.

Fig3.1.png
Рис. 3.1. Парная корреляционная функция g(r) для жидкого аргона в приведённых единицах (σ= 1, ε=1).


Fig3.2.png
Рис. 3.2. Автокорреляционная функция скорости для жидкого аргона при разных термодинамических условиях. Для газа спадает экспоненциально, для жидкости имеет минимум, который становится глубже для переохлажденной жидкости, а в застеклованном состоянии наблюдаются затухающие осцилляции.


Помимо структурной, в жидкости имеется динамическая корреляция. Движение атома в жидкости ограничено ближайшими соседями, которые образуют «клетку», где атом вынужден находиться некоторое время. Обычно это время не превышает пикосекунды, что видно из автокорреляционной функции скорости для аргона, см. Рис. 3.2. Напомним, что для газа эта функция спадает экспоненциально, а в твердой фазе имеются осцилляции.

Однако на более длинных временах, превышающих время «сидения» в клетке, поведение атома в жидкости хорошо описывается диффузионными законами, Рис. 3.3. Средний квадрат смещения атомов растёт линейно со временем согласно формуле Эйнштейна (r2=6Dt). Угол наклона этой прямой определяет коэффициент самодиффузии атомов.

Fig3.3.jpg
Рис. 3.3. Средний квадрат смещения атомов от времени в различных жидкостях.


Указания к выполнению работы

Создание начальной конфигурации для моделирования

Работаем в директории 03_argon/

cd 03_argon/

В ней лежит файл топологии topol.top и файл молекулярно-динамических параметров grompp.mdp. Файл структуры conf.gro мы приготовим самостоятельно, см. ниже.

В отличие от предыдущей работы, где мы моделировали леннард-джонсовскую жидкость в приведённых единицах (σ=1, ε=1, m=1), здесь мы будем моделировать «реальный» аргон, т.е. работать с размерными величинами. Напомним, что в пакете Gromacs расстояние измеряется в нанометрах, а время в пикосекундах. Используем стандартное поле сил OPLS (optimized potential for liquid simulations), в котором для аргона приняты следующие значения параметров потенциал Леннарда-Джонса: σ=0,3401 нм, ε=0,979 кДж/моль.

Посмотрим файл топологии topol.top:

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               3               no             0.5     0.5

[ atomtypes ]
; from opls forcefield opls_097 atom type
; name  bond_type    mass    charge   ptype     sigma      epsilon
AR   AR  18   39.94800     0.000       A    3.40100e-01  9.78638e-01

[ moleculetype ]
; molname       nrexcl
Argon	1

[ atoms ]
;   nr   type  resnr residue  atom   cgnr     charge       mass
     1	   AR	   1	 AR     AR      1     0            39.948

[ system ]
; Name
Argon

[ molecules ]
; Compound        #mols
Argon		216

Масса атома аргона задается в единицах г/моль, в данном случае это 39,948 а.е.м., а глубина потенциальной ямы в КДж/моль: 9.79e-01. Обратите внимание на отличия параметров этого файла от файла топологии предыдущего задания (№ 2).

Посмотрим файл МД-параметров grompp.mdp:

;       Input file
;
title               =  Argon                    ; a string
dt                  =  0.002                    ; time step
nsteps              =  100000                   ; number of steps (200ps)
nstcomm             =  5                        ; reset c.o.m. motion
nstxout             =  500                      ; write coords 1 ps
nstvout             =  500                      ; write velocities 1ps
nstlog              =  1000                     ; print to logfile 1ps
nstcalcenergy       =  5                        ; print energies 0.1ps
nstenergy           =  100                      ; print energies 0.1ps
nstlist             =  6                        ; update pairlist
ns_type             =  grid                     ; pairlist method
coulombtype         =  cut-off
DispCorr            =  EnerPres
rlist               =  1.5                      ; cut-off for ns
rvdw                =  1.5                      ; cut-off for vdw
rcoulomb            =  1.5                      ; cut-off for coulomb
Tcoupl              =  v-rescale                ; temperature coupling
ref_t               =  90.0                     ; reference temperature,K
tc-grps             =  System
tau_t               =  0.5
Pcoupl              =  berendsen                ; pressure bath
Pcoupltype          =  isotropic                ; pressure geometry
tau_p               =  1.0                      ; p-coupling time
compressibility     =  4.5e-5
ref_p               =  1.0                      ; reference pressure, bar
gen_vel             =  yes                      ; gener. init. velocit.
gen-temp            =  90
gen-seed            =  173529

Заметим, что в файле параметров могут быть заданы не все величины, а лишь те, которые отличаются от параметров, задаваемых «по умолчанию».

Обратим внимание на следующие параметры:

  • gen_vel = yes означает, что в первый момент времени программа сгенерирует скорости для всех молекул. Они будут соответствовать заданной температуре и иметь максвелловское рапределение.
  • Pcoupl = berendsen включает баростат Берендсена. Это означает, что мы используем NPT-ансамбль (в предыдущей работе был NVT).
  • dt = 0.002 указывает, что шаг интегрирования равен 0.002 пикосекундам. Обратим внимание на то, что в работе № 2 время t* было в безразмерных единицах. Не трудно убедиться. что оно связано с реальным временем формулой . Подставляя соответствующие значения для аргона, найдем [пс] ≈ 2.17 t*. Таким образом, используемое здесь число 0.002 в качестве шага интегрирования для аргона, составляет примерно 0.001 в безразмерных единицах. Из предыдущей работы мы знаем, что это разумная величина для молекулярно-динамического моделирования.

Теперь займемся файлом структуры. Предлагается создать его на основе конечной конфигурации из предыдущей работы. Копируем в нашу текущую директорию файл конечной конфигурации confout.gro, полученный в предыдущей работе. Например, из директории 002, расположенной в timestep/, находящейся в 02_timestep/. Скопированный файл назовём стандартным именем conf.gro:

cp ../02_timestep/timestep/002/confout.gro conf.gro

Эта конфигурация содержит 216 атомов, координаты которых записаны в приведенных единицах. Сделаем из неё более крупную модель и запишем координаты атомов в размерных единицах. Новую исходную конфигурацию удобно создать в специальной директории. Создадим поддиректорию preparing, перейдем в неё и перепишем туда файл conf.gro из директории 03_argon/:

mkdir preparing
cd preparing
mv ../conf.gro .

Увеличить модель можно с помощью программы genconf, которая транслирует (размножает) исходный бокс вдоль его осей, Рис. 3.4. Сделаем это по всем трем направлениям дважды:

genconf –nbox 2 2 2

Здесь по умолчанию в качестве входного файла используется файл со стандартным именем conf.gro из текущей директории, числа после ключа –nbox означают количество трансляций вдоль каждой оси. Новая модель запишется в файл со стандартным именем out.gro. (Если бы мы хотели, чтобы выходной файл имел другое имя, то его нужно было бы указать явно, поставив перед ним ключ -o ).

Fig3.4.png
Рис. 3.4. Трансляция исходного бокса по трем направлениям дважды

Теперь отмасштабируем модельный бокс. В исходном файле величина σ равнялась единице, а теперь она 0,3401 нм. Координаты атомов модели достаточно умножить на коэффициент, равный к этому значению. Для этого в пакете Gromacs есть команда editconf:

editconf –f out.gro –o scaled.gro –scale 0.34 0.34 0.34

В данном примере мы используем ключ –f , который означает, что в качестве исходного файл будет взят файл с именем out.gro. После ключа -о указываем желаемое имя выходного файла, в данном случае это scaled.gro. Числа после ключа –scale являются множителями для масштабирования вдоль каждой оси. Для простоты, мы используем числа 0,34 вместо 0.3401, поскольку здесь мы задаем всего лишь исходные положения атомов. В процессе релаксации атомы начнут двигаться, подчиняясь заданным термодинамическим условиям. При этом исходные значения координат не так важны, главное, чтобы в начальной конфигурации не было сильных перекрываний между атомами.

В нашей новой модели другое число атомов. Это надо учесть в файле топологии. Напомним, он находится в директории 03_argon/. В последней строке этого файла содержится количество атомов.

 [ molecules ]
; Compound        #mols
Argon           216

Теперь вместо числа 216 надо поставить число молекул, которое получилось в нашей новой модели. Для этого надо войти в текстовый редактор и исправить это число.Чему равно это число? (подсказка: 2х2х2х216).

Напомним, что в операционной системе Linux имеется несколько текстовых редакторов:

  • редактор файлового менеджера mc (midnight commander). Файловый менеджер запускается командой mc, он напоминает всем известный far – менеджер. Режим редактирования запускается выбором нужного файла и нажатием клавиши F4.
  • редактор vi. Его любят те, кто привык пользоваться сочетаниями горячих клавиш и предпочитает интерфейс командной строки (для выхода из редактора надо нажать q). Редактирование данного файла (например, topol.top) в этом редакторе запускает команда:
vi topol.top
  • редактор kate, kwrite, gedit или другие им подобные. Они имеют оконный интерфейс и более простой, интуитивно понятный набор команд. В этом случае запускаем команду:
kate topol.top

Релаксация модели

Наша большая модель составлена из восьми одинаковых частей, полученных трансляцией. Необходимо провести предварительную релаксацию, чтобы жидкость стала структурно-однородной. Эту работу также удобно проводить в отдельной директории. Создадим поддиректорию equi1 (equilibration) в директории 03_argon/ и перейдем в неё:

mkdir equi1
cd equi1

Работая в этой директории, удобно сделать ссылку на файл структуры scaled.gro (который находится в директории preparing) как на файл со стандартным именем conf.gro:

ln –s ../preparing/scaled.gro conf.gro

Файл параметров grompp.mdp удобно скопировать в equi1/из директории 03_argon/:

cp ../grompp.mdp .

(Напомним, что точка в конце этой команды означает, что файл имеет то же имя, что и исходный файл, т.е. grompp.mdp). Откроем этот файл в текстовом редакторе, чтобы убедиться, что в разделах, отвечающих за термостат и баростат, прописаны термостат v-rescale и баростат Берендсена:

tcoupl = v-rescale,
Pcoupl = Berendsen

Их использование позволяет системе быстро отрелаксировать к заданным температуре и давлению.

Число шагов моделирования:

nsteps= 100000

Это соответствует 200 пикосекундам, поскольку шаг интегрирования равен 0.002. Можно оставить эти параметры без изменения.

Создадим ссылку на файл топологии из данной директории. Мы не будем его больше менять, поэтому нет смысла его копировать.

ln –s ../topol.top .

После этого запускаем препроцессинг:

grompp

Система из 1728 частиц уже достаточно большая, чтобы использовать параллельные вычисления. По умолчанию Gromacs (в последних версиях) использует все доступные ядра процессора на данном компьютере, т.е. для запуска моделирования можно просто набрать команду

mdrun

Можно добавить ключ -v , если нужен вывод текущей инфрормации о ходе вычислений. При необходимости можно задать конкретное число ядер, которые будут использоваться при расчёте:

mdrun –nt 8 –dd 2 2 2 -v

Ключ –nt 8 означает, что будут использоваться восемь ядер. Ключ –dd 2 2 2 определяет, каким образом модель будет поделена между ядрами компьютера. Числа 2 2 2 показывают, что каждое ребро модельного бокса будет поделено на 2 части. Таким образом, наша модель разбивается на 8 частей, каждая из которых будет обрабатываться своим ядром.

После выполнения работы убедимся, что релаксация прошла успешно. Для этого можно воспользоваться командой g_energy, чтобы посмотреть, например, на изменение плотности (density) модели. Она должна выйти на стационар, Рис. 3.5.

g_energy –o density.xvg -w
Fig3.5.png
Рис. 3.5. Первый этап релаксации модели. Изменение плотности исходной конфигурации со временем. Видно, что плотность быстро выходит на стационарное значение.

Результатом релаксации является файл конечной конфигурации confout.gro, который мы будем использовать в дальнейшем.

Второй этап релаксации

Последняя конфигурация после первого этапа релаксации (файл confout.gro) ещё не является равновесной. Нужно провести дополнительную релаксацию, на этот раз с параметрами основного моделирования, где обычно используется термостат Nose-Hoover’a и баростат Parrinello-Rahman’а. Напомним, на первом шаге релаксации использовались термостат v-rescale и баростат Берендзена. Они быстро и надежно приводят систему к заданным значениям температуры и давления, однако генерируемый ансамбль конфигураций получается не вполне правильным (характер флуктуаций энергии и других характеристик отличается от флуктуаций равновесной системы). Смысл второго этапа релаксации – дать системе прийти к равновесию с новыми термостатом и баростатом. Это предлагается сделать в отдельной поддиректории equi2. Итак, создаем директорию equi2 и переходим в нее:

cd ../
mkdir equi2
cd equi2

Делаем ссылку на конечный файл предыдущего этапа релаксации, чтобы использовать его как файл структуры:

ln –s ../equi1/confout.gro conf.gro

Делаем ссылку на файл топологии (напомним, он лежит в 03_argon/):

ln –s ../topol.top .

Копируем в нашу директорию файл параметров моделирования grompp.mdp:

cp ../grompp.mdp .

и проведём в нём изменения параметров:

Tcoupl = Nose-Hoover
Pcoupl = Parrinello-Rahman.

Кроме того, здесь нужно отказаться от генерации скоростей, т.е задать gen_vel = no , поскольку скорости атомов записаны в конечном файле предыдущего этапа релаксации.

После этого запускается препроцессинг и моделирование:

grompp
mdrun -v

После завершения вычислений можно посмотреть флуктуации плотности системы так же, как после первого шага релаксации. Можно нарисовать обе кривые на одном графике.

Fig3.6.png
Рис. 3.6. Флуктуации плотности модели в процессе первого этапа релаксации (красная линия, см. Рис. 3.4) и во время второго этапа (синяя).

Вспомним, что первая кривая у нас записана в директории equi1/ в файле под именем density.xvg. Извлечём новую кривую и запишем её также под именем density.xvg в текущей директории (equi2):

g_energy –o density.xvg

Указав эти файлы для программы xmgrace, мы получим обе кривые на одном графике, см. Рис. 3.6:

xmgrace density.xvg ../equi1/density.xvg

Видно, что термостат Нозе-Хувера и баростат Парринелло-Рамана дают более хаотические флуктуации, что лучше соответствует реальности. Однако если бы мы применили их на первом шаге релаксации, где начальная плотность системы заметно отличалась от равновесной, то возникли бы большие флуктуации, и выход на стационар занял бы продолжительное время.

Основное моделирование (production run)

Наконец, создадим модель, которую будем использовать для анализа. Это моделирование проведём в основной директории 03_argon/. Перейдём в неё и создадим ссылку на конечный файл второго этапа релаксации confout.gro, который будет файлом структуры conf.gro, представляющим начальную конфигурацию основного моделирования:

cd ../
ln –s ../equi2/confout.gro conf.gro

Скопируем сюда также файл параметров grompp.mdp, который мы использовали в папке equi2:

cp equi2/grompp.mdp grompp.mdp

Число шагов моделирования в нём задано как nsteps = 100000. При необходимости можно его увеличить, чтобы основная молекулярно-динамическая траектория соответствовала большему времени для улучшения статистики. Однако в данной работе можно не менять этот параметр, ограничиться 200 пикосекундами. Заметим, что в нашем файле задано значение nstvout=500, которое указывает, через сколько шагов записываются скорости атомов. Поскольку мы планируем рассчитывать автокорреляционную функцию скорости, то их следует записывать чаще. Зададим nstvout=5, т.е. через каждые 0.01 пс.

Теперь проведём препроцессинг и запустим основное моделирование:

grompp
mdrun -v

После завершения вычислений в директории 03_argon/ появятся файлы, содержащие молекулярно-динамическую траекторию нашей леннард-джонсовской системы: traj.trr, ener.edr и др. Теперь становится понятно, зачем мы заводили дополнительные поддиректории: preparing/, equi1/, qui2/. Все ненужные для дальнейшей работы файлы остались в них и не загромождают основную директорию 03_argon/.

Анализ модели

Функция радиального распределения

Программа g_rdf рассчитывает парную корреляционную функцию g(r) между атомами модели. По умолчанию она усредняется по всем конфигурациям из файла траектории traj.trr:

g_rdf -s –w

Ключ –s сообщает, что будет также использоваться информация о топологии из файла topol.tpr. Для программы g_rdf необходимо, чтобы был указан или файл топологии, или индекс-файл, с которым мы познакомимся позднее. В данном случае удобнее воспользоваться первым вариантом. В процессе работы программа g_rdf спрашивает, между какими группами атомов рассчитывать парную корреляционную функцию. Предлагается на выбор три группы. В нашем случае все они содержат одну и ту же информацию, каждая из них определяет все атомы системы. Поэтому можно выбрать два раза любую из них, например, группу 0 ‘System’:

Select a reference group and 1 group
Group 0 ( System) has 1728 elements
Group 1 ( Other) has 1728 elements
Group 2 ( AR) has 1728 elements
Select a group:

Ключ –w вызывает автоматическое обращение к программе xmgrace, которая сразу же рисует график, Рис. 3.7.

Fig3.7.png
Рис. 3.7. Парная корреляционная функция для жидкого аргона.

Этот график также записывается в текстовом виде в файле rdf.xvg.

Автокорреляционная функция скорости

Для расчёта автокорреляционной функции скорости существует программа g_velacc. Здесь достаточно выполнить команду:

g_velacc –w

Вся необходимая информация возьмётся из файла траектории со стандартным именем traj.trr. Результат расчёта записывает в файл vac.xvg, а ключ -w обеспечивает автоматический вывод на экран графика этой функции, Рис. 3.8.

Fig3.8.png
Рис. 3.8. Автокорреляционная функция скорости для модели, показанной на Рис. 3.6


Обратим внимание на то, что автокорреляционная функция скорости будет рассчитана по умолчанию на всём возможном интервале времён для нашей траектории, что составляет половину времени эволюции модели, которое у нас 200 пс. Однако характерные изменения эта функция имеет только на коротких временах, порядка пикосекунды. Поэтому надо изменить масштаб рисунка по оси времени.

Среднеквадратичное смещение атомов со временем. Коэффициент самодиффузии

Расчёт среднеквадратичного смещения атомов производится программой g_msd. (Напомним, что подробную информацию об использовании программы можно найти по ключу -h , например, в данном случае получить подсказку можно по команде g_msd -h). Все необходимые файлы у нас имеют стандартные имена, поэтому достаточно написать:

g_msd -s -w

Как и в случае расчёта парной корреляционной функции, координаты атомов берутся по умолчанию из файла траектории traj.trr. Также используется файл топологии topol.tpr. Программа также запросит информацию, для какой группы атомов следует проводить расчёт. В нашем случае все они содержат одну и ту же информацию, каждая из них определяет все атомы системы. Можно снова выбрать любую группу, например, 0.

Рассчитанная зависимость среднеквадратичного смещения от времени записывается по умолчанию в файл msd.xvg и автоматически изображается на графике, Рис. 3.9.

Fig3.9.png
Рис. 3.9. Среднеквадратичное смещение атомов как функция от времени.


Получается линейная зависимость, по наклону которой определяется коэффициент самодиффузии D. Он рассчитывается автоматически и выводится на экран после завершения работы программы:

Last frame 200 time 200.000
Used 21 restart points spaced 10 ps over 200 ps
Fitting from 20 to 180 ps
D[ System] 2.3555 (+/– 0.0836) 1e-5 cm^2/s
Executing 'xmgrace msd.xvg &'