Taopackage 用の原子番号追記ソフト

Gaussian で ONIOM 法により計算を行う際に Tao Package を利用して input ファイルを作ると非常に楽です。

Tao Package は非常に便利なプログラムなのですが、何箇所か手作業でやらなくてはならない不便な部分がありました。

今回は、prep ファイルの右側に原子番号を書いていく作業を簡略化するためのプログラムをご紹介します。

簡単な背景

Tao Package を使用するためには、下図のように prep ファイルの右端に原子番号を書いていかなければなりませんでした。

数原子程度の prep ファイルならば手動で行っても良いのですが、間違えが起こる可能性もありますし、原子数や化合物数が増えたりすると非常に面倒です。

このようなルーティンワークは、コンピューターにやらせるのが一番です!

    0    0    2

4-hydroxybenzoate anion -2 RESP charges
PHA.res
PHA    INT  0
CORRECT     OMIT DU   BEG                                             
   0.00000                                                      
   1  DUMM  DU    M    0  -1  -2     0.000      .0        .0      .00000
   2  DUMM  DU    M    1   0  -1     1.449      .0        .0      .00000
   3  DUMM  DU    M    2   1   0     1.523   111.21       .0      .00000
   4  O1'   o     M    3   2   1     1.540   111.208  -180.000 -0.862070 8
   5  C1'   c     M    4   3   2     1.269    63.213     0.000  1.017970 6
   6  O2'   o     E    5   4   3     1.269   126.426   180.000 -0.862060 8
   7  C1    cc    M    5   4   3     1.534   116.787    -0.000 -0.457170 6
   8  C2    cc    M    7   5   4     1.409   122.212     0.000  0.034310 6
   9  H1    ha    E    8   7   5     1.090   116.518     0.000  0.031060 1
  10  C3    cd    M    8   7   5     1.393   122.747   180.000 -0.431580 6
  11  H2    ha    E   10   8   7     1.095   120.448   180.000  0.070070 1
  12  C4    c     M   10   8   7     1.447   123.111     0.047  0.543870 6
  13  O4    o     E   12  10   8     1.280   123.646  -179.946 -0.788260 8
  14  C5    cd    M   12  10   8     1.448   112.675    -0.085 -0.431580 6
  15  H3    ha    E   14  12  10     1.095   116.362  -179.960  0.070070 1
  16  C6    cc    M   14  12  10     1.392   123.144     0.085  0.034310 6
  17  H4    ha    E   16  14  12     1.090   120.735   179.953  0.031060 1

LOOP
   C6   C1

IMPROPER
   C1   O1'  C1'  O2'
   C1'  C6   C1   C2
   C1   C3   C2   H1
   C4   C2   C3   H2
   C5   C3   C4   O4
   C4   C6   C5   H3
   C1   C5   C6   H4

DONE
STOP

prep ファイルの構成

このようなテキスト処理は perl などで書いたほうが楽だと思いますが、今回は C++ を選択しました。

今回は、例として以下に示す prep ファイルを使用します。

    0    0    2

4-hydroxybenzoate anion -2 RESP charges
PHA.res
PHA    INT  0
CORRECT     OMIT DU   BEG                                             
   0.00000                                                      
   1  DUMM  DU    M    0  -1  -2     0.000      .0        .0      .00000
   2  DUMM  DU    M    1   0  -1     1.449      .0        .0      .00000
   3  DUMM  DU    M    2   1   0     1.523   111.21       .0      .00000
   4  O1'   o     M    3   2   1     1.540   111.208  -180.000 -0.862070
   5  C1'   c     M    4   3   2     1.269    63.213     0.000  1.017970
   6  O2'   o     E    5   4   3     1.269   126.426   180.000 -0.862060
   7  C1    cc    M    5   4   3     1.534   116.787    -0.000 -0.457170
   8  C2    cc    M    7   5   4     1.409   122.212     0.000  0.034310
   9  H1    ha    E    8   7   5     1.090   116.518     0.000  0.031060
  10  C3    cd    M    8   7   5     1.393   122.747   180.000 -0.431580
  11  H2    ha    E   10   8   7     1.095   120.448   180.000  0.070070
  12  C4    c     M   10   8   7     1.447   123.111     0.047  0.543870
  13  O4    o     E   12  10   8     1.280   123.646  -179.946 -0.788260
  14  C5    cd    M   12  10   8     1.448   112.675    -0.085 -0.431580
  15  H3    ha    E   14  12  10     1.095   116.362  -179.960  0.070070
  16  C6    cc    M   14  12  10     1.392   123.144     0.085  0.034310
  17  H4    ha    E   16  14  12     1.090   120.735   179.953  0.031060


LOOP
   C6   C1

IMPROPER
   C1   O1'  C1'  O2'
   C1'  C6   C1   C2
   C1   C3   C2   H1
   C4   C2   C3   H2
   C5   C3   C4   O4
   C4   C6   C5   H3
   C1   C5   C6   H4

DONE
STOP

この prep ファイルを下のように変換します。

コード

prep ファイルは、

CORRECT OMIT DU BEG

と書かれた行の 2 つ下の行から、次の空行までが座標情報です。そのため、まず CORRECT の行と空行を探します。
次に、座標情報の 2 列目の元素記号を原子番号に変換し、各行の末尾に追加していきます。

実際のコードを下に示します。

/*   Copyright (c) 2017 計算化学.com  All Rights Reserved.
This software contains proprietary and confidential information, including trade secrets, belonging to 計算化学.com.
*/
#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
#include <cstring>
using namespace std;
int element_convert(string atom);    //function to convert elements
int main (int argc, char *argv[]){
if(argc == 1){
cout << "\nERROR -- no argumanet --\nprep_converter [prep_file_name]\n" << endl;
exit(1);
}
ifstream ifs(argv[1]);
string line;
int file_length (0);
for(int i (0); argv[1][i] !='\0'; ++i)file_length +=1;
char filename_without_prep[50];
for(int i = 0; i <= file_length-5; i++) filename_without_prep[i] = argv[1][i];
filename_without_prep[file_length-5] = '\0';
int i(0), k(0), CORRECT_line_number(0), Empty_line_number(0);
stringstream outputfile_name;
outputfile_name << filename_without_prep << "_convert.prep";
ofstream ofs(outputfile_name.str());
while (getline(ifs,line)){
i++;
std::string::size_type CORRECT_Find = line.find("CORRECT     OMIT DU   BEG");
if (CORRECT_line_number != 0 && Empty_line_number !=0){
ofs << line << endl;
}if (CORRECT_line_number != 0 && Empty_line_number ==0){
if (i > CORRECT_line_number+4){
if (line.substr(0,1) == "\0"){
Empty_line_number = i;
}if (line.substr(0,1) != "\0"){
ofs << line << " " << element_convert(line.substr(6,1)) << endl;
}
}if (i <= CORRECT_line_number+4){
ofs << line << endl;
}
}if (CORRECT_line_number == 0){
if(CORRECT_Find != std::string::npos ){
CORRECT_line_number = i;
}
ofs << line << endl;
}
}
ifs.close();
ofs.close();
}
//Function to Convert Elements
int element_convert (string atom){
string elements[128]={"Error","H","He","Li","Be","B","C","N","O","F","Ne","Na","Mg","Al","Si","P","S","Cl","Al","K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb","Sr","Y","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I","Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W","Re","Os","Ir","Pt","Au","Hg","Ti","Pb","Bi","Po","At","Rn","Fr","Ra","Ac","Th","Pa","U","Np","Pu","Am","Cm","Bk","Cf","Es","Md","Bo","Lr","Rf","Db","Sg","Bh","Hs","Mt","Unn","Uuu","Unb","Uut","Uuq","Uup","Uuh","Uus","Uuo"};
int j = 1;
if(atom.substr(1,1)== " "){
for(int i(1); i<=118; i++){
if(strcmp(elements[i].c_str(),(atom.substr(0,1)).c_str())==0){
j = i;break;
}
else continue;}
}else{
for(int i(1); i<=118; i++){
if(strcmp(elements[i].c_str(),(atom.substr(0,2)).c_str())==0){
j = i;break;
}
else continue;}}
return j;
}

上記のコードは仕事の空き時間に書いたため、まだまだダメなところが何箇所もあります。まず、このプログラムは 1 文字表記の元素にしか対応していません。少し修正すれば、すぐに 2 文字表記元素にも対応できますが、、、
2列目の元素名を読み込む際に正規表現を使用するべきだとは思いますが、、、

いずれにせよ、自分が仕事で使う分には十分だと思いました。

管理人は、プログラムの専門家ではありませんので、より良い方法がある場合などはコメントまたはメールで教えていただければ幸いです。
このようなコードを公開するのは、みなさんからのアドバイスを受けたいという気持ちが強いからです。
よろしくお願い致します。

関連する記事

汎関数一覧に戻る

計算手法に戻る

コメントを残す(投稿者名のみ必須)