Необходимо получить файл с топологией этана. Для этого нужно взять заготовку файла et.top из данного задания и заменить типы атомов и значения углов связей в соответствии с данными из файла atomtypes.atp на kodomo.
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 os
import subprocess
import matplotlib.pyplot as plt
# создадим этан
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
# и придумаем циклы
## связи
print('Bonds')
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
print(' ',b.GetBeginAtomIdx()+1 ,' ',b.GetEndAtomIdx()+1,' 1')
## углы
print('Angles')
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')
## dihedrals
print('Dihedrals')
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
for b3 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and b1.GetIdx() != b2.GetIdx() and b2.GetEndAtomIdx()==b3.GetBeginAtomIdx() and b3.GetIdx() != b2.GetIdx():
print(' ',b1.GetEndAtomIdx()+1,' ',b1.GetBeginAtomIdx()+1,' ',b2.GetEndAtomIdx()+1,' ',b3.GetEndAtomIdx()+1,' ','3')
В качестве типов атомов были взяты opls_139 (alkane C) и opls_140 (alkane H)
#запишем результаты в файл топологии
with open('et.top', 'w') as f:
f.write('''#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_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 ]
; 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
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 ]
; 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 ]
; список атомов 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''')
Даны 5 конфигурационных файлов с разными алгоритмами контроля температуры:
be.mdp - метод Берендсена для контроля температуры.
vr.mdp - метод "Velocity rescale" для контроля температуры.
nh.mdp - метод Нуза-Хувера для контроля температуры.
an.mdp - метод Андерсена для контроля температуры.
sd.mdp - метод стохастической молекулярной динамики.
Все файлы поместила с рабочую директорию.
for program in ['be','vr','nh','an','sd']:
#построение входных файлов
cmd = 'grompp -f %s.mdp -c ethane.gro -p et.top -o et_%s.tpr >& log.grompp' %(program,program)
subprocess.call(cmd,shell=True)
#запуск молекулярно-динамического движка mdrun
cmd = 'mdrun -deffnm et_%s.tpr -v -nt 1' %program
subprocess.call(cmd,shell=True)
#конвертация в pdb
cmd = 'echo 0 | trjconv -f et_%s.tpr.trr -s et_%s.tpr -o et_%s.pdb >& log.trjconv' %(program,program,program)
subprocess.call(cmd,shell=True)
Посмотрим на итоговые pdb в PyMol
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
from pymol import cmd,stored
cmd.do('''
delete all
load et_an.pdb
mset 1-250
set cache_frames=0
mclear
viewport 300,300
mpng et_an/mov
''')
import imageio
import os
folder = 'et_an'
files = [f"{folder}/{file}" for file in os.listdir(folder)]
images = [imageio.imread(file) for file in files]
imageio.mimsave('et_an.gif', images, fps=24)
from IPython.display import Image
display('Метод Берендсена')
display(Image("et_be.gif", format='png'))
display('Mетод "Velocity rescale"')
display(Image("et_vr.gif", format='png'))
display('метод стохастической молекулярной динамики"')
display(Image("et_sd.gif", format='png'))
display('метод Нуза-Хувера"')
display(Image("et_nh.gif", format='png'))
display('метод Андерсена"')
display(Image("et_an.gif", format='png'))