%matplotlib inline
from os import system
from IPython.display import Image,display
import numpy as np
import pylab as plt
import seaborn as sns
! wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
! mv etane.gro ethane.gro
Напишем топологию
cat ethane.top
Индекс-файл для считывания длины связи
! echo [ cc ] > cc.ndx
! echo 1 2 >> cc.ndx
Произведем расчеты
for mdparam in ('be','vr','nh','an','sd'):
# get MD parameters
system('wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/%s.mdp' % mdparam)
# run MD
system('grompp -f %s.mdp -c ethane.gro -p ethane.top -o ethane_%s.tpr' % (mdparam,mdparam))
system('mdrun -deffnm ethane_%s -nt 1' % mdparam)
# dump trajectory
system('echo 0 | trjconv -f ethane_%s.trr -s ethane_%s.tpr -o ethane_%s.pdb' % (mdparam, mdparam, mdparam))
# calculate energy. U - potential, T - kinetic
system('echo 10 | g_energy -f ethane_%s.edr -o ethane_%s_U.xvg -xvg none' % (mdparam,mdparam))
system('echo 11 | g_energy -f ethane_%s.edr -o ethane_%s_T.xvg -xvg none' % (mdparam,mdparam))
# dump bond length
system('g_bond -f ethane_%s.trr -s ethane_%s.tpr -o ethane_%s_bond.xvg -n cc.ndx -xvg none' % (mdparam,mdparam,mdparam))
Посмотрим на траектории расчетов
Термостат Андерсона
Не происходит никаких вращений. Связь C-C вибрирует с постоянной частотой и амплитудой около 0.04 ангстрем.
Термостат Нозе-Хувера
Происходят врящения метильных групп, наблюдаемых изменений длины связи нет. Целиком как твердое тело молекула не движется.
Velocity Rescale
Происходит вращение метильных групп и небольшой поворот молекулы как целого.
Термостат Берендсена
Происходит вращение метильных групп и поворот молекулы как целого, выраженный лучше, чем в случае Velocity Rescale
Стохастическая динамика
Движения всей системы стохастичны
Rolling-window function for analysis
def rolling(array, window, func):
if (window-1) % 2 == 0:
margin_l = margin_r = (window-1) // 2
else:
margin_l = (window-1) // 2
margin_r = margin_l + 1
tmp_array = np.zeros(len(array) + window - 1)
tmp_array[margin_l:-margin_r] = array
res_array = [func(tmp_array[j:j+window]) for j in range(len(tmp_array)-window+1)]
return res_array
Energy-plot function
def en_plot():
for mdparam in ('be','vr','nh','an','sd'):
plt.figure(figsize=(14,4))
print mdparam.upper(), 'Kinetic', 'Potential'
for var in (('T',121),('U',122)):
X = np.loadtxt('ethane_%s_%s.xvg' % (mdparam, var[0]))
plt.subplot(var[1])
plt.plot(X[:,0],X[:,1],linewidth = 0.5)
plt.plot(X[25:-25,0],rolling(X[:,1], 50, np.mean)[25:-25], linewidth = 2.0)
plt.show()
en_plot()
Лучше всего сходится Берендсен. Нозе-Хувер тоже, если закрыть глаза на выбросы.
Построим распределение длин связи C-C. Предполагаем форму, похожую на распределение Максвелла-Больцмана. По этому параметру неплохо выглядят все варианты, кроме Андерсона, хотя тяжелого правого хвоста не наблюдается ни в одном случае.
for mdparam in ('be','vr','nh','an','sd'):
print 'C-C bond length in MD using %s temperature control' % mdparam.upper()
X = np.loadtxt('ethane_%s_bond.xvg' % mdparam)
plt.bar(X[:,0],X[:,1],X[1,0]-X[0,0])
plt.show()
! grep Performance *log > perf.txt
Рассчитаем быстродействие
with open('perf.txt','r') as perf:
print 'method','\t','ns/day'
for line in perf:
method = line.split(':')[0][7:9].upper()
ns_day = line.split()[3]
print method,'\t',ns_day
Хорошие результаты по сходимости, форме распределения длин связи C-C и быстродействию показал термостат Берендсона