Работа 2

Материал из md
Перейти к: навигация, поиск
Изучение влияния величины шага ∆t и радиуса обрезания потенциала на точность моделирования

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

Цель

На примере моделирования системы атомов с потенциалом Леннарда-Джонса изучить влияние шага интегрирования ∆t и величины радиуса обрезания потенциала rcut на точность молекулярно-динамического расчёта.

Задания

∆t:

  1. Провести МД-моделирование леннард-джонсовской жидкости при значениях шага интегрирования ∆t = 0.001, 0.002, 0.005, 0.01, 0.015, 0.02. Для каждого случая использовать свое число шагов моделирования nsteps, чтобы время эволюции было одинаковым для каждой модели и составляло 50 единиц времени.
  2. Найти время Т, затраченное на расчёт каждой модели для указанных значения ∆t.
  3. Найти отношения флуктуаций энергии R = Δ Etot / Δ Ekin.
  4. Построить зависимости полученных T и R от шага интегрирования ∆t.

rcut:

  1. Провести моделирование леннард-джонсовской системы при разных значениях радиуса обрезания rcut = 1.5, 2.0, 2.5, 3.0, 3.5, 4. Для каждой модели использовать шаг интегрирования ∆t = 0.002 и продолжительность моделирования, соответствующую 50 единицам времени.
  2. Найти время T, затраченное на расчёт моделей для каждого значения rcut .
  3. Найти отношения флуктуаций энергии R = Δ Etot/Δ Ekin.
  4. Построить зависимость полученных T и R от величины rcut.

Скрипт (дополнительно):

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

Вопросы

  1. Оценить оптимальный шаг интегрирования ∆t.
  2. Оценить оптимальную величину rcut.
  3. Как изменяется полная энергия системы, если шаг интегрирования выбран большим?

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

Метод молекулярной динамики используется для моделирования поведения большого числа взаимодействующих частиц. С точки зрения механики это обычная консервативная система, в которой должны сохраняться полная энергия (Etot), импульс и момент импульса для всей системы. С термодинамической точки зрения она соответствует микроканоническому ансамблю (NVE), т.е. системе с постоянным числом частиц, постоянным объёмом и постоянной полной энергией. Полная энергия складывается из потенциальной и кинетической: Etot = Ekin + Epot , где обе составляющих флуктуируют со временем. Кинетическая энергия Ekin, связывается с температурой системы. Масштаб этих флуктуаций будем обозначать как ΔEkin.

Уравнения движения для многочастичной системы решаются численно, поэтому очень важен правильный выбор величины шага интегрирования ∆t. Основным критерием выбора шага в данном случае является сохранение полной энергии. Если шаг излишне большой, то поведение системы становится нефизическим, при этом, как правило, наблюдается возрастание полной энергии системы Etot. Для более точного решения уравнений движения всегда можно уменьшить шаг, но тогда увеличится продолжительность расчёта. Поэтому приходится находить компромисс между точностью расчёта и временем, затраченным на моделирование.

Вместо "механистического" ансамбля NVE в молекулярной динамике используют также "статистические" ансамбли – NVT и NPT. В первом поддерживаются постоянными объём и температура, во втором – давление и температура. Полная энергия в этом случае не является интегралом движения, а флуктуирует около некоторого среднего значения. Понятно, что для описания движения атомов в таких ансамблях классических уравнений аналитической механики не достаточно. Нужны дополнительные ухищрения, имитирующие наличие термостата и/или баростата. Однако в любом случае, большой шаг ∆t будет приводить к нефизическому поведению среднего значения полной энергии и масштаба флуктуаций ΔEtot и ΔEkin .

Для проверки стабильности решения уравнений движения можно использовать простой приём, общий для разных ансамблей. Предлагается следить за флуктуациями полной и кинетической энергий и вычислять их отношение R=ΔEtot /ΔEkin. Если расчёт идёт без существенных ошибок, то это отношение должно быть постоянным, и сохраняться при небольших вариациях шага интегрирования ∆t. Если отношение R меняется, то это убедительное свидетельство того, что траектория рассчитывается с ошибкой.

Многие задачи требуют долгого моделирования (продолжительной "жизни" системы на больших временах). Это нужно, например, при изучении фазовых переходов, конформаций макромолекул, моделировании самоорганизации молекул в растворах. Один из способов ускорения моделирования – сокращение времени на расчёт межатомных взаимодействий. Если учитывать все взаимодействия, то для системы из N атомов трудоемкость оценивается как О(N2), т.е. возрастает пропорционально квадрату числа атомов. Это означает, что увеличение системы в 10 раз с неизбежностью замедлит моделирование в 100 раз. Однако если потенциал достаточно короткодействующий, как в случае ван-дер ваальсовских взаимодействий (которые обычно описываются потенциалом Леннарда-Джонса), то взаимодействие становится пренебрежимо малым уже на расстояниях два-три диаметра атома. Для «отсечения» далеких атомов используется параметр rcut , т.е. взаимодействие рассчитывается только между атомами, которые расположены не далее, чем на указанном расстоянии rcut друг от друга. Число таких атомов Ncut сравнительно небольшое (десятки-сотня) и не зависит от полного числа атомов N в системе. Поэтому с увеличением размера модели время расчёта будет возрастать не так быстро. Выбор величины rcut также является результатом компромисса. С одной стороны, хотелось бы ограничиться как можно меньшим числом соседей и не тратить времени на расчёт взаимодействия с далекими атомами, но с другой стороны, надо быть уверенным, что мы учитываем всех соседей, оказывающих влияние на данный атом, поскольку недоучет влияния окружения также приведёт к неправильному поведению моделируемой системы, даже при малом значении шага ∆t. Для контроля правильности расчёта здесь также можно использовать параметр R, введённый выше.

Заметим, что использование параметра rcut ускоряет моделирование, но не обеспечивает, как может показаться, линейной зависимости времени расчёта от числа N. Оно также возрастает пропорционально N2, хотя намного медленнее. Здесь мы по-прежнему перебираем все пары атомов и вычисляем расстояние между ними. Выигрыш только в том, что теперь взаимодействие (силы) рассчитываются лишь для атомов, расстояние между которыми удовлетворяет критерию rcut. Заметим, что в настоящее время существуют алгоритмы молекулярной динамики, имеющие эффективность линейную с числом N. Без них было бы невозможно работать с системами, содержащими десятки тысяч атомов и более. Однако для этого используются дополнительные приемы.

Данная работа основана на демонстрационной разработке David van der Spoel, одного из соавторов пакета Gromacs:

https://extras.csc.fi/chem/courses/gmx2004/exercises/

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

Данная работа содержится в директории 02_timestep/, где для изучения влияния шага моделирования и радиуса обрезания выделены отдельные поддиректории: timestep/ и rcutoff/.

Изучение влияния ∆t

Перейти в папку timestep/:

cd 02_timestep/timestep/

Посмотрим содержимое папки:

ls

В ней лежат три файла, необходимые для моделирования: conf.gro, grompp.mdp, и topol.top. Они соответствуют системе из 216 леннард-джонсовских атомов. Посмотрим файл структуры conf.gro с помощью команды less (для выхода из программы less достаточно нажать клавишу q):

less conf.gro

В данной работе мы используем этот файл без изменений.

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

less topol.top

Его содержимое:

[ defaults ]
; nbfunc	comb-rule
1		3

[ atomtypes ]
; full atom descriptions are available in ffoplsaa.atp
; name   bond_type  mass   charge   ptype      sigma    epsilon
  AR 	   AR	   1	     0		A	1	1

[ molecule_type ]
Argon   1

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

[ system ]
; Name
Argon

[ molecules ]
; Compound #mols

Argon 216

Леннард-джонсовский атом здесь именуется аргоном. Однако это может быть любой другой благородный газ, так как в данной работе используются приведённые единицы. В них считается, что глубина потенциальной ямы LJ-потенциала ε(epsilon), «диаметр атома» σ (sigma) и масса (mass) равны единице. Время также измеряется в приведенных единицах. (Подчеркнем, приведенные параметры леннард-джонсовского потенциала удобны для теоретических исследований, но не имеют особого значения для численных расчётов. Поэтому в пакете Gromacs, как правило, используются размерные единицы: расстояние измеряется в нанометрах, а время – в пикосекундах. Поэтому не удивляйтесь, если на стандартных графиках по оси времени будут подписаны пикосекунды. В данной работе время выводится в приведённых единицах. Ниже мы обсудим связь приведенного времени с пикосекундами подробнее).

Посмотрите файл с параметрами моделирования:

less grompp.mdp

Показаны лишь некоторые его строки. Жирным выделены те, в которые будем вносить изменения:

; RUN CONTROL PARAMETERS
dt = 0.002
nsteps = 25000
;
.................................................................
..................................................................
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 10
; ns algorithm (simple or grid)
ns-type = grid
; Periodic boundary conditions: xyz, no, xy
pbc = xyz
; nblist cut-off 
rlist = 2.5
...................................................................
...................................................................
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = Cut-off
coulomb-modifier         = Potential-shift-Verlet
rcoulomb-switch = 0
rcoulomb = 2.5
;  Relative dielectric constant for the medium and the reaction field
epsilon-r = 1
epsilon-rf               = 0 
; Method for doing Van der Waals
vdw-type = Cut-off
vdw-modifier             = Potential-shift-Verlet
; cut-off lengths 
rvdw-switch = 0
rvdw = 2.5

Прежде всего будем менять значения dt и nsteps. Предлагается для каждого значения dt создавать свою директорию и в ней работать. Например, для dt=0.001 действуем следующим образом:

mkdir 001– создаём поддиректорию 001 в директории timestep.

cp * 001/– копируем все 3 файла из рабочей директории timestep в директорию 001.

cd 001/ переходим в директорию 001:

Теперь в файле параметров grompp.mdp задаём значения параметров: dt=0.001; nsteps = 50000 (чтобы эволюция модели соответствовала 50-ти единицам времени для данного шага). Для этого можно использовать любой текстовый редактор.

Запускаем препроцессинг:

gmx grompp

Этой командой мы собираем файлы со стандартными именами conf.gro, grompp.mdp и topol.top в стартовый файл topol.tpr. Убедимся, что нужный файл появился, просмотрев содержимое директории:

ls

Файл topol.tpr – бинарный, чтобы посмотреть его, нужно использовать программу gmxdump:

gmx dump -s topol.tpr | less(для выхода использовать клавишу «q»)

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

gmx mdrun –nt 1 -v

Программа по умолчанию использует полученный стартовый файл со стандартным именем topol.tpr. Ключ –nt с параметром 1 означает, что мы запускаем расчёт на одном ядре (без распараллеливания), использование ключа -v позволяет выводить на экран текущий номер шага интегрирования и оставшееся время расчёта. В конце работы программа выдаёт на экран краткий отчёт о работе, в том числе, полное время T, затраченное на расчёт, также это время сохраняется в файле md.log.

После завершения работы программы появляются четыре файла:

traj.xtc– траектория системы, содержащая координаты атомов в последовательные моменты времени (бинарный файл),

confout.gro– конечная конфигурация,

ener.edr– файл «энергии» (содержит кинетическую, потенциальную и полную энергии, а также ряд других характеристик, рассчитываемых по ходу моделирования),

md.log– log-файл. Для нас сейчас важно, что кроме параметров расчёта в нем записано время расчёта Т.

С помощью программы gmx energy извлечём полную и кинетическую энергии, а также их флуктуации (среднеквадратичный разброс) RMSD (root mean square deviation):

gmx energy -w

При запуске программы нужно будет выбрать (в интерактивном режиме) интересующие нас характеристики. Можно выбрать несколько характеристик. После ввода первой мы нажимаем клавишу Enter, затем вводим название или номер второй характеристики, и т.д. Дважды нажатый Enter говорит о том, что мы выбрали всё, что нам хотелось. Таким образом, в нашем случае выполняем: Total (Enter) Kinetic (Enter Enter). После этого программа генерирует файл и рисует график. Можно набирать не слово, а соответствующую ему цифру: 6 (для полной энергии) и 5 (для кинетической). Ключ -w указывает, что график автоматически откроется в графической программе xmgrace. Программа выдаёт на экран средние значения потенциальной и кинетической и энергий, а также их RMSD, которые мы используем в качестве ΔEtot и ΔEkin. Найденные величины флуктуаций (вместе с выше-найденной величиной T) следует записать на листе бумаги, или в отдельный файл, или сразу в таблицу программы qtiplot, с помощью которой мы рассчитаем параметр R =ΔEtot /ΔEkin и нарисуем графики (о программе qtiplot см. Вводную часть). После этого возвращаемся в директорию timestep и переходим к созданию модели с другим значением шага интегрирования.

Создаем директорию 002/ и, как было сделано ранее, переходим в неё, предварительно скопировав в неё файлы из timestep . Повторяем всю работу, проделанную с предыдущим значением шага и получаем соответствующие значения R и T.

Аналогично выполняем работу с остальными заданными значением dt (0.005, 0.01, 0.015, 0.02), создав для каждого свою директорию.

Собрав значения параметра R для разных шагов интегрирования, нарисуем его зависимость от величины шага. Должен получиться график, подобный показанному на Рис. 2.1. Аналогично надо построить зависимость от dt для времени расчёта T, Рис. 2.2.

Fig2.1.png

Рис. 2.1. Зависимость параметра R=ΔEtot/ΔEkin от величины шага интегрирования Δt. Флуктуации полной энергии ΔEtot резко возрастают по сравнению с флуктуациями кинетической энергии с увеличением шага Δt.

Fig2.2.png

Рис. 2.2. Зависимость времени моделирования от величины шага интегрирования Δt. Чем меньше шаг, тем медленнее расчёт.

Для работы с графиками удобно использовать программу qtiplot – аналог программы Origin. qtiplot входит в стандартную комплектацию Linux. Как альтернативу, можно использовать программу из Openoffice, аналог Excel.

Изучение влияния rcut

Перейти в папку rcutoff/. В ней лежат те же три файла, необходимые для моделирования: conf.gro, grompp.mdp и topol.top. Работа выполняется аналогичным образом, как с шагом по времени Δt.

Для значения rcut=1.5 создаем директорию, например, с именем 15, копируем в неё файлы из директории rcutoff и переходим в нее:

mkdir 15; cp conf.gro grompp.mdp topol.top 15/; cd 15.

(Напомним, что, с использованием символа ";" (точка с запятой), несколько команд может быть записано в одной строке).

Параметр rcut следует менять в файле grompp.mdp (см выше) в трёх местах: как для ван-дер ваальсовского (rvdw), так и для электростатического (rcoulomb) взаимодействий (не смотря на то, что в данной задаче нет зарядов), а также для списка соседей (rlist). Построение списка соседей Верле (Verlet list) всегда используется в пакете Gromacs с целью ускорения расчётов. В остальном надо следовать указаниям, описанным для работы с ∆t. На Рис. 2.3 и Рис. 2.4 показано, что приблизительно должно получиться.

Fig2.3.png

Рис. 2.3. Зависимость параметра R=ΔEtot/ΔEkin от величины радиуса обрезания rcut. Флуктуации полной энергии возрастают по сравнению с флуктуациями кинетической энергии при уменьшении rcut.

Fig2.4.png

Рис. 2.4. Зависимость времени моделирования от величины радиуса обрезания rcut. Чем больше радиус обрезания, тем медленнее расчёт.

Скрипт (дополнительно):

Можно не создавать каждую модель вручную, а написать скрипт, автоматизирующий эту работу. Ниже приведён пример такого скрипта для разных значений шага интегрирования и записи их в разные директории с указанными именами. Скрипт снабжён комментариями, которые приводятся после знака решётки #. Необходимая информация об используемых командах приведена во Введении

timestep.sh:
#!/bin/bash
# this script performs simulation of argon with different timesteps. It is needed to prepare template for grompp.mdp to substitute dt and n. This template is named grompp.mdb_tmpl here.
for dt in 0.001 0.002 0.005 0.01 0.015 0.02                 # dt – integration timestep
do
  n=$(echo 50 / $dt | bc)                                # use bc-calculator for computing the number of timesteps, corresponding 50ps of total time
  cat grompp.mdp_tmpl |  sed -e "s/@dt@/$dt/; s/@n@/$n/" > grompp.mdp      # making mdp-file with different timesteps (first you should prepare template for grompp.mdp). Here comamd sed replases @dt@ with the current value of dt. And replases @n@ with value of n.
  mkdir $dt
  cd $dt
     cp ../conf.gro .
     cp ../topol.top .
     cp ../grompp.mdp .
     gmx grompp -v
     gmx mdrun -nt 1 -v
     echo '5 6' | gmx energy -f ener.edr -o etot_ekin.xvg -b 1 > output.ener      # the output report is sent to file output.ener where we can find desired fluctuations later
  cd ../
done

# this section is just optional. Here we'll calculate R-factor and collect all numbers in file "r_factor.txt" and times in "count_time.txt"
echo "Dtot   Dkin   R=Dtot/Dkin" > r_factor.txt
echo "delta_t   count_time" > count_time.txt
for dt in 0.001 0.002 0.005 0.01 0.015 0.02                 # dt – integration timestep
do
  cd $dt                # awk is a handy language for manipulating strings
     cat output.ener | awk '/^Total Energy/ {Dtot=$5}; /^Kinetic En/ {Dkin=$5}; END {print Dtot, Dkin, Dtot/Dkin}'>> ../r_factor.txt    # output energies
     cat md.log | awk '/delta_t/ {dt=$3}; /Time:/ {count_t=$3}; END {print dt, count_t}' >> ../count_time.txt                           # output count-time
  cd ../
done

Приводим также пример скрипта для создания моделей с разными значениями rcut.

cutoff.sh:
for r in 1.5 2.0 2.5                 # r – cutoff
do
  cat grompp.mdp_tmpl |  sed -e "s/@r@/$r/" > grompp.mdp      # making mdp-file with different timesteps
  mkdir cutoff_$r
  cd cutoff_$r/
  cp ../conf.gro .
  cp ../topol.top .
  mv ../grompp.mdp .
  gmx grompp -v >& output.grompp
  gmx mdrun -v
  cd ../
done

echo "Dtot   Dkin   R=Dtot/Dkin" > r_factor.txt
echo "delta_t   count_time" > count_time.txt
 for r in 1.5 2.0 2.5                 # r – cutoff
 do
   cd cutoff_$r/
     echo '5 6' | gmx energy -f ener.edr -o etot_ekin.xvg >& output.ener
     cat output.ener | awk '/^Total Energy/ {Dtot=$5}; /^Kinetic En/ {Dkin=$5}; END {print Dtot, Dkin, Dtot/Dkin}'>> ../r_factor.txt    # output energies
     cat md.log | awk '/delta_t/ {dt=$3}; /Time:/ {count_t=$3}; END {print dt, count_t}' >> ../count_time.txt                           # output count-time
   cd ../
 done