Изучение работы методов контроля температуры в молекулярной динамике¶
import numpy as np
import matplotlib.pyplot as plt
import time
import os
%matplotlib inline
from IPython.display import Image,display
Используем потенциал 3 для торсионных углов, типы атомов opls_139 - углерод и opls_140 - водород. Достроим файл: добавим все углы
os.system('wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro')
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('et.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
os.system('grompp -f {0}.mdp -c etane.gro -p et.top -o et_{0}.tpr'.format(m))
os.system('mdrun -deffnm et_{0} -v -nt 1'.format(m))
Построим анимацию (pdb) молекулярной динамики, извлечем потенциальную и кинетическую энергию для шагов симуляции с каждым методом контроля, а также оценим длину C-C связи
for m in methods:
os.system('echo 1 | trjconv -f et_{0}.trr -s et_{0}.tpr -o et_{0}.pdb'.format(m))
os.system('echo 10 11 0 | g_energy -f et_{0}.edr -o et_{0}_en.xvg -xvg none'.format(m))
os.system('g_bond -f et_{0}.trr -s et_{0}.tpr -o bond_{0}.xvg -n b.ndx -xvg none'.format(m))
display(Image(filename='all.png',retina=True))
Выше представлены полученные визуализации конформаций этана при различных методах контроля температуры. Если проранжировать методы от минимума различий между конформациями до максимума, то будет так: an < nh < vr < be < sd.
- в методе an молекула слабо колеблется по длине связи
- в методе be молекула быстро вращается в плоскости связи
- в методе nh молекула вращается вокруг связи, но слабо задерживается в некоторых позициях
- в методе sd молекула быстро вращается вокруг центра масс
- в методе vr молекула поддерживает заторможенную конформацию молекулы, иногда происходит поворот на 120 градусов вокруг связи
Построим графики изменения энергий¶
fig, ax = plt.subplots(nrows=1, ncols=5, figsize=(24,5))
c = 0
for i, name in zip(methods, ['Метод Берендсена', 'Метод Velocity rescale',
'Метод Нуза-Хувера',
'Метод Андерсена', 'Метод стохастической молекулярной динамики']):
etan = np.loadtxt('et_%s_en.xvg' % i)
x = etan[:,0]
pot = etan[:,1]
kin = etan[:,2]
ax[c].plot(x, pot, alpha=0.7, label = 'Potential energy', linewidth=0.5)
ax[c].plot(x, kin, alpha=0.7, label = 'Kinetic energy', linewidth=0.5)
ax[c].legend()
ax[c].title.set_text(name)
c += 1
fig.text(0.5, 0.04, 'Time (ps)', ha='center')
fig.text(0.04, 0.5, 'Energy (kJ/mol)', va='center', rotation='vertical')
plt.show()
Как можно заметить метод Velocity rescale и стохастической молекулярной динамики показывает сходные результаты на всем промежутке времени: в одном и том же диапазоне энергий происходят постоянные колебания. В случае метода Нуза-Хавера видно, что сильным колебаниям подвержена кинетическая энергия, в то время как потенциальная еолеблется с меньшей амплитудой. В методе Берендсена сперва наблюдаются колебания с большой амплитудой, однаков затем она уменьшается. Метод Андерсена характеризуется очень низкими значениями энергий и амплитуд.
Постройте графики распределения длинн связей С-С¶
fig, ax = plt.subplots(nrows=1, ncols=5, figsize=(24,5))
c = 0
for i, name in zip(methods, ['Метод Берендсена', 'Метод Velocity rescale',
'Метод Нуза-Хувера',
'Метод Андерсена', 'Метод стохастической молекулярной динамики']):
bond = np.loadtxt('bond_%s.xvg' % i)
x = bond[:,0]
y = bond[:,1]
ax[c].bar(x, y, width = 0.0005, label = 'length of C-C bond', color = 'purple', linewidth = 0)
ax[c].title.set_text(name)
c += 1
fig.text(0.5, 0.04, 'Bond length (nm)', ha='center')
fig.text(0.04, 0.5, 'Probability density', va='center', rotation='vertical')
plt.show()
У всех методов максимум распределений соответствует длине связи ~0.153. Наименьшая дисперсия у методов Берендсена и Нуза-Хувера и они больше всего похожи на распределение Больцмана, а следовательно лучше всего поддерживают температуру в системе.
Определим время, потраченное на работу GROMACS для различных симуляций¶
for m in methods:
with open("et_%s.log" % m) as f:
print(m, f.readlines()[-4].split()[2])
Итого, самый быстрый - метод Андерсена, а самый длительный - метод стохастической молекулярной динамики. На мой взгляд, по всем параметрам лучше всего подходит метод Берендсена