Работа 1

Знакомство с пакетом Gromacs

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

Цель

Познакомиться с пакетом 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/ или же /opt/gromacs-2016.4/share/gromacs/top/, а исполняемые файлы в директории /usr/bin/ или /opt/gromacs-2016.4/bin/.

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

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

/home/student/

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

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

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

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

Для создания стартового файла (topol.tpr) нужно провести «препроцессинг» (preprocessing) Это делает программа gmx 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, то с помощью программы gmx pdb2gmx можно сформировать для этой молекулы файл топологии (topol.top), используя готовое поле сил из библиотеки пакета Gromacs.

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

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

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

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

gmx rdf – расчёт парной корреляционной функции g(r),

gmx velacc – расчёт автокорреляционной функции скорости,

gmx msd – расчёт среднеквадратичного смещения и коэффициента самодиффузии,

gmx hbond – расчёт числа водородных связей,

gmx density – расчёт плотности,

gmx 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 "gromos53a6.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 лежит в директории gromos53a6.ff/ внутри стандартной директории Gromacs /usr/share/Gromacs/top/ или /opt/gromacs-2016.4/share/gromacs/top/:

[ defaults ]
 nbfunc        comb-rule        gen-pairs        fudgeLJ        fudgeQQ
 1             1                no               1.0            1.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), прописаны длины связей (в нанометрах) и жёсткость (кДж/моль/нм²).
• В разделе [ pairs ] прописаны параметры 1-4 взаимодействия, в данном случае С6 и С12.
• В разделе [ angles ] для каждой тройки атомов указано, что деформация угла также описывается гармоническим потенциалом (funct 1), указаны параметры связей: С0 – угол (в данном случае 120 градусов.), С1 – жесткость (кДж/моль/рад²).
• В разделе [ dihedrals ] перечислены четвёрки атомов, составляющих двугранные углы стандартного типа 1 (funct 1), задаваемого периодическим потенциалом , где – текущий угол, – равновесный угол (град.), – жёсткость (кДж/моль) и – мультиплетность.
• Вторая секция двугранных углов [ dihedrals ] отличается от первой типом функции (funct 2) и призвана к тому, чтобы удерживать некоторые атомы в одной плоскости (0 градусов), например, ароматические кольца. Данный тип функции соответствует гармоническому потенциалу.
• В разделе [ position_restraints ] можно указать, если это требуется, следует ли удерживать какие-то атомы в данной точке пространства – как например первый атом (первая строчка в этом разделе): после его номера следует тип функции – гармоническая (funct 1), и три жёсткости гармонического потенциала k_x, k_y, k_z. Если молекулу требуется привязать к линии вдоль оси y, ограничив движение по двум другим осям, то подойдёт пример из второй строчки: для атома с номером 2 применяется, опять же, гармонический потенциал (funct 1) и его жёсткости k_x=1000, k_y=0, k_z=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
integrator               = md
dt                       = 0.002
nsteps                   = 500000

nstlog                   = 5000
nstenergy                = 5000
nstxout-compressed       = 5000

continuation             = yes
constraints              = all-bonds
constraint-algorithm     = lincs

cutoff-scheme            = Verlet

coulombtype              = PME
rcoulomb                 = 1.0

vdwtype                  = Cut-off
rvdw                     = 1.0
DispCorr                 = EnerPres

tcoupl                   = V-rescale
tc-grps                  = Protein  SOL
tau-t                    = 0.1      0.1
ref-t                    = 300      300

pcoupl                   = Berendsen
tau-p                    = 2.0
compressibility          = 4.5e-5
ref-p                    = 1.0

gen_vel = yes
gen_temp = 300
gen_seed = 173529

Здесь указано, что шаг интегрирования равен 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 надо набрать команду: gmx grompp -h. Более подробно с особенностями работы программ мы будем знакомиться постепенно, по мере выполнения работ данного практикума.

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

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

grompp.mdp – файл параметров моделирования,

conf.gro – файл молекулярной структуры,

topol.top – файл топологии,

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

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

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

gmx grompp

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

gmx 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 – файл сопровождающей информации.

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

gmx mdrun

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

gmx mdrun -s sys1.tpr

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

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

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

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

gmx mdrun -deffnm name

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

genconf –nbox 2 3 2

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

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

protein.gro – файл структуры белка (ключ -cp),

water.gro – файл структуры растворителя (ключ -cs),

topol.top – файл топологии (ключ -p).

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

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

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

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

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

gmx genbox -cp methanol_box.gro -cs h2o.gro

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

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

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

topol.top – файл топологии (ключ -p),

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

ngmx

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

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

Проведём моделирование воды, содержащей 216 молекул с потенциалом SPC.

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

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

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

less conf.gro

less topol.top

less grompp.mdp

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

vmd conf.gro

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

gmx grompp

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

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

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

gmx mdrun -nt 2 -v

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

ll -t

где ключ -t означает сортировку по времени. Среди новых файлов есть файл траектории traj_comp.xtc . Прежде чем его визуализировать, необходимо восстановить молекулы, разрезанные периодическими границами. Это делает команда gmx trjconv с флагом -pbc mol:

gmx trjconv -f traj_comp.xtc -pbc mol

Данная программа выдаёт на выходе файл trajout.xtc , его можно посмотреть (визуализовать) с помощью программы VMD. Для этого следует набрать команду:

vmd conf.gro trajout.xtc

В программе VMD можно включить воспроизведение траектории, нажав на "play", а также настроить внешний вид молекул, например, на шарики и палочки Drawing method -> CPK.

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

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

gmx rdf -s -n

где ключ -s означает, что также будет использован файл топологии, по умолчанию называющийся topol.tpr. -n означает подключение индекс-файла со стандартным именем index.ndx. После запуска требуется дважды выбрать группу 0, затем нажать комбинацию клавиш Ctrl+D. В результате работы программы будет создан файл rdf.xvg , содержимое которого можно посмотреть командой less:

less rdf.xvg

Файл содержит искомую функцию g(r), которую можно посмотреть с помощью программы xmgrace:

xmgrace rdf.xvg

.