Работа 5

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

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

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

Цель

Провести моделирование жидкого метанола. Исследовать его и свойства.

Задания

1. Получить МД модель жидкого метанола.

(создать начальный файл, провести релаксацию, провести основное моделирование).

2. Рассчитать функции радиального распределения (полноатомные и парциальные):

  • g(r) для всех атомов, учитывая внутримолекулярные расстояния.
  • g(r) для всех атомов без учета внутримолекулярных расстояний.
  • gO-O(r) только для атомов кислорода.
  • gС-С(r) только для атомов углерода.
  • gO-HO(r) для кислородов и гидроксильных водородов.
  • gO-HС(r) для кислородов и метильных водородов.

3. Найти коэффициент самодиффузии.

4. Найти среднее число водородных связей на одну молекулу.

Вопросы

  1. Дать интерпретацию внутримолекулярным пикам g(r), исходя из структурной формулы молекулы метанола.
  2. Объяснить пики межмолекулярной g(r) жидкого метанола.
  3. Показать проявление водородных связей между молекулами метанола на парциальных g(r).

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

Молекула метилового спирта является одной из простейших, но важных органических молекул. Жидкий метанол относится к основным продуктам промышленного органического синтеза. Молекулы спиртов, благодаря наличию группы ОН, образуют друг с другом водородные связи, которые существенным образом влияют на физические свойства таких жидкостей. Например, алканы, у которых нет водородных связей, являются газообразными при нормальных условиях (от метана до бутана), тогда как все соответствующие спирты – жидкости, причём имеющие достаточно высокую температуру кипения. Группа ОН спирта может образовать не более трёх водородных связей (одну через водород и две – через неспаренные электроны кислорода), Рис. 5.21. Однако углеводородные части молекул мешают им сблизиться нужным образом, поэтому обычно реализуется одна – две (или даже ни одной) водородные связи. В отличие от воды, в спиртах нет «трехмерной» (в топологическом смысле) сетки связей, в них имеются линейные или ветвистые агрегаты водородносвязанных молекул.

Напомним, что водородная связь не задаётся каким-либо специальным потенциалом взаимодействия, а образуется сама собой как следствие кулоновского взаимодействия атомов кислорода и водорода. Атомам кислорода выгодно располагаться на некотором расстоянии друг от друга, не ближе, не дальше, и находиться примерно на одной прямой с водородом. По этим признакам водородные связи и находят в модели: если расстояние между кислородами ≤ 0,35 нм, а угол между направлением OH и линией кислород-кислород ≤ 30º, то считается, что молекулы связаны водородной связью.

Для жидких метанола и этанола практически все водороды гидроксильной группы участвуют в водородных связях. В спиртах с большим углеводородным радикалом остаётся заметная доля OH групп, которые не могут образовать водородную связь. Это хорошо видно из экспериментов по колебательной спектроскопии, так как частота ОН-колебания сильно зависит от того, участвует данный гидроксильный водород в водородной связи или нет. В таких спиртах виден чёткий пик, соответствующий несвязанным ОН-осцилляторам.

Fig5.1.png
Рис. 5.1. Водородные связи между молекулами метилового спирта. Теоретически возможны три связи на каждую молекулу, но из-за стерических препятствий не все из них могут возникнуть.

В молекулярной динамике компактная атомная группа (например, метильная группа) часто рассматривается как единая частица – «объединенный атом», Рис. 5.22. Это разумное физическое приближение, которое позволяет уменьшить время вычислений по сравнению с «полноатомным» описанием, где учитывается каждый атом. Оно хорошо работает, если целью моделирования является изучение структурных и динамических свойств самой жидкости. Если позволяют вычислительные мощности, то используют полноатомное описание, которое содержит больше деталей. В данной работе мы используем полноатомное представление молекулы метанола.

Fig5.2.jpg
Рис. 5.2. Молекула метилового спирта. Cлева: «полноатомное» представление (all-atoms). Справа: метильная группа представлена «объединенным атомом» (united atom).

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

Моделирование жидкого метанола

Работа выполняется в директории 05_methanol/. Там находятся файлы topol.top и поддиректории, необходимые для работы.

Рассмотрим файл topol.top:

#include "ffoplsaa.itp"
#include "oplsaa.ff/methanol.itp"

[ system ]
Pure Methanol – Yummie!

[ molecules ]
MET	1

Include-файл ffoplsaa.itp находится в стандартной директории GROMACS (/usr/share/gromacs/top/). Он подключает поле сил OPLS (optimized potential for liquid simulations). Все файлы этого поля лежат в папке oplsaa.ff/. Там же находится файл топологии для метанола.

Так будет выглядеть полностью собранный файл топологии, после подстановки необходимых разделов из поля сил ffG43a1.itp и methanol.itp:

;information from ffoplsaa.itp, oplsaa.ff/forcefield.itp:
#define _FF_OPLS
#define _FF_OPLSAA
[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               3               yes             0.5     0.5

;information from /usr/share/gromacs/top/oplsaa/ffnonbonded.itp:
[ atomtypes ]
; name   bond_type      mass     charge   ptype    sigma      ep-silon
 opls_157   CT    6     12.01100     0.145       A    3.50000e-01  2.76144e-01
 opls_156   HC    1      1.00800     0.040       A    2.50000e-01  1.25520e-01
 opls_154   OH    8     15.99940    -0.683       A    3.12000e-01  7.11280e-01
 opls_155   HO    1      1.00800     0.418       A    0.00000e+00  0.00000e+00

;information from /usr/share/gromacs/top/oplsaa/ffbonded.itp:
[ bondtypes ]
; i     j      func     b0       kb
  CT    HC      1    0.10900   284512.0   ; CHARMM 22 parameter file
  CT    OH      1    0.14100   267776.0   ;
  HO    OH      1    0.09450   462750.4   ; SUG(OL) wlj mod 0.96-> 0.945
[ angletypes ]
;  i     j      k    func       th0       cth
  HC     CT     OH      1   109.500    292.880   ;
  CT     OH     HO      1   108.500    460.240   ;
[ dihedraltypes ]
;  i     j      k      l    func      coefficients
  HC     CT     OH     HO      3      0.94140   2.82420   0.00000  -3.76560   0.00000   0.00000 ; alcohols AA

;information from oplsaa.ff/methanol.itp
;
; Methanol, Jorgensen et al. JACS 118 pp. 11225 (1996)
; 
[ moleculetype ]
; name  nrexcl
MET     3

[ atoms ]
;    nr      type        resnr   residue   atom  cgnr  charge
     1       opls_157       1     MET        C     1   0.145
     2       opls_156       1     MET        H     1   0.04
     3       opls_156       1     MET        H     1   0.04
     4       opls_156       1     MET        H     1   0.04
     5       opls_154       1     MET       OA     1   -0.683
     6       opls_155       1     MET       HO     1   0.418

[ bonds ]
; ai    aj  funct           c0           c1
1       2       1
1       3       1
1       4       1
1       5       1
5       6       1

[ pairs ]
; i     j       funct
2       6
3       6
4       6

[ angles ]
;  ai    aj    ak     funct           c0        c1
; H3
 2      1       5       1
 3      1       5       1
 4      1       5       1
 4      1       3       1
 4      1       2       1
 3      1       2       1
 1      5       6       1

;
[ dihedrals ]
; ai    aj     ak       al    funct      c0        c1        c3...
2       1       5       6       3
3       1       5       6       3
4       1       5       6       3

[ system ]
Pure Methanol – Yummie!

[ molecules ]
Methanol	1

Итак, наша молекула метанола описывается атомами четырех типов:

  • C – углерод метильной группы (имя атома в поле сил: opls_157, тип атома CT),
  • H – водород метильной группы (имя атома в поле сил: opls_156, тип атома HC),
  • OA – кислород (имя атома opls_154, тип атома OH)
  • HO – водород гидроксильной группы (имя атома opls_155, тип атома HO).

Это указано в разделах [ atomtypes ] и [ atoms ]. Там также указаны параметры потенциалов, в частности заданы значения sigma, характеризующие ван-дер ваальсовские диаметры атомов.

В разделах [ bonds ] указана связность атомов (первого со вторым, первого с третьим, и т.д.), но не указаны параметры, характеризующие ковалентные взаимодействия (длина связи и константа жёсткости). Однако эти параметры программа подставляет из раздела [ bondtypes ], где заданы длины связей:

В разделе [ angles ] указаны тройки атомов, составляющих углы, в то время как параметры, характеризующие деформационное взаимодействие, прописаны в разделе [angletypes ]: равновесный угол H-C-O равен 109.5º, C-O-H – 108.5º.

В разделе [dihedrals] указаны четвёрки атомов, составляющих двугранные углы, которые описываются периодической функцией Рикаэрта-Бэллмана (func 3), параметры прописаны в разделе [dihedraltypes ].

Приготовление исходной конфигурации

Перейдём в папку preparing. Координаты атомов молекулы находятся в файле methanol.gro. Его содержимое:

Methanol
    6
    1MeO      C    1   4.119   1.207   3.971
    1MeO      H    2   4.130   1.315   3.963
    1MeO      H    3   4.208   1.151   3.942
    1MeO      H    4   4.044   1.169   3.902
    1MeO     OA    5   4.065   1.176   4.098
    1MeO     HO    6   4.129   1.123   4.144

Обратим внимание, что координаты пишутся в нанометрах в некоторой заданной системе координат. Свяжем с нашей молекулой модельный бокс, который потом размножим нужное число раз. Размер этого бокса разумно выбрать таким, чтобы он давал желаемую плотность, т.е. чтобы его объём был примерно равен объёму, который приходится на одну молекулу в жидкости. Плотность метанола равна 791 кг/м3, откуда можно найти Vна молекулу (Получите эту оценку, используя известные значения атомных весов атомов). После того, как объём на молекулу будет найден, найдем длину ребра бокса, которая, очевидно, есть корень кубический из объёма. Проведя вычисления, получим, что искомое ребро равно 0.394 нм. Чтобы поместить нашу молекулу в такой бокс, можно воспользоваться командой editconf с ключом –box:

editconf -f methanol.gro -box 0.394

За ключом, в общем случае, следует указать размеры бокса по трём осям. Однако в данном случае можно написать одно число, программа воспримет его для всех трех осей. Полученный кубический бокс с молекулой метанола будет записан в файл out.gro.

Теперь получим бокс, содержащий 512 молекул метана. Для этого размножим наш бокс с одной молекулой метанола, транслируя его по восемь раз вдоль каждой оси:

genconf -f out.gro -nbox 8 8 8 -o methanol_512.gro

Получившийся большой бокс с ребром 3,152 нм. (0.394*8) запишется в файл methanol_512.gro. Посмотрите содержимое этого файла. Молекулы в нём расположены регулярно (в узлах простой кубической решетки), что хорошо видно, если визуализировать полученную конфигурацию с помощью VMD.

Релаксация

Проведем первый этап релаксации в директории equi1/. В ней находится файл grompp.mdp с параметрами первого этапа релаксации.

grompp.mdp:

integrator               = md		; RUN CONTROL PARAMETERS
dt                       = 0.002	; timestep in ps
nsteps                   = 10000
emtol                    = 100		; ENERGY MINIMIZATION OPTIONS ; Force tolerance and initial step-size

; OUTPUT CONTROL OPTIONS
nstxout                  = 0		; Output frequency for coords (x), velocities (v) and forces (f)
nstvout                  = 0
nstlog                   = 100		; Output frequency for ener-gies to log file and energy file
nstenergy                = 100
nstxtcout                = 100		; Output frequency and pre-cision for xtc file
xtc-precision            = 1000

; NEIGHBORSEARCHING PARAMETERS
nstlist                  = 5		; nblist update frequency
ns_type                  = grid		; ns algorithm (simple or grid)
rlist                    = 1.2		; nblist cut-off    

; OPTIONS FOR ELECTROSTATICS AND VDW
coulombtype              = pme		; Method for doing electro-statics
rcoulomb                 = 1.2
vdw-type                 = Cut-off	; Method for doing Van der Waals
rvdw                     = 1.2
DispCorr                 = EnerPres	; Apply long range disper-sion corrections for Energy and Pressure

; OPTIONS FOR WEAK COUPLING ALGORITHMS
Tcoupl                   = v-rescale	; Temperature coupling
tc-grps                  = System	; Groups to couple separate-ly
tau_t                    = 0.2		; Time constant (ps) should be at least 20 times larger than nsttcouple*dt (0.01)
ref_t                    = 300		; reference temperature (K)
Pcoupl                   = berendsen	; Pressure coupling 
Pcoupltype               = isotropic
tau_p                    = 1.0		; Time constant (ps), com-pressibility (1/bar) and reference P (bar)
compressibility          = 4.5e-5
ref_p                    = 1.0

annealing                = no		; SIMULATED ANNEALING  ; Type of annealing for each temperature group (no/single/periodic)

; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel                  = yes
gen_temp                 = 300
gen_seed                 = 2011

; OPTIONS FOR BONDS
constraints              = all-bonds
constraint-algorithm     = Lincs	; Type of constraint algo-rithm
continuation             = no		; Do not constrain the start

В нём задан термостат v-rescale и баростат Берендсена, которые обычно используются для релаксации, в случаях, когда исходное состояние заметно отличаются от равновесного. Параметры gen_vel=yes, gen_temp=300 означают, что в начальный момент молекулам будут заданы случайные (по направлению) скорости, имеющие максвелловское распределение по величине, и их значения соответствуют температуре 300 К. Указано, что рассчитанную траекторию необходимо записывать в формате .xtc с периодичностью каждые 100 шагов (nstxtcout=100). Файл .trr останется пустым, так как для него параметры записи координат и скоростей заданы нулевыми: nstxout=0, nstvout=0.

Чтобы далее работать со стандартными именами, делаем нужные ссылки. Исходный файл молекулярной структуры связываем с именем conf.gro:

ln -s ../preparing/methanol_512.gro conf.gro

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

ln -s ../topol.top .

Создаем стартовый файл topol.tpr при помощи команды grompp:

grompp

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

Запускаем МД-моделирование для первого этапа релаксации:

mdrun –nt 8 -v -nice 11

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

Чтобы убедиться, что релаксация прошла успешно, построим (при помощи g_energy) зависимость потенциальной энергии и плотности от времени. Видно, что энергия стартовой конфигурации была далека от оптимальной, Рис. 5.3 (слева). На Рис. 5.3 (справа) показано изменение плотности системы. Хотя плотность была задана равной экспериментальной (см. выше), исходная структура модели не соответствовала реальной жидкости. Видно, что система вначале расширилась (плотность упала), а потом начала постепенно выходить на равновесное значение. В нашем моделировании (за 20 пс) мы не достигли асимптотического значения для плотности, но уже близко подошли к нему. Можно было бы несколько увеличить время релаксации, однако в данном случае это не так важно, так как предстоит ещё один шаг релаксации. Целью первого этапа было перейти от искусственной исходной «кристаллической» конфигурации к разупорядоченной, более разумной для жидкости.

Fig5.3.png
Рис. 5.3. Изменение от времени потенциальной энергии (слева) и плотности (справа) на первом этапе релаксации.

Второй этап релаксации проведём в директории equi2/. Напомним, что параметры в файле grompp.mdp (кроме числа шагов моделирования) должны соответствовать параметрам, которые будут использоваться при основном моделировании. Сравним файлы параметров из разных директорий (текущей и предыдущей) и при необходимости внесем нужные изменения:

vimdiff grompp.mdp ../grompp.mdp.

После этого делаем ссылки на нужные файлы, запускаем процессинг и моделирование:

ln -s ../equi1/confout.gro conf.gro
ln -s ../topol.top .
grompp
mdrun –nt 8 -v -nice 11

Основное моделирование

Основное моделирование проводим в папке 05_methanol/. В файле параметров grompp.mdp укажем число шагов, соответствующее эволюции модели в течение 100 пс. В качестве начальной конфигурации используем файл confout.gro, полученный в директории equi2/:

ln -s ../equi2/confout.gro conf.gro

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

grompp
mdrun –nt 8 -v -nice 11

После завершения моделирования в текущей директории появятся файлы, содержащие молекулярно-динамическую траекторию нашей модели.

Анализ жидкого метанола

Функции радиального распределения

Прежде всего, рассчитаем g(r) между всеми атомами системы, включая внутримолекулярные расстояния. Для этого достаточно выполнить команду:

g_rdf -s conf.gro -o rdf-all-intramol.xvg -w

Программа по умолчанию использует файл траектории traj.xtc. Для каждой конфигурации модели, записанной в этом файле, рассчитываются расстояния между всеми атомами и строится их гистограмма. После этого производится усреднение по всем конфигурациям и нормировка на общую плотность системы. Здесь используется файл conf.gro, подаваемый после ключа -s. Результатом этой работы, в данном случае, будет файл rdf-all-intramol.xvg, в котором записана парная корреляционная функция для всех атомов, Рис. 5.4. Напомним, расширение .xvg выходного файла используется по умолчанию встроенным графическим редактором xmgrace.

Fig5.4.png
Рис. 5.4. Полноатомная парная корреляционная функция для жидкого метанола с учетом внутримолекулярных расстояний. На вставке показана область первых пиков.

Зная структурную формулу молекулы, можно понять, каким атомам соответствуют наблюдаемые пики ФРР. Однако сначала имеет смысл выяснить, какие из пиков являются внутримолекулярными, а какие межмолекулярными. Рассчитать функцию g(r) только для межмолекулярных расстояний можно также программой g_rdf , однако для этого вместо файла conf.gro нужно использовать файл topol.tpr. Если он имеет стандартное имя, то после ключа -s можно ничего не писать (однако сам ключ -s необходим, он указывает, что нужно взять файл topol.tpr.). Рассчитанную таким образом функцию запишем в файл с именем rdf-all.xvg:

g_rdf -s -o rdf-all.xvg -w

Сравним обе наши функции на одном графике, Рис. 5.5:

xmgrace rdf-all* -legend load

Здесь запись rdf-all* означает, что будут использованы все файлы, начинающиеся на rdf-all. Разумеется, можно было бы написать полностью имена нужных файлов через пробел. Ключ -legend load сообщает программе, что мы хотим, чтобы на графике были написаны названия кривых. (Подчеркнем, что обе функции мы называем "полноатомными", поскольку при расчёте расстояний используются все атомы модели. Ниже мы будем рассматривать расстояния между заданными атомами, например, только между кислородами или только между углеродами. Такие "частичные" функции называются "парциальными").

Fig5.5.png
Рис. 5.5. Полноатомные парные корреляционные функции для жидкого метанола: включая внутримолекулярные расстояния, показанная на Рис. 5.4 (кривая all-total), и для межмолекулярных расстояний (кривая all).

Понятно, что внутримолекулярные пики должны быть более острыми, чем межмолекулярные. Действительно, самые чёткие пики соответствуют ковалентно-связанным парам атомов. Очевидно, что три самых первых пика (см. подробнее вставку на Рис. 5.4) соответствуют расстояниям: (1) — между атомами углерода и его водородами, (2) — между кислородом и его водородом, (3) — между углеродом и кислородом нашей молекулы. Вспомним, что эти расстояния были заданы в файле топологии молекулы, см. раздел [bondtypes]. Последующие внутримолекулярные пики (на интервале от 0.15 до 0.30 нм) соответствуют другим парам атомов молекулы: между углеродом и гидроксильным водородом, между кислородом и метильными водородами и, наконец, между водородами. Самый дальний внутримолекулярный пик на 0.3 нм. соответствует, очевидно, расстояниям между гидроксильным и метильными водородами. (Дополнительные вопросы: Где пик между метильными водородами? Почему дальние внутримолекулярные пики более размыты?).

Нетривиальным является наличие достаточно узкого пика на 0.2 нм на кривой для межмолекулярных расстояний. Казалось бы, что все они должны быть широкими. Однако этот пик соответствует расстоянию между кислородом и гидроксильным водородом другой молекулы, с которыми наш кислород образует водородную связь.

Парциальные функции радиального распределения

Для более детальной интерпретации структуры молекулярной жидкости рассчитывают парциальные ФРР, т.е. между выделенными группами атомов. Для работы с парциальными функциями необходимо создать индекс-файл с помощью команды:

make_ndx -f

Ключ –f подгружает файл conf.gro, из которого программа считывает данные об именах атомов. Программа работает в интерактивном режиме и приглашает к диалогу. На экране появляется следующая информация:

Reading structure file
Going to read 0 old index file(s)
Analysing residue names:
There are:   512    Protein residues
Analysing Protein...

  0 System             :  3072 atoms
  1 Protein            :  3072 atoms
  2 Protein-H          :  1024 atoms
  3 C-alpha            :     0 atoms
  4 Backbone           :   512 atoms
  5 MainChain          :   512 atoms
  6 MainChain+Cb       :   512 atoms
  7 MainChain+H        :  2048 atoms
  8 SideChain          :  1024 atoms
  9 SideChain-H        :   512 atoms

 nr: group       !   'name' nr name   'splitch' nr    Enter: list groups
 'a': atom        &   'del' nr         'splitres' nr   'l': list residues
 't': atom type   |   'keep' nr        'splitat' nr    'h': help
 'r': residue         'res' nr         'chain' char
 "name": group        'case': case sensitive           'q': save and quit
 'ri': residue index

Gromacs исторически ориентирован на работу с биологическими молекулами и по умолчанию воспринимает любую модель как белок, для которого создает стандартные группы атомов. Например, группа 1 (Protein) объединяет все атомы модели, а группа 2 (Protein-H) – это все «тяжелые» атомы, т.е водороды исключены. Имеется также список команд, с помощью которых можно, в частности, создавать и удалять группы.

Чтобы создать группу, объединяющую заданные атомы, надо знать имена этих атомов. Они записаны в файле молекулярной структуры. В нашем случае, если мы хотим создать группу из атомов кислородов, надо использовать имя "OA", которое указано в файле methanol.gro. Таким образом, достаточно выполнить команду:

> a OA

После нажатия Enter программа сообщит, что создана новая группа с номером 10 с именем ‘OA’, содержащая 512 атомов:

Found 512 atoms with name OA
 10 OA                 :   512 atoms

Аналогично, создадим группу углеродов:

> a C
Found 512 atoms with name C
 11 C                  :   512 atoms

группу метильных водородов:

> a H
Found 1536 atoms with name H
 12 H                  :  1536 atoms

и группу гидроксильных водородов:

> a HO
Found 512 atoms with name HO
 13 HO                 :   512 atoms

Для удобства работы можно удалить ненужные нам группы, созданные автоматически, начиная с Protein, заканчивая SideChain-H. (Группа 0 System, содержащая все атомы системы, нам понадобится в дальнейшем). Удаление групп с первой по девятую производится с помощью команды:

del 1-9

После удаления оставшиеся группы перенумеровываются.

Работа с программой make_ndx заканчивается нажатием клавиши q. После чего возникает файл с именем index.ndx, где для каждой группы перечислены номера соответствующих атомов согласно нумерации, заданной в файле conf.gro. В нашем случае индекс-файл должен выглядеть следующим образом:

[ 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

Для расчёта парциальной парной корреляционной функции для кислородов (gO-O(r)) выполняется команда:

g_rdf -n -o rdf-oo.xvg -w

Ключ -n подключает индекс-файл, выходной файл будет иметь имя rdf-oo.xvg. В интерактивном режиме выбираем дважды группу, относящуюся к кислородам OA:

Select a reference group and 1 group
Group     0 (         System) has  3072 elements
Group     1 (             OA) has   512 elements
Group     2 (              C) has   512 elements
Group     3 (              H) has  1536 elements
Group     4 (             HO) has   512 elements
Select a group: 1
Selected 1: 'OA'
Select a group: 1
Selected 1: 'OA'
Last frame        500 time  100.000

Аналогичным образом построим и запишем в файл rdf-cc.xvg функцию gC-C(r) для углеродов:

g_rdf -n -o rdf-cc.xvg

Нарисуем рассчитанные функции на одном графике, Рис. 5.6:

xmgrace rdf-all.xvg rdf-oo.xvg rdf-cc.xvg -legend load
Fig5.6.png
Рис. 5.6. Парциальные функции радиального распределения для атомов кислорода gO-O(r) (кривая О-О) и углерода gC-C(r) (кривая C-C) в сравнении с полноатомной межмолекулярной (кривая all, из Рис. 5.5) для жидкого метанола.

Обратим внимание, что первый пик gO-O(r) лежит примерно на 0.28 нм, т.е. примерно так же, как в воде, и соответствует водородно связанным атомам кислорода. Его положение определяется кулоновским взаимодействием соседних кислородов с атомом водорода, располагающимся между ними. В случае чистого леннард-джонсовского взаимодействия (см. работу №3 по моделированию жидкого аргона) положение первого пика в жидкости примерно соответствует величине параметра sigma, который иногда трактуют как "диаметр твердой сердцевины" атома. В нашем случае для кислорода метанола sigma = 0.35 нм, см. файл топологии, раздел [atomtypes]. Примерно такая же величина используется для кислорода воды, см. работу №4. Однако из-за дополнительного кулоновского взаимодействия центры соседних кислородов оказываются, в среднем, ближе друг к другу. Это связано с образованием водородной связи между соседними молекулами. Тот факт, что обсуждаемый пик узкий и мономодальный, показывает, что в метаноле практически все молекулы участвуют в образовании водородных связей.

Рассмотрим теперь первый пик функции gC-C(r). Он лежит примерно на 0.4 нм. Это больше, чем значение sigma для метанольного углерода, которое равно 0.35 нм, см. выше файл топологии. Однако это можно понять, в нашем случае "упаковываются" не углероды, а метильные группы, эффективный размер которых несколько больше, чем размер атома углерода.

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

g_rdf -n -o rdf-o_ho.xvg -s

Для сравнения рассчитаем также функцию gO-HС(r) между кислородами и водородами при углероде разных молекул:

g_rdf -n -o rdf-o_hc.xvg -s

Нарисуем эти функции вместе с полноатомной межмолекулярной g(r), Рис. 5.27:

xmgrace rdf-all.xvg rdf-o_ho.xvg rdf-o_hc.xvg –legend load
Fig5.7.png
Рис. 5.7. Парциальные функции для атомов кислорода с гидроксильными водородами gO-HO(r) (кривая О-Н(О)) и для кислорода с метильными водородами gO-HC(r) (кривая О-Н(С))) для жидкого метанола. Для сравнения показана полноатомная межмолекулярная g(r) (из Рис. 5.5)

На парциальной функции gO-HO(r) чётко виден пик при 0.2 нм. Очевидно, он соответствует расстоянию от кислорода до водорода соседней молекулы, с котором образовалась водородная связь, см. красные пунктирные линии на Рис. 5.8. Таким образом, с помощью расчёта парциальной функции мы убедительно подтвердили наше предположение о природе маленького пика при 0.2 нм на полноатомной g(r). (Заметим, что разная интенсивность одного и того же пика на данных кривых связана с нормировкой. В перовом случае функция нормируется на число всех пар кислород – гидроксильный водород другой молекулы, а в последнем – на полное число всех пар атомов в системе).

Fig5.8.png
Рис. 5.8. Водородносвязанные молекулы метилового спирта. Красным пунктиром показаны водородные связи. Тонкая пунктирная линия показывает расстояние между кислородом и гидроксильным водородом соседней молекулы.

Второй четкий пик функции gO-HO(r) с максимумом на 3.4 нм соответствует расстояниям между кислородом и гидроксильным водородом другой соседней молекулы, см. тонкий пунктир на Рис. 5.8. Направленность водородных связей обеспечивает преимущественную ориентацию молекул, в результате чего второе расстояние получается тоже довольно хорошо определённым.

Функция gO-HС(r) на Рис. 5.7 имеет совершенно иной вид. Между кислородом и метильными водородами другой молекулы нет специфических взаимодействий. Кроме того, здесь надо учитывать различные взаимные ориентации соседних молекул. Во-первых, в случае, показанном на Рис. 5.8, положения метильных протонов задаются геометрией водородных связей. Такие молекулы определяют относительно острый максимум при 4.3 нм. С другой стороны, молекула метана может быть повернута своей метильной группой к кислороду соседней молекулы. Именно такие ситуации дают вклад в широкую полосу от 2.5 до 4 нм.

Самодиффузия

Среднеквадратичное смещение атомов рассчитывается командой:

g_msd -n -w

При этом по умолчанию используется траектория traj.xtc и файл топологии topol.tpr. Ключ -n указывает, что будет использоваться индекс-файл. При работе программа выдаёт следующую информацию:

Reading file topol.tpr, VERSION 4.5.5 (single precision)
Select a group to calculate mean squared displacement for:
Group     0 (         System) has  3072 elements
Group     1 (             OA) has   512 elements
Group     2 (              C) has   512 elements
Group     3 (              H) has  1536 elements
Group     4 (             HO) has   512 elements

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

Select a group: 1
Selected 1: 'OA'

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

Last frame        500 time  100.000
Used 11 restart points spaced 10 ps over 100 ps
Fitting from 10 to 90 ps
D[        OA] 2.8672 (+/– 0.1371) 1e-5 cm^2/s
Executing 'xmgrace  msd.xvg &'

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

График среднеквадратического смещения сохраняется в файле msd.xvg. Ключ -w запускает автоматически программу xmgrace для отображения полученной кривой, Рис. 5.9.

Fig5.9.png
Рис. 5.9. Среднеквадратичное смещение молекул в жидком метаноле в зависимости от времени.

Водородные связи

Для установления водородных связей Gromacs использует геометрический критерий (как для воды). Напомним, водородная связь между водородом OH группы одной молекулы и кислородом O другой молекулы существует, если расстояние между кислородами rO-O ≤ 0,35 нм, и угол между направлением OH связи и линией кислород-кислород  ≤ 30º.

Водородные связи рассчитываются командой:

g_hbond

По умолчанию подгружается файл траектории traj.xtc и топологии topol.tpr. Программа работает в интерактивном режиме и выводит следующую информацию:

Reading file topol.tpr, VERSION 4.5.5 (single precision)
Specify 2 groups to analyze:
Group     0 (    
Group     1 (  ......

Необходимо выбрать, между какими группами мы хотим подсчитать количество водородных связей. В данном случае нас интересует вся система и связи внутри неё, то есть, надо два раза выбрать 0):

Select a group: 0
Selected 0: 'System'
Select a group: 0
Selected 0: 'System'

После завершения работы программа выдаёт следующую информацию:

Calculating hydrogen bonds in System (3072 atoms)
Found 512 donors and 512 acceptors
Reading frame       0 time    0.000
Will do grid-seach on 7x7x7 grid, rcut=0.35
Last frame        500 time  100.000
Average number of hbonds per timeframe 462.257 out of 131072 pos-sible

В последней строчке указано полное число водородных связей в системе. В данном случае программа нашла в среднем 462.257 водородных связей между 512 молекулами метанола (напомним, что усреднение идёт по всем конфигурациям траектории). Представляет интерес количество водородных связей, приходящихся на одну молекулу. Из полученных данных это легко найти, используя встроенный калькулятор bc (binary calculator):

echo "462.257 /512 * 2" | bc -l (это латинская буква l ).

(Множитель 2 здесь также означает тот факт, что каждая из найденных связей даёт вклад в две молекулы).

С помощью пакета VMD можно нарисовать водородные связи, Рис. 5.30. Сначала загрузим нашу модель.

vmd conf.gro

В окне Graphical Representations выберите способ отображения (Drawing Method): HBonds, в этом же окне появятся настройки критерия водородной связи – увеличьте критерий расстояния (Distance Cutoff) до 3,5 (Ǻ), угла (Angle Cutoff) до 30º. Также выберите подходящий для вас способ окраски в поле Coloring Method: режим ColorID позволяет задать конкретный цвет, например, синий. На экране появятся водородные связи. Чтобы отобразились ещё и сами молекулы, необходимо создать «копию» модели (нажать Create Rep), и для неё выставить Drawing Method – например, Licorice (лакричные палочки), а способ окраски Coloring Method – по имени Name.

Fig5.10.jpg
Рис. 5.10. Слева: фрагмент модели жидкого метанола. Водородные связи показаны пунктиром. Справа – соответствующее окно меню пакета VMD

В VMD можно нарисовать не всю толщу воды, а какой-то конкретный её участок. Например, чтобы отобразить молекулы в пределах 5 Å вокруг молекулы номер 10, необходимо в окне Graphical Representations в строке Selected Atoms написать: same resid as within 5 of resid 10 . Если у вас при этом ещё загружена траектория, то во вкладке Trajectory отметьте Update Selection every frame , тогда на каждом новом кадре соседние молекулы будут выбираться заново.