Практикум 6. Изучение методов контроля температуры в МД¶
import pymol
from pymol import cmd
from IPython.display import Image, Video, display
import numpy as np
from matplotlib import pyplot as plt
Подготовим файлы координат и топологии для молекулы этана:
cat ethane.gro
ethane
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
cat topol.top
#include "/usr/local/gromacs/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.243 12.01
2 opls_135 1 ETH C2 2 -0.243 12.01
3 opls_140 1 ETH H1 3 0.081 1.008
4 opls_140 1 ETH H2 4 0.081 1.008
5 opls_140 1 ETH H3 5 0.081 1.008
6 opls_140 1 ETH H4 6 0.081 1.008
7 opls_140 1 ETH H5 7 0.081 1.008
8 opls_140 1 ETH H6 8 0.081 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 7 1
7 2 8 1
6 2 8 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
ethane
[ molecules ]
;Name count
et 1
Для каждого из предложенных для сравнения методов контроля температуры (Берендсена, масштабирования скоростей, Нуза-Хувера, Андерсена, стохастической динамики) напишем скрипт для запуска молдинамики и анализа результатов:
for temp in ['be', 'vr', 'nh', 'an', 'sd']:
text = '''gmx grompp -f {0}.mdp -c ethane.gro -p topol.top -o ethane_{0}.tpr &&
gmx mdrun -deffnm ethane_{0} -v -nt 1 &&
gmx trjconv -f ethane_{0}.trr -s ethane_{0}.tpr -o ethane_{0}.pdb &&
gmx energy -f ethane_{0}.edr -o ethane_{0}_en.xvg -xvg none &&
gmx distance -f ethane_{0}.trr -s ethane_{0}.tpr -oh bond_{0}.xvg -n b.ndx -xvg none -select 0
'''.format(temp)
with open('run_{}.bash'.format(temp), 'w') as f:
f.write(text)
Визуализируем траектории, построим графики кинетической и потенциальной энергии и распределения длины связи углерод-углерод:
pymol.finish_launching(['pymol', '-xQ'])
def render_movie(temp):
cmd.reinitialize()
cmd.bg_color('white')
cmd.set('stick_radius', 0.35)
cmd.set('stick_h_scale', 1)
cmd.mclear()
cmd.load('ethane_{}.pdb'.format(temp))
cmd.show_as('licorice', 'all')
cmd.mset('1-{}'.format(cmd.count_states()))
cmd.movie.produce('{}.mp4'.format(temp))
def plot_energy(name):
data = np.loadtxt(f'ethane_{name}_en.xvg')
time = data[:, 0]
pot = data[:, 1]
kin = data[:, 2]
total = data[:, 3]
fig, ax = plt.subplots(figsize=(12, 5))
plt.plot(time, pot, time, kin, time, total, alpha=0.8)
plt.xlabel('Time, ps')
plt.ylabel('Energy')
plt.legend(['Potential energy', 'Kinetic energy', 'Total energy'])
plt.show()
def plot_bond(name):
data = np.loadtxt(f'bond_{name}.xvg')
length = data[:, 0]
num = data[:, 1]
plt.bar(length, num, width=0.001)
plt.xlabel('C-C bond length, nm')
plt.xlim([0.11, 0.17])
plt.show()
render_movie('be')
render_movie('vr')
render_movie('nh')
render_movie('an')
render_movie('sd')
Метод Берендсена¶
Из графика видно, как система под действием этого алгоритма постепенно пришла к состоянию почти постоянной кинетической энергии с крайне незначительными флуктуациями. Значение длины С-С связи имеет узкое распределение; отличные от среднего значения недопредставлены. Термостат Берендсена не позволяет сгенерировать корректный канонический ансамбль.
Video('be.mp4')
plot_energy('be')
plot_bond('be')
Метод Нуза-Хувера¶
Данный результат, пожалуй, можно назвать самым нереалистичным. На графике наблюдаются гигантские периодические флуктуации энергии, при этом распределение длины связи асимметрично: состояния с длиной связи больше исходной совсем не были сэмплированы. Вообще термостат Нуза-Хувера способен сэмплировать корректное распределение, но только в системе с большим числом степеней свободы (т. е. для данного случая с одной молекулой он не подходит).
Video('nh.mp4')
plot_energy('nh')
plot_bond('nh')
Метод velocity rescale¶
Video('vr.mp4')
plot_energy('vr')
plot_bond('vr')
Метод Андерсена¶
Video('an.mp4')
plot_energy('an')
plot_bond('an')
Метод стохастической динамики¶
Video('sd.mp4')
plot_energy('sd')
plot_bond('sd')
Результаты симуляции с алгоритмами velocity rescale, Андерсена и стохастической динамики достаточно похожи друг на друга: энергия с примерно постоянной ампилтудой флуктуирует вокруг среднего значения, а распределение длины связи имеет колоколообразную форму.