В этом практикуме сравнивается поведение этана при различных термостатах:
# загрузим модули
import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem import Descriptors
from rdkit import RDConfig
import numpy as np
В качестве типов атомов взяты углерод и водород из CH3 группы алканов (opls_135 = CT = alkane CH3 и opls_140 = HC = alkane H).
Параметры атомов и связей оставлены по умолчанию из поля.
# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
0
## связи
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
print('', b.GetBeginAtomIdx() + 1, b.GetEndAtomIdx() + 1, 1, sep=' ')
1 2 1 1 3 1 1 4 1 1 5 1 2 6 1 2 7 1 2 8 1
# углы
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and\
b1.GetIdx() != b2.GetIdx() and\
b1.GetEndAtomIdx()<b2.GetEndAtomIdx():
print('', b1.GetEndAtomIdx()+1, b1.GetBeginAtomIdx()+1, b2.GetEndAtomIdx()+1, 1, sep=' ')
2 1 3 1 2 1 4 1 2 1 5 1 3 1 4 1 3 1 5 1 4 1 5 1 6 2 7 1 6 2 8 1 7 2 8 1
# dihedrals
for Hi in [3, 4, 5]:
for Hj in [6, 7, 8]:
print('', Hi, 1, 2, Hj, 3, sep=' ')
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
for Hi in [3, 4, 5]:
for Hj in [6, 7, 8]:
print('', Hi, Hj, 1, sep=' ')
3 6 1 3 7 1 3 8 1 4 6 1 4 7 1 4 8 1 5 6 1 5 7 1 5 8 1
!cat et.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_135 1 ETH C1 1 -0.180 12.011
2 opls_135 1 ETH C2 2 -0.180 12.011
3 opls_140 1 ETH H1 3 0.060 1.008
4 opls_140 1 ETH H2 4 0.060 1.008
5 opls_140 1 ETH H3 5 0.060 1.008
6 opls_140 1 ETH H4 6 0.060 1.008
7 opls_140 1 ETH H5 7 0.060 1.008
8 opls_140 1 ETH H6 8 0.060 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
; C2-C1-H
2 1 3 1
2 1 4 1
2 1 5 1
; C1-C2-H
1 2 6 1
1 2 7 1
1 2 8 1
; H-C1-H
3 1 4 1
3 1 5 1
4 1 5 1
; H-C2-H
6 2 7 1
6 2 8 1
7 2 8 1
[ dihedrals ]
; ai aj ak al funct
; H-C-C-H
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
ethane
[ molecules ]
;Name count
et 1
Приводятся команды только для одного из термостатов, для прочих аналогичны.
Для термостата Андерсона было изменено в файле параметров значение integrator с md на md-vv, поскольку метод Андерсена требует интегратора Верле, и nstcalcenergy (число шагов между подсчётами энергии) со 100 на 1, поскольку иначе grompp автоматически подставляет значение nstcalcenergy в nstcomm (число шагов между коррекцией движения центра масс. Если коррекция центра масс осуществляется чаще, чем подсчёт энергии, то этот подсчёт лишается смысла), а nstcomm не равное 1 несовместимо с термостатом Андерсона (в данном варианте на каждом шаге перевыбираются скорости для каких-то частиц, поэтому центр масс нужно корректировать также каждый шаг).
gmx grompp -f an.mdp -c et.gro -p et.top -o et_an.tpr
gmx mdrun -deffnm et_be -v -nt 1
Конвертация в pdb (в интерактивной части выбрала 0 — всю систему):
gmx trjconv -f et_be.trr -s et_be.tpr -o et_be.pdb
gmx energy -f et_be.edr -o et_be_en.xvg -xvg none
gmx distance -f et_be.trr -s et_be.tpr -oh bond_be.xvg -n b.ndx -xvg none
!cat b.ndx
[ b ] 1 2
gmx angle -f et_be.trr -od tors_be.xvg -n t.ndx -xvg none -type dihedral
!cat t.ndx
[ t ] 3 1 2 6
Визуальный анализ в PyMOL показывает следующие преобладающие движения, "осевое" в данном случае = округ оси связи СС:
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', size=13)
colors = plt.cm.tab10.colors
methods = ['be', 'vr', 'an', 'nh', 'sd']
fig, axes = plt.subplots(5, 1, figsize=(20, 3*5),
sharex=False, sharey=True
)
for s, ax, c in zip(methods, axes, colors):
a = np.loadtxt(f'et_{s}_en.xvg')
t = a[:,0]
e_pot = a[:,1]
e_kin = a[:,2]
ax.plot(t, e_pot, c=c)
ax.set_title(s)
ax.set_xlabel('t')
ax.set_ylabel('Potential en.')
fig.suptitle('Потенциальная энергия', size=25)
plt.tight_layout()
fig, axes = plt.subplots(5, 1, figsize=(20, 3*5), sharex=False, sharey=True)
for s, ax, c in zip(methods, axes, colors):
a = np.loadtxt(f'et_{s}_en.xvg')
t = a[:,0]
e_pot = a[:,1]
e_kin = a[:,2]
ax.plot(t, e_kin, c=c)
ax.set_title(s)
ax.set_xlabel('t')
ax.set_ylabel('Kinetic en.')
fig.suptitle('Кинетическая энергия', size=25)
plt.tight_layout()
И для потенциальной, и для кинетической энергии наблюдается сходая картина: в термостате Берендсена отклонений спустя небольшое время практически нет, в Velocity rescale и стохастической колебания схожей амплитуды, в Андерсоне также схожая амплитуда, но горазда большая частота изменений (ведь скорости меняются на каждом шаге, а не раз в 100). В методе Нуза-Хувера же наблюдаются периодичные сильные отклонения.
fig, axes = plt.subplots(3, 2, figsize=(20, 3*5), sharex=False, sharey=True)
for s, ax, c in zip(methods, axes.flatten(), colors):
a = np.loadtxt(f'bond_{s}.xvg')
l = a[:,0]
n = a[:,1]
ax.bar(l, n, width=0.001, color=c)
ax.set_title(s)
ax.set_xlim(0.1, 0.2)
ax.set_xlabel('CC length')
ax.set_ylabel('count')
plt.tight_layout()
plt.savefig('len.png')
Энергия пропорциональна квадрату длины (E = kb * (l - b0)^2), подставив это в формулу распределения Больцмана p ~ exp(-E/kT), получим зависимость p ~ exp(-kb * (l - b0)^2)/kT), т.е. C1 * exp( (l - l0)^2 / C2 ), т.е. форму гауссиана. Распределение длин должно иметь форму колокола со средним в оптимальной длине связи (примерно 0.153).
В данном случае примерно одинаковое колоколообразное распределение мы получаем для Velocity rescale, Андерсена и стохастического. Берендсен тоже околоколокол, но гораздо более узкий — допускает меньше вариации длины, чем должно быть, ровно эту проблему решает vr, добавляя стохастический член. В методе Нуза-Хувера распределение смещено вправо — много маленьких значений, а большие довольно резко обрубаются на 0.154 — оптимальном значении с учётом ширины бина в 0.001. Таким образом методе Нуза-Хувера почему-то не семплирует значения ниже исходных.
fig, axes = plt.subplots(3, 2, figsize=(20, 3*5), sharex=False, sharey=True)
for s, ax, c in zip(methods, axes.flatten(), colors):
a = np.loadtxt(f'tors_{s}.xvg')
t = a[:,0]
n = a[:,1]
ax.bar(t, n, width=1, color=c)
ax.set_title(s)
ax.set_xlim(-180, +180)
ax.set_ylim(0, 0.1)
ax.set_xlabel('Torsion')
ax.set_ylabel('Count')
plt.tight_layout()
plt.savefig('tors.png')
Соответсвует выводам по визуальному анализу.
Три заторможенные конформации: +60 (исходная), +180 (= -180) и -60.
В методе Нуза-Хувера торсионный угол практически не меняется, в методе Берендсена меняется в рамках одной заторможенной конформации, для прочих методов семплирует разные конформации.
| Method | Performance (ns/day) |
|---|---|
| be | 11034.211 |
| vr | 10689.147 |
| an | 4300.348 |
| nh | 11327.913 |
| sd | 9973.649 |
Самые быстрые — be и nh, чуть отстаёт vr, затем sd и в два раза медленнее прочих — an.
Метод Нуза-Хувера ведёт себя странно — вся кинетическая энергия воплощается в колебании по связям без вращательной, имеет периодичные выбросы энергии, некорректное распределение длин связи и торсиона. Пусть он при этом и является самым быстрым, его достоверность (если при запуске не было ошибок) очень сомнительна.
Метод Бередсена лишь ненамного отстаёт по скорости. Он удерживает энергию в узком диапазоне, не воспроизодит корректное распределение длин (завышает долю близких к оптимальному значений) и торсинов (застревает в одном заторможенном состоянии и не семплирует остальные). Пусть в нём молекула и вращается, в основном вращение реализуется только вокруг оси молекулы, что выглядит странно.
Т.е. самые быстрые методы — самые неточные.
Методы Velocity rescale, Андерсона и стохастической МД дают схожие картинки по энергиям и распределению длин связи и торсионов. При этом Андерсон более чем в два раза медленнее прочих методов, зато в торсионах все три заторможенные конформации у него разнозначны, в то время как в остальных двух конфомация на -60 недопредставлена относительно +60 и 180. Таким образом Velocity rescale и стохастическая МД — "средние" по скорости и точности, а Андерсон самый точный и самый медленный.