自動でAMBER形式のトポロジーファイルからGROMACS形式のトポロジーに変換するときのスクリプト

使用条件

macOSの場合、Homebrewをインストールした状態でbrew install gnu-sed gawkを入力してコマンドgsedとgawkコマンドが使える状態にしておく。
Linuxの場合はデフォルトのsedとawkがそのまま使えるはずなので特に問題なし。

やり方

下に載せている3つのスクリプトを同じディレクトリに用意する
AmberToolsのtleapでAMBER形式のトポロジーファイルleap.parm7,初期座標ファイルleap.rst7を用意する。
groconvert.sh内のオプション-i, -oに適切な文字列を入力して実行する。必要に応じて-rもつける（通常つけるべし）。
だいたいの場合はこれでGROMACSのminimizationを実行する準備が整ったことになる。

スクリプト

acpype.py

~~公式からダウンロードする。~~今はオリジナルがなくなったみたいなので、そのクローンであるhttps://github.com/t-/acpype/blob/master/acpype.py からダウンロードする。

gen_posre.pl

GROMACSの座標ファイルをインプットとして、リガンドも含めた水・イオン以外の重原子に対してGROMACSのNVT, NPT平衡化用のposition restraintファイルを作り出してくれるスクリプト。使い方はperl ~/apps/gen_posre.pl hogehoge.groといった感じ。

#!/usr/bin/perl -w

use strict;
use warnings;

my $fx=4.184;   # in kJ nm^-2
my $fy=4.184;
my $fz=4.184;
my @sca=(1000,500,200,100,50,20,10,0);  # 1000 -> 10 kcal A^-2
#my @sca=(10000);       # 1000 -> 10 kcal A^-2

foreach my $s (@sca) {
  open(IN,"$ARGV[0]");
  open(OUT,">posre$s.itp");
  printf(OUT "[ position_restraints ]\n");
  printf(OUT "; atom  type      fx      fy      fz\n");

  my $title=<IN>;
  my $natom=<IN>;
  chomp($natom);
  for(my $i=0;$i<$natom;$i++) {
    $_=<IN>;
    my @data;
    $data[0]=substr($_,0,5);  # residue number (5 positions, integer)
    $data[1]=substr($_,5,5);  # residue name (5 characters)
    $data[2]=substr($_,10,5); # atom name (5 characters)
    $data[3]=substr($_,15,5); # atom number (5 positions, integer)
    for (my $j=0;$j<3;$j++) {
      $data[$j]=trim($data[$j]);
    }
    my $resnum=$data[0];
    my $resnam=$data[1];
    if($resnam ne "Na+" && $resnam ne "Cl-" && $resnam ne "WAT") {
      if($data[2] =~ /^[^H]/) {
        printf(OUT "%6d%6d%6d%6d%6d\n",$data[3],1,$fx*$s,$fy*$s,$fz*$s);
      }
    }
  }
}
close(IN);
close(OUT);

sub trim {
  my $val = shift;
  $val =~ s/^ *(.*?) *$/$1/;
  return $val;
}

groconvert.sh

AmberTools 18のtleapで生成されたleap.parm7, leap.rst7ファイルをインプットとして、上記acpype.py とgen_posre.plをこのgroconvert.shと同じディレクトリに置いておく。（私は~/apps/scriptsディレクトリにおいて呼び出しています。）

使用例は/path/to/groconvert.sh -i leap -o outputname-rみたいな感じ。-rをつけるとgen_posre.plを起動させる。
処理の流れは、

acpype.pyでAMBER用.parm7と.rst7ファイルをGROMACS用に変換する
その中の陽イオン・陰イオンの記述がIP, IM表記だった場合、それらをNa+, Cl-にそれぞれ置き換える
acpype.pyによる変換後、GROMACS用.topファイルの各原子の部分電荷の小数点第6位になぜか1がつくことがあるので、小数点第6位を削って0に揃える
-rをつけた場合、GROMACS用.topファイルの適切な位置にposition restraintの作動用記述を付け加え、さらにgen_posre.plでposition restraintの.itpファイルを作成してくれる。

#!/bin/bash

PROGNAME=$(basename $0)
VERSION="1.0"

usage() {
    echo "Usage: $PROGNAME [OPTIONS] FILE"
    echo "  This script is ~."
    echo
    echo "Options:"
    echo "  -h, --help"
    echo "      --version"
    echo "  -i, --input ARG"
    echo "  -o, --output [ARG]"
    echo "  -c, --long-c"
    echo
    exit 1
}

for OPT in "$@"
do
    case "$OPT" in
        '-h'|'--help' )
            usage
            exit 1
            ;;
        '--version' )
            echo $VERSION
            exit 1
            ;;
        '-i'|'--input' )
            if [[ -z "$2" ]] || [[ "$2" =~ ^-+ ]]; then
                echo "$PROGNAME: option requires an argument -- $1" 1>&2
                exit 1
            fi
            input="$2"
            shift 2
            ;;
        '-o'|'--output' )
            if [[ -z "$2" ]] || [[ "$2" =~ ^-+ ]]; then
                echo "$PROGNAME: option requires an argument -- $1" 1>&2
                exit 1
            fi
            output="$2"
            shift 2
            ;;
        '-r'|'--restraint' )
        R_flag=1
            shift 1
            ;;
        '--'|'-' )
            shift 1
            param+=( "$@" )
            break
            ;;
        -*)
            echo "$PROGNAME: illegal option -- '$(echo $1 | sed 's/^-*//')'" 1>&2
            exit 1
            ;;
        *)
            if [[ ! -z "$1" ]] && [[ ! "$1" =~ ^-+ ]]; then
                #param=( ${param[@]} "$1" )
                param+=( "$1" )
                shift 1
            fi
            ;;
    esac
done

#gsed, gawkコマンドが存在する場合はそちらを優先する（Mac OSとLinuxとの互換性のため）
#Mac OSで使う場合には先にgawk, gsed(gnu-sed)をHomebrewか何かで入れておいてください。
if type gsed 2>/dev/null 1>/dev/null;then
  SED=gsed
else
  SED=sed
fi
if type gawk 2>/dev/null 1>/dev/null;then
  AWK=gawk
else
  AWK=awk
fi

# if [ -z $param ]; then
#     echo "$PROGNAME: too few arguments" 1>&2
#     echo "Try '$PROGNAME --help' for more information." 1>&2
#     exit 1
# fi

##################
python ~/apps/acpype.py -p ${input}.parm7 -x ${input}.rst7
rm md.mdp em.mdp
if [ ! -e ${input}_GMX.gro ];then
    echo "No ${input}_GMX.gro file. Exit."
    exit 1
fi
mv ${input}_GMX.top ${output}.top
mv ${input}_GMX.gro ${output}.gro
${SED} -i -e "s/  IP  /  Na+ /g" ${output}.top
${SED} -i -e "s/  IM  /  Cl- /g" ${output}.top
${SED} -i -e "s/NA+/Na+/g" ${output}.top
${SED} -i -e "s/NA+/Na+/g" ${output}.gro
${SED} -i -e "s/CL-/Cl-/g" ${output}.top
${SED} -i -e "s/CL-/Cl-/g" ${output}.gro

##insert position restraint files into the top file##
if [ $R_flag ]; then
    echo "insert position restraint files into the top file..."
    #search the second [ moleculetype ] block, or [ system ] block."
    moleculetypearray=()
    moleculetypearray+=(`grep "\[ moleculetype \]" ${output}.top -n | sed -e "s/:\[ moleculetype \]//g"`)
    if [ ${#moleculetypearray[@]} -gt 1 ]; then
      i=${moleculetypearray[1]}
    elif [ ${#moleculetypearray[@]} -eq 1 ]; then
      i=`grep "\[ system \]" ${output}.top -n | sed -e "s/:\[ system \]//g"`
    else
      echo "Not found [ moleculetype ] in ${output}.top"
      exit 1
    fi
    j=`expr ${i} - 1`
    ${SED} -i -e "${j}a ; Include Position restraint file\n#ifdef POSRES1000\n#include \"posre1000.itp\"\n#endif\n#ifdef POSRES500\n#include \"posre500.itp\"\n#endif\n#ifdef POSRES200\n#include \"posre200.itp\"\n#endif\n#ifdef POSRES100\n#include \"posre100.itp\"\n#endif\n#ifdef POSRES50\n#include \"posre50.itp\"\n#endif\n#ifdef POSRES20\n#include \"posre20.itp\"\n#endif\n#ifdef POSRES10\n#include \"posre10.itp\"\n#endif\n#ifdef POSRES0\n#include \"posre0.itp\"\n#endif\n" ${output}.top
    ##end of insertion##
    echo "Done!"

    echo "making position restraints file..."
    perl ~/apps/gen_posre.pl ${output}.gro
    echo "Done!"
fi

##round charges##
echo "Round off to 5 decimal place..."
num=`grep "\[ atoms \]" ${output}.top -n | sed -e 's/:.*//g' | head -1`
bondnum=`grep "\[ bonds \]" ${output}.top -n | sed -e 's/:.*//g' | head -1`
lines=$(( bondnum - num + 1 ))
${AWK} -v num=${num} -v bondnum=${bondnum} '{if (NR > num && NR < bondnum && $1 !~ ";" && $7 ~ /(^[0-9\.]+$|^-[0-9\.]+$)/){ chg=sprintf("%.5f",$7) ; gsub($7, chg);print } else{ print }}' ${output}.top > ${output}.top.out
mv ${output}.top.out ${output}.top
echo "Done!"