Изучение работы методов контроля температуры в молекулярной динамике¶
Построим файл топологии для молекулы этана
# загрузим модули
import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem import Descriptors
from rdkit import RDConfig
from IPython.display import Image,display
import numpy as np
import os
Используем потенциал 3 для торсионных углов, типы атомов opls_139 и opls_140 для атомов углерода и водорода ссответственно. Достроим файл: добавим все углы и двугранные углы.
top = '''
#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_139 1 ETH C1 1 -0.189 12.01
2 opls_139 1 ETH C2 2 -0.155 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
1 3 1
1 4 1
1 5 1
2 6 1
2 7 1
2 8 1
[ angles ]
; ai aj ak funct phi0 kphi
;around c1
3 1 4 1
4 1 5 1
3 1 5 1
2 1 3 1
2 1 4 1
2 1 5 1
;around c2
6 2 8 1
7 2 8 1
6 2 7 1
1 2 6 1
1 2 7 1
1 2 8 1
[ dihedrals ]
; ai aj ak al funct
3 1 2 6 3
3 1 2 7 3
3 1 2 8 3
4 1 2 6 3
4 1 2 7 3
4 1 2 8 3
5 1 2 6 3
5 1 2 7 3
5 1 2 8 3
[ pairs ]
; список атомов 1-4
; ai aj funct
3 6
3 7
3 8
4 6
4 7
4 8
5 6
5 7
5 8
[ System ]
; any text here
first one
[ molecules ]
;Name count
et 1
'''
with open('eth.top', 'w') as f:
f.write(top)
Использованы следующие методы контроля температуры:
be - метод Берендсена для контроля температуры.
vr - метод "Velocity rescale" для контроля температуры.
nh - метод Нуза-Хувера для контроля температуры.
an - метод Андерсена для контроля температуры.
sd - метод стохастической молекулярной динамики.
methods = ['be', 'vr', 'nh', 'an', 'sd']
for m in methods:
os.system('wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/{0}.mdp'.format(m))
Построим систему для симуляции МД и запустим GROMACS:
for m in methods:
os.system('grompp -f {0}.mdp -c et.gro -p et.top -o et_{0}.tpr'.format(m))
os.system('mdrun -deffnm et_{0} -v -nt 1'.format(m))
Построим анимацию (pdb) молекулярной динамики, извлечем потенциальную и кинетическую энергию для шагов симуляции с каждым методом контроля.
for m in methods:
os.system('echo 2 | trjconv -f et_{0}.trr -s et_{0}.tpr -o et_{0}.pdb'.format(m))
os.system('echo "10 11 \n \n" |g_energy -f et_{0}.edr -o et_{0}_en.xvg -xvg none'.format(m))
Опишем результат наблюдения динамики (визуализация pdb в PyMOL). В методе an молекула слабо колеблется по длине связи, в be -- быстро вращается вокруг связи, nh -- вращается вокруг связи, но слабо задерживается в некоторых позициях, sd -- быстро вращается вокруг центра масс, vr -- поддерживает заторможенную конформацию молекулы, иногда происходит поворот на 120 градусов вокруг связи. Предположительно, картина из моделирования с методом поддержания температуры vr лучше всего приближает реальную динамикку молекулы.
Построим зависимость потенциальной и кинетической энергии от времени:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
#t = {}
kin = {}
pot = {}
for m in methods:
a= np.loadtxt('et_%s_en.xvg' % m)
t=a[:,0]
pot[m]=a[:,1]
kin[m] =a[:,2]
methods_s = sorted(methods, key = lambda x : np.std(kin[x]), reverse = True)
plt.set_cmap('Accent')
for n, m in enumerate(methods_s):
plt.plot(t, pot[m], linewidth=0.5, label=m)
plt.legend(methods_s)
plt.ylim(0, 100)
plt.show()
plt.set_cmap('Accent')
for n, m in enumerate(methods_s):
plt.plot(t, kin[m], linewidth=0.5)
plt.legend(methods_s)
plt.ylim(0, 100)
plt.show()
np.mean(pot['an'])
Метод Андерсена работает некорректно и выдает низкие значения кинетической и потенциальной энергии (фиолетовая линия): далее не будем его рассматривать.
Остальные методы упорядочиваются по разбросу значений так: nh > (sd, vr) > be.
Извлечем значения длины C-C связи и построим гистограмму:
for m in methods:
os.system('g_bond -f et_{0}.trr -s et_{0}.tpr -o bond_{0}.xvg -n b.ndx -xvg none'.format(m))
import matplotlib.pyplot as plt
bond = {}
for m in methods:
a= np.loadtxt('bond_%s.xvg' % m)
bond[m]=a[:,1]
for m in methods:
plt.hist(bond[m], histtype = 'step', bins=10)
plt.legend(methods)
По полученным распределениям методы можно упорядочить по "вытянутости" распределения вправо (то есть "мягкости" ограничений, накладываемых на систему): an > vr > sd > be > nh.
Получим время (real) расчета GROMACS для этих алгоритмов контроля:
for m in methods:
with open("et_%s.log" % m) as f:
print(m, f.readlines()[-4].split()[2])
Отметим, что производительность методов слабо отличается (10-20%).
Подведем итоги. В рассмотренной системе производительность методов контроля температуры практически не отличается. Методы выдают разные результаты для распределений состояний системы во времени по кинетической\потенциальной энергии и по длине C-C связи, различие этих результатов легко описать как различие дисперсии распределения.
Метод Андерсена выдает некорректные (в ~50 меньше прочих методов) значения энергии системы.
Методы Нуза-Хувера и Берендсена дают сильно отличающиеся от прочих методов значения дисперсии: это объясняется тем, что эти методы не предназначены для описания систем с малым числом степеней свободы.
Методы стохастической молекулярной динамики и velocity rescale дают близкие значения распределений, средние среди прочих методов: вероятно, для рассматриваемой системы стоит использовать один из них: мне показалось, что результат симуляции velocity rescale ближе к ожидаемому движению молекулы этана.