Изучение работы методов контроля температуры в молекулярной динамике

In [1]:
from subprocess import run
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
In [2]:
! export GMX_MAXBACKUP=-1

Скачаем все необходимые файлы:

  • etane.gro - файл с координатами этана
  • name.mdp - 5 конфигурационных файлов с разными алгоритмами контроля температуры
In [3]:
! rm -rf *.mdp etane.gro
! wget -q http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp http://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp

Подготовим файл координат etane.gro и файл топологии et.top.

В поле OPLS тип функционала для торсионных углов не "1". Посмотрим файл /usr/share/gromacs/top/oplsaa.ff/forcefield.itp Там есть ссылка на файл с параметрами для ковалентных взаимодействий ffbonded.itp, где указано, что нужное значение для типа dihedrals = 3. Типы атомов: opls_135 для углерода и opls_140 для водорода согласно /usr/share/gromacs/top/oplsaa.ff/atomtypes.atp

In [4]:
def run_and_check(command):
    print(f'Running the command: {command}')
    result = run(command, shell=True, capture_output=True, text=True)
    if result.returncode != 0:
        print(result.stderr)
t_stats = ['be', 'vr', 'nh', 'sd']

Молекулярно-динамический анализ

In [5]:
for thermo in t_stats:
    # В более поздних версиях Gromacs'а требуется, чтобы rvdw = rcoulomb. Сделаем их равными 0.70.
    run_and_check(f'sed -ie "s/0.74/0.70/" {thermo}.mdp')
    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')
    run_and_check(f'echo 0 | gmx trjconv -f et_{thermo}.trr -s et_{thermo}.tpr -o et_{thermo}.pdb')

Running the command: sed -ie "s/0.74/0.70/" be.mdp
Running the command: gmx grompp -f be.mdp -c et.gro -p et.top -o et_be.tpr
Running the command: gmx mdrun -deffnm et_be -v -nt 1
Running the command: echo 0 | gmx trjconv -f et_be.trr -s et_be.tpr -o et_be.pdb
Running the command: sed -ie "s/0.74/0.70/" vr.mdp
Running the command: gmx grompp -f vr.mdp -c et.gro -p et.top -o et_vr.tpr
Running the command: gmx mdrun -deffnm et_vr -v -nt 1
Running the command: echo 0 | gmx trjconv -f et_vr.trr -s et_vr.tpr -o et_vr.pdb
Running the command: sed -ie "s/0.74/0.70/" nh.mdp
Running the command: gmx grompp -f nh.mdp -c et.gro -p et.top -o et_nh.tpr
Running the command: gmx mdrun -deffnm et_nh -v -nt 1
Running the command: echo 0 | gmx trjconv -f et_nh.trr -s et_nh.tpr -o et_nh.pdb
Running the command: sed -ie "s/0.74/0.70/" sd.mdp
Running the command: gmx grompp -f sd.mdp -c et.gro -p et.top -o et_sd.tpr
Running the command: gmx mdrun -deffnm et_sd -v -nt 1
Running the command: echo 0 | gmx trjconv -f et_sd.trr -s et_sd.tpr -o et_sd.pdb

Анализ потенциальной и кинетической энергии системы

In [6]:
for thermo in t_stats:
    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)

Running the command: echo 8 9 | gmx energy -f et_be.edr -o et_be.xvg
Running the command: echo 8 9 | gmx energy -f et_vr.edr -o et_vr.xvg
Running the command: echo 8 9 | gmx energy -f et_nh.edr -o et_nh.xvg
Running the command: echo 8 9 | gmx energy -f et_sd.edr -o et_sd.xvg

Распределение длины связи С-С за время моделирования

In [7]:
for thermo in t_stats:
    print(f'Thermostat {thermo}')
    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.show()

Thermostat be
Running the command: echo 0 | gmx distance -f et_be.trr -s et_be.tpr -n b.ndx -oh bond_be.xvg -xvg none 

Thermostat vr
Running the command: echo 0 | gmx distance -f et_vr.trr -s et_vr.tpr -n b.ndx -oh bond_vr.xvg -xvg none 

Thermostat nh
Running the command: echo 0 | gmx distance -f et_nh.trr -s et_nh.tpr -n b.ndx -oh bond_nh.xvg -xvg none 

Thermostat sd
Running the command: echo 0 | gmx distance -f et_sd.trr -s et_sd.tpr -n b.ndx -oh bond_sd.xvg -xvg none