Основы МД: Шаг интегрирования, радиусы обрезания, сохранение энергии

PDF version

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

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

Цель:

изучение работы базовых алгоритмов молекулярной динамики на примере моделирования инертного газа (аргона).

Расчеты осуществляются при помощи программного пакета Gromacs.

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

Для выполнения задачи необходим компьютер, работающий под OS Linux, либо родственными ОС. Расчеты молекулярной динамики рекомендуется выполнять на высокопроизводительных компьютерных кластерах.

Необходимое ПО :

• Gromacs - осуществляет расчеты молекулярной динамики (в данной статье использована версия 4.5.4) как в обычной, так и в MPI версиях (для распределенных вычислений на кластере)
• Gnuplot - используется для построения графиков (в статье рассмотрена версия 5)
• LAM-MPI - драйвер, необходимый для распределенных вычислений

Для удобства работы рекомендуем использовать файловые менеджеры, например, Midnigt Commander.

Система представляет собой 1000 атомов инертного газа Аргона в коробке 10х10х10 нм.

Описание и подготовка подробно описаны в другой статье.

Cкачайте архив в интересующую вас директорию

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

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

tar -xvf tutor01.tar.gz

Влияние шага интегрирования

Перейдите в директорию istep

 cd tutorial/istep 

Команда ls возвращает следующий список файлов

ls
do.sh       energy.sh  eps_to_img.sh       topol1000.top   
em1000.gro  enter.txt  gnuplotenergy1.scr  

Где

• do.sh - основной скрипт, варьирующий шаг интегрирования.
• energy.sh, eps_to_img.sh, enter.txt, gnuplotenergy1.scr - файлы, необходимые для обработки данных и построения графиков.
• em1000.gro и topol1000.top - файлы, описывающие систему из тысячи атомов аргона после минимизации энергии.

Рассмотрим подробно основной скрипт

do.sh

#!/bin/bash


l="1 2 4 6 8";     #перечисляем интересующие нас шаги интегрирования


for a in $l;
do
let "b = 500000 * 8 / a"        #вычисляет количество шагов для равенства времен
cat << _EOF_ > NVE${a}ps.mdp    #создает файл конфигурации громакс с маской имени
title           = OPLS ARGON NVE equilibration
; Run parameters
integrator      = md            ; leap-frog интегратор
nsteps          = $b            ; количество шагов
dt              = 0.00$a        ; шаг в пикосекундах
; 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		= 1.0		; радиус обрезания Ван-дер-Ваальсовых взаимодействий
;vdwtype	= Switch        ; способ выключения взаимодействий Ван-дер-Ваальса
;rvdw_switch	= 0.7           ; расстояние до выключателя задаваемого в Switch
; Electrostatics
coulombtype     = cut-off       ; способ выключения электростатических взаимодействий
; Temperature coupling is off
tcoupl          = no            ; не термостатирует в NVE ансамбле
; Pressure coupling is off
pcoupl          = no            ; не поддерживает давление в NVE ансамбле
; Periodic boundary conditions
pbc             = xyz           ; задает трехмерную периодическую систему ячеек
; Dispersion correction
DispCorr        = EnerPres      ; если rvdw_switch закомментирован 
; Velocity generation
gen_vel         = yes           ; задать скорости движения атомов по Максвеллу
gen_temp        = 300           ; задать начальную температуру системы
gen_seed        = 7             ; использовать этот ключ для рандомайзера скоростей
_EOF_                           # закрывает файл
grompp_4.5.4 -f NVE${a}ps.mdp -c em1000.gro -p topol1000.top -o NVE${a}ps.tpr # см ниже
cat << _EOF_ > NVE${a}ps.sh     # создает скрипт для системы постановки задач
#!/bin/sh
#PBS -N ARGON${a}
#PBS -l nodes=1:ppn=4           ; один узел, четыре процессора
#PBS -q day
#PBS -V

lamboot -v -ssi boot tm
cd ${PWD}
mpirun C  mdrun_4.5.4_MPI -v -deffnm NVE${a}ps ;запуск молекулярной динамики 
lamhalt -v $PBS_NODEFILE        
_EOF_                           # завершает файл
qsub NVE${a}ps.sh               # ставит задачу в очередь
done

Точка с запятой обозначают комментарий в файле конфигурации, решетка - комментарий в скрипте, таким образом, параметры vdwtype и rvdw_switch изначально закомментированы.

У команды grompp_4.5.4, создающей бинарные файлы для запуска молекулярной динамики (mdrun), есть несколько ключей

-f файл на вход с параметрами МД (в данном случае файл конфигурации)
-c файл структуры системы
-p файл топологии
-o создаваемый бинарный файл

Если вы не используете систему постановки задач, закомментируйте с помощью # все после grompp и до done, далее запускайте каждую обработку по отдельности с помощью mdrun.

Запустите скрип do.sh

./do.sh

Теперь в директории появится россыпь файлов вида NVE(цифра)ps.расширение.

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

Например, для шага интегрирования 2:

mdrun_4.5.4 -v -deffnm NVE2ps

В ином случае - проверьте с помощью команды qstat, запустились ли ваши задачи.

Если все хорошо - можно подождать или приступить к выполнению второй части задачи.

После завершения рассчетов запустите скрипт energy.sh

./energy.sh

Рассмотрим его подробно:

energy.sh

#!/bin/bash

l="1 2 4 6 8";  #для тех же шагов интегрирования


for a in $l;
do
g_energy_4.5.4 -f NVE${a}ps.edr -o energyNVE${a}ps.xvg < enter.txt # выполняет
#используя значения из файла enter.txt
done

Файл enter.txt содержит в себе одну цифру и два перехода на новую строку, что задает расчет полной энергии, если удалить файл или закомментировать его прикрепление, то можно задавать их вручную.

Скрипт energy.sh создаст файлы вида energyNVE(шаг интегр)ps.xvg, их можно просмотреть в текстовом редакторе или с помощью xmgrace.
Например, для шага 8:

 xmgrace energyNVE8ps.xvg

Постройте один общий график, используя gnuplot

gnuplot5 gnuplotenergy1.scr

Скрипт gnuplotenergy1.scr обрабатывает только файлы вышеупомянутого вида с заданными шагами интегрирования и на выходе дает файл energy1.eps

Далее, используя скрипт eps_to_img.sh, конвертируем eps файл в tif

./eps_to_img.sh

(Так же можно изменить разрешение и формат файла, изменив eps_to_img.sh)

Теперь можно просмотреть получившийся график

если вы работаете в Linux или в Windows с запущенным X сервером:

display energy1.tif

можно скачать файл, например, для нашего кластера команда имеет такой вид:

scp student@biosim.moldyn.org:/путь до файла/energy1.tif /куда скачать/

Должен получиться график такого вида:

Откуда видно, что с ростом шага увеличивается и быстрее накапливается ошибка интегрирования.

Но накопление ошибки и изменение энергии системы в симуляции микроканонического ансамбля связано не только с ошибкой интегрирования, но и с заведомо неверными параметрами конфигурации. Мы использовали простое обрезание Ван-дер-Ваальсовых взаимодействий, что вызывает резкий скачок потенциала и бесконечный скачок силы взаимодействия на границе обрезания. Исправим этот недостаток, изменив файл do.sh и раскомментировав параметры vdwtype и rvdw_switch.

 vi do.sh
...
vdwtype	= Switch       
rvdw_switch = 0.7
...

Для того чтобы выйти с сохранением из редактора vi, надо в командном режиме ввести :x

Далее проделайте все операции аналогично первому случаю.

В итоге вы получите график такого вида:

Откуда видно, что в целом энергия системы сохраняется при всех значениях шага интегрирования, но с увеличением шага увеличивается дисперсия энергии.

Влияние радиуса обрезания взаимодействий Ван-дер-Ваальса

Перейдем в директорию cutoffrange

cd ../cutoffrange

Вторая часть задачи выполняется так же, как и первая, исключая изменение типа обрезания взаимодействий Ван-дер-Ваальса. По умолчанию используется Cut-off. В нашем случае это не критично, так как мы хотим наблюдать изменения в расчетах только за счет изменения радиуса обрезания.

Рассмотрим скрипт do.sh

do.sh


#!/bin/bash

l="0.4 0.8 1 1.5 2";#   перечисляет интересующие нас шаги интегрирования


for a in $l;
do
cat << _EOF_ > NVE${a}coff.mdp    #создает файл конфигурации громакс с маской имени
title           = OPLS ARGON NVE equilibration
; Run parameters
integrator      = md            ; leap-frog интегратор
nsteps          = 500000        ; количество шагов
dt              = 0.002         ; шаг в пикосекундах
; Output control
nstxout         = 100           ; сохранять координаты каждые 100 шагов
nstvout         = 0              ; не сохранять скорости
nstenergy       = 100           ; сохранять энергии каждые 100 шагов
nstlog          = 100           ; обновлять лог каждые 100 шагов
; Neighborsearching
ns_type         = grid          ; искать соседей по сетке
nstlist         = 5             ; количество шагов
rlist		= $a		; радиус обрезания построения списка соседей
rcoulomb	= 2.0		; радиус обрезания электростатических взаимодействий
rvdw		= $a		; радиус обрезания Ван-дер-Ваальсовых взаимодействий
; Electrostatics
coulombtype     = cut-off       ; способ выключения электростатических взаимодействий
; Temperature coupling is off
tcoupl          = no            ; не термостатирует в NVE ансамбле
; Pressure coupling is off
pcoupl          = no            ; не поддерживает давление в NVE ансамбле
; Periodic boundary conditions
pbc             = xyz           ; задает трехмерную периодическую систему ячеек
; Dispersion correction
DispCorr        = EnerPres      ; если rvdw_switch закомментирован 
; Velocity generation
gen_vel         = yes           ; задать скорости движения атомов по Максвеллу
gen_temp        = 300           ; задать начальную температуру системы
gen_seed        = 7             ; использовать этот ключ для рандомайзера скоростей
_EOF_                           # закрывает файл
grompp_4.5.4 -f NVE${a}coff.mdp -c em1000.gro -p topol1000.top -o NVE${a}coff.tpr -maxwarn 2 
# см ниже
cat << _EOF_ > NVE${a}coff.sh # создает скрипт для системы постановки задач
#!/bin/sh
#PBS -N ARGON_${a}_co
#PBS -l nodes=1:ppn=4    ; один узел, четыре процессора
#PBS -q day
#PBS -V

lamboot -v -ssi boot tm
cd ${PWD}
mpirun C  mdrun_4.5.4_MPI -v -deffnm NVE${a}coff ;запуск молекулярной динамики 
lamhalt -v $PBS_NODEFILE
_EOF_                           # завершает файл
qsub NVE${a}coff.sh             # ставит задачу в очередь
done

К команде grompp_4.5.4 добавился ключ -maxwarn 2, который позволяет программе игнорировать до 2х предупреждений.

После завершения вычислений, выполнять задачу так же как и ранее:

./energy.sh
...
gnuplot5 gnuplotenergy1.scr
...
./eps_to_img.sh

Должен получиться график вида:

Из которого видно, что с уменьшением радиуса обрезания взаимодействий энергия сильно изменяется. Это легко объяснить тем, что взаимодействия между атомами на меньших расстояниях довольно сильны, но мы ими пренебрегаем. При радиусе обрезания больше одного нанометра энергии неизменны и практически совпадают, это свидетельствует о том, что взаимодействия на этих расстояниях пренебрежительно малы.