import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem import Descriptors
from rdkit import RDConfig
import numpy as np
Для понимания.
be — Берендсена. Oтклонение от нужной температуры исправляется корректировкой скоростей экспоненциально по времени.
vr — Velocity rescale. То же, но стохастичнее.
an — Андерсена. Скорости перевыбираются из Больцмана раз в несколько шагов.
nh — Нуза-Хувера. Добавляет в гамильтониан трение.
sd — стохастической МД. Добавляет к гамильтонианн трение и шум
from subprocess import run
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
В качестве типов атомов - opls_135 = CT = alkane CH3 и opls_140 = HC = alkane H).
mol=Chem.MolFromSmiles('CC')
AllChem.Compute2DCoords(mol)
m3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
0
# связи
bonds = m3d.GetBonds()
for i,b in enumerate(bonds):
print('', b.GetBeginAtomIdx() + 1, b.GetEndAtomIdx() + 1, 1, sep=' ')
1 2 1 1 3 1 1 4 1 1 5 1 2 6 1 2 7 1 2 8 1
# углы
for b1 in m3d.GetBonds():
for b2 in m3d.GetBonds():
if b1.GetBeginAtomIdx() == b2.GetBeginAtomIdx() and\
b1.GetIdx() != b2.GetIdx() and\
b1.GetEndAtomIdx()<b2.GetEndAtomIdx():
print('', b1.GetEndAtomIdx()+1, b1.GetBeginAtomIdx()+1, b2.GetEndAtomIdx()+1, 1, sep=' ')
2 1 3 1 2 1 4 1 2 1 5 1 3 1 4 1 3 1 5 1 4 1 5 1 6 2 7 1 6 2 8 1 7 2 8 1
# dihedrals
for Hi in [3, 4, 5]:
for Hj in [6, 7, 8]:
print('', Hi, 1, 2, Hj, 3, sep=' ')
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
#Записываем в файл топологии
!cat et.top
#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.180 12.011 2 opls_135 1 ETH C2 2 -0.180 12.011 3 opls_140 1 ETH H1 3 0.060 1.008 4 opls_140 1 ETH H2 4 0.060 1.008 5 opls_140 1 ETH H3 5 0.060 1.008 6 opls_140 1 ETH H4 6 0.060 1.008 7 opls_140 1 ETH H5 7 0.060 1.008 8 opls_140 1 ETH H6 8 0.060 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 ; C2-C1-H 2 1 3 1 2 1 4 1 2 1 5 1 ; C1-C2-H 1 2 6 1 1 2 7 1 1 2 8 1 ; H-C1-H 3 1 4 1 3 1 5 1 4 1 5 1 ; H-C2-H 6 2 7 1 6 2 8 1 7 2 8 1 [ dihedrals ] ; ai aj ak al funct ; H-C-C-H 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 ethane [ molecules ] ;Name count et 1
#сделаем функцию для запуска
def run_and_check(command):
result = run(command, shell=True, capture_output=True, text=True)
if result.returncode != 0:
print(result.stderr)
thermostats = ['be', 'vr', 'an', 'nh', 'sd']
integrator must be md-vv, not md for Andersen nstcomm -> nstcalcenergy must be 1, not 100 for Andersen
gmx grompp -f %s.mdp -c et.gro -p et.top -o et_%s.tpr gmx mdrun -deffnm et_%s -v -nt 1 gmx trjconv -f et_%s.trr -s et_%s.tpr -o et_%s.pdb gmx energy -f et_%s.edr -o et_%s_en.xvg -xvg none
thermostats = ['be']
thermostats = ['be', 'vr', 'an', 'nh', 'sd']
import time
times = []
for thermo in thermostats:
start = time.perf_counter()
run_and_check(f'gmx grompp -f {thermo}.mdp -c et.gro -p et.top -o et_{thermo}.tpr')
run_and_check(f'gmx mdrun -deffnm et_{thermo} -v -nt 1')
stop = time.perf_counter()
times.append(stop-start)
print(f'Runned **{thermo}** in {stop-start} seconds')
## run_and_check(f'gmx trjconv -f et_{thermo}.trr -s et_{thermo}.tpr -o et_{thermo}.pdb')
Runned **be** in 2.388442035764456 seconds Runned **vr** in 2.300229623913765 seconds Runned **an** in 5.434225313365459 seconds Runned **nh** in 2.32022038474679 seconds Runned **sd** in 2.420413237065077 seconds
an очень медленный, остальные примерно наравне с небольшим выигрышем vr.
#Energy visualisation
for thermo in thermostats:
run_and_check(f'echo 8 9 | gmx energy -f et_{thermo}.edr -o et_{thermo}.xvg')
df = pd.read_table(f'et_{thermo}.xvg', skiprows=25, delim_whitespace=True, names=['Time, ps','Potential E, kJ/mol','Kinetic E, kJ/mol'])
df.plot(x='Time, ps', title=thermo, alpha = 0.7)
Для обеих энергий результат похож: в be отклонения минимизируются со временем, в vr, sd и an случайные колебания около равновесного значения с похожей амплитудой, но в an резкие изменения чаще из-за пересчетов. В методе nh появляются частые выбросы.
for thermo in thermostats:
run_and_check(f'echo 0 | gmx distance -f et_{thermo}.trr -s et_{thermo}.tpr -n b.ndx -oh bond_{thermo}.xvg -xvg none ')
df_bar = pd.read_fwf(f'bond_{thermo}.xvg', names=['Bond Length, nm', 'Density'])
plt.bar(df_bar['Bond Length, nm'], df_bar['Density'], width=0.005)
plt.title(f'Thermostat {thermo}')
plt.xlabel('Bond Length, nm')
plt.ylabel('Count')
plt.show()
#Не длинами одними...
for thermo in thermostats:
run_and_check(f'echo 0 | gmx angle -f et_{thermo}.trr -od tors_{thermo}.xvg -n a.ndx -xvg none -type angle')
df_bar = pd.read_fwf(f'tors_{thermo}.xvg', names=['Angle', 'Density'])
plt.bar(df_bar['Angle'], df_bar['Density'], width=3)
plt.title(f'Thermostat {thermo}')
plt.xlabel('Angle')
plt.ylabel('Count')
plt.show()
for thermo in thermostats:
run_and_check(f'echo 0 | gmx angle -f et_{thermo}.trr -od tors_{thermo}.xvg -n tors.ndx -xvg none -type dihedral')
df_bar = pd.read_fwf(f'tors_{thermo}.xvg', names=['Angle', 'Density'])
plt.bar(df_bar['Angle'], df_bar['Density'], width=3)
plt.title(f'Thermostat {thermo}')
plt.xlabel('Angle')
plt.ylabel('Count')
plt.show()
ВЫВОДЫ
Метод Нуза-Хувера - отсутствуют вращения, много выбросов, почти равномерное распределение СС расстояний. Неправдоподобно, зато быстро))
Метод Бередсена лишь ненамного отстаёт по скорости. Он удерживает энергию в узком диапазоне, но распределение длинн крайне несимметрично, нет торсионного вращения. Но тоже шустрый малый))
Методы Velocity rescale, Андерсона и стохастической МД уже довольно похожи. Неплохие распределения энергий и расстояний. Андерсон гораздо медленнее остальных, но умеет в адекватное распределение конформаций.