In [ ]:
Задание 7. Изучение работы методов контроля температуры в молекулярной динамике
В этом задании на примере молекулы этана нужно было изучить, как работает контроль температуры в молекулярной динамике на примере GROMACS.
In [4]:
%%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
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
In [5]:
# загрузим модули
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
In [ ]:
Файл с топологией этана et.top был создан вручную из заготовки в задании.
In [6]:
import subprocess
for method in ['be','vr','nh','an','sd']:
cmd = 'grompp -f %s.mdp -c etane.gro -p et.top -o et_%s.tpr >& log.grompp' %(method, method)
subprocess.call(cmd, shell=True)
In [7]:
for filename in ['be','vr','nh','an','sd']:
cmd = 'mdrun -deffnm et_%s.tpr -v -nt 1' %filename
subprocess.call(cmd,shell=True)
In [9]:
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)
In [15]:
import nglview
import ipywidgets
w1 = nglview.show_structure_file('et_be.pdb')
w1
In [ ]:
метод Берендсена
In [16]:
from IPython.display import Image
Image(url='et_be.gif')
Out[16]:
In [ ]:
метод Андерсена
In [19]:
from IPython.display import Image
Image(url='et_an.gif')
Out[19]:
In [ ]:
метод Нуза-Хувера
In [18]:
from IPython.display import Image
Image(url='et_nh.gif')
Out[18]:
In [ ]:
Метод стохастической молекулярной динамики
In [21]:
При визуальном анализе движения молекулы этана видно, что при методе Берендсена колебания заключаются в довольно быстрых вращениях вокруг С-С связи вместе с небольшими изменениями длины С-С связи, метод Андерсена
в минимальных колебаниях только длины С-С связи, метод Нуза-Хувера похож на Берендсена, при методе стохастической мол динамикитепловые колебания быстрые и хаотичные, Velocity rescale похож на nh и be.
Теперь сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем с помощью g_energy.
In [20]:
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)
In [22]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
for method in ["be","vr","nh","an","sd"]:
a_pot = np.loadtxt('et_%s_pot_en.xvg' % method)
t_pot = a_pot[:,0]
y_pot = a_pot[:,1]
a_kin = np.loadtxt('et_%s_kin_en.xvg' % method)
t_kin = a_kin[:,0]
y_kin = a_kin[:,1]
plt.plot(t_pot, y_pot, "bo", t_kin, y_kin, "ro", alpha=0.8)
red_patch = mpatches.Patch(color='red', label='Kinetic energy')
blue_patch = mpatches.Patch(color='blue', label='Potential energy')
plt.legend(handles=[red_patch,blue_patch])
plt.title('Energies (%s method)' % method)
plt.show()
In [ ]:
Из рисунков видно, что кинетическая и потенциальная энергия примерно совпадают, хотя кинетическая в an и sd немного выше потенциальной, а для nh довольно сильно выше. Ну, оно и видно на гифке, скачет молекула там ну очень быстро и хаотично.
In [ ]: