Изучение работы методов контроля температуры в молекулярной динамике¶
Будем изучать, как реализован контроль температуры в GROMACS, и сравнивать различные модели на примере молекулы этана.
%matplotlib inline
(Некоторые ячейки не запущены в отчетном варианте блокнота, чтобы не было большой выдачи)
%%bash
cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp
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 matplotlib.pyplot as plt
# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
#посмотрим, какие нужно задать углы и связи
# и придумаем циклы
## связи
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
print b.GetBeginAtomIdx()+1 , b.GetEndAtomIdx()+1
## углы
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():
print b1.GetEndAtomIdx()+1 , b1.GetBeginAtomIdx()+1, b2.GetEndAtomIdx()+1
## dihedrals
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
for b3 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b2.GetEndAtomIdx() == b3.GetBeginAtomIdx()and b1.GetIdx() != b2.GetIdx():
print b1.GetEndAtomIdx()+1 , b1.GetBeginAtomIdx()+1,b2.GetBeginAtomIdx()+1, b2.GetEndAtomIdx()+1
t = """#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
7 2 8 1
8 2 6 1
7 2 6 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 topology:
topology.write(t)
et_gro = """etane
8
1ETH C1 1 0.577 0.217 0.574
1ETH C2 2 0.680 0.252 0.467
1ETH H1 3 0.478 0.241 0.538
1ETH H2 4 0.597 0.274 0.664
1ETH H3 5 0.583 0.111 0.597
1ETH H4 6 0.676 0.358 0.445
1ETH H5 7 0.660 0.195 0.377
1ETH H6 8 0.780 0.228 0.504
1.50000 1.50000 1.50000
"""
with open("et.gro","w") as grot:
grot.write(et_gro)
Переходим к моделированию:
1)построим входные файлы для молекулярно-динамического движка mdrun
2)для каждого из них запустим mdrun
3)конвертируем в pdb файлы, которые нужно визуализировать
4)сравниваем потенциальную энергию и кинетическую энергию для каждой из 5 систем
%%bash
export i='be vr nh an sd'
for k in $i; do
grompp -f $k.mdp -c et.gro -p et.top -o et_$k.tpr
mdrun -deffnm et_$k -v -nt 1
echo 0 | trjconv -f et_$k.trr -s et_$k.tpr -o et_$k.pdb
echo 10 11 | g_energy -f et_$k.edr -o et_en_$k.xvg -xvg none;
done
Файлы для визуального анализа:
Andersen.pdb
Berendsen.pdb
Velocity rescale.pdb
Nose-Hoover.pdb
Stochastic.pdb
Визуальный анализ показывает, что метод Андерсена чрезмерно ограничивает движение молекулы, а в методе стохастической молекулярной динамики, напротив, слишком велико отклонени молекулы от первоначального положения.
for i in ['be','vr','nh','an','sd']:
a= np.loadtxt("et_en_"+i+'.xvg')
t=a[:,0]
y1=a[:,1]
y2=a[:,2]
plt.plot(t, y1, 'g-',label='potential')
plt.plot(t, y2,'y-',label='kinetic' )
plt.xlabel('time')
plt.ylabel('energy, kJ/mol')
plt.title('method '+i)
plt.legend()
plt.show()
Графики зависимости энергий от времени.
Метод Берендсена: амплитуда изменения энергии уменьшается со временем до небольшого значения - 3-5 kJ/mol.
Метод "Velocity rescale": изменение энергии в диапазоне 10-50 kJ/mol, амплитуда в среднем не меняется.
Метод Нуза-Хувера: не понятно, каково было начальное значение температуры, разброс очень большой, амплитуда примерно сохраняется.
Метод Андерсена: амплитуда не меняется, ее значение мало по сравнению с другими методами.
Метод стохастической молекулярной динамики: зависимость схожа с графиком для "Velocity rescale", но потенциальная энергия здесь чуть меньше кинетической.
%%bash
export i='be vr nh an sd'
for k in $i; do
g_bond -f et_$k.trr -s et_$k.tpr -o bond_$k.xvg -d distance_$k -n b.ndx -xvg none;
done
from scipy.stats import gaussian_kde
for i in ['be','vr','nh','an','sd']:
a= np.loadtxt("distance_"+i+'.xvg')
y=a[:,1]
plt.hist(y, bins=17, color="green", range=(0.143,0.167))
plt.xlabel("C-C bond's length")
plt.title('method '+i)
plt.show()
Распределение Больцмана показывает распределение молекул по скоростям.
Метод Берендсена дает неплохие результаты: гистограмма отдаленно напоминает распределение Больцмана (можно заметить небольшое смещение пика влево).
Метод "Velocity rescale": наибольшее сходство с распределением Больцмана по сравнению со всеми остальными гистограммами.
Метод Нуза-Хувера и метод стохастической молекулярной динамики: распределение слишком симметрично и скорее напоминает нормальное.
Метод Андерсена: узкое распределение - разброс длины связи маленький.
Метод "Velocity rescale" дает самую реалистичную картину поведения молекулы этана при поддержании температуры. Также не плохи методы Берендсена и Нуза-Хувера. Однако первый накладывает значительные ограничения на динамику молекулы, как видно на графике изменения энергии. Второй также показывает странное изменение энергии от времени - слишком большой разброс (хотя длина связи изменяется не значительно по сравнению с другими). Метод Андерсена плохо подходит для описания естественного поведения молекулы при постоянной температуре: молекула слишком мало двигается. Видео, гистограмма и график для энергий подтверждают это. Метод стохастической молекулярной динамики плохо подходит: молекула слишком подвижна, меняет не только длину связи, торсионные углы, но и положение в пространстве.
Лучшим методом для поддержания температуры является "Velocity rescale". Время работы не сильно различается между методами. Заметно дольше работает только метод стохастической молекулярной динамики.
from IPython.display import HTML, display
import tabulate
table = [["method", "time, s"],["Andersen",3.361],["Berendsen",3.418],["Velocity rescale",3.482],["Nose-Hoover",3.517],["Stochastic",4.442]]
display(HTML(tabulate.tabulate(table, tablefmt='html')))