Практикум №6. Изучение работы методов контроля температуры в GROMACS¶
1. Подготовка файлов с координатами и топологией молекулы этана¶
Из gro-файла для 38 молекул этана в жидкости был создан ndx-файл, а затем gro-файл для одной молекулы этана:
In [ ]:
make_ndx -f box_38.gro -o 1.ndx
# r1
# q
editconf -f box_38.gro -o et1.gro -n 1.ndx
# 3
Зададим ячейку так, чтобы молекула этана располагалась по центру и от нее был отступ 2 нанометра:
In [ ]:
editconf -f et1.gro -o et.gro -d 2 -c
Создадим файл топологии для одной молекулы этана:
In [ ]:
#include "/usr/share/gromacs/top/oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
et 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 opls_135 1 ETH C1 1 -0.0189 12.01
2 opls_135 1 ETH C2 2 -0.0155 12.01
3 opls_140 1 ETH H1 3 0.0059 1.008
4 opls_140 1 ETH H2 4 0.0059 1.008
5 opls_140 1 ETH H3 5 0.0059 1.008
6 opls_140 1 ETH H4 6 0.0056 1.008
7 opls_140 1 ETH H5 7 0.0056 1.008
8 opls_140 1 ETH H6 8 0.0056 1.008
[ bonds ]
; ai aj funct b0 kb
1 2 1 0.1554 250000.0
1 3 1 0.1085 300000.0
1 4 1 0.1085 300000.0
1 5 1 0.1085 300000.0
2 6 1 0.1085 300000.0
2 7 1 0.1085 300000.0
2 8 1 0.1085 300000.0
[ angles ]
; ai aj ak funct phi0 kphi
;around c1
2 1 3 1 111.291 400.400
2 1 4 1 111.291 400.400
2 1 5 1 111.291 400.400
3 1 4 1 111.291 200.400
3 1 5 1 111.291 200.400
4 1 5 1 111.291 200.400
;around c2
1 2 6 1 111.291 400.400
1 2 7 1 111.291 400.400
1 2 8 1 111.291 400.400
6 2 7 1 111.291 200.400
6 2 8 1 111.291 200.400
7 2 8 1 111.291 200.400
[ dihedrals ]
; ai aj ak al funct t0 kt mult
3 1 2 6 1 0.0 0.62760 3
3 1 2 7 1 0.0 0.62760 3
3 1 2 8 1 0.0 0.62760 3
4 1 2 6 1 0.0 0.62760 3
4 1 2 7 1 0.0 0.62760 3
4 1 2 8 1 0.0 0.62760 3
5 1 2 6 1 0.0 0.62760 3
5 1 2 7 1 0.0 0.62760 3
5 1 2 8 1 0.0 0.62760 3
[ pairs ]
; ai aj funct
3 6
3 7
3 8
4 6
4 7
4 8
5 6
5 7
5 8
[ System ]
first one
[ molecules ]
;Name count
et 1
2-5. Визуальный анализ¶
Файлы с параметрами контроля температуры разными методами:
Создадим файлы для молекулярно-механического движка, оптимизируем геометрию молекул этана и конвертируем соответствующие файлы в формат pdb:
In [ ]:
%%bash
for i in $(echo an be nh vr sd); do
grompp -f ${i}.mdp -c et.gro -p et.top -o et_${i}.tpr
mdrun -deffnm et_${i} -v -nt 1
trjconv -f et_${i}.trr -s et_${i}.tpr -o et_${i}.pdb
done
| Метод Андерсена для контроля температуры | Метод Берендсена для контроля температуры |
|---|---|
| Метод Нуза-Хувера для контроля температуры | Метод "Velocity rescale" для контроля температуры |
| Метод стохастической молекулярной динамики | |
- При контроле температуры методом Андерсена молекула этана совершает только небольшие колебательные движения.
- При контроле температуры методом Нуза-Хувера молекула этана совершает вращательные движения вокруг СС-связи, при этом остается на месте.
- При контроле температуры методом Берендсена молекула этана быстро вращается вокруг СС-связи.
- При контроле температуры методом "Velocity rescale" молекула этана как бы разгоняется и вращается не только вокруг СС-связи, но и вокруг других осей.
- При контроле температуры методом стохастической молекулярной динамики молекула этана быстро вращается вокруг различных "случайных" осей.
6-7. Сравнение потенциальной и кинетической энергий. Распределение длины СС-связи¶
Вычислим зависимость потенциальной и кинетической энергии каждой молекулы этана (в кДж/моль) от времени (в пикосекундах) и распределение длины СС-связи каждой молекулы этана во время моделирования (в нанометрах):
In [ ]:
%%bash
for i in $(echo an be nh vr sd); do
g_energy -f et_${i}.edr -o et_${i}_en.xvg
# 10
# 11
g_bond -f et_${i}.trr -s et_${i}.tpr -o bond_${i}.xvg -n b.ndx
done
# b.ndx:
# [ b ]
# 1 2
Построим соответствующие графики:
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import re
def read_et_en(filename):
xs_time = list()
ys_epot = list()
ys_ekin = list()
xvg = open(filename, 'r')
for line in xvg:
if line.startswith('#') == False and line.startswith('@') == False:
(time, epot, ekin) = re.split(' +', line.strip())
xs_time.append(time)
ys_epot.append(epot)
ys_ekin.append(ekin)
xvg.close()
xs_time = map(float, xs_time)
ys_epot = map(float, ys_epot)
ys_ekin = map(float, ys_ekin)
return xs_time, ys_epot, ys_ekin
def read_bond(filename):
xs = list()
ys = list()
xvg = open(filename, 'r')
for line in xvg:
if line.startswith('#') == False and line.startswith('@') == False:
(x, y) = re.split(' +', line.strip())
xs.append(x)
ys.append(y)
xvg.close()
xs = map(float, xs)
ys = map(float, ys)
return xs, ys
titles = {
'an': 'Anderson method of temperature control',
'be': 'Berendsen method of temperature control',
'nh': 'Nose-Hoover method of temperature control',
'vr': '"Velocity rescale" method of temperature control',
'sd': 'Stochastic molecular dynamics method'}
i = 1
plt.figure(figsize=(15, 25))
for method in ('an', 'be', 'nh', 'vr', 'sd'):
xs_time, ys_epot, ys_ekin = read_et_en('methods/et_%s_en.xvg'%method)
xs, ys = read_bond('methods/bond_%s.xvg'%method)
plt.subplot(5, 2, i)
plt.plot(xs_time, ys_epot, color='blue', marker='o', ls='*', label='Potential')
plt.plot(xs_time, ys_ekin, color='lightgreen', marker='o', ls='*', label='Kinetic')
plt.legend(loc='upper center')
plt.xlabel("Time (ps)")
plt.ylabel("Energy (kJ/mol)")
plt.title(titles[method])
plt.subplot(5, 2, i+1)
plt.bar(xs, ys, 0.0002, color='grey', edgecolor='none')
plt.xlabel("CC-bond length (nm)")
plt.title(titles[method])
i += 2
На приведенных выше графиках видно, что при контроле темтературы методом Андерсена поведение молекулы не описывается распределением Больцмана, а внутренняя энергия системы близка к нулю. При контроле температуры методом Нуза-Хувера внутренняя энергия системы в целом выше, чем при контроле температуры другими методами, а также наблюдаются аномально высокие значения внутренней энергии системы (~400 кДж/моль) и аномально большие длины СС-связи (~2 А). При контроле температуры методом Берендсена потенциальная и кинетическая энергии находятся в узком диапазоне и длина СС-связи меняется незначительно, что кажется не очень правдоподобным. Оставшиеся метод "Velocity rescale" и метод стохастической молекулярной динамики, на мой взгляд, наиболее подходят для контроля температуры.
8. Выводы¶
- Контроль температуры методом Андерсена предполагает случайные столкновения молекул, в результате которых молекулы приобретают новые скорости, при этом температура системы остается постоянной за счет обмена энергией с термостатом (канонический ансамбль); данный метод не годится для системы, состоящей из одной молекулы, т.к. одной молекуле не с кем сталкиваться.
- При контроле температуры методом "Velocity rescale" температура системы остается постоянной за счет пересчета скоростей молекул на каждом шаге, что дает неплохие результаты моделирования при визуальном анализе и при анализе постоенных графиков.
- При контроле температуры методом Берендсена пересчет скоростей молекул производится не на каждом шаге, а реже, что допускает небольшие колебания температуры в системе, при этом энергия системы близка к постоянной (микроканонический ансамбль).
- При контроле температуры методом Нуза-Хувера изменения энергии системы (канонический ансамбль) детерменированы "эффективной массой тепловой бани"; данный метод не годится для системы, состоящей из одной молекулы, т.к. наблюдаются аномально высокие значения энергии системы.
- Контроль температуры методом стохастической молекулярной динамики предполагает, что на систему действуют случайные силы, зависящие от температуры и аппроксимируемые Гауссовой функцией. По-видимому, данный метод является наиболее универсальным для моделирования.
- Таким образом, случайный фактор допускается при контроле температуры методом Андерсена и методом стохастической молекулярной динамики, в отличие от трех остальных рассмотренных методов.