Работа 1

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

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

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

Цель

Познакомиться с пакетом Gromacs, приобрести необходимые навыки для работы с ним.

Задания

  1. Изучить основные команды и файлы пакета Gromacs.
  2. Выполнить стандартное демонстрационное задание.

Вопросы

  1. Назначение основных программ, реализующих молекулярно-динамическое моделирование в пакете Gromacs. Порядок запуска программ, «препроцессинг».
  2. Основные файлы, их стандартные имена, форматы записи.

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

Gromacs – это пакет программ, используемый для молекулярно-динамического моделирования и анализа полученных моделей. Он разрабатывался для работы с биологическими молекулами (белки, липиды, нуклеиновые кислоты), в которых много атомов, химических связей и разных взаимодействий. Однако его можно использовать и для более простых систем, для моделирования жидкостей, растворов, а также твердой фазы. Популярность Gromacs’а для научно-исследовательской работы объясняется тем, что он работает быстрее многих других известных пакетов, имеет широкие возможности для анализа получаемых МД-моделей, содержит подробное описание и обратную связь с разработчиками. Все это можно найти по ссылкам:

http://www.gromacs.org

http://manual.Gromacs.org/

http://manual.Gromacs.org/online/getting_started.html

Данный практикум проводится на многоядерных персональных компьютерах с установленным пакетом Gromacs в операционной системе Linux. Файлы с параметрами используемых полей сил расположены, как правило, в директории /usr/share/Gromacs/top/, а исполняемые файлы в директории /usr/bin/.

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

Чтобы начать работу, нужно войти в систему под именем student, набрав его в окне-приглашении операционной системы. После входа в систему необходимо запустить «терминал», программу эмулирующую работу с командной строкой. В результате мы оказываемся в директории

/home/student/

В которой расположены директории для выполнения заданий. Номера директорий соответствуют номеру работы. Для выполнения данной работы нужно перейти в директорию 01_demo/. Выполнить демонстрационную работу по моделированию воды, см. ниже раздел "Демонстрационная работа", где описано, как это делается. Однако сначала следует познакомиться с общей схемой работы Gromacs'а, используемыми файлами и основными программами.

Общая схема работы Gromacs

Моделирование

Основной командой, запускающей молекулярно-динамическое (МД) моделирование является mdrun . Она читает «стартовый файл» (run input file), который по умолчанию имеет имя topol.tpr. В нём собрана вся информация, необходимая для проведения моделирования. По завершению молекулярно-динамического расчёта создается несколько выходных файлов. Главный из них – «файл траектории» (trajectory file), который содержит мгновенные конфигурации ("мгновенные снимки", snapshots) системы в последовательные моменты времени. Файл траектории может задаваться двумя способами, т.е. могут быть два файла, с именами traj.trr и traj.xtc. Первый содержит координаты атомов и скорости, а второй имеет более компактный формат (бинарный) и не содержит скоростей. Полезным является «файл энергии» (energy file) с именем ener.edr, содержащий энергию, температуру, давление и ряд других характеристик системы, рассчитываемых в процессе моделирования. Эти имена файлов являются "стандартными", т.е. могут использоваться "по умолчанию". Во многих случаях это позволяет существенно упростить работу.

Для создания стартового файла (topol.tpr) нужно провести «препроцессинг» (preprocessing) Это делает программа grompp, она объединяет все необходимые для работы данные в один файл:

  • координаты и скорости атомов моделируемой системы, находящиеся обычно в файле «молекулярной структуры» (molecular structure file) со стандартным именем conf.gro. Для краткости мы будем его также называть «файл струкуры».
  • информация о топологии молекулы и используемых взаимодействиях, находящаяся в файле «молекулярной топологии» (molecular topology file, для краткости — «файл топологии») с именем topol.top. Обратим внимание, что это совсем не то же самое, что стартовый файла topol.tpr, а отличается от него, хотя бы, своим расширением.
  • общие параметры моделирования, находящиеся в файле «молекулярно-динамических параметров» (molecular dynamics parameter file, для краткости, «файл параметров») с именем grompp.mdp.

Создание файлов для препроцессинга (conf.gro, topol.top, grompp.mdp) является отдельной работой. Они могут быть созданы независимо, с помощью каких-либо специальных программ, или же вручную. В некоторых случаях удобно пользоваться пограммами Gromacs для подготовки этих файлов. Например, если координаты молекулы записаны в формате PDB, то с помощью программы pdb2gmx можно сформировать для этой молекулы файл топологии (topol.top), используя готовое поле сил из библиотеки пакета Gromacs.

Работа с молекулярно-динамической моделью

Полученная в результате моделирования МД-траектория (файлы traj.trr или traj.xtc) представляет собой последовательные мгновенные конфигурации модели, т.е. координаты (и скорости) всех атомов в последовательные моменты времени. Их можно визуализировать, наблюдая поведение системы со временем, программой ngmx, входящей в пакет Gromacs. С её помощью можно запустить анимацию (последовательное отображение кадров на экране монитора), вращать модельную систему относительно осей координат. Однако больше возможностей для наблюдения и анализа модели даёт пакет VMD[1] (см. Error: Reference source not found), разработанный специально для визуализации молекулярно-динамических расчётов.

Основные характеристики модели, рассчитываемые по МД-траектории, записаны в файле энергии ener.edr. Посмотреть и извлечь эти характеристики можно c помощью программы g_energy. Она работает в интерактивном режиме, позволяя выбрать интересующую нас характеристику модели, а также изобразить её поведение во времени в виде графика с помощью программы xmgrace.

Для структурного и динамического анализа МД-траектории пакет Gromacs имеет обширную библиотеку программ различного назначения. В частности, будут использоваться следующие программы:

g_rdf – расчёт парной корреляционной функции g(r),
g_velacc – расчёт автокорреляционной функции скорости,
g_msd – расчёт среднеквадратичного смещения и коэффициента самодиффузии,
g_hbond – расчёт числа водородных связей,
g_density – расчёт плотности,
g_gyrate – расчёт радиуса гирации молекулы.

По мере выполнения работ данного практикума мы научимся работать со всеми этими программами.

Следующие два раздела («Основные файлы» и «Основные программы») следует рассматривать как справочные материалы, к которым можно возвращаться по мере необходимости при выполнении всех последующих заданий.

Основные файлы

Рассмотрим назначение и структуру основных файлов пакета Gromacs. В отличие от Linux’a, здесь существенно используются расширения файлов, кроме того, существуют стандартные имена, используемые по умолчанию.

Файл молекулярной структуры (molecular structure file) (.gro,.pdb)

Содержит исходную конфигурацию моделируемой системы, которой может быть как одна молекула, например, молекула воды или белка, так и вся система целиком, например, макромолекула и окружающие её молекулы растворителя. Файл содержит имена всех атомов (центров взаимодействия) в данной системе, их координаты и скорости, а также размер модельного бокса. По умолчанию для него используется имя conf.gro. Однако важным здесь является расширение, а имя может быть любым, например, water.gro или protein.gro. Пакет Gromacs может использовать в качестве файла молекулярной структуры файлы формата .pdb (Protein Data Bank), который широко используется для записи структуры белков. Однако файл с расширением .gro может содержать скорости атомов, в отличие от файлов .pdb, в которых скоростей нет.

Ниже показано содержимое файла формата .gro, для системы из двух молекул воды. (Заметим, что в данном примере используется трехцентровая модель воды).

MD of 2 waters, t= 0.0
6
1WATER OW1 1 0.126 1.624 1.679  0.1227 -0.0580  0.0434
1WATER HW2 2 0.190 1.661 1.747  0.8085  0.3191 -0.7791
1WATER HW3 3 0.177 1.568 1.613 -0.9045 -2.6469  1.3180
2WATER OW1 4 1.275 0.053 0.622  0.2519  0.3140 -0.1734
2WATER HW2 5 1.337 0.002 0.680 -1.0641 -1.1349  0.0257
2WATER HW3 6 1.326 0.120 0.568  1.9427 -0.8216 -0.0244
1.82060 1.82060 1.82060

Строки содержат следующую информацию (сверху вниз):

1 – информационная строка (текстовый формат)
2 – полное число атомов (точнее, центров взаимодействия) (формат integer),
3-8 – информация об используемых атомах, их координаты и скорости. Для каждого атома своя строка (см. ниже),
9 – длины ребер модельного бокса в нанометрах (свободный формат).

Столбцы для строк 3-8 (слева направо):

1 – номер "аминокислотного остатка". Он занимает 5 позиций, начиная с первой (формат integer). Термин «остаток» (residue) используется исторически, поскольку изначальной целью было моделирование белков. В общем случае он нумерует заданную группу (центров взаимодействия). Здесь их две, 1 и 2 гуппы, относящиеся к первой и второй молекулам воды. В каждой молекуле по три центра взаимодействия, координаты и скорости которых записаны в отдельных строках. Поэтому числа 1 и 2 повторяются в первом столбце по три раза.
2 – имя остатка (5 букв). Наши группы одинаковы и обозначаются именем WATER.
3 – имя атома (5 букв). Молекула воды здесь имеет три центра взаимодействия, соответствующих атому кислорода и двум атомам водорода, которые обозначены как OW1, HW1, HW2. По этим именам устанавливается специфика атомов, т.е. соответсвующие им параметры взаимодействия, записанные в других файлах, см. ниже.
4 – номер атома (центра взаимодействия) (5 позиций, integer). Все атомы пронумерованы в том порядке, как они выписаны в файле.
5–7 координаты x y z (в нанометрах) каждая содержит по 8 цифр с тремя знаками после запятой.
8–10 скорости (в нм/пс (или км/с), каждая по 8 цифр с 4 знаками после запятой.

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

Файл топологии молекулы (molecular topology file) (.top)

Файл топологии содержит полное описание всех взаимодействий в моделируемой молекуле: как, собственно, топологию (связность атомов) молекулы, так и используемые взаимодействия. По умолчанию для него используется имя topol.top. Он является текстовым и содержит следующие основные разделы:

[defaults] – здесь указывается тип взаимодействий между «несвязанными» атомами (атомами разных молекул и относительно далекими атомами внутри одной молекулы), например, потенциал Леннард-Джонса или Букингема, а также правила комбинирования параметров потенциала при взаимодействии атомов разного типа.
[atomtypes]– содержит имена атомов, их заряд, атомный вес, тип центра взаимодействия (атом обозначается буквой A, оболочка – S, виртуальный центр – V), параметры леннард-джонсовского потенциала С6 и С12 (или же σ, ε, в зависимости от указанного в разделе [defaults] типа взаимодействия и правил комбинирования).
[moleculetype] указывается имя молекулы, описание которой приводится ниже.
[atoms] – указывается, как имена данных атомов соотносятся с типами атомов данного поля сил, а также приписанные им заряды.
[bonds] – выписаны номера атомов, связанных ковалентной связью, а также тип связи, длина и жёсткость (если они отличаются от стандартных для данного поля сил).
[pairs] –выписаны пары атомов, разделённые тремя последовательными связями. Они расположены достаточно далеко друг от друга, поэтому взаимодействие между ними считается как между несвязанными атомами, но учитывается частично (так называемое 1 – 4 взаимодействие).
[angles] – выписаны тройки атомов, определяющие угол между соседними ковалентными связями (деформационное взаимодействие), а также тип потенциала, величина и жёсткость угла.
[dihedrals] – четвёрки атомов, составляющих двугранный угол (торсионное взаимодействие), тип потенциала, величина и жёсткость.
[system] – имя системы.
[molecules] – имя молекулы и количество таких молекул.

Файл также может содержать оператор #include , который подключает, например, файл поля сил или файл топологии отдельного остатка ( .itp), или даже целые разделы.

Ниже показан пример файла топологии для молекулы мочевины в воде, взятый из описания Gromacs. (Точка с запятой в первой позиции означает, что данная строка является комментарием).

;Example topology file

; The force field files to be included
#include "gmx.ff/forcefield.itp"

[ moleculetype ]
; name nrexcl
Urea 3

[ atoms ]
; nr type resnr residu atom cgnr  charge
  1   C    1     UREA   C1   1    0.683
  2   O    1     UREA   O2   1   -0.683
  3  NT    1     UREA   N3   2   -0.622
  4   H    1     UREA   H4   2    0.346
  5   H    1     UREA   H5   2    0.276
  6  NT    1     UREA   N6   3   -0.622
  7   H    1     UREA   H7   3    0.346
  8   H    1     UREA   H8   3    0.276

[ bonds ]

; ai aj funct    c0           c1
  3  4    1   1.000000e-01 3.744680e+05
  3  5    1   1.000000e-01 3.744680e+05
  6  7    1   1.000000e-01 3.744680e+05
  6  8    1   1.000000e-01 3.744680e+05
  1  2    1   1.230000e-01 5.020800e+05
  1  3    1   1.330000e-01 3.765600e+05
  1  6    1   1.330000e-01 3.765600e+05

[ pairs ]

; ai aj funct   c0            c1
  2  4    1   0.000000e+00 0.000000e+00
  2  5    1   0.000000e+00 0.000000e+00
  2  7    1   0.000000e+00 0.000000e+00
  2  8    1   0.000000e+00 0.000000e+00
  3  7    1   0.000000e+00 0.000000e+00
  3  8    1   0.000000e+00 0.000000e+00
  4  6    1   0.000000e+00 0.000000e+00
  5  6    1   0.000000e+00 0.000000e+00

[ angles ]
; ai aj ak funct   c0            c1
  1  3  4   1   1.200000e+02 2.928800e+02
  1  3  5   1   1.200000e+02 2.928800e+02
  4  3  5   1   1.200000e+02 3.347200e+02
  1  6  7   1   1.200000e+02 2.928800e+02
  1  6  8   1   1.200000e+02 2.928800e+02
  7  6  8   1   1.200000e+02 3.347200e+02
  2  1  3   1   1.215000e+02 5.020800e+02
  2  1  6   1   1.215000e+02 5.020800e+02
  3  1  6   1   1.170000e+02 5.020800e+02

[ dihedrals ]

; ai aj ak al funct   c0            c1            c2
  2  1  3  4   1   1.800000e+02 3.347200e+01 2.000000e+00
  6  1  3  4   1   1.800000e+02 3.347200e+01 2.000000e+00
  2  1  3  5   1   1.800000e+02 3.347200e+01 2.000000e+00
  6  1  3  5   1   1.800000e+02 3.347200e+01 2.000000e+00
  2  1  6  7   1   1.800000e+02 3.347200e+01 2.000000e+00
  3  1  6  7   1   1.800000e+02 3.347200e+01 2.000000e+00
  2  1  6  8   1   1.800000e+02 3.347200e+01 2.000000e+00
  3  1  6  8   1   1.800000e+02 3.347200e+01 2.000000e+00

[ dihedrals ]
; ai aj ak al funct   c0            c1
  3  4  5  1   2   0.000000e+00 1.673600e+02
  6  7  8  1   2   0.000000e+00 1.673600e+02
  1  3  6  2   2   0.000000e+00 1.673600e+02

[ position_restraints ]
; you wouldn’t normally use this for a molecule like Urea,
; but we include it here for didactic purposes
; ai funct fc
  1   1    1000 1000 1000 ; Restrain to a point
  2   1    1000    0 1000 ; Restrain to a line (Y-axis)
  3   1    1000    0    0 ; Restrain to a plane (Y-Z-plane)

; Include SPC water topology
#include "spc.itp"

[ system ]
Urea in Water

[ molecules ]
Urea1
SOL1000

Обратим внимание, что в данном примере не видно разделов [ defaults ] и [atomtypes ]. В данном случае они содержатся в файле forcefield.itp, который включается в наш файл топологии с помощью оператора #include. Соответствующий файл forcefield.itp лежит в директории gmx.ff/ внутри стандартной директории Gromacs /usr/share/Gromacs/top/:

[ defaults ]
; nbfunccomb-rulegen-pairsfudgeLJfudgeQQ
11no1.01.0

#include "ffnonbonded.itp"
#include "ffbonded.itp"

Видим, что он содержит раздел [defaults], в котором находятся необходимые параметры для расчёта взаимодействий, и подключает ещё два файла ffnonbonded.itp и ffbonded.itp .

Первому параметру (nbfunc) сопоставлено значение 1. Это означает, что для описания ван-дер-Ваальсовских взаимодействий будет использоваться леннард-джонсовский потенциал. Значение 1 для следующего параметра (comb-rule) означает, что потенциал задается параметрами С6 и С12 (а не σ и ε), и параметры взаимодействия между атомами разного типа рассчитываются как среднее геометрическое: , . Параметру gen-pairs соответствует no. Это означает, что параметры 1-4 – взаимодействия (между атомами, разделёнными тремя связями) не будут генерироваться при препроцессинге. Здесь этого и не требуется, так как все параметры взаимодействий указаны в явном виде.

В данном случае в файле ffnonbonded.itp содержится информация о взаимодействиях между несвязанными атомами, а в файле ffbonded.itp – между ковалентно-связанными. Эти два файла содержат параметры для большого набора атомов, относящихся к разным молекулам. Мы приводим только те строки этих файлов, которые относятся к атомам нашей молекулы.

Файл ffnonbonded.itp:

[ atomtypes ]
;name at.num mass     charge ptype     c6          c12
  O     8    15.99940  0.000  A     0.22617E-02 0.74158E-06
  NT    7    14.00670  0.000  A     0.24362E-02 0.16924E-05
  C     6    12.01100  0.000  A     0.23402E-02 0.33740E-05
  H     1     1.00800  0.000  A     0.00000E+00 0.00000E+00

[ nonbond_params ]
; i j func    c6           c12
  O  O 1   0.22617E-02 0.74158E-06
  O NT 1   0.23473E-02 0.25313E-05
 NT NT 1   0.24362E-02 0.16924E-05
  C  C 1   0.23402E-02 0.33740E-05

[ pairtypes ]
; i j func    cs6         cs12
  O  O 1   0.22617E-02 0.74158E-06
  O NT 1   0.23473E-02 0.11203E-05
  O  C 1   0.23006E-02 0.15818E-05
  O  H 1   0.00000E+00 0.00000E+00
 NT NT 1   0.24362E-02 0.16924E-05
 NT  C 1   0.23877E-02 0.23896E-05
 NT  H 1   0.00000E+00 0.00000E+00
  C  C 1   0.23402E-02 0.33740E-05
  C  H 1   0.00000E+00 0.00000E+00
  H  H 1   0.00000E+00 0.00000E+00

Файл ffbonded.itp:

[ bondtypes ]
; i j func b0       kb
  C  O 1  0.12300 502080.
  C NT 1  0.13300 376560.
  H NT 1  0.10000 374468.

[ angletypes ]
; i j  k func   th0     cth
  H NT C  1   120.000 292.880
  H NT H  1   120.000 334.720
  NT C O  1   124.000 502.080

[ dihedraltypes ]
; i l func q0     cq
  H NT 2  0.000 167.360

В разделе [ atomtypes ] перечислены типы атомов, которые встречаются в данной модели, их масса, заряд и параметры С6, С12. Разделы [ nonbond_params ], [pairtypes] и все разделы из файла ffbonded.itp. не являются обязательными. Они нужны, если в последующих разделах файла топологии ([ bonds ], [ pairs ], [ angles ], [dihedrals ]) отсутствуют заданные параметры. В этом случае они могут быть рассчитаны из общих установок данного поля сил.

Вернёмся к содержимому файла топологии молекулы мочевины в воде.

  • В разделе [ moleculetype ] пишется имя молекулы, описание которой приводится далее.
  • В разделе [ atoms ] прописано, что, например, атом с номером 1 имеет тип (type) C, относится к первому (resnr=1) остатку с именем UREA и именуется (atom) C1. Именно это имя C1 должно фигурировать в файле молекулярной структуры для мочевины (.gro). Также в последнем столбце этого раздела указан частичный заряд (charge) на данном атоме (в целом молекула, естественно, электронейтральна).
  • В разделе [ bonds ] для пары ковалентно-связанных атомов указано, что связь удерживается гармоническим потенциалом (funct 1), прописаны длины связей (в нанометрах) и жёсткость (кДж/моль/нм2).
  • В разделе [ pairs ] прописаны параметры 1-4 взаимодействия, в данном случае С6 и С12.
  • В разделе [ angles ] для каждой тройки атомов указано, что деформация угла также описывается гармоническим потенциалом (funct 1), указаны параметры связей: С0 – угол (в данном случае 120 градусов.), С1 – жесткость (кДж/моль/рад2).
  • В разделе [ dihedrals ] перечислены четвёрки атомов, составляющих двугранные углы стандартного типа 1 (funct 1), задаваемого периодическим потенциалом , где – текущий угол, – равновесный угол (град.), – жёсткость (кДж/моль) и – мультиплетность.
  • Вторая секция двугранных углов [ dihedrals ] отличается от первой типом функции (funct 2) и призвана к тому, чтобы удерживать некоторые атомы в одной плоскости (0 градусов), например, ароматические кольца. Данный тип функции соответствует гармоническому потенциалу.
  • В разделе [ position_restraints ] можно указать, если это требуется, следует ли удерживать какие-то атомы в данной точке пространства – как например первый атом (первая строчка в этом разделе): после его номера следует тип функции – гармоническая (funct 1), и три жёсткости гармонического потенциала kx, ky, kz. Если молекулу требуется привязать к линии вдоль оси y, ограничив движение по двум другим осям, то подойдёт пример из второй строчки: для атома с номером 2 применяется, опять же, гармонический потенциал (funct 1) и его жёсткости kx=1000, ky=0, kz=1000. Третья строка соответствует удерживанию атома номер 3 в плоскости x.
  • Далее, #include подключает файл топологии воды. В данном примере для воды используется трёхцентровая модель SPC, параметры которой читаются из подключаемого файла spc.itp, который программа найдёт в стандартной директории usr/share/Gromacs/top/
  • Затем в разделе [ system ] написано имя системы Urea in Water
  • В разделе [ molecules ] указано, что модель состоит из одной молекулы мочевины и 1000 молекул воды.

Файл молекулярно-динамических параметров (molecular dynamics parameter file) (.mdp)

Он содержит общую информацию о моделировании, в частности, шаг интегрирования dt, число шагов моделирования nsteps, температуру T, давление P, параметры термостата и баростата и т.п. Если не указать значения каких-то параметров, то будут использоваться параметры по умолчанию. Ниже показан пример такого файла для моделирования молекулы белка в воде:

; parameter file
title = Yo
cpp = /lib/cpp
include = -I../top
define =
integrator = md
dt = 0.002
nsteps = 500000
nstxout = 5000
nstvout = 5000
nstlog = 5000
nstenergy = 250
nstxtcout = 250
xtc_grps = Protein
energygrps = Protein SOL
nstlist = 10
ns_type = grid
rlist = 0.8
coulombtype = cut-off
rcoulomb = 1.4
rvdw = 0.8
tcoupl = v-rescale
tc-grps = Protein SOL
tau_t = 0.1 0.1
ref_t = 300 300
Pcoupl = Berendsen
tau_p = 1.0
compressibility = 4.5e-5
ref_p = 1.0
gen_vel = yes
gen_temp = 300
gen_seed = 173529
constraints = all-bonds

Здесь указано, что шаг интегрирования равен 0.002 пикосекунды (dt=0.002), будет проведено 500000 шагов (nsteps=500000), температура моделирования 300К (ref_t=300), давление 1 Бар (tau_p=1), будет использован термостат v-rescale, в котором масштабируются скорости (tcoupl=v-rescale) и баростат Берендсена (Pcoupl=Berendsen), в котором масштабируется размер периодического бокса. Строка (gen_vel=yes) даёт команду сгенерировать скорости атомов в момент запуска моделирования. Направления и величины скоростей задаются случайным образом, но с требованием, чтобы распределение по скоростям соответствовало Максвелловскому, а средняя величина соответствовала заданной температуре 300 К (gen_temp = 300). По умолчанию файл параметров имеет имя grompp.mdp.

Индекс-файл (Index file) (.ndx)

Иногда бывает необходимо выделить группы атомов, например, кислороды воды или кислороды спирта. В Gromacs предусмотрена такая возможность. Индекс-файл служит для хранения списков номеров интересующих нас атомов. Создать индекс-файл можно при помощи программы make_ndx.

Ниже показан пример содержимого индекс-файла, который мы получим при моделировании жидкого метанола в работе №5:

 System ]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
...
3061 3062 3063 3064 3065 3066 3067 3068 3069 3070 3071 3072
[ OA ]
5 11 17 23 29 35 41 47 53 59 65 71 77 83 89
...
3065 3071
[ C ]
1 7 13 19 25 31 37 43 49 55 61 67 73 79 85
...
3061 3067
[ H ]
2 3 4 8 9 10 14 15 16 20 21 22 26 27 28
...
3062 3063 3064 3068 3069 3070
[ HO ]
6 12 18 24 30 36 42 48 54 60 66 72 78 84 90
...
3066 3072

Система содержит 512 молекул метанола, каждая из которых состоит из шести атомов. Все атомы пронумерованы от 1 до 3072 в своем файле conf.gro. Раздел [ System ] представляет группу, соответствующую всем атомам системы, раздел [ OA ]содержит номера кислородов, раздел [ C ] — углеродов, а разделы [ H ] и [ HO ] — номера протонов, относящихся к метильной группе и к кислороду, соответственно.

Стартовый файл (Run input file) (.tpr)

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

Файл траектории (Trajectory file) (.trr)

Файл с указанным расширением содержит молекулярно-динамическую траекторию, где записаны координаты и скорости атомов в моменты времени, заданные в файле МД-параметров grompp.mdp (переменные nstxout и nstvout). Файл бинарный, но может быть прочитан программой gmxdump. Результат обработки файла .trr программой gmxdump (запущенной командой gmxdump –f traj.trr | less) выглядит так:

traj.trr frame 0:
atoms= 216 step= 0 time=0.0000000e+00 lambda= 0
box (3x3):
box[ 0]={ 6.65637e+00, 0.00000e+00, 0.00000e+00}
box[ 1]={ 0.00000e+00, 6.65637e+00, 0.00000e+00}
box[ 2]={ 0.00000e+00, 0.00000e+00, 6.65637e+00}
x (216x3):
x[ 0]={ 1.03600e+00, 3.94800e+00, 4.91100e+00}
x[ 1]={ 2.16800e+00, 7.73000e-01, 5.70500e+00}
x[ 2]={ 4.50000e+00, 2.11600e+00, 1.26800e+00}
x[ 3]={ 3.59700e+00, 5.46700e+00, 5.55900e+00}
x[ 4]={ 4.61200e+00, 2.09400e+00, 5.76700e+00}
...
x[ 215]={ 4.24600e+00, 4.30000e-01, 1.12500e+00}
v (216x3):
v[ 0]={ 1.39500e+00, 1.05700e-01, -8.68800e-01}
v[ 1]={-6.68400e-01, -8.70200e-01, 1.88340e+00}
v[ 2]={ 1.48900e-01, 1.48600e-01, -5.82600e-01}
v[ 3]={ 1.29790e+00, -7.02400e-01, 6.19300e-01}
v[ 4]={ 7.21400e-01, 8.43800e-01, -1.75460e+00}
...
v[ 214]={-2.10000e-03, -3.67400e-01, 6.02800e-01}
v[ 215]={-1.29860e+00, -6.38200e-01, 2.03590e+00}
traj.trr frame 1:
natoms= 216 step= 2000 time=4.0000000e+00 lambda= 0
box (3x3):
box[ 0]={ 6.65637e+00, 0.00000e+00, 0.00000e+00}
box[ 1]={ 0.00000e+00, 6.65637e+00, 0.00000e+00}
box[ 2]={ 0.00000e+00, 0.00000e+00, 6.65637e+00}
x (216x3):
x[ 0]={ 1.36948e+00, 2.75648e+00, 5.24567e+00}
x[ 1]={ 4.66562e+00, 6.50124e+00, 3.67500e-01}
x[ 2]={ 5.20250e+00, 2.24663e+00, 9.63009e-01}
...

Показаны начала двух конфигураций (мгновенных снимков, «фрэймов»): frame 0 и frame 1. Каждая из них начинается с информационной строки, где указывается число атомов и номер шага моделирования (номер данной конфигурации). После этого идут три строки, задающие модельный бокс, и затем – строки с координатами атомов (начиная с 0-ого по 215-й, в данном случае). После координат аналогичным способом записаны значения скоростей.

Файл траектории (Trajectory file) (.xtс)

Файл с указанным расширением хранит траекторию в более компактном виде по сравнению с файлом .trr . Он также бинарный, содержит координаты, но не содержит скоростей. Можно задать точность, с которой требуется сохранять данные. По умолчанию хранится только три знака после запятой, этого обычно бывает достаточно для анализа структуры и динамики. В файле параметров grompp.mdp за xtc-формат отвечает ряд переменных: nstxtcout, xtc-precision, xtc-grps). Заглянуть внутрь .xtc-файла поможет также программа gmxdump.

Файл энергии (Energy file) (.edr)

Бинарный файл, для работы с которым предназначена программа g_energy. В интерактивном режиме появляется список имеющихся характеристик, рассчитанных в процессе МД-моделирования и хранящихся в .edr-файле:

1 LJ-(SR)     2 Disper.-corr. 3 Coulomb-(SR) 4 Potential 
5 Kinetic-En. 6 Total-Energy  7 Temperature  8 Pres.-DC 
9 Pressure   10 Vir-XX       11 Vir-XY      12 Vir-XZ 
13 Vir-YX    14 Vir-YY       15 Vir-YZ      16 Vir-ZX 
17 Vir-ZY    18 Vir-ZZ       19 Pres-XX     20 Pres-XY 
21 Pres-XZ   22 Pres-YX      23 Pres-YY     24 Pres-YZ 
25 Pres-ZX   26 Pres-ZY      27 Pres-ZZ     28 #Surf*SurfTen
29 Mu-X      30 Mu-Y         31 Mu-Z        32 T-System 

Для выбора характеристики достаточно набрать её номер или же само название характеристики из списка.

Файл сопровождающей информации (Log-file) (.log)

При работе программы mdrun формируется также файл сопровождающей информации (log-файл), имеющий расширение (.log), содержание которого помогает контролировать расчёт. Это текстовый файл, частота записи задаётся в файле параметров grompp.mdp (nstlog – число МД-шагов, через которое записывается информация в log-файл).

Основные программы

Отметим кратко основные программы Gromacs'a, с которыми будем работать, их назначение и параметры. Полную информацию по каждой программе можно найти в описании или во встроенном справочнике Gromacs'a с помощью ключа -h. Например, для получения описания программы grompp надо набрать команду: grompp -h. Более подробно с особенностями работы программ мы будем знакомиться постепенно, по мере выполнения работ данного практикума.

grompp – (препроцессинг) – создает стартововый файл .tpr для начала моделирования.

Входные файлы (стандартные имена):

grompp.mdp – файл параметров моделирования,
conf.gro – файл молекулярной структуры,
topol.top – файл топологии,

Выходной файл:

topol.tpr – стартовый файл для последующего запуска МД-моделирования.

Запуск с использованием стандартных имен:

grompp

в общем случае, для нестандартных имен, эта программа использует ключи -f , -c, -p, -o, например:

grompp -f name1.mdp -c name2.gro -p name3.top -o name4.tpr

mdrun – запуск моделирования. Создает молекулярно-динамическую траекторию.

Входные файлы (стандартные имена):

topol.tpr – стартовый файл.

Выходные файлы:

traj.trr – файл траектории (полная, вместе со скоростями),
traj.xtc – файл траектории (сжатой),
confout.gro – файл конечной конфигурации,
ener.edr – файл энергии,
md.log – файл сопровождающей информации.

Запуск с использованием стандартных имен:

mdrun

Запуск моделирования с указанием нестандартного имени входного файла (ключ -s):

mdrun -s sys1.tpr

Все выходные файлы в данном примере будут имеет стандартные имена.

Запуск моделирования с заданием нестандартных имен выходных файлов производится с помощью ключей -o, -x, -c, -e, -g:

mdrun -o name.trr -x name.xtc -c name.gro -e name.edr -g name.log

Здесь молекулярно-динамическая траектория и сопутствующие ей файлы будут иметь общее имя (name), а различаться расширением. Это бывает удобно при дальнейшей работе с системой. Входной файл в данном примере, предполагается, имеет стандартное имя (topol.tpr).

С помощью дополнительных ключей можно управлять параметрами счёта, например, устанавливать приоритет задачи (-nice), задавать параметры распараллеливания вычислений, (–nt, -dd), включать режим вывода на экран текущей информации в процессе счёта (-v) и т.п.

make_ndx – создание индекс-файла (определение групп атомов). Работает в интерактивном режиме, в котором нужно выбрать имена атомов, объединяющихся в группу.

Входные файлы (стандартные имена):

conf.gro – файл структуры (ключ -f),

возможно также использование файла формата .pdb.

Выходные файлы:

index.ndx – индекс-файл (ключ -o).

editconf — подготовка исходной конфигурации, масштабирование, задание размеров и формы бокса.

Входные файлы (стандартные имена):

conf.gro – файл структуры (ключ -f).

возможно также использование файла траектории (.trr).

index.ndx – индекс-файл (ключ -n).

Выходные файлы:

out.gro – файл с полученной конфигурацией (ключ -o).

Возможно задать размер бокса (ключ -box), форму бокса (-bt), центрирование молекулы в боксе (-c), задать минимальное расстояние от молекулы до края бокса (-d) и др.

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

Входные файлы (стандартные имена):

conf.gro – файл структуры (ключ -f),

возможно также использование файла траектории (.trr).

Выходные файлы:

out.gro – файл с полученной конфигурацией (ключ -o)

Количество трансляций вдоль каждой из трех осей задаётся числами после ключа –nbox , например:

genconf –nbox 2 3 2

genbox – создает модель раствора, заполняет свободное пространство вокруг заданной молекулы белка (системы молекул) растворителем (водой).

Входные файлы (имена указываются по ключу):

protein.gro – файл структуры белка (ключ -cp),
water.gro – файл структуры растворителя (ключ -cs),
topol.top – файл топологии (ключ -p).

Выходные файлы:

out.gro – файл с полученной конфигурацией (ключ -o).

Пример "растворения" белка в воде:

genbox -cp 1AKI_box.gro -cs spc216.gro -o 1AKI_solv.gro -p

Пример создания раствора метанол-вода:

genbox -cp methanol_box.gro -cs h2o.gro

genion – добавляет в раствор небольшие ионы (обычно это необходимо, чтобы скомпенсировать заряд белка).

Входные файлы (стандартные имена):

topol.tpr – стартовый файл для раствора (ключ -s),
topol.top – файл топологии (ключ -p),

указывается по ключу:

ions.mdp – файл параметров для иона (ключ -f).

Выходные файлы:

out.gro – файл с полученной конфигурацией (ключ -o).

Ионы заменяют случайно выбранные молекулы воды. Тип иона задаётся ключом –nname (для отрицательного) или –pname (для положительного). Количество ионов определяется ключом –nn. Стартовый файл topol.tpr для программы genion создается специально с помощью программы grompp с учетом файла параметров для ионов, например, ions.mdp. По завершению программа исправляет файл топологии topol.top, добавляя информацию о вставленных ионах.

gmxdump – читает бинарный файл (стартовый файл, файлы траектории и энергии), представляет информацию в текстовом виде на экране компьютера. Для разных файлов используется разные ключи (-s,-f,-e). Примеры:

gmxdump -s topol.tpr– читает бинарный стартовый файл,

gmxdump -f traj.xtc– читает бинарный файл траектории,

gmxdump -е ener.edr– читает бинарный файл энергии.

xmgrace – рисует график функции из файлов с расширением .xvg. Пример:

xmgrace potential.xvg– рисует график из указанного файла.

xmgrace rdf*.xvg– рисует вместе все графики из файлов текущей директории, имена которых начинаются с букв rdf.

ngmx – изображает на экране последовательные конфигурации молекулярно-динамической траектории (анимация).

Входные файлы (стандартные имена):

traj.xtc – файл МД-траектории, можно также файл traj.xrr,
topol.tpr– стартовый файл.

Данную программу имеет смысл использовать для оперативного просмотра полученной траектории. Пример:

ngmx

Можно выбрать участок траектории, указав после ключей -b и -е время (в пикосекундах), соответствующее начальной и конечной конфигурациям. Можно смотреть не каждую конфигурацию, а через шаг (в пикосекундах), указанный после ключа -dt.

Демонстрационная работа

На сайте http://manual.gromacs.org/online/water.html находится демонстрация моделирования воды, содержащей 216 молекул с потенциалом SPC. Эта простая задача реализована также на вашем компьютере.

Получение модели

В директории 01_demo/water/ имеются готовые: файл молекулярной структуры conf.gro с координатами всех молекул, файл топологии topol.top , а также файл с МД-параметрами grompp.mdp.

Все эти файлы можно посмотреть с помощью команд Linux:

less conf.gro
less topol.top
less grompp.mdp

Файл структуры можно визуализировать с помощью программы VMD, посмотрев исходное расположение молекул в системе:

vmd conf.gro

С помощью программы grompp создадим из заданных файлов стартовый файл topol.tpr:

grompp

(Не путать имя программы grompp с именем файла grompp.mdp) Здесь не указаны имена используемых файлов, так как все они имеют стандартные имена и расположены в текущей рабочей директории. После работы этой программы в рабочей директории появляется файл topol.tpr. Он записан в бинарном формате, однако его содержимое можно просмотреть с помощью программы gmxdump:

gmxdump -s topol.tpr | less(выход – по нажатию клавиши q)

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

mdrun -v

Ключ -v включает режим расширенного вывода информации, т.е. в процессе работы программы на экран выводится текущая информация о её работе. После выполнения этой программы в рабочей директории появляется файл траектории traj.xtc. Можно посмотреть (визуализовать) полученную траекторию с помощью собственной программы Gromacs:

ngmx

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

vmd conf.gro traj.xtc
Fig1.1.jpg

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

Рассчитаем функцию радиального распределения для атомов кислорода (ФРР), (radial distribution function), в данном случае, парную корреляционную функцию g(r). Программа g_rdf по умолчанию вычисляет парную корреляционную функцию для всех атомов (кислородов и водородов в данном случае). Однако интерес представляет функция, учитывающая только кислороды. Для выделения атомов заданного типа используется индекс-файл (index.ndx). В нашем примере он уже создан и в нём выделена одна группа атомов, а именно – кислороды. Поэтому достаточно выполнить команду:

g_rdf -n

где ключ -n означает подключение индекс-файла со стандартным именем index.ndx. В результате работы программы будет создан файл rdf.xvg. Он содержит искомую функцию g(r), которую можно посмотреть с помощью программы xmgrace:

xmgrace rdf.xvg

Для удобства можно было ещё на предыдущем этапе добавить к команде g_rdf ключ –w.

g_rdf -n -w
Этот ключ означает, что после завершения расчёта автоматически запустится программа xmgrace, которая сразу нарисует график рассчитанной функции.
  1. http://www.ks.uiuc.edu/Research/vmd/