PDF versionОписание задачи
Данная задача посвящена рассмотрению поведения системы при различных режимах термостатирования и баростатирования.
Цель:
изучение работы базовых алгоритмов молекулярной динамики на примере моделирования инертного газа (аргона).
Расчеты осуществляются при помощи программного пакета Gromacs.
Подготовка к работе
Для выполнения задачи необходим компьютер, работающий под OS Linux, либо родственными ОС. Расчеты молекулярной динамики рекомендуется выполнять на высокопроизводительных компьютерных кластерах.
Необходимое ПО :
- Gromacs - осуществляет расчеты молекулярной динамики (в данной статье использована версия 4.5.4) как в обычной, так и в MPI версиях (для распределенных вычислений на кластере)
- Gnuplot - используется для построения графиков (в статье рассмотрена версия 5)
- LAM-MPI - драйвер, необходимый для распределенных вычислений
Для удобства работы рекомендуем использовать файловые менеджеры, например, Midnigt Commander.
Система представляет собой 1000 атомов инертного газа Аргона в коробке 10х10х10 нм.
Описание и подготовка подробно описаны в другой статье.
Cкачайте архив в интересующую вас директорию
wget АРХИВ НЕ ГОТОВ
Распакуйте его
tar -xvf file.tar.gz
Режимы термостатирования
Мы будем использовать 4 различных термостата
- Berendsen
- Nose-Hoover
- Nose-Hoover with chains
- Velocity rescaling thermostat
Термостат Берендсена.
В данном термостате температура системы корректируется так, что отклонение экспоненциально затухает с некоторой постоянной времени τ.
Термостат Берендсена подавляет колебания кинетической энергии системы и, следовательно, не может воспроизводить траектории в соответствии с каноническим ансамблем.
Влияние времени релаксации на поведение системы при нагревании с нуля Кельвинов
Для этого необходимо рассмотреть скрипт do.sh
do.sh
#!/bin/bash l="1 10 100 1000 5000"; #Задает необходимые времена релаксации for a in $l; do cat << _EOF_ > NVTBER${a}ps.mdp #создает файл параметров МД title = OPLS Ar NVT equilibration in berendsen thermostat ; Run parameters integrator = md ; leap-frog интегратор nsteps = 5000000 ; количество шагов dt = 0.002 ; шаг в пс ; Output control nstxout = 1000 ; сохранять координаты раз в 1000 шагов nstvout = 0 ; не записывать скорости nstenergy = 1000 ; сохранять координаты раз в 1000 шагов nstlog = 1000 ; изменять лог раз в 1000 шагов ; Neighborsearching ns_type = grid ; искать соседей по сетке nstlist = 5 ; количество шагов до обновления списка соседей rlist = 1.0 ; радиус обрезания построения списка соседей rcoulomb = 1.0 ; радиус обрезания электростатических взаимодействий rvdw = 2.0 ; радиус обрезания Ван-дер-Ваальсовых взаимодействий ; Electrostatics coulombtype = cut-off ; способ выключения электростатических взаимодействий ; Temperature coupling is on tcoupl = berendsen ; Berendsen thermostat tc-grps = AR ; группа для термостатирования tau_t = $a ; время релаксации ref_t = 300 ; температура термостата ; Pressure coupling is off pcoupl = no ; не термостатировать в NVT ; Periodic boundary conditions pbc = xyz ; задает трехмерную периодическую систему ячеек ; Dispersion correction DispCorr = EnerPres ; вносит поправку в дальние взаимодействия ; Velocity generation gen_vel = no ; не задает скорости по Максвеллу gen_temp = 0 ; стартовая температура gen_seed = 6 ; этот ключ для рандомайзера _EOF_ #закрывает файл grompp_4.5.4 -f NVTBER${a}ps.mdp -c em.gro -p topol.top -n indx.ndx -o NVTBER${a}ps.tpr #см ниже cat << _EOF_ > NVTBER${a}ps.sh #создает скрипт системы постановки задач #!/bin/sh #PBS -N ARGON_BE #PBS -l nodes=1:ppn=2 #один узел, 2 процессора #PBS -q day #PBS -V lamboot -v -ssi boot tm cd ${PWD} mpirun C mdrun_4.5.4_MPI -v -deffnm NVTBER${a}ps lamhalt -v $PBS_NODEFILE _EOF_ #закрывает файл qsub NVTBER${a}ps.sh #ставит задачу в очередь done
Точка с запятой обозначают комментарий в файле конфигурации, решетка - комментарий в скрипте, таким образом, параметры vdwtype и rvdw_switch изначально закомментированы.
У команды grompp_4.5.4, создающей бинарные файлы для запуска молекулярной динамики (mdrun), есть несколько ключей
- -f файл на вход с параметрами МД (в данном случае файл конфигурации)
-c файл структуры системы
-p файл топологии
-o создаваемый бинарный файл
Если вы не используете систему постановки задач, закомментируйте с помощью # все после grompp и до done, далее запускайте каждую обработку по отдельности с помощью mdrun.
Запустите скрипт do.sh
./do.sh
Теперь в директории появится россыпь файлов вида NVTBER(цифра)ps.расширение.
Если вы не использовали систему постановки задач, то необходимо запустить их обработку вручную по очереди.
Например, для времени релаксации 100:
mdrun_4.5.4 -v -deffnm NVTBER100ps
В ином случае - проверьте с помощью команды qstat, запустились ли ваши задачи.
Если все хорошо - можно подождать или приступить к дальнейшему выполнению задачи.
После завершения рассчетов запустите скрипт energy.sh
./energy.sh
Рассмотрим его подробно:
energy.sh
#!/bin/bash l="1 10 100 1000 5000"; #берет те же времена релаксации for a in $l; do g_energy_4.5.4 -f NVTBER${a}ps.edr -o energyNVT${a}ps.xvg < enter.txt #запускает расчет температуры, файл txt содержит номер параметра done gnuplot5 gnuplotenergy1.scr #создает eps файл, scr файл содержит параметры
Файл enter.txt содержит в себе одну цифру и два перехода на новую строку, что задает расчет температуры, если удалить файл или закомментировать его прикрепление, то можно задавать их вручную.
Скрипт energy.sh создаст файлы вида energyNVT(время релаксации)ps.xvg, их можно просмотреть в текстовом редакторе или с помощью xmgrace.
Например, для времени1000:
xmgrace energyNVT1000ps.xvg
Скрипт gnuplotenergy1.scr обрабатывает только файлы вышеупомянутого вида с заданными шагами интегрирования и на выходе дает файл energy1.eps
Скрипт можно изменить, например, толщина линий задается ключом lw 2 после каждого файла. Параметр set key задает положение легенды, например, set key bottom rigth поместит легенду в верхний правый угол.
Параметр set xrange задает масштаб оси абсцисс.
Далее, используя скрипт eps_to_img.sh, конвертируем eps файл в png
./eps_to_img.sh
(Так же можно изменить разрешение и формат файла, изменив eps_to_img.sh)
Теперь можно просмотреть получившийся график
если вы работаете в Linux или в Windows с запущенным X сервером: display energy1.tif можно скачать файл, например, для нашего кластера команда имеет такой вид: scp student@biosim.moldyn.org:/путь до файла/energy1.tif /куда скачать/
Должны получиться графики такого вида:
_______________________________
Как и ожидалось, температура возрастает экспоненциально, практически отсутствуют флуктуации.
Влияние времени релаксации на поведение системы при поддержании температуры 300 К
Для этого измените файл do.sh следующим образом:
do.sh
#!/bin/bash ... ... ... ; Velocity generation gen_vel = yes ; задает скорости по Максвеллу gen_temp = 300 ; стартовая температура gen_seed = 6 ; этот ключ для рандомайзера ... ... ... done
Постройте графики так же как и в случае разогрева системы.
Должны получиться аналогичные графики:
_______________________________
Гистограмма распределения температур:
Температура держится в очень узком интервале и флуктуации незначительны. Этот термостат хорошо поддерживает необходимую температуру, но не задает распределение энергии, нарушая канонический ансамбль.
Nose-Hoover
Алгоритм Берендсена удобен для приведения системы к равновесной температуре, но после ее достижения гораздо важнее моделировать корректный NVT ансамбль. Это невозможно при использовании термостата Берендсена.
Для корректного моделирования канонического ансамбля можно использовать алгоритм Nose-Hoover'а. Его суть заключается в том, что гамильтониан системы изменяется путем добавления к нему члена, описывающего тепловой резервуар, а сообщение с ним добавляется в уравнения движения.
В связи с тем что происходит постоянный периодический обмен энергии, возникают ее колебания, что лучше описывает канонический ансамбль
Перейдите в директорию /Nose-Hoover
cd Nose-Hoover
Влияние времени релаксации на поведение системы при нагревании с нуля Кельвинов
Для запуска расчетов так же запускается скрипт do.sh, от предыдущих он отличается следующим:
do.sh
#!/bin/bash l="1 10 100 1000 5000"; #Задает необходимые времена релаксации for a in $l; do cat << _EOF_ > NVTNH${a}ps.mdp #создает файл параметров МД ... ... ... ; Temperature coupling is on tcoupl = nose-hoover ; Nose-Hovver thermostat tc-grps = AR ; группа для термостатирования ... ... ... grompp_4.5.4 -f NVTNH${a}ps.mdp -c em.gro -p topol.top -n indx.ndx -o NVTNH${a}ps.tpr #см ниже cat << _EOF_ > NVTNH${a}ps.sh #создает скрипт системы постановки задач #!/bin/sh #PBS -N ARGON_NH ... ... ... mpirun C mdrun_4.5.4_MPI -v -deffnm NVTNH${a}ps ... ... ... done
То есть изменился лишь режим термостатирования и названия файлов.
Выше описанным путем получим графики:
_______________________________
Из графиков видно, что вместо того, чтобы поддерживать температуру на нужной отметке, алгоритм скачкообразно изменяет ее с 0К до некоторых очень высоких значений температуры. Средняя температура таких колебаний близка к 300К так что формально термостат работает. По-видимому, проблема заключается в большой разнице температур при малой теплоемкости исследуемой системы.
Влияние времени релаксации на поведение системы при поддержании температуры 300 К
Как и в случае термостата Берендсена, необходимо изменить скрипт do.sh:
do.sh
#!/bin/bash ... ... ... ; Velocity generation gen_vel = yes ; задает скорости по Максвеллу gen_temp = 300 ; стартовая температура gen_seed = 6 ; этот ключ для рандомайзера ... ... ... done
Получите графики:
_______________________________
Гистограмма распределения температур:
Как видно из графиков, на всем временном промежутке температура близка к 300К, но ее изменение явно не соответствует хаотическому поведению.
Nose-Hoover with chains
Динамика Nose-Hoover'а неприменима для несложных систем, так как в этом случае с определенной вероятностью в энергообмене будет задействован только некоторый участок фазового пространства.
Для того, чтобы исправить этот недостаток применяются цепи Nose-Hoover'a, которые представляют собой набор термостатов, каждый из которых имеет свою температуру и способен обмениваться энергией с системой и другими термостатами. Если устремить количество термостатов к бесконечности, то система гарантировано будет эргодичной, так как каждый атом будет иметь одинаковую вероятность энергообмена с термостатом. На практике же было показано, что включение даже небольшого количества цепей сильно повышает качество термостатирования.
Перейдите в директорию /Nose-Hoover_with_chains
cd Nose-Hoover_with_chains
Влияние времени релаксации на поведение системы при нагревании с нуля Кельвинов
Для запуска расчетов так же запускается скрипт do.sh, от предыдущих он отличается следующим:
do.sh
#!/bin/bash l="1 10 100 1000 5000"; #Задает необходимые времена релаксации for a in $l; do cat << _EOF_ > NVTNH${a}ps.mdp #создает файл параметров МД ... integrator = md-vv ; интегратор velocity Verle ... ; Temperature coupling is on tcoupl = nose-hoover ; Nose-Hovver thermostat tc-grps = AR ; группа для термостатирования ... ; Velocity generation gen_vel = no ; не задает скорости по Максвеллу gen_temp = 0 ; стартовая температура ... grompp_4.5.4 -f NVTNH${a}ps.mdp -c em.gro -p topol.top -n indx.ndx -o NVTNH${a}ps.tpr #см ниже cat << _EOF_ > NVTNH${a}ps.sh #создает скрипт системы постановки задач #!/bin/sh #PBS -N ARGON_NH ... ... ... mpirun C mdrun_4.5.4_MPI -v -deffnm NVTNH${a}ps ... ... ... done
То есть изменился режим термостатирования и интегратор, так как leap-frog не поддерживает режим цепей Nose-Hoover'a.
Выше описанным путем получим графики:
_______________________________
В отличие от обычного термостата Nose-Hoover'a, температура поднимается и колеблется вблизи необходимого уровня.400K?
Влияние времени релаксации на поведение системы при поддержании температуры 300 К
Как и в случае термостата Берендсена, необходимо изменить скрипт do.sh:
do.sh
#!/bin/bash ... ... ... ; Velocity generation gen_vel = yes ; задает скорости по Максвеллу gen_temp = 300 ; стартовая температура ... ... ... done
Получите графики:
_______________________________
Гистограмма распределения температур:
Velocity rescaling thermostat
Данный термостат является модифицированным термостатом Берендсена со стохастической прибавкой к кинетической энергии, что обеспечивает ее верное распределение.
Перейдите в директорию /v-resc
cd v-resc
Влияние времени релаксации на поведение системы при нагревании с нуля Кельвинов
Для запуска расчетов так же запускается скрипт do.sh, от предыдущих он отличается следующим:
do.sh
#!/bin/bash l="1 10 100 1000 5000"; #Задает необходимые времена релаксации for a in $l; do cat << _EOF_ > NVTVR${a}ps.mdp #создает файл параметров МД ... ... ... ; Temperature coupling is on tcoupl = V-rescale ; Nose-Hovver thermostat tc-grps = AR ; группа для термостатирования ... ... ... grompp_4.5.4 -f NVTVR${a}ps.mdp -c em.gro -p topol.top -n indx.ndx -o NVTVR${a}ps.tpr #см ниже cat << _EOF_ > NVTVR${a}ps.sh #создает скрипт системы постановки задач #!/bin/sh #PBS -N ARGON_VR ... ... ... mpirun C mdrun_4.5.4_MPI -v -deffnm NVTVR${a}ps ... ... ... done
То есть изменился лишь режим термостатирования и названия файлов.
Выше описанным путем получим графики:
_______________________________
Из графиков видно, что температура поднимается до необходимого уровня по закону, близкому к экспоненциальному, и со временем изменяется хаотично. Таким образом данный термостат позволяет симулировать флуктуации температуры.
Влияние времени релаксации на поведение системы при поддержании температуры 300 К
Как и в случае термостата Берендсена, необходимо изменить скрипт do.sh:
do.sh
#!/bin/bash ... ... ... ; Velocity generation gen_vel = yes ; задает скорости по Максвеллу gen_temp = 300 ; стартовая температура gen_seed = 6 ; этот ключ для рандомайзера ... ... ... done
Получите графики:
_______________________________
Гистограмма распределения температур:
