Работа 4

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

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

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

Цель:

Провести моделирование жидкой воды. Исследовать её свойства. Сравнить результаты для разных моделей воды.

Задания:

  1. Получить SPC-модель воды.
  2. Получить Tip4p-модель воды.

Для обеих моделей:

  1. создать исходный файл конфигурации с 1000 молекулами воды.
  2. провести релаксацию.
  3. получить равновесную МД-траекторию.
  4. рассчитать функцию радиального распределения.
  5. рассчитать коэффициент самодиффузии молекулы воды.
  6. рассчитать число водородных связей, приходящихся на молекулу.
  7. Сравнить ФРР, коэффициент самодиффузии и среднее число водородных связей на молекулу для обеих моделей.
  8. Рассчитать плотности обеих моделей для разных температур и найти положение максимума плотности.

Вопросы

  1. С какой моделью воды (трехцентровой или четырехцентровой) расчёт быстрее?
  2. В чем отличие ФРР воды от ФРР простой жидкости?
  3. Почему у воды имеется максимум плотности?

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

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

Fig4.1.jpg
Рис. 4.1. Иллюстрация трехцентровой (SPC) и четырехцентровой (Tip4p) моделей воды. (Рисунки с сайта http://www.sklogwiki.org/SklogWiki/).


Для моделирования жидкой воды существуют различные эмпирические потенциалы взаимодействия между молекулами воды – различные модели воды. Наиболее распространёнными в настоящее время в молекулярной динамике являются «трёхчастичные» (Tip3p, SPC) и «четырёхчастичные» (Tip4p) модели воды, см. Табл. 4.1.

Модель SPC (simple point charge model) представляет собой жёсткий треугольник с зарядами на каждом из атомов, Рис. 4.14(слева). Заряд qO располагается на центре атома кислорода, а на протонах находятся положительные заряды величиной qO / 2. Кроме кулоновского взаимодействия имеется леннард-джонсовское с центром на атоме кислорода. Таким образом, здесь имеется только три центра взаимодействия: центр кислорода и два протона.

В четырехцентровых моделях Tip4p (4-point-transferable-intermolecular-potential). В ней положение отрицательного заряда не совпадает с центром атома кислорода, Рис. 4.14(справа), отрицательный заряд (точка M) лежит на заданном расстоянии rOM от центра атома кислорода, располагаясь в одной плоскости с кислородом и протонами.

Табл. 4.1. Значения параметров некоторых моделей воды:

σ Å
ε kJ mol-1
rOH Å
HOH°
qO (e)
qH (e)
qM (e)
rOM Å
SPC
3.166
0.650
1.0000
109.47
-0.8200
+0.410
0
Tip4p
3.154
0.648
0.9572
104.52
0
+0.52
-1.04
0.1546
Tip4p-Ew
3.154
0.648
0.9572
104.52
0
+0.52422
−1.04844
0.125
Tip4p-2005
3.154
0.648
0.9572
104.52
0
+0.5564
−1.1128
0.1546

Тетраэдрическая координация молекул воды при моделировании обеспечивается кулоновскими взаимодействиями. Это соответствует реальной физике. Квантово-механические расчёты показывают, что водородная связь имеет, в значительной мере, электростатическую природу. С другой стороны, такой простой учёт ориентационно-зависимого взаимодействия между молекулами воды позволяет существенно ускорить МД-расчёты.

Поскольку водородная связь между молекулами воды не определена явно, то для выявления водородносвязанных атомов используют геометрический критерий. Известно, что при оптимальной водородной связи атом водорода, через который образуется связь, расположен на линии, соединяющей центры кислородов соседних молекул воды (как во льду). Однако в жидкости возможны большие отклонения от этого угла. Считается, что молекулы образуют водородную связь, если угол отклонения атома водорода от линии, соединяющей кислороды, не превышает 30 градусов. Кроме того, молекулы должны быть расположены достаточно близко, т.е. нужен ещё один параметр – расстояние между атомами кислорода. Считается, что оно не должно превышать 0.35 nm, Рис. 4.2.

Fig4.2.png
Рис. 4.2. Геометрический критерий водородной связи, используемый в Gromacs: r<0.35 nm, α < 30°.


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

Для выполнения данной работы в директории 04_water/ подготовлено следующее дерево поддиректорий, Рис. 4.3.

Fig4.3.jpg
Рис. 4.3. Дерево директорий для работы по моделированию и анализу воды.


Предлагается работать с каждой моделью воды в своей директории spc/ и tip4p-2005/. В каждой из них лежит файл топологии для данной модели воды (topol.top) и файл с общими параметрами моделирования (grompp.mdp).

Поддиректории для разных моделей имеют одинаковые имена. Директория preparing/ служит для подготовки стартовой конфигурации, в ней лежит файл молекулярной структуры для одной молекулы воды (h2o.gro). Директории equi1/ и equi2/ предназначены для двух последовательных релаксаций. В каждой из них лежит свой файл grompp.mdp с параметрами моделирования для данного этапа релаксации.


SPC-вода

Подготовка начальной конфигурации модели

Координаты атомов одной молекулы воды заданы в файле h2o.gro в директории preparing. Содержимое этого файла (less preparing/h2o.gro):

216H2O,WATJP01,SPC216,SPC-MODEL,300K,BOX(M)=1.86206NM,WFVG,MAR. 1984
    3
    1SOL     OW    1   0.186   0.169   0.174
    1SOL    HW1    2   0.093   0.167   0.211
    1SOL    HW2    3   0.187   0.130   0.082
   0.31143   0.31143   0.31143

Файл топологии для SPC-молекулы воды topol.top находится в директории spc/.

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

#include "ffgmx.itp"
#include "spc.itp"

[ system ]
Pure Water

[ molecules ]
SOL     1000

Оператор include подключает файлы ffgmx.itp и spc.itp, которые находятся в библиотеке Gromacs (директория /usr/share/gromacs/top/).

В первом файле ffgmx.itp и в файлах, включённых в него, записаны типы атомов и параметры взаимодействий. Во втором файле spc.itp указаны имена атомов (используемые в файле .gro) и типы атомов, которым они соответствуют. Эти файлы содержат информацию о взаимодействии между атомами.

/usr/share/gromacs/top/ffgmx.itp
; This file is for backward compatibility only.
; Please directly include the file below in your .top file.
#include "gmx.ff/forcefield.itp"
/usr/share/gromacs/top/gmx.ff/forcefield.itp
#define _FF_GROMACS
#define _FF_GROMACS1

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
  1             1               no              1.0     1.0

#include "ffnonbonded.itp"
#include "ffbonded.itp"
/usr/share/gromacs/top/spc.itp
; This file is for backward compatibility only.
; Please directly include the spc.itp file from your force field directo-ry.
#ifdef _FF_GROMACS
#include "gmx.ff/spc.itp"
#endif
#ifdef _FF_GROMOS96
#include "gromos43a1.ff/spc.itp"
#endif
#ifdef _FF_OPLS
#include "oplsaa.ff/spc.itp"
#endif
/usr/share/gromacs/top/gmx.ff/spc.itp
[ moleculetype ]
; molname       nrexcl
SOL             2

[ atoms ]
;   nr   type  resnr residue  atom   cgnr     charge       mass
     1     OW      1    SOL     OW      1      -0.82
     2     HW      1    SOL    HW1      1       0.41
     3     HW      1    SOL    HW2      1       0.41

#ifndef FLEXIBLE
[ settles ]
; OW    funct   doh     dhh
1       1       0.1     0.16330

[ exclusions ]
1       2       3
2       1       3
3       1       2
#else
[ bonds ]
; i     j       funct   length  force.c.
1       2       1       0.1     345000  0.1     345000
1       3       1       0.1     345000  0.1     345000

[ angles ]
; i     j       k       funct   angle   force.c.
2       1       3       1       109.47  383     109.47  383
#endif

Размножим молекулу внутри модельного бокса. Для этого переходим в директорию preparing:

cd preparing

и сделаем из неё ссылку на файл топологии:

ln -s ../topol.top .

Используем команду genconf для создания модельного бокса с 1000 молекулами воды:

genconf -f h2o.gro -nbox 10 10 10 -o h2o_1000.gro

ключ -f указывает, что используется молекула из файла h2o.gro, ключ -nbox с параметрами 10 10 10 задаёт, что модельный кубик с одной молекулой воды оттранслируется по 10 раз вдоль каждого ребра, т.е. мы будем иметь большой бокс, содержащий 1000 молекул воды, ключ -o определяет имя выходного файла, назовём его h2o_1000.gro.

Можно посмотреть получившуюся конфигурацию с помощью VMD:

vmd h2o_1000.gro

При выполнении команды genconf происходит автоматическое изменение файла топологии, т.е. файл topol.top теперь соответствует конфигурации из 1000 молекул воды. Для этого последняя строка SOL 1 заменилась на SOL 1000.

Релаксация модели. Первый этап

В полученной начальной конфигурации молекул воды расположены регулярно, но не оптимально. На первом этапе релаксации сделаем её более оптимальной. При этом зададим, чтобы плотность и температура соответствовали планируемым значениям, в данном случае – давлению P=1бар и температуре T=300K. Фактически, мы проводим обычное МД-моделирование в NPT ансамбле. Цель данного шага — быстрее перевести систему в заданное состояние. Для этого обычно используют термостат v-rescale и баростат Берендсена.

Переходим в директорию equi1:

cd ../equi1/

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

Посмотрим его содержимое:

;       Input file
;
title		= water	; a string
; Run parameters
integrator	= md		; leap-frog integrator
nsteps	= 5000	; 10 ps
dt		= 0.002	; 2 fs
; Output control
nstxout	= 0		; save coordinates every 0 ps
nstvout	= 0		; save velocities every 0 ps
nstxtcout	= 100		; save coord-s in compressed format every 0.2 ps
nstenergy	= 100		; save energies every 0.2 ps
nstlog	= 100		; update log file every 0.2 ps
; Neighborsearching
ns_type	= grid	; search neighboring grid cels
nstlist	= 5		; 10 fs
rlist		= 1.2		; short-range neighborlist cutoff (in nm)
rcoulomb	= 1.2		; short-range electrostatic cutoff (in nm)
rvdw		= 1.2		; short-range van der Waals cutoff (in nm)
DispCorr	= EnerPres	; long range disp. corr. for Energy and Pressure
; Electrostatics
coulombtype = PME	; Particle Mesh Ewald for long-range elec-trostatics
pme_order	= 4		; cubic interpolation
fourierspacing	= 0.12	; grid spacing for FFT
optimize_fft	= yes
ewald_rtol	= 1.0e-5
; Temperature coupling is on
tcoupl	= V-rescale	; modified Berendsen thermostat
tc-grps	= System	; two coupling groups – more accurate
tau_t		= 0.5		; time constant, in ps
ref_t		= 300		; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl	= Berendsen	; Pressure coupling on in NPT
pcoupltype	= isotropic	; uniform scaling of box vectors
tau_p		= 2.0		; time constant, in ps
ref_p		= 1.0		; reference pressure, in bar
compressibility = 4.5e-5	; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc		= xyz		; 3-D PBC
; Velocity generation
gen_vel	= yes		; Velocity generation is off
gen_temp	= 300		; initial temperature
gen_seed	= 2011 	; random seed

Удобно работать со стандартными именами файлов (использовать имена по умолчанию). Для этого создадим соответствующие ссылки на файлы, подготовленные в директории preparing:

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

а также на файл топологии

ln -s ../topol.top .

Для создания стартового файла проведём «препроцессинг» с помощью команды grompp , которую теперь можно использовать без параметров:

grompp

В результате получается стартовый файл с именем topol.tpr, в который собраны все необходимые данные для МД-моделирования. Это бинарный файл, который можно посмотреть с помощью команды gmxdump:

gmxdump -s topol.tpr | less

Наконец, запускаем МД-моделирование, осуществляющее первый этап релаксации. Полученный стартовый файл topol.tpr подключается по умолчанию:

mdrun –nt 8 –dd 2 2 2 -nice 11 -v

Ключи -nt и -dd указывают, что мы для вычислений используем 8 ядер и разбиваем модельный бокс на 8 частей (2х2х2). Однако при обычной работе их можно не использовать. Напомним, что современные версии Gromacs'a умеют сами оптимальным образом выбирать ядра. Приоритет работы программы моделирования можно установить чуть ниже, чем у других программ, загруженных на компьютере (графической среды или текстового редактора), чтобы не замедлять работу с ними. Стандартный приоритет соответствует числу 10 параметра ключа -nice. Ключ -v включает режим вывода на экран текущей информации о работе программы (оставшееся время до окончания расчёта и предупреждения).

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

Для контроля проведённой релаксации можно посмотреть, как изменилась потенциальная энергия системы. Эти данные записаны в файле ener.edr. Используем команду g_energy, которая позволяет выбрать желаемую характеристику в интерактивном режиме:

g_energy -o potential.xvg -w

В файл potential.xvg будет записана выбранная характеристика (потенциальная энергия) Расширение .xvg ассоциирует этот файл со стандартной графической программой xmgrace. Ключ -w рисует автоматически график выбранной характеристики с помощью программы xmgrace.

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

Релаксация модели. Второй этап

В файле confout.gro мы имеем конфигурацию модели, близкую к равновесной. Рекомендуется продолжить релаксацию, но уже с параметрами основного моделирования. Включаем термостат Нозе-Хувера и баростат Паринелло-Рамана, вместо баростата и термостата Берендсена (или его разновидности – термостата v-rescale). Напомним, что термостат Нозе-Хувера и баростат Паринелло-Рамана дают более правильные флуктуации системы в равновесном состоянии, но плохо работают с неравновесной системой.

Переходим в директорию equi2. Там находится свой файл grompp.mdp. Различия с одноименным файлом в директории equi1 можно посмотреть с помощью команды vimdiff:

vimdiff grompp.mdp ../equi1/grompp.mdp

(выход:q ) Посмотрите, в чем эти различия.

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

ln -s ../equi1/confout.gro conf.gro
ln -s ../topol.top .
grompp

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

mdrun -nice 11 -v

Полученный файл конечной конфигурации confout.gro будет начальной конфигурацией для расчёта равновесной траектории.

Основное моделирование (production run)

Рассчитаем теперь равновесную траекторию для SPC-воды. Возвращаемся в папку spc/:

cd ../

Там находится файл параметров моделирования. Убеждаемся, что параметры совпадают с параметрами, используемыми на втором этапе релаксации. Задаем время моделирования 100 ps. Задаем ссылку на стартовый файл структуры:

ln -s equi2/confout.gro conf.gro

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

grompp
mdrun –nt 8 -nice 11 -v

После выполнения программы в директории spc будут записаны файлы равновесной молекулярно-динамической траектории: traj.xtc, ener.edr и др.

Tip4p -вода

Перейдем к работе с четырёхцентровой моделью воды. Существует несколько вариаций этой модели, предложенных в разное время. Они незначительно различаются значениями зарядов и величиной rOM — сдвиг отрицательного заряда от центра кислорода, см. Табл. 4.1. Наилучшей в настоящее время считается модель tip4p-2005, хотя и предыдущие версии тоже неплохо описывают разные свойства воды, в частности, разумно воспроизводят фазовые диаграммы воды (льдов) и положение максимума плотности воды.

Получение МД-модели для Tip4p-воды производится в точности так же, как было показано для SPC-воды. Отличие в том, что теперь работа идёт в поддиректориях директории tip4p-2005/, где находятся соответствующие файлы topol.top и h2o.gro , отличающиеся от SPC воды. Ниже показано содержимое файла координат h2o.gro в директории tip4p-2005/preparing и файла topol.top, находящегося в директории tip4p-2005.

h2o.gro
;tip4p water model
    4
    1SOL     OW    1   1.736   0.839   0.257
    1SOL    HW1    2   1.777   0.781   0.322
    1SOL    HW2    3   1.643   0.831   0.247
    1SOL     MW    4   1.730   0.831   0.267
   0.31143   0.31143   0.31143
topol.top
;
; TIP4P/2005 Water
;

[ defaults ]
; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ
  1		2		no		1.0	1.0 

#include "tip4p2005.itp"

[ system ]
water tip4p/2005

[ molecules ]
SOL 1000

Обратите внимание, что в файле топологии присутствует "включённый" файл tip4p2005.itp, в котором, собственно, содержится вся топология модели TIP4P-2005. Для поиска включённых файлов препроцессор grompp обращается сначала к стандартной директории /usr/share/gromacs/top/, если же там их не обнаруживает (как будет в данном случае), то ищет в текущей директории. Поэтому необходимо, чтобы файл tip4p2005.itp лежал в текущей директории.

Содержимое этого файла:

tip4p2005.itp
; J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005).

[ atomtypes ]
; name      at.num  mass     charge ptype  sigma      epsilon
OW_tip4p2005 8     15.9994  0.0000  A   3.15890e-01  7.749e-01
HW_tip4p2005 1     1.008    0.0000  A   0.00000e+00  0.000e+00
MW_tip4p2005 0     0.000    0.0000  A   0.00000e+00  0.000e+00

[ moleculetype ]
; molname	nrexcl
SOL		1

[ atoms ]
; 	at_type      res_nr  res_name	at_name   cg_nr	charge     mass
1     OW_tip4p2005 	1    SOL	  OW		1	0         15.9994
2	HW_tip4p2005  	1    SOL	 HW1		1	0.5564    1.00800
3	HW_tip4p2005	1    SOL	 HW2		1	0.5564    1.00800
4	MW_tip4p2005	1    SOL	  MW		1    -1.1128    0.00000

[ settles ]
; i	funct	doh	dhh ---> same as standard tip4p
1	1	0.09572	0.15139

[ virtual_sites3 ]
; Dummy from			funct	a		b
4	1	2	3	1	0.131937768	0.131937768

[ exclusions ]
1	2	3	4
2	1	3	4
3	1	2	4
4	1	2	3

; The position of the dummy is computed as follows:
;
;		O
;  	  
;	    	D
;	
;	H		H
;
; tip4p/2005:
; const = distance (OD) / [ cos (angle(DOH)) 	* distance (OH) ]
;	  0.01546 nm	/ [ cos (52.26 deg)	* 0.09572 nm	]
;	then a = b = 0.5 * const
;
; Dummy pos x4 = x1 + a*(x2-x1) + b*(x3-x1)

Отличие от предыдущей SPC-модели в том, что здесь, кроме атомов кислорода и двух протонов, имеется дополнительный виртуальный центр с именем MW. Видно, что отрицательный заряд находится не на кислороде, а на этом виртуальном центре.

Если все Ваши предыдущие действия, совершаемые с моделью SPC, Вы сохраняли в файл с именем, скажем, run.sh, то теперь его можно сделать исполняемым:

chmod +x run.sh

и запустить в директории tip4p-2005:

./run.sh

Если же не сохраняли, то историю команд, запускаемых в командной строке, можно посмотреть в файле less ~/.bash_history.

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

Анализ производить для каждой модели в своей директории (т.е в spc/, и в tip4p-2005/).

Расчёт функции радиального распределения

Между молекулами воды существует три типа расстояний: кислород-кислород, кислород-водород и водород-водород. При выполнении команды g_rdf по умолчанию рассчитываются все эти расстояния. Однако пространственное расположение молекул воды будет более понятным, если рассматривать корреляции только между атомами кислорода. В этом случае мы должны создать индекс-файл, определяющий группу атомов кислорода (см. ниже). Однако есть возможность рассчитать функцию радиального распределения между центрами масс молекул. Для этого достаточно воспользоваться ключом -rdf mol_com. Как показывают расчёты, для молекул воды оба способа дают примерно одинаковые функции, поскольку кислород является здесь доминирующим атомом. Рассчитаем, для примера, парную корреляционную функцию воды g(r) между центрами масс молекул воды:

g_rdf -s -rdf mol_com –w

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


Fig4.4.png
Рис. 4.4. Функции g(r) для SPC tip4p-2005 воды при T= 300К и P = 1 бар, рассчитанные относительно центров масс молекул.


На Рис. 4.4 показаны функции g(r), рассчитанные для разных моделей воды (SPC и Tip4p-2005). Заметим, что Tip4p-2005 лучше описывает экспериментальную g(r), чем SPC, у которой излишне размыты последующие осцилляции (второй и третий максимумы).

Расчёт коэффициента самодиффузии молекулы воды

Рассчитаем среднеквадратичное смещение молекул воды в зависимости от времени. Но прежде нужно создать индекс-файл, содержащий номера атомов кислорода. Для этого воспользуемся программой make_ndx:

make_ndx -f

Здесь ключ -f означает, что мы используем файл координат, называющийся по умолчанию conf.gro, из которого берётся информация об имени атомов и их порядковом номере. У программы make_ndx предусмотрен интерактивный режим. На экране возникает информация:

Reading structure file
Going to read 0 old index file(s)
Analysing residue names:
There are:  1000      Water residues
  0 System             :  3000 atoms
  1 Water              :  3000 atoms
  2 SOL                :  3000 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

Мы должны создать группу из атомов кислородов воды, которые у нас называются OW. Для этого используем команду:

> a OW

после чего получим сообщение:

Found 1000 atoms with name OW
3 OW : 1000 atoms
Fig4.5.png
Рис. 4.5. Зависимость среднеквадратичного смещения молекул воды от времени. (Для SPC и tip4p-2005 моделей, Т=300К, давление 1 бар).


После нажатия клавиши q программа запишет индекс-файл и выйдет из интерактивного режима. В директории появится файл index.ndx .

Расчёт среднеквадратичного смещения и коэффициента самодиффузии производится командой g_msd:

g_msd -n –w

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

Среднеквадратичное смещение с хорошей точностью описывается прямой линией, Рис. 4.5. Коэффициент самодиффузии определяется из наклона этой прямой. Видно, что для SPC он почти в два раза больше, чем для Tip4p-2005. Его значение рассчитывается автоматически программой g_msd и отображается на экране компьютера в последней строке сообщения:

Select a group to calculate mean squared displacement for:
Group 0 ( OW) has 1000 elements
There is one group in the index
Last frame 500 time 100.000
Used 11 restart points spaced 10 ps over 100 ps
Fitting from 10 to 90 ps
D[ O*] 4.1406 (+/– 0.0902) 1e-5 cm^2/s

Аналогичным образом для модели tip4p/2005 получается:

D[ OW] = 2.2842 (+/– 0.0152) (1e-5 cm^2/s)

Экспериментальное значение коэффициента самодиффузии чистой воды равно 2.30·10-5 см2/с при 25°С. Таким образом, модель tip4p/2005 лучше описывает диффузию молекул воды, чем SPC модель.

Расчёт среднего числа водородных связей на молекулу

Команда g_hbonds выдаёт полное число водородных связей в системе (см. подсказку g_hbond -h). Для расчёта по умолчанию используется траектория traj.xtc, а также файл топологии topol.tpr. В интерактивном режиме программа запрашивает, между какими группами требуется рассчитать водородные связи. В нашем случае все предлагаемые группы идентичны и соответствуют молекулам воды. Можно выбирать любую, указав её дважды, например, дважды выбрать 0.

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

echo "0 0" | g_hbond

После выполнения команды получим отчёт, содержащий в последней строке среднее число водородных связей в модели:

Reading file topol.tpr, VERSION 4.5.5 (single precision)
Specify 2 groups to analyze:
Group 0 ( System) has 3000 elements
Group 1 ( Water) has 3000 elements
Group 2 ( SOL) has 3000 elements
Select a group: 0
Selected 0: 'System'
Select a group: 0
Selected 0: 'System'
Calculating hydrogen bonds in System (3000 atoms)
Found 1000 donors and 1000 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 1728.880 out of 500000 possible

Химики обычно говорят о среднем числе водородных связей на молекулу воды. Зная число молекул в модели, это легко вычислить. Для простых вычислений можно использовать встроенный калькулятор bc (binary calculator). Например:

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

Отличаются ли наши модели воды по количеству водородных связей на молекулу?

Расчёт плотности

Мгновенная плотность (плотность данной конфигурации) рассчитывается как отношение массы молекул воды к объёму модельного бокса. Поскольку мы работаем в NPT-ансамбле, где объём бокса меняется для разных конфигураций, то плотность системы флуктуирует, см. Рис. 4.6. Средняя плотность системы рассчитывается усреднением по всей МД-траектории. Все эти данные находятся в файле ener.edr и доступны по команде:

g_energy -o density -w > output.density

После ключа -o стоит имя файла, в который будет записан график. В данном случае это density.xvg – расширение .xvg подставляется по умолчанию. Значок > направляет выходную информацию в указанный файл output.density. Там запишется табличка с кратким отчётом о средней плотности и её разбросе:

Statistics over 50001 steps [ 0.0000 through 100.0000 ps ], 1 data sets
All statistics are over 10001 points

Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
Density 995.314 1.2 7.25942 -1.92807 (kg/m^3)
Fig4.6.png
Рис. 4.6. Зависимость плотности воды для моделей SPC и tip4-2005 от времени в процессе моделирования в NPT ансамбле. Экспериментальное значение плотности воды при температуре 300 К равно 996,56 кг/м3. (см. сайт webbook.nist.gov/chemistry/fluid).


Изучение поведения плотности воды в зависимости от температуры

Основная особенность воды – максимум плотности при 4ºС. Воспроизводят ли это свойство воды наши модели? Для этого надо построить зависимость плотности воды от температуры для обеих моделей, SPC и Tip4p-2005. Проведём моделирование и получим молекулярно-динамические траектории воды для разных значений температур (от 260 до 300 К с шагом 10 К) и постоянным давлением 1 bar. Найти плотность и сравнить рассчитанные зависимости с экспериментальными данными.

Для ускорения работы предлагается автоматизировать расчёты, используя скрипт. С его помощью можно перебрать все температуры с помощью цикла:

for i in $(seq 260 10 300)
do
mkdir $i
cd $i
......
......
cd ../
done

В цикле необходимо менять значение температуры в файлах grompp.mdp – для этого создадим файл-шаблон (например, с именем equi1_grompp.mdp_tmpl):

cp grompp.mdp equi1_grompp.mdp_tmpl

В нём (при помощи текстового редактора) запишем вместо значения температуры любую узнаваемую последовательность символов, например @i@. Затем при помощи команды построчной обработки файлов, sed, будем заменять @i@ на нужную нам температуру $i:

cat ../equi1_grompp.mdp_tmpl | sed -e “/@i@/$i/” > grompp.mdp

Здесь программа cat подаёт строчки из файла equi1_grompp.mdp_tmpl (находящегося на директорию выше ../), на вход команды sed. А результат направляется в файл grompp.mdp. (sed обрабатывает файл построчно, ищет в строке шаблон @i@ и заменяет его на значение переменной $i )

Итак, необходимо создать 3 шаблона – к двум этапам релаксации и к основному моделированию. Обратите внимание, что первый этап релаксации лучше сделать длиннее, например, 30 пс (так как системе будет требоваться больше времени, чтобы отрелаксировать с температуры 300 до 260 К). Во всех шаблонах число, соответствующее температуре термостата ref_t необходимо заменить на @i@. Также на первом этапе релаксации необходимо вставить ту же узнаваемую последовательность в поле gen_temp, соответствующее генерации скоростей в первый момент времени.

Экспериментальные значения зависимости плотности воды от температуры при нормальном давлении можно найти, например, в справочнике nist.gov:

http://webbook.nist.gov/chemistry/fluid/

или Wagner,Pruss,J.Phys.Chem.Ref.Data,2002,31,2,387-535

Для построения графиков можно использовать qtiplot. На Рис. 4.7 собраны рассчитанные зависимости плотности от температуры.

Fig4.7.png
Рис. 4.7. Зависимость плотности воды от температуры при нормальном давлении. Красные точки – эксперимент, зеленые — SPC-модель, синие – tip4p-2005 модель.


Таким образом, сравнивая наши модели воды, можно опять сделать вывод в пользу четырехцентровой модели. Тем не менее, модель SPC широко используется в современных исследованиях, главным образом, при моделировании растворов макромолекул в воде. Для этого есть разумные причины. Во-первых, работа с трёхцентровой моделью требует меньше вычислительного времени, что важно, поскольку в окружение макромолекулы входит большое число молекул воды (нередко это десятки тысяч). Кроме того, в таких расчётах основной интерес представляет сама исследуемая молекула, а тонкости описания воды здесь не так важны.

Пример скрипта

density.sh:
for i in $(seq 260 5 280)
do
  mkdir $i
  cd $i
    mkdir equi1
    cd equi1
        ln -s ../../confout.gro conf.gro
        ln -s ../../topol.top .
        ln -s ../../tip4p2005.itp .
        cat ../../equi1/grompp.mdp_tmpl | sed -e "s/@i@/$i/" > grompp.mdp

        grompp
        mdrun -nt 8 -dd 2 2 2  -nice 11 -v
    cd ../

    mkdir equi2
    cd equi2
        ln -s ../equi1/confout.gro conf.gro
        ln -s ../../topol.top .
        ln -s ../../tip4p2005.itp .
        cat ../../equi2/grompp.mdp_tmpl | sed -e "s/@i@/$i/" > grompp.mdp
        grompp
        mdrun -nt 8 -dd 2 2 2  -nice 11 -v
    cd ../
    ln -s equi2/confout.gro conf.gro
    cat ../grompp.mdp_tmpl | sed -e "s/@i@/$i/" > grompp.mdp
    ln -s ../topol.top .
    ln -s ../tip4p2005.itp .
    grompp
    mdrun -nt 8 -dd 2 2 2  -nice 11 -v
    echo 'Density' | g_energy -o density.xvg > output.density

cd ../
done