自動で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 -ooutputname-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!"