【GROMACS】gmx genrestrで原子の位置を拘束する

#はじめに
MD計算では平衡化の過程でタンパク質等の位置を拘束したいときがあります。
この記事ではgmx genrestrを使ってMD計算中の原子の位置を拘束する方法を説明します。

#今回扱う系

![vmdscene (コピー).jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/6abd522c-492e-c7c6-ceef-f843f0b78a02.jpeg)

こちらの系を扱います。

* タンパク質
* リガンド
* 水
* イオン

が含まれています。

topファイルの最後の[molecules]の欄は

```
[ molecules ]
; Compound        #mols
Protein_chain_A     1
ligand              1
SOL             10126
CL                  8
```
となります。

ここでは「Protein_chain_A」と「ligand」を拘束するとして進めます。

#タンパク質だけのgroファイル、リガンドのみのgroファイルを作成

タンパク質（Protein_chain_A）だけのgroファイル、リガンド（ligand）のみのgroファイルを作成する必要があります。

※タンパク質が複数種類ある場合（トポロジーファイル内の記述が別れている場合）にはそれぞれのgroファイルを作成する必要があります。

まずタンパク質だけのgroファイルを作成します。

手動で作成しても良いですが、ここではtrjconvを使いましょう。

全系が含まれるwhole.groがあるとすると、

```bash
$gmx trjconv -s whole.gro -f whole.gro -o protein.gro
                     :-) GROMACS - gmx trjconv, 2021.2 (-:

                            GROMACS is written by:
     Andrey Alekseenko              Emile Apol              Rossen Apostolov     
         Paul Bauer           Herman J.C. Berendsen           Par Bjelkmar       
       Christian Blau           Viacheslav Bolnykh             Kevin Boyd        
     Aldert van Buuren           Rudi van Drunen             Anton Feenstra      
    Gilles Gouaillardet             Alan Gray               Gerrit Groenhof      
       Anca Hamuraru            Vincent Hindriksen          M. Eric Irrgang      
      Aleksei Iupinov           Christoph Junghans             Joe Jordan        
    Dimitrios Karkoulis            Peter Kasson                Jiri Kraus        
      Carsten Kutzner              Per Larsson              Justin A. Lemkul     
       Viveca Lindahl            Magnus Lundborg             Erik Marklund       
        Pascal Merz             Pieter Meulenhoff            Teemu Murtola       
        Szilard Pall               Sander Pronk              Roland Schulz       
       Michael Shirts            Alexey Shvetsov             Alfons Sijbers      
       Peter Tieleman              Jon Vincent              Teemu Virolainen     
     Christian Wennberg            Maarten Wolf              Artem Zhmurov       
                           and the project leaders:
        Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel

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

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

GROMACS:      gmx trjconv, version 2021.2
Executable:   /usr/local/gromacs/2021.2_GPU/bin/gmx
Data prefix:  /usr/local/gromacs/2021.2_GPU
Working dir:  /home/tacoma/work/genrestr
Command line:
  gmx trjconv -s whole.gro -f whole.gro -o protein.gro

Note that major changes are planned in future for trjconv, to improve usability and utility.
Will write gro: Coordinate file in Gromos-87 format
Select group for output
Group     0 (         System) has 33005 elements
Group     1 (        Protein) has  2601 elements
Group     2 (      Protein-H) has  1291 elements
Group     3 (        C-alpha) has   162 elements
Group     4 (       Backbone) has   486 elements
Group     5 (      MainChain) has   649 elements
Group     6 (   MainChain+Cb) has   800 elements
Group     7 (    MainChain+H) has   810 elements
Group     8 (      SideChain) has  1791 elements
Group     9 (    SideChain-H) has   642 elements
Group    10 (    Prot-Masses) has  2601 elements
Group    11 (    non-Protein) has 30404 elements
Group    12 (          Other) has    18 elements
Group    13 (            LIG) has    18 elements
Group    14 (             CL) has     8 elements
Group    15 (          Water) has 30378 elements
Group    16 (            SOL) has 30378 elements
Group    17 (      non-Water) has  2627 elements
Group    18 (            Ion) has     8 elements
Group    19 (            LIG) has    18 elements
Group    20 (             CL) has     8 elements
Group    21 ( Water_and_ions) has 30386 elements
Select a group: 1
Selected 1: 'Protein'
Reading frames from gro file 'Complex', 33005 atoms.
Reading frame       0 time    0.000   
Precision of whole.gro is 0.001 (nm)

Back Off! I just backed up protein.gro to ./#protein.gro.1#
Last frame          0 time    0.000   


GROMACS reminds you: "Good judgement is the result of experience; experience is the result of bad judgement." (Mark Twain)
```

続いてリガンドです。

```bash
$ gmx trjconv -s whole.gro -f whole.gro -o ligand.gro
                     :-) GROMACS - gmx trjconv, 2021.2 (-:

                            GROMACS is written by:
     Andrey Alekseenko              Emile Apol              Rossen Apostolov     
         Paul Bauer           Herman J.C. Berendsen           Par Bjelkmar       
       Christian Blau           Viacheslav Bolnykh             Kevin Boyd        
     Aldert van Buuren           Rudi van Drunen             Anton Feenstra      
    Gilles Gouaillardet             Alan Gray               Gerrit Groenhof      
       Anca Hamuraru            Vincent Hindriksen          M. Eric Irrgang      
      Aleksei Iupinov           Christoph Junghans             Joe Jordan        
    Dimitrios Karkoulis            Peter Kasson                Jiri Kraus        
      Carsten Kutzner              Per Larsson              Justin A. Lemkul     
       Viveca Lindahl            Magnus Lundborg             Erik Marklund       
        Pascal Merz             Pieter Meulenhoff            Teemu Murtola       
        Szilard Pall               Sander Pronk              Roland Schulz       
       Michael Shirts            Alexey Shvetsov             Alfons Sijbers      
       Peter Tieleman              Jon Vincent              Teemu Virolainen     
     Christian Wennberg            Maarten Wolf              Artem Zhmurov       
                           and the project leaders:
        Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel

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

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

GROMACS:      gmx trjconv, version 2021.2
Executable:   /usr/local/gromacs/2021.2_GPU/bin/gmx
Data prefix:  /usr/local/gromacs/2021.2_GPU
Working dir:  /home/tacoma/work/genrestr
Command line:
  gmx trjconv -s whole.gro -f whole.gro -o ligand.gro

Note that major changes are planned in future for trjconv, to improve usability and utility.
Will write gro: Coordinate file in Gromos-87 format
Select group for output
Group     0 (         System) has 33005 elements
Group     1 (        Protein) has  2601 elements
Group     2 (      Protein-H) has  1291 elements
Group     3 (        C-alpha) has   162 elements
Group     4 (       Backbone) has   486 elements
Group     5 (      MainChain) has   649 elements
Group     6 (   MainChain+Cb) has   800 elements
Group     7 (    MainChain+H) has   810 elements
Group     8 (      SideChain) has  1791 elements
Group     9 (    SideChain-H) has   642 elements
Group    10 (    Prot-Masses) has  2601 elements
Group    11 (    non-Protein) has 30404 elements
Group    12 (          Other) has    18 elements
Group    13 (            LIG) has    18 elements
Group    14 (             CL) has     8 elements
Group    15 (          Water) has 30378 elements
Group    16 (            SOL) has 30378 elements
Group    17 (      non-Water) has  2627 elements
Group    18 (            Ion) has     8 elements
Group    19 (            LIG) has    18 elements
Group    20 (             CL) has     8 elements
Group    21 ( Water_and_ions) has 30386 elements
Select a group: 13
Selected 13: 'LIG'
Reading frames from gro file 'Complex', 33005 atoms.
Reading frame       0 time    0.000   
Precision of whole.gro is 0.001 (nm)

Back Off! I just backed up ligand.gro to ./#ligand.gro.1#
Last frame          0 time    0.000   


GROMACS reminds you: "Live for Liposuction" (Robbie Williams)

```

#gmx genrestrで拘束条件が記されたitpファイルを作成する

タンパク質に対して実行してposre_protein.itpを生成します。

```bash
$ gmx genrestr -f protein.gro -o posre_protein.itp
                     :-) GROMACS - gmx genrestr, 2021.2 (-:

                            GROMACS is written by:
     Andrey Alekseenko              Emile Apol              Rossen Apostolov     
         Paul Bauer           Herman J.C. Berendsen           Par Bjelkmar       
       Christian Blau           Viacheslav Bolnykh             Kevin Boyd        
     Aldert van Buuren           Rudi van Drunen             Anton Feenstra      
    Gilles Gouaillardet             Alan Gray               Gerrit Groenhof      
       Anca Hamuraru            Vincent Hindriksen          M. Eric Irrgang      
      Aleksei Iupinov           Christoph Junghans             Joe Jordan        
    Dimitrios Karkoulis            Peter Kasson                Jiri Kraus        
      Carsten Kutzner              Per Larsson              Justin A. Lemkul     
       Viveca Lindahl            Magnus Lundborg             Erik Marklund       
        Pascal Merz             Pieter Meulenhoff            Teemu Murtola       
        Szilard Pall               Sander Pronk              Roland Schulz       
       Michael Shirts            Alexey Shvetsov             Alfons Sijbers      
       Peter Tieleman              Jon Vincent              Teemu Virolainen     
     Christian Wennberg            Maarten Wolf              Artem Zhmurov       
                           and the project leaders:
        Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel

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

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

GROMACS:      gmx genrestr, version 2021.2
Executable:   /usr/local/gromacs/2021.2_GPU/bin/gmx
Data prefix:  /usr/local/gromacs/2021.2_GPU
Working dir:  /home/tacoma/work/genrestr
Command line:
  gmx genrestr -f protein.gro -o posre_protein.itp


Reading structure file
Select group to position restrain
Group     0 (         System) has  2601 elements
Group     1 (        Protein) has  2601 elements
Group     2 (      Protein-H) has  1291 elements
Group     3 (        C-alpha) has   162 elements
Group     4 (       Backbone) has   486 elements
Group     5 (      MainChain) has   649 elements
Group     6 (   MainChain+Cb) has   800 elements
Group     7 (    MainChain+H) has   810 elements
Group     8 (      SideChain) has  1791 elements
Group     9 (    SideChain-H) has   642 elements
Select a group: 0
Selected 0: 'System'

Back Off! I just backed up posre_protein.itp to ./#posre_protein.itp.1#

GROMACS reminds you: "This May Come As a Shock" (F. Black)
```

リガンドに対して実行してposre_ligand.itpを生成します。

```bash
$ gmx genrestr -f ligand.gro -o posre_ligand.itp
                     :-) GROMACS - gmx genrestr, 2021.2 (-:

                            GROMACS is written by:
     Andrey Alekseenko              Emile Apol              Rossen Apostolov     
         Paul Bauer           Herman J.C. Berendsen           Par Bjelkmar       
       Christian Blau           Viacheslav Bolnykh             Kevin Boyd        
     Aldert van Buuren           Rudi van Drunen             Anton Feenstra      
    Gilles Gouaillardet             Alan Gray               Gerrit Groenhof      
       Anca Hamuraru            Vincent Hindriksen          M. Eric Irrgang      
      Aleksei Iupinov           Christoph Junghans             Joe Jordan        
    Dimitrios Karkoulis            Peter Kasson                Jiri Kraus        
      Carsten Kutzner              Per Larsson              Justin A. Lemkul     
       Viveca Lindahl            Magnus Lundborg             Erik Marklund       
        Pascal Merz             Pieter Meulenhoff            Teemu Murtola       
        Szilard Pall               Sander Pronk              Roland Schulz       
       Michael Shirts            Alexey Shvetsov             Alfons Sijbers      
       Peter Tieleman              Jon Vincent              Teemu Virolainen     
     Christian Wennberg            Maarten Wolf              Artem Zhmurov       
                           and the project leaders:
        Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel

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

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

GROMACS:      gmx genrestr, version 2021.2
Executable:   /usr/local/gromacs/2021.2_GPU/bin/gmx
Data prefix:  /usr/local/gromacs/2021.2_GPU
Working dir:  /home/tacoma/work/genrestr
Command line:
  gmx genrestr -f ligand.gro -o posre_ligand.itp


Reading structure file
Select group to position restrain
Group     0 (         System) has    18 elements
Group     1 (          Other) has    18 elements
Group     2 (            LIG) has    18 elements
Select a group: 0
Selected 0: 'System'

Back Off! I just backed up posre_ligand.itp to ./#posre_ligand.itp.1#

GROMACS reminds you: "And It Goes a Little Something Like This" (Tag Team)
```


もし指定したいGroupがないときはmake_ndxで拘束する原子のインデックスファイルを作成して使いましょう。

https://qiita.com/tacoma/items/083ba58aabe56fc66525

タンパク質の全原子のインデックスを作成する

参考[2]では水素原子を除外しているのでGroup 2を選択しても良いかもしれません。

#トポロジーファイルを修正する

タンパク質、リガンドそれぞれ修正します。
たとえば、タンパク質の場合ではdihedralの節の後に

```
#ifdef POSRES
#include "posre_protein.itp"
#endif
```

リガンドに対しても同様に行います。

```
#ifdef POSRES
#include "posre_ligand.itp"
#endif
```

出力されたitpファイルを見ると、

```
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
   2    1       1000       1000       1000
   3    1       1000       1000       1000
   4    1       1000       1000       1000
   5    1       1000       1000       1000
   6    1       1000       1000       1000
   7    1       1000       1000       1000
   8    1       1000       1000       1000
   9    1       1000       1000       1000
  10    1       1000       1000       1000
  11    1       1000       1000       1000
  12    1       1000       1000       1000
  13    1       1000       1000       1000
  14    1       1000       1000       1000
  15    1       1000       1000       1000
  16    1       1000       1000       1000
  17    1       1000       1000       1000
  18    1       1000       1000       1000
```

1列目は原子の番号、2列目は関数系、3~5列目がxyz方向の力の定数です。
つまり、参照構造から各原子がバネで繋がれているようなことを行って拘束するわけです。


#mdpファイルの修正

mdpファイルの頭に

```
define                  = -DPOSRES
```

と記述しましょう。

これにより「POSRES」が指定され、トポロジーファイルに書いた「#ifdef POSRES」から「#endif」までに書いた処理が実行されます。

つまり、posre_protein.itp、posre_ligand.itpがincludeされます。

もしリガンドだけ指定したいときは、

```
#ifdef POSRES_LIG
#include "posre_ligand.itp"
#endif
```

とし、mdpファイルで

```
define                  = -DPOSRES_LIG
```

とすればOKです。


#gromppの実行

gromppを実行してtprファイルを作成します。

```
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr
```

通常と違うのは拘束する参照構造を「-r」オプションで指定していることです。
「-c」オプションと同じファイルで大丈夫です。

あとはmdrunでMD計算を流せばOKです。
やや揺らぎますが、大きくは動かないはずです。

#補足
タンパク質がgroファイルの最初にある場合はインデックスがずれないので、whole.groを処理してしまっても構いません。

#参考

https://manual.gromacs.org/documentation/2019/onlinehelp/gmx-genrestr.html

http://www.mdtutorials.com/gmx/complex/06_equil.html