Методы контроля температуры в молекулярной динамике (на примере молекулы этана)
In [10]:
%%bash
# файл координат этана
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
In [7]:
%%bash
cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp
# нужные типы атомов для файла топологии:
# opls_135 12.01100 ; alkane CH3
# opls_140 1.00800 ; alkane H.
In [31]:
%%bash
# подготовленный файл топологии этана
cat ethane.top
In [11]:
%%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 # метод "Velocity rescale"
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 # метод стохастической молекулярной динамики
In [45]:
%%writefile b.ndx
[ b ]
1 2
In [13]:
%%bash
for i in $(ls *.mdp | sed 's/.mdp//')
do
grompp -f $i.mdp -c etane.gro -p ethane.top -o et_$i.tpr
mdrun -deffnm et_$i -v -nt 1
echo 2 | trjconv -f et_$i.trr -s et_$i.tpr -o et_$i.pdb
echo '10 11' | g_energy -f et_$i.edr -o et_$i_en.xvg -xvg none
g_bond -f et_$i.trr -s et_$i.tpr -o bond_$i.xvg -n b.ndx -xvg none
done
In [1]:
# Launch PyMOL
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
from pymol import cmd
pymol.finish_launching()
# Set image
import time
from IPython.display import Image, display
defaultImage = 'pymolimg.png'
def prepareImage(width=500, height=500, sleep=2, filename=defaultImage):
cmd.ray(width, height)
cmd.png(filename)
time.sleep(sleep)
In [2]:
# готовим .png из состояний этана в .pdb файле
import glob
for pdb in glob.glob('*.pdb'):
cmd.load(pdb)
cmd.show_as('sticks')
cmd.mset('1-'+str(cmd.count_states('all')))
cmd.mpng(pdb.rsplit('.')[0]+'/'+pdb.rsplit('.')[0])
cmd.delete('all')
In [3]:
# видео для каждого из 5 методов контроля температуры
# %%bash
# ffmpeg -pattern_type glob -i "*.png" -c:v libx264 -pix_fmt yuv420p -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2" ~/et_*.mp4
In [7]:
from IPython.display import HTML
def show_video(file):
video = open(file, "rb").read().encode("base64")
video_tag = '<video controls src="data:video/mp4;base64, %s">' % video
display(HTML(data=video_tag))
In [11]:
print ('метод Андерсена')
show_video('et_an.mp4')
print ('метод Берендсена ')
show_video('et_be.mp4')
print ('метод Нуза-Хувера')
show_video('et_nh.mp4')
print ('метод "Velocity rescale" ')
show_video('et_vr.mp4')
print ('метод стохастической молекулярной динамики ')
show_video('et_sd.mp4')
метод Андерсена: колебания длины связи без вращения
метод Нуза-Хувера: вращение вокруг С-С связи
методы Берендсена и "Velocity rescale": вращение вокруг С-С связи и в прострастве
метод стохастической молекулярной динамики: наиболее интенсивное вращение молекулы в прострастве
In [8]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import glob
import sys
reload(sys)
sys.setdefaultencoding('utf-8')
Сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем.
In [19]:
dict_method = {'an': 'Метод Андерсена', 'be':'Метод Берендсена', 'nh':'Метод Нуза-Хувера', 'vr':'Velocity rescale',
'sd':'Метод стохастической молекулярной динамики'}
for file in sorted(glob.glob('energy*.xvg'), key=str):
data = np.loadtxt(file)
time = data[:,0]
u = data[:,1] # потеницальная энергия
k = data[:,2] # кинетическая энергия
plt.plot(time, u, color='red', linewidth=0.2)
plt.plot(time, k, color='black', linewidth=0.2)
plt.title('%s' % dict_method[file.split('_')[1].split('.')[0]])
u_lab = mpatches.Patch(color='red', label='Potential energy')
k_lab = mpatches.Patch(color='black', label='Kinetic energy')
plt.legend(handles=[u_lab,k_lab])
plt.show()
Построим распределение длины связи С-С за время моделирования
In [28]:
for file in sorted(glob.glob('bond*.xvg'), key=str):
data = np.loadtxt(file)
bond = data[:,1]
plt.hist(bond, bins=10, color='red')
plt.title('%s' % dict_method[file.split('_')[1].split('.')[0]])
plt.show()
In [1]:
print 'runtime \n an 3.010 \n be 2.970 \n nh 3.140 \n sd 3.540 \n vr 3.030'
Наибольшее сходство с распределением Больцмана имеет метод Берендсена; кроме того, он показал наименьшее время работы из 5 алгоритмов. Вероятно, из рассмотренных методов поддержания температуры в системе он является оптимальным.
In [12]:
%%bash
jupyter nbconvert --to html hw8_filippova.ipynb