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列目の元素名を読み込む際に正規表現を使用するべきだとは思いますが、、、
いずれにせよ、自分が仕事で使う分には十分だと思いました。
管理人は、プログラムの専門家ではありませんので、より良い方法がある場合などはコメントまたはメールで教えていただければ幸いです。
このようなコードを公開するのは、みなさんからのアドバイスを受けたいという気持ちが強いからです。
よろしくお願い致します。
関連する記事
- ファイル名と原子数を読み取るプログラム（複数ファイル）
- GRRM の振動計算の結果を Gauss View で解析する方法
- SCFの収束について
- 構造最適化の閾値は、何を意味しているのか？
- 量子化学計算で a.u. を使う理由 〜Why Atomic Unit?〜
