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

In [1]:
!wget http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
--2021-05-27 22:11:54--  http://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 192.168.180.1
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:80... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro [following]
--2021-05-27 22:11:54--  https://kodomo.fbb.msu.ru/FBB/year_08/term6/etane.gro
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 399
Saving to: `etane.gro'

100%[======================================>] 399         --.-K/s   in 0s      

2021-05-27 22:11:55 (5,39 MB/s) - `etane.gro' saved [399/399]

Файл et.top был модифицирован: добавлена информация про связи и углы, тип функционала для торсионных углов заменен с 1 на 3.

In [3]:
print(open('et.top', 'r').read())
#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      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 ]
; list of atoms 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 конфигурационных файлов с алгоритмами контроля температуры:

  1. be.mdp - метод Берендсена для контроля температуры.
  2. vr.mdp - метод "Velocity rescale" для контроля температуры.
  3. nh.mdp - метод Нуза-Хувера для контроля температуры.
  4. an.mdp - метод Андерсена для контроля температуры.
  5. sd.mdp - метод стохастической молекулярной динамики.
In [14]:
methods = ['be', 'vr', 'nh', 'an', 'sd']
method_names = ['Берендсена', '"Velocity rescale"', 'Нуза-Хувера',
           'Андерсена', 'стохастической молекулярной динамики']
In [4]:
!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
--2021-05-27 22:15:10--  http://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp
Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 192.168.180.1
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:80... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp [following]
--2021-05-27 22:15:10--  https://kodomo.fbb.msu.ru/FBB/year_08/term6/be.mdp
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1356 (1,3K)
Saving to: `be.mdp'

100%[======================================>] 1 356       --.-K/s   in 0s      

2021-05-27 22:15:11 (13,7 MB/s) - `be.mdp' saved [1356/1356]

--2021-05-27 22:15:11--  http://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp
Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 192.168.180.1
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:80... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp [following]
--2021-05-27 22:15:11--  https://kodomo.fbb.msu.ru/FBB/year_08/term6/vr.mdp
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1427 (1,4K)
Saving to: `vr.mdp'

100%[======================================>] 1 427       --.-K/s   in 0s      

2021-05-27 22:15:12 (15,1 MB/s) - `vr.mdp' saved [1427/1427]

--2021-05-27 22:15:12--  http://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp
Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 192.168.180.1
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:80... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp [following]
--2021-05-27 22:15:12--  https://kodomo.fbb.msu.ru/FBB/year_08/term6/nh.mdp
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1429 (1,4K)
Saving to: `nh.mdp'

100%[======================================>] 1 429       --.-K/s   in 0s      

2021-05-27 22:15:13 (14,1 MB/s) - `nh.mdp' saved [1429/1429]

--2021-05-27 22:15:14--  http://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp
Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 192.168.180.1
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:80... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp [following]
--2021-05-27 22:15:14--  https://kodomo.fbb.msu.ru/FBB/year_08/term6/an.mdp
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1426 (1,4K)
Saving to: `an.mdp'

100%[======================================>] 1 426       --.-K/s   in 0s      

2021-05-27 22:15:14 (15,1 MB/s) - `an.mdp' saved [1426/1426]

--2021-05-27 22:15:15--  http://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp
Resolving kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)... 192.168.180.1
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:80... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp [following]
--2021-05-27 22:15:15--  https://kodomo.fbb.msu.ru/FBB/year_08/term6/sd.mdp
Connecting to kodomo.fbb.msu.ru (kodomo.fbb.msu.ru)|192.168.180.1|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1441 (1,4K)
Saving to: `sd.mdp'

100%[======================================>] 1 441       --.-K/s   in 0s      

2021-05-27 22:15:15 (15,2 MB/s) - `sd.mdp' saved [1441/1441]

In [5]:
!ls
an.mdp	be.mdp	et.top	etane.gro  nh.mdp  prac7.ipynb	sd.mdp	vr.mdp

Получим команды для запуска расчетов с различными типами алгоритмов контроля температуры. Из-за интерактивного ввода на некоторых стадиях их пришлось запускать вручную.

In [8]:
for method in methods:
    mdp_file = method + '.mdp'
    tpr_file = method + '.tpr'
    trr_file = method + '.trr'
    pdb_file = method + '.pdb'
    
    print(f"grompp -f {mdp_file} -c etane.gro -p et.top -o {tpr_file}")
    print(f"mdrun -deffnm {method} -v -nt 1")
    print(f"trjconv -f {trr_file} -s {tpr_file} -o {pdb_file}")
    print("")
grompp -f be.mdp -c etane.gro -p et.top -o be.tpr
mdrun -deffnm be -v -nt 1
trjconv -f be.trr -s be.tpr -o be.pdb

grompp -f vr.mdp -c etane.gro -p et.top -o vr.tpr
mdrun -deffnm vr -v -nt 1
trjconv -f vr.trr -s vr.tpr -o vr.pdb

grompp -f nh.mdp -c etane.gro -p et.top -o nh.tpr
mdrun -deffnm nh -v -nt 1
trjconv -f nh.trr -s nh.tpr -o nh.pdb

grompp -f an.mdp -c etane.gro -p et.top -o an.tpr
mdrun -deffnm an -v -nt 1
trjconv -f an.trr -s an.tpr -o an.pdb

grompp -f sd.mdp -c etane.gro -p et.top -o sd.tpr
mdrun -deffnm sd -v -nt 1
trjconv -f sd.trr -s sd.tpr -o sd.pdb

In [9]:
!ls
#mdout.mdp.1#  an.gro  be.cpt  be.tpr	  nh.edr  nh.trr       sd.mdp  vr.gro
#mdout.mdp.2#  an.log  be.edr  be.trr	  nh.gro  prac7.ipynb  sd.pdb  vr.log
#mdout.mdp.3#  an.mdp  be.gro  et.top	  nh.log  sd.cpt       sd.tpr  vr.mdp
#mdout.mdp.4#  an.pdb  be.log  etane.gro  nh.mdp  sd.edr       sd.trr  vr.pdb
an.cpt	       an.tpr  be.mdp  mdout.mdp  nh.pdb  sd.gro       vr.cpt  vr.tpr
an.edr	       an.trr  be.pdb  nh.cpt	  nh.tpr  sd.log       vr.edr  vr.trr
In [10]:
from IPython.display import Image

Получим гифки с молекулярными динамиками по каждому методу.

In [17]:
for method, name in zip(methods, method_names):
    print("Метод", name)
    display(Image(f"{method}.gif", format='png'))
Метод Берендсена
Метод "Velocity rescale"
Метод Нуза-Хувера
Метод Андерсена
Метод стохастической молекулярной динамики

Вывод: Для разных методов наблюдается разная подвижность. Наименьшая для контроля температуры методом Андерсена, наибольшая для стохастической молекулярной динамики.

Расчет потенциальной и кинетической энергии для каждого метода

Получим команды для ввода на kodomo. Необходимо дополнительно будет ввести числа 10 и 11, что соответствует расчету потенциальной и кинетической энергии.

In [19]:
for method in methods:
    print(f"g_energy -f {method}.edr -o {method}.xvg -xvg none")
g_energy -f be.edr -o be.xvg -xvg none
g_energy -f vr.edr -o vr.xvg -xvg none
g_energy -f nh.edr -o nh.xvg -xvg none
g_energy -f an.edr -o an.xvg -xvg none
g_energy -f sd.edr -o sd.xvg -xvg none
In [26]:
import numpy as np
import matplotlib.pyplot as plt
def plot_energy(method, name):
    data_in = np.genfromtxt(method+'.xvg')
    x, pot, kin = data_in[:,0], data_in[:,1], data_in[:,2]
    
    fig, ax = plt.subplots(figsize=(15,5))
    ax.plot(x, pot, label='Потенциальная E')
    ax.plot(x, kin, label='Кинетическая E')
    ax.set_title(f'Метод {name}')
    ax.set_xlabel('Время (пс)')
    ax.set_ylabel('Энергия (кДж/моль)')
    ax.legend()

Изменение кинетической и потенциальной энергии от времени:

In [27]:
for method, name in zip(methods, method_names):
    plot_energy(method, name)

Вывод: Метод Бендерсона - наименьший диапазон колебаний энергии,значения не изменяются. Разброс наблюдается в методе Нуза-Хувера. Есть пики свыше 350 кДж/моль. Наименьшей энергия оказалась в методе Андерсона - не превосходит 1.4 кДж/моль.

Распределение длины C-C связи

Cоздадим индекс файл с одной связью.

In [32]:
ndx = '''[ bonds ]
1 2'''
with open('b.ndx', 'w') as f:
    f.write(ndx)
In [33]:
for method in methods:  
    command = f"g_bond -f {method}.trr -s {method}.tpr -o bond_{method}.xvg -n b.ndx -xvg none"
    !$command
                         :-)  G  R  O  M  A  C  S  (-:

                     Gnomes, ROck Monsters And Chili Sauce

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                                :-)  g_bond  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f         be.trr  Input        Trajectory: xtc trr trj gro g96 pdb cpt
  -n          b.ndx  Input        Index file
  -s         be.tpr  Input, Opt!  Structure+mass(db): tpr tpb tpa gro g96 pdb
  -o    bond_be.xvg  Output       xvgr/xmgr file
  -l      bonds.log  Output, Opt. Log file
  -d   distance.xvg  Output, Opt. xvgr/xmgr file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    19      Set the nicelevel
-b           time   0       First frame (ps) to read from trajectory
-e           time   0       Last frame (ps) to read from trajectory
-dt          time   0       Only use frame when t MOD dt = first time (ps)
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-xvg         enum   none    xvg plot formatting: xmgrace, xmgr or none
-blen        real   -1      Bond length. By default length of first bond
-tol         real   0.1     Half width of distribution as fraction of -blen
-[no]aver    bool   yes     Average bond length distributions
-[no]averdist  bool yes     Average distances (turns on -d)

Warning: file does not end with a newline, last line:
1 2
Group     0 (          bonds) has     2 elements
There is one group in the index
Will gather information on 1 bonds
trn version: GMX_trn_file (single precision)
Reading frame     200 time  200.000   

Total number of samples               : 251
Mean                                  : 0.153135
Standard deviation of the distribution: 0.00177711
Standard deviation of the mean        : 0.00011217

gcq#162: "Confirmed" (Star Trek)

                         :-)  G  R  O  M  A  C  S  (-:

                     Gnomes, ROck Monsters And Chili Sauce

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                                :-)  g_bond  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f         vr.trr  Input        Trajectory: xtc trr trj gro g96 pdb cpt
  -n          b.ndx  Input        Index file
  -s         vr.tpr  Input, Opt!  Structure+mass(db): tpr tpb tpa gro g96 pdb
  -o    bond_vr.xvg  Output       xvgr/xmgr file
  -l      bonds.log  Output, Opt. Log file
  -d   distance.xvg  Output, Opt. xvgr/xmgr file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    19      Set the nicelevel
-b           time   0       First frame (ps) to read from trajectory
-e           time   0       Last frame (ps) to read from trajectory
-dt          time   0       Only use frame when t MOD dt = first time (ps)
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-xvg         enum   none    xvg plot formatting: xmgrace, xmgr or none
-blen        real   -1      Bond length. By default length of first bond
-tol         real   0.1     Half width of distribution as fraction of -blen
-[no]aver    bool   yes     Average bond length distributions
-[no]averdist  bool yes     Average distances (turns on -d)

Warning: file does not end with a newline, last line:
1 2
Group     0 (          bonds) has     2 elements
There is one group in the index
Will gather information on 1 bonds
trn version: GMX_trn_file (single precision)
Reading frame       0 time    0.000   
Back Off! I just backed up distance.xvg to ./#distance.xvg.1#
Reading frame     200 time  200.000   

Total number of samples               : 251
Mean                                  : 0.153297
Standard deviation of the distribution: 0.00371384
Standard deviation of the mean        : 0.000234415

gcq#162: "Confirmed" (Star Trek)

                         :-)  G  R  O  M  A  C  S  (-:

                God Rules Over Mankind, Animals, Cosmos and Such

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                                :-)  g_bond  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f         nh.trr  Input        Trajectory: xtc trr trj gro g96 pdb cpt
  -n          b.ndx  Input        Index file
  -s         nh.tpr  Input, Opt!  Structure+mass(db): tpr tpb tpa gro g96 pdb
  -o    bond_nh.xvg  Output       xvgr/xmgr file
  -l      bonds.log  Output, Opt. Log file
  -d   distance.xvg  Output, Opt. xvgr/xmgr file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    19      Set the nicelevel
-b           time   0       First frame (ps) to read from trajectory
-e           time   0       Last frame (ps) to read from trajectory
-dt          time   0       Only use frame when t MOD dt = first time (ps)
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-xvg         enum   none    xvg plot formatting: xmgrace, xmgr or none
-blen        real   -1      Bond length. By default length of first bond
-tol         real   0.1     Half width of distribution as fraction of -blen
-[no]aver    bool   yes     Average bond length distributions
-[no]averdist  bool yes     Average distances (turns on -d)

Warning: file does not end with a newline, last line:
1 2
Group     0 (          bonds) has     2 elements
There is one group in the index
Will gather information on 1 bonds
trn version: GMX_trn_file (single precision)
Reading frame       0 time    0.000   
Back Off! I just backed up distance.xvg to ./#distance.xvg.2#
Reading frame     200 time  200.000   

Total number of samples               : 251
Mean                                  : 0.153033
Standard deviation of the distribution: 0.00293297
Standard deviation of the mean        : 0.000185128

gcq#277: "Lets Kill the Fucking Band"  (Tom Savini - From Dusk til Dawn) 

                         :-)  G  R  O  M  A  C  S  (-:

                God Rules Over Mankind, Animals, Cosmos and Such

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                                :-)  g_bond  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f         an.trr  Input        Trajectory: xtc trr trj gro g96 pdb cpt
  -n          b.ndx  Input        Index file
  -s         an.tpr  Input, Opt!  Structure+mass(db): tpr tpb tpa gro g96 pdb
  -o    bond_an.xvg  Output       xvgr/xmgr file
  -l      bonds.log  Output, Opt. Log file
  -d   distance.xvg  Output, Opt. xvgr/xmgr file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    19      Set the nicelevel
-b           time   0       First frame (ps) to read from trajectory
-e           time   0       Last frame (ps) to read from trajectory
-dt          time   0       Only use frame when t MOD dt = first time (ps)
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-xvg         enum   none    xvg plot formatting: xmgrace, xmgr or none
-blen        real   -1      Bond length. By default length of first bond
-tol         real   0.1     Half width of distribution as fraction of -blen
-[no]aver    bool   yes     Average bond length distributions
-[no]averdist  bool yes     Average distances (turns on -d)

Warning: file does not end with a newline, last line:
1 2
Group     0 (          bonds) has     2 elements
There is one group in the index
Will gather information on 1 bonds
trn version: GMX_trn_file (single precision)
Reading frame       0 time    0.000   
Back Off! I just backed up distance.xvg to ./#distance.xvg.3#
Reading frame     200 time  200.000   

Total number of samples               : 251
Mean                                  : 0.153024
Standard deviation of the distribution: 0.00143875
Standard deviation of the mean        : 9.08129e-05

gcq#277: "Lets Kill the Fucking Band"  (Tom Savini - From Dusk til Dawn) 

                         :-)  G  R  O  M  A  C  S  (-:

                God Rules Over Mankind, Animals, Cosmos and Such

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                                :-)  g_bond  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f         sd.trr  Input        Trajectory: xtc trr trj gro g96 pdb cpt
  -n          b.ndx  Input        Index file
  -s         sd.tpr  Input, Opt!  Structure+mass(db): tpr tpb tpa gro g96 pdb
  -o    bond_sd.xvg  Output       xvgr/xmgr file
  -l      bonds.log  Output, Opt. Log file
  -d   distance.xvg  Output, Opt. xvgr/xmgr file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    19      Set the nicelevel
-b           time   0       First frame (ps) to read from trajectory
-e           time   0       Last frame (ps) to read from trajectory
-dt          time   0       Only use frame when t MOD dt = first time (ps)
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-xvg         enum   none    xvg plot formatting: xmgrace, xmgr or none
-blen        real   -1      Bond length. By default length of first bond
-tol         real   0.1     Half width of distribution as fraction of -blen
-[no]aver    bool   yes     Average bond length distributions
-[no]averdist  bool yes     Average distances (turns on -d)

Warning: file does not end with a newline, last line:
1 2
Group     0 (          bonds) has     2 elements
There is one group in the index
Will gather information on 1 bonds
trn version: GMX_trn_file (single precision)
Reading frame       0 time    0.000   
Back Off! I just backed up distance.xvg to ./#distance.xvg.4#
Reading frame     200 time  200.000   

Total number of samples               : 251
Mean                                  : 0.153279
Standard deviation of the distribution: 0.0031202
Standard deviation of the mean        : 0.000196945

gcq#277: "Lets Kill the Fucking Band"  (Tom Savini - From Dusk til Dawn) 

In [36]:
def plot_bonds(method, name):
    data_in = np.genfromtxt('bond_'+method+'.xvg')
    x, l = data_in[:,0], data_in[:,1]
    
    fig, ax = plt.subplots(figsize=(15,3))
    ax.bar(x, l, width=0.0005)
    ax.set_title(f'Метод {name}')
    ax.set_xlabel('Длина связи, нм')
In [37]:
for method, name in zip(methods, method_names):
    plot_bonds(method, name)

Вывод: Метода Андерсена - диапазон колебаний длин связи оказался наименьшим. Это согласуется с наименьшей энергией и общей заторможенностью структуры. Также заметно, что распределение для данного метода имеет 2 заметных пика, помимо центрального. У метода Берендсена наблюдается узкий пик для длины связи, что объясняется узким диапазоном энергий. У остальных методов распределения по связям примерно сопоставимы. В среднем же длины связей оказываются сопоставимыми в диапазое 1.52-1.53 ангстрем, что похоже на правду.

Распределение по потенциальной энергии

In [40]:
def plot_energy_hist(method, name):
    data_in = np.genfromtxt(method+'.xvg')
    l = data_in[:,1]
    
    fig, ax = plt.subplots(figsize=(15,5))
    ax.hist(l, bins=250)
    ax.set_title(f'Метод {name}')
    ax.set_xlabel('Потенциальная энергия, кДж/моль')
In [41]:
for method, name in zip(methods, method_names):
    plot_energy_hist(method, name)

Вывод: Форма распределений соответствует распределению Больцмана у всех методов, кроме Андерсена и Берендсена. Лучше всего выглядит распределени в методе Нуза-Хувера, но большинство молекул имеют низкое значение энергии.

Время работы mdrun

Астрономическое время работы алгоритмов по данным из логов:

  1. Метод Берендсена: 4.118
  2. Метод "Velocity rescale": 4.077
  3. Метод Нуза-Хувера: 4.405
  4. Метод Андерсена: 4.056
  5. Метод стохастической молекулярной динамики: 4.913

Вывод: Времена работы отличаются мало. Быстрее всего работает метод Андерсена. Медленне всего - метод стохастической молекулярной динамики. Колебания молекулы проще вычислить при малых изменениях энергии, что характерно для Андерсена.

Лучшие результаты дал метод Нуза-Хувера в общем случае. Для отдельных случаев с особыми условями могут подойти другие методы.