Изучение работы методов контроля температуры в молекулярной динамике
In [33]:
from os import system
In [8]:
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 [1]:
! wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
In [11]:
less /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp
In [ ]:
# подберем наиболее подходящие типы атомов:
# opls_139 12.01100 ; alkane C
# opls_140 1.00800 ; alkane H.
In [ ]:
# Соберем файл ethane.top:
In [9]:
# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200)
Out[9]:
In [24]:
# получим значения углов...
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetEndAtomIdx() < b2.GetEndAtomIdx():
print (' ' + str(b1.GetEndAtomIdx() + 1) +
' ' + str(b1.GetBeginAtomIdx() + 1) +
' ' + str(b2.GetEndAtomIdx() + 1) +
' ' + str(1))
In [23]:
# и торсионых углов
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
for b3 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetEndAtomIdx() < b2.GetEndAtomIdx():
if b1.GetEndAtomIdx() == b3.GetBeginAtomIdx():
print (' ' + str(b2.GetEndAtomIdx() + 1) +
' ' + str(b1.GetBeginAtomIdx() + 1) +
' ' + str(b1.GetEndAtomIdx() + 1) +
' ' + str(b3.GetEndAtomIdx() + 1) +
' ' + str(3))
In [21]:
top = open('ethane.top', 'w')
#head
top.write('''#include "/usr/share/gromacs/top/oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
et 3
''')
#atoms
top.write('''[ 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
top.write('''[ bonds ]
; ai aj funct b0 kb
''')
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
top.write(' ' + str(b.GetBeginAtomIdx() + 1) + ' ' + str(b.GetEndAtomIdx() + 1) + ' ' + str(1) + '\n')
#angles
top.write('''\n[ 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
top.write('''
[ 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
top.write('''[ 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''')
#done!
top.close()
In [22]:
! cat ethane.top
In [28]:
# создадим индекс-файл с одной с-с связью:
with open('cc.ndx', 'w') as f:
f.write('''[ cc ]
1 2''')
In [32]:
configs = ['be','vr','nh','an','sd']
In [34]:
for config in configs:
# скачаем конфгурацонный файл
system('wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/%s.mdp' % config)
# построим входные файлы для молекулярно-динамического движка
system('grompp -f %s.mdp -c etane.gro -p ethane.top -o ethane_%s.tpr' % (config,config))
# запустим mdrun
system('mdrun -deffnm ethane_%s -v -nt 1' % config)
# проведем конвертацию в pdb
system('echo 0 | trjconv -f ethane_%s.trr -s ethane_%s.tpr -o ethane_%s.pdb' % (config,config,config))
# Сравним потенциальную энергию и кинетическую энергию для каждой из 5 систем
system('echo 10 | g_energy -f ethane_%s.edr -o ethane_%s_pot.xvg -xvg none' % (config,config))
system('echo 11 | g_energy -f ethane_%s.edr -o ethane_%s_kin.xvg -xvg none' % (config,config))
#запустим утилиту по анализу связей g_bond:
system('g_bond -f ethane_%s.trr -s ethane_%s.tpr -o ethane_%s_bond.xvg -n cc.ndx -xvg none' % (config,config,config))
#checkpoint
print '%s done!' % config
In [35]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
In [40]:
for config in ['be','vr','nh','an','sd']:
pot = np.loadtxt('ethane_%s_pot.xvg' % config)
kin = np.loadtxt('ethane_%s_kin.xvg' % config)
t1 = pot[:,0]
y1 = pot[:,1]
t2 = kin[:,0]
y2 = kin[:,1]
plt.plot(t1, y1, "g-", t2, y2, "b-" )
blue_patch = mpatches.Patch(color='blue', label='Kinetic')
green_patch = mpatches.Patch(color='green', label='Potential')
plt.legend(handles=[blue_patch,green_patch])
plt.title('Energies in %s system' % config)
plt.show()
Лучше всего сходится метод Берендсена - колебания начинают затухать уже в первой трети всего процесса. Чуть-чуть похуже ведет себя метод Нуза-Хувера - слишком много выбросов (особенно кинетическо составляющей). Остальные методы не сходятся от слова совсем.
Теперь рассмотрим распределение длинны связи С-С за время моделирования. Для этого построим гистограммы для каждого метода:
In [43]:
for config in ['be','vr','nh','an','sd']:
data = np.loadtxt('ethane_%s_bond.xvg' % config)
x = data[:,0]
y = data[:,1]
width = 0.0001
plt.bar(x,y,width, color="grey")
plt.title(' C-C bond length distribution in %s system' % config)
plt.show()
На распределение Больцмана похожи только методы Берендсена и Нуза-Хувера, что согласуется с предыдущим наблюдением. Максимумы распределений находятся как раз в районе цифры 1.53, что опять-таки говорит в пользу их правдоподобности. Остальные распределения представляют из себя какую-то маловразумительную кашу (как и в случае с энергиями для этих методов).
Теперь сравним методы по быстродействию.
In [57]:
from timeit import default_timer as timer
from subprocess import Popen, PIPE
In [58]:
def count_times(names):
results = []
for config in names:
start = timer()
proc = Popen(
'grompp -f %s.mdp -c etane.gro -p ethane.top -o ethane_%s.tpr' % (config,config),
shell=True, stdout=PIPE, stderr=PIPE)
_ = proc.communicate()
proc = Popen('mdrun -deffnm ethane_%s -v -nt 1' % config, shell=True, stdout=PIPE, stderr=PIPE)
results.append(timer() - start)
return results
In [48]:
import pandas as pd
In [60]:
states = pd.DataFrame(data={'Methods':configs, 'Times':count_times(configs)})
In [61]:
states
Out[61]:
И вновь метрод Берендсена показывает отличный результат - в этот раз наравне с методом "Velocity rescale", который был далеко не самым лучшим при измерении энергий или колебания длины связи. Таким образом, можно с увененностью сказать, что сегодня у нас выиграл метод be.