Изучение работы методов контроля температуры в молекулярной динамике¶
Скачаем файл с координатами этана и создадим файл с топологией
> http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
> 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_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 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
first one
[ molecules ]
;Name count
et 1
Загрузим файлы с методами контроля температуры
> 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
Построим входные файлы для mdrun
> grompp -f be.mdp -c et.gro -p et.top -o et_be.tpr
> grompp -f vr.mdp -c et.gro -p et.top -o et_vr.tpr
> grompp -f nh.mdp -c et.gro -p et.top -o et_nh.tpr
> grompp -f an.mdp -c et.gro -p et.top -o et_an.tpr
> grompp -f sd.mdp -c et.gro -p et.top -o et_sd.tpr
Запустим молекулярную динамику для каждого метода
> mdrun -deffnm et_be -v -nt 1
> mdrun -deffnm et_vr -v -nt 1
> mdrun -deffnm et_nh -v -nt 1
> mdrun -deffnm et_an -v -nt 1
> mdrun -deffnm et_sd -v -nt 1
Проанализируем результат визуально: переконвертируем в pdb и посмотрим в pymol
> for i in be vr nh an sd; do trjconv -f et_$i.trr -s et_$i.tpr -o et_$i.pdb; done
Select a group: ETH
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
from pymol import cmd
cmd.reinitialize()
cmd.load('et_sd.pdb')
метод Берендсена
Наблюдается в основном вращение вокруг С-С связи и, в меньшей степени, вокруг перпендикулярной ей оси
метод "Velocity rescale"
Наблюдается вращение вокруг С-С связи
метод Нуза-Хувера
Наблюдается вращение вокруг С-С связи и конформационные изменения
метод Андерсена
Наблюдается колебание длин связей
метод стохастической молекулярной динамики
Наблюдается хаотичное движение молекулы в пространстве
Теперь сравним потенциальную и кинетическую энергии
> for i in be vr nh an sd; do g_energy -f et_$i.edr -o et_$i_en.xvg -xvg none; done
import numpy as np
import matplotlib.pyplot as plt
a= np.loadtxt('et_be_en.xvg')
t=a[:,0]
pot=a[:,1]
kin=a[:,2]
tot=a[:,3]
print('Berendsen')
plt.figure(figsize=[16,4.8])
n=plt.plot(t, pot,c='blue', linewidth=0.5)
n=plt.plot(t, kin,c='red', linewidth=0.5)
n=plt.plot(t, tot,c='black', linewidth=0.5)
Потенциальная и кинетическая энергии примерно одинаковы, и амплитуда колебаний уменьшается со временем от 40 до 5 кДж/моль и затем остается постоянной. При этм полная энергия системы остается постоянной практически на всем отрезке времени
a= np.loadtxt('et_vr_en.xvg')
t=a[:,0]
pot=a[:,1]
kin=a[:,2]
tot=a[:,3]
print('velocity rescale')
plt.figure(figsize=[16,4.8])
n=plt.plot(t, pot,c='blue', linewidth=0.5)
n=plt.plot(t, kin,c='red', linewidth=0.5)
n=plt.plot(t, tot,c='black', linewidth=0.5)
Потенциальная и кинетическая энергии сравнимы, их колебания хаотичны на всем отрезке времени.
a= np.loadtxt('et_nh_en.xvg')
t=a[:,0]
pot=a[:,1]
kin=a[:,2]
tot=a[:,3]
print('nose-hoover')
plt.figure(figsize=[16,4.8])
n=plt.plot(t, tot,c='black', linewidth=0.5)
n=plt.plot(t, pot,c='blue', linewidth=0.5)
n=plt.plot(t, kin,c='red', linewidth=0.5)
Кинетическая энергия несколько выше потенциальной, колебания сначала довольно большие (порядка 200 кДж/моль), затем
резко уменьшаются до 50 кДж/моль.
a= np.loadtxt('et_an_en.xvg')
t=a[:,0]
pot=a[:,1]
kin=a[:,2]
tot=a[:,3]
print('andersen')
plt.figure(figsize=[16,4.8])
n=plt.plot(t, pot,c='blue', linewidth=0.5)
n=plt.plot(t, kin,c='red', linewidth=0.5)
n=plt.plot(t, tot,c='black', linewidth=0.5)
Кинетическая энергия выше потенциальной, их амплитуды довольно малы (1.4 кДж/моль) и постоянны на всем отрезке времени.
Полная энергия не меняется
a= np.loadtxt('et_sd_en.xvg')
t=a[:,0]
pot=a[:,1]
kin=a[:,2]
tot=a[:,3]
print('stochastic')
plt.figure(figsize=[16,4.8])
n=plt.plot(t, pot,c='blue', linewidth=0.5)
n=plt.plot(t, kin,c='red', linewidth=0.5)
n=plt.plot(t, tot,c='black', linewidth=0.5)
Потенциальная и кинетическая энергии сравнимы и демонстрируют хаотичные колебания
Рассмотрим распределение полной энергии системы
a= np.loadtxt('et_be_en.xvg')
tot=a[:,3]
print('Berendsen')
plt.figure(figsize=[16,4])
n=plt.hist(tot,color='green', bins=500)
Распределение энергий симметрично, много выбросов, дисперсия мала
a= np.loadtxt('et_vr_en.xvg')
tot=a[:,3]
plt.figure(figsize=[16,4])
n=plt.hist(tot,color='green', bins=500)
Распределение симметрично, с большой дисперсией
a= np.loadtxt('et_nh_en.xvg')
tot=a[:,3]
print('nose-hoover')
plt.figure(figsize=[16,4])
n=plt.hist(tot,color='green', bins=500)
Асимметричное распределение, похоже на расределение Больцмана, но почему-то имеет два пика
a= np.loadtxt('et_an_en.xvg')
tot=a[:,3]
print('andersen')
plt.figure(figsize=[16,4])
n=plt.hist(tot,color='green', bins=500)
Распределение обладает наименьше дисперсией, также имеется небольшая асимметрия, напоминающая распределение Больцмана
a= np.loadtxt('et_sd_en.xvg')
tot=a[:,3]
print('stochastic')
plt.figure(figsize=[16,4])
n=plt.hist(tot,color='green', bins=500)
Распределение со средней дисперсией, и некоторой асимметрией.
Распределение энергий в системе должно соответствовать распределению Больцмана. Из всех методов больше ему соответствует метод Нуза-Хувера.
Рассчитаем распределение длины С-С связи в каждом случае
> cat b.ndx
[ b ]
1 2
> for i in be vr nh an sd; do g_bond -f et_$i.trr -s et_$i.tpr -o bond_$i.xvg -n b.ndx -xvg none; done
a= np.loadtxt('bond_be.xvg')
length=a[:,0]
dens=a[:,1]
print('Berendsen')
plt.figure(figsize=[16,4])
n=plt.hist(length,weights=dens, bins=len(dens), color='green')
Длина связи колеблется в пределах 0.151-0.156 нм, однако есть отдельные выбросы до 0.023 нм в обе стороны
a= np.loadtxt('bond_vr.xvg')
length=a[:,0]
dens=a[:,1]
plt.figure(figsize=[16,4])
n=plt.hist(length,weights=dens, bins=len(dens), color='green')
Длина связи колеблется от 0.145 до 0.16 нм
a= np.loadtxt('bond_nh.xvg')
length=a[:,0]
dens=a[:,1]
print('nose-hoover')
plt.figure(figsize=[16,4])
n=plt.hist(length,weights=dens, bins=len(dens), color='green')
Длина связи колеблется от 0.145 до 0.16 нм, более четко выражен пик в области 0.153 нм
a= np.loadtxt('bond_an.xvg')
length=a[:,0]
dens=a[:,1]
print('andersen')
plt.figure(figsize=[16,4])
n=plt.hist(length,weights=dens, bins=len(dens), color='green')
Длина связи колеблется в пределах 0.150-0.156 нм
a= np.loadtxt('bond_sd.xvg')
length=a[:,0]
dens=a[:,1]
print('stochastic')
plt.figure(figsize=[16,4])
n=plt.hist(length,weights=dens, bins=len(dens), color='green')
Длина колеблется от 0.145 до 0.16 нм
Рассмотрим время работы всех методов.
> for i in be vr nh an sd; do echo $i; cat et_$i.log | tail -10 | grep Time | awk '{print $3}' ; done
be
4.088
vr
4.139
nh
4.184
an
4.003
sd
4.793
Как можно заметить, самый быстрый метод - Андерсена, а самый медленный - метод стохастической молекулярной динамики
Можно сделать вывод, что метод Нуза-Хувера позволяет наиболее реалистично поддерживать температуру в системе. Однако у других методов есть свои преимущества: например, метод Андерсена самый быстрый, также он позволяет моделировать колебательные движения молекул. Метод стохастической молекулярной динамики позволяет моделировать перемещения молекулы в пространстве, таким образом, его можно использовать для докинга.