Методы контроля температуры в молекулярной динамике
В этом практикуме будем смотреть, как реализован контроль температуры в молекулярной динамике на примере GROMACS. Объект исследования - одна молекула этана.
(1) Подготовим файл координат и файл топологии.
У нас уже есть gro файл с координатами этана. Нужно сделать файл топологии et.top. У нас уже есть "скелет" файла, надо его заполнить:
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
# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
В файле /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp нашли подходящие типа атомов: opls_139 (alkane C) и opls_140 (alkane H). Записали их в файл.
# связи
bonds = m3d.GetBonds()
print "ai aj funct"
for i,b in enumerate(bonds):
print ' ' +str(b.GetBeginAtomIdx()+1) +' ' +str(b.GetEndAtomIdx()+1)+ ' ' +"1"
# углы
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetEndAtomIdx() < b2.GetEndAtomIdx():
print ' ' + str(b1.GetEndAtomIdx() + 1) + ' ' + str(b1.GetBeginAtomIdx() + 1) + ' ' +str(b2.GetEndAtomIdx() + 1)+' '+ "1"
# dihedrals
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
for b3 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetEndAtomIdx() < b2.GetEndAtomIdx():
if b1.GetEndAtomIdx() == b3.GetBeginAtomIdx():
print ' '+str(b3.GetEndAtomIdx()+1)+' ' + str(b1.GetEndAtomIdx() + 1) + ' ' + str(b1.GetBeginAtomIdx() + 1) + ' ' +str(b2.GetEndAtomIdx() + 1)+' '+ "3"
В итоге получили такой файл топологии et.top.
(2) Даны пять файлов с разными параметрами контроля температуры:
- be.mdp - метод Берендсена для контроля температуры.
- vr.mdp - метод "Velocity rescale" для контроля температуры.
- nh.mdp - метод Нуза-Хувера для контроля температуры.
- an.mdp - метод Андерсена для контроля температуры.
- sd.mdp - метод стохастической молекулярной динамики.</ul>
Посмотрим, какой метод работает лучше для этана. Делаем входные файлы для mdrun с помощью grompp, получаем 5 .tpr файлов (по сути это объединение gro, top и mdp). Запускаем mdrun
%%bash
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp
import subprocess
for filename in ['be','vr','nh','an','sd']:
cmd = 'grompp -f %s.mdp -c etane.gro -p et.top -o et_%s.tpr >& log.grompp' %(filename,filename)
subprocess.call(cmd,shell=True)
for filename in ['be','vr','nh','an','sd']:
cmd = 'mdrun -deffnm et_%s.tpr -v -nt 1' %filename
subprocess.call(cmd,shell=True)
Конвертируем файлы в pdb-формат.
for filename in ['be','vr','nh','an','sd']:
cmd = 'echo 0 | trjconv -f et_%s.tpr.trr -s et_%s.tpr -o et_%s.pdb >& log.trjconv' %(filename,filename,filename)
subprocess.call(cmd,shell=True)
Сравним файлы визуально. В et_sd движения самые хаотичные, амплитуда колебаний самая большая, центр молекулы сильно меняется. et_be, et_nh, et_vr, крутятся только вокруг С-С связи и вокруг центра молекулы. У et_an колебаний почти нет.
Теперь сравним колебания в структурах более точно. Посмотрим потенциальную энергию и кинетическую энергию для каждой из 5 систем.
for filename in ['be','vr','nh','an','sd']:
cmd1 = 'echo 10 | g_energy -f et_%s.tpr.edr -o et_%s_pot_en.xvg -xvg none >& log.g_energy' %(filename,filename)
cmd2 = 'echo 11 | g_energy -f et_%s.tpr.edr -o et_%s_kin_en.xvg -xvg none >& log.g_energy' %(filename,filename)
subprocess.call(cmd1,shell=True)
subprocess.call(cmd2,shell=True)
В выходных файлах две колонки: первая - время, вторая - энергия.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
for filename in ['be','vr','nh','an','sd']:
open(filename + '.png', 'w').close() #rewrite
a= np.loadtxt('et_%s_pot_en.xvg' % filename)
b= np.loadtxt('et_%s_kin_en.xvg' % filename)
t1=a[:,0]
y1=a[:,1]
t2=b[:,0]
y2=b[:,1]
plt.plot(t1, y1, "b-", t2, y2, "r-" )
red_patch = mpatches.Patch(color='red', label='Kinetic')
blue_patch = mpatches.Patch(color='blue', label='Potential')
plt.legend(handles=[red_patch,blue_patch])
plt.title('Energies of %s system' % filename)
name = filename + '.png'
plt.savefig(name)
plt.show()
be отличается от остальных тем, что амплитуда колебаний энергии уменьшается со временем. В структуре an значения энергии самые низкие: в среднем на 100 порядков ниже, чем в остальных структурах (это согласуется с незначительным колебанием, которые мы видели в PyMOL). Удивительно, что колебание в sd было визуально намного более активным и хаотичным, чем в остальных структурах, а по графику энергии этого и не скажешь.
Посмотрим на колебание связи С-С. Для этого создаем индекс-файл с одной связью b.ndx с индексами атомов, образующих связь. А затем строим распределения длин.
import subprocess
for filename in ['be','vr','nh','an','sd']:
cmd = 'g_bond -f et_%s.tpr.trr -s et_%s.tpr -o bond_%s.xvg -n b.ndx -xvg none' % (filename,filename,filename)
subprocess.call(cmd,shell=True)
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
for filename in ['be','vr','nh','an','sd']:
a= np.loadtxt('bond_%s.xvg' % filename)
t1=a[:,0]
y1=a[:,1]
width = 0.0001 #эмпирически
plt.bar(t1,y1,width, color="blue")
plt.title('Length of C-C bond in %s system' % filename)
plt.show()
По форме на распределение Больцмана похожи распределения длин С-С связей в структурах be и nh. sd тоже похожа, но пик не так сильно выражен. Средняя длина С-С связи в этане 1.535 ангстрем. В nh и be пики распределений совпадают с этим значением, что хорошо. Если выбирать из двух распределений: hn и be, - то в nh правый хвост длиннее левого, как в распределении Больцмана. Можно сделать вывод, что метод nh позволяет наиболее реалистично поддерживать температуру.
Сравним время работы всех этих методов. Его можно посмотреть в log-файлах mdrun (поле time:).
Algorithm | Real time, s |
|---|---|
be | 4.046 |
vr | 4.166 |
nh | 4.261 |
an | 3.996 |
sd | 4.665 |
Самый быстрый метод - an, самый медленный - sd. Это согласуется с колебаниями молекул, которые мы видели: координаты атомов после an почти не менялись, а в sd был самый большой разброс координат, поэтому и понадобилось больше времени на расчеты этих разных структур.