Директория на ShadBox
/hse/amironov/Task9
Изучение работы методов контроля температуры в молекулярной динамике
In [1]:
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
import matplotlib
from matplotlib import pyplot as plt
%matplotlib inline
In [ ]:
_Bash_
wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
cat /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp > tmp
Какие соображения могут нам помочь при поиске нужных типов атомов:
Этан - второй член гомологического ряда алканов, его формула $H_3C-CH_3$
Ethane can be viewed as two methyl groups joined, that is, a dimer of methyl groups.
Посмотрим на вывод функции cat:
In [6]:
Image('alkanes.png')
Out[6]:
In [ ]:
_Bash_
grep 'methyl' tmp
Смотрим на результат:
In [7]:
Image('grep methyl.png')
Out[7]:
Ничего особо полезного здесь не нашли.
Возьмем тип alkane $CH_3$
Для остальных типов атомов не ясно, относятся ли они к этану.
Посмотрим, как устроен файл .gro по этану:
In [ ]:
# Файл ethane.gro устроен следующим образом:
etane
8
1ETH C1 1 0.577 0.217 0.574
1ETH C2 2 0.680 0.252 0.467
1ETH H1 3 0.478 0.241 0.538
1ETH H2 4 0.597 0.274 0.664
1ETH H3 5 0.583 0.111 0.597
1ETH H4 6 0.676 0.358 0.445
1ETH H5 7 0.660 0.195 0.377
1ETH H6 8 0.780 0.228 0.504
1.50000 1.50000 1.50000
Попробуем угадать типы атомов по
In [2]:
# загрузим модули
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 [3]:
# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
Out[3]:
In [4]:
# и придумаем циклы
## связи
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
print b.GetBeginAtomIdx() , b.GetEndAtomIdx()
## углы
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():
print b1.GetEndAtomIdx() , b1.GetBeginAtomIdx(), b2.GetEndAtomIdx()
## dihedrals
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx():
print b1.GetEndAtomIdx()+1 , b1.GetBeginAtomIdx()+1, b2.GetEndAtomIdx()+1, 1
Не понятно, как нам это дальше помогает...
Используем тип alkane $CH_3$
Контроль температуры
Создаем файл et.top в рабочей директории средствами jupyter notebook -> new -> Text File
In [ ]:
#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_135 1 ETH C1 1 -0.189 12.01
2 opls_135 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
1 2 6 1
1 2 7 1
1 2 8 1
6 2 7 1
6 2 8 1
7 2 8 1
[ dihedrals ]
; ai aj ak al funct
3 1 2 6 1
3 1 2 7 1
3 1 2 8 1
4 1 2 6 1
4 1 2 6 1
4 1 2 7 1
4 1 2 8 1
5 1 2 6 1
5 1 2 7 1
5 1 2 8 1
[ 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
In [ ]:
# Загружаем 5 файлов с разными параметрами для процедуры контроля температуры:
_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
# Строим входные файлы для молекулярно-динамического движка mdrun с помощью grompp:
for i in *.mdp; do
grompp -f $i -c etane.gro -p et.top -o et_$i.tpr
done
In [8]:
Image('error.png')
Out[8]:
In [ ]:
# Удаляем ненужные файлы
_bash_
find . -name '#mdout.mdp*' -delete
Попробуем вообще удалить часть "dihedrals" из файла et.top и перезапустить расчет.
Ура, получили файлы .tpr !!!
In [ ]:
_bash_
# запускаем для каждого .tpr mdrun
for i in *.tpr; do
mdrun -deffnm $i -v -nt 1
done
# конвертируем в .pdb
for i in *.mdp; do
trjconv -f et_$i.trr -s et_$i.tpr -o et_$i.pdb
done