GRRM は非常に便利なソフトですが、log ファイルを GaussView で開くことができないという面倒な点があります。
構造最適化計算の log ファイルであれば、末尾の座標情報を切り出してきて、input ファイルを作ることで Gauss View で開くことができます。
しかし、振動計算の log ファイルの内容は複雑なため、視覚化しないことには分子がどの方向に振動しているかが分かりません。また、GaussView 形式のファイルに変換するのも手間がかかります。
本記事では、GRRM による振動計算の Gauss View での解析方法について紹介します。
まずは動画
まずは、実際にプログラムで GRRM での振動計算の log ファイルを変換している様子をご覧ください。今回は、sample.log というファイルを処理して sample_convert.log というファイルを生成しています。
目次
- GRRM と Gaussian の組み合わせ
- Gaussian の振動計算の結果の読み方
- GRRM の振動計算の log ファイルの読み方
- Gauss View で視覚化するには?
- プログラムで変換しよう!
- GRRM での振動計算はやめた方が良い?
- 関連する記事
GRRM と Gaussian の組み合わせ
GRRM は外部の量子化学計算ソフトを利用して計算を進めていきます。計算は GRRM 独自のアルゴリズムで進められますが、各ポイントのエネルギー計算は外部ソフトが行っています。外部の量子化学計算ソフトとして Gaussian が使われる例が多いため、本記事では gaussian と組み合わせた場合について説明します。
Gaussian の振動計算の結果の読み方
GRRM の振動計算の log ファイルを Gauss View で読み取れる形式に変換するためには、まず Gaussian での振動計算の log ファイルの形式を理解しておく必要があります。
振動計算の log ファイル内での各振動の値は以下のように記述されています。”and normal coordinates:” と検索すると該当箇所が見つかります。
振動数順に 1 番から順に並んでいます。 横ならびで 3 つずつ記述されています。
Frequencies – と書かれている行にあるのが振動数です。その下に各原子ごとの振動方向のベクトルが記述されています。例えば 1 番の振動では 、-39.6625 という虚振動があり、48 番目の原子が最も大きく動いていることが分かります(記事を短縮するために省略しました)。
and normal coordinates: 1 2 3 A A A Frequencies -- -39.6625 20.3414 29.8710 Red. masses -- 3.1469 3.7085 3.7319 Frc consts -- 0.0029 0.0009 0.0020 IR Inten -- 0.0709 0.1613 0.5312 Atom AN X Y Z X Y Z X Y Z 1 6 0.03 0.01 0.03 -0.02 -0.04 0.01 -0.01 0.04 0.03 2 6 -0.01 0.02 -0.01 0.04 -0.01 0.00 0.00 0.03 0.03 3 6 -0.03 0.03 0.03 0.02 -0.01 0.05 0.00 0.03 0.02 4 6 0.03 0.05 0.00 0.00 -0.05 0.05 0.03 0.05 0.05 (中略) 72 1 0.01 -0.03 0.03 0.01 -0.04 0.08 0.11 -0.03 -0.03 73 1 0.01 0.00 0.02 0.01 -0.02 0.06 0.09 0.01 0.02 4 5 6 A A A Frequencies -- 38.1194 47.4094 54.7824 Red. masses -- 3.9525 4.1334 4.5046 Frc consts -- 0.0034 0.0055 0.0080 IR Inten -- 0.6458 0.0977 1.2680 Atom AN X Y Z X Y Z X Y Z 1 6 -0.02 -0.01 0.05 0.04 0.01 0.05 -0.04 -0.03 0.04 2 6 0.00 -0.03 0.03 -0.03 -0.03 -0.02 0.08 0.00 0.05 3 6 0.01 -0.04 0.01 -0.03 -0.04 0.01 0.07 0.03 -0.02 4 6 0.01 0.01 0.06 0.06 0.02 0.06 -0.04 0.01 0.00 (中略) 72 1 0.03 0.00 -0.02 -0.13 0.09 -0.13 0.03 0.06 -0.06 73 1 0.00 -0.01 0.06 -0.13 -0.03 -0.08 0.02 0.00 -0.01 7 8 9 A A A Frequencies -- 59.9279 76.5129 84.9341 Red. masses -- 4.4125 4.5072 2.7520 Frc consts -- 0.0093 0.0155 0.0117 IR Inten -- 0.1743 0.7956 1.4194 Atom AN X Y Z X Y Z X Y Z 1 6 0.10 -0.03 -0.05 0.01 0.05 0.06 0.01 -0.04 0.06 2 6 0.01 -0.02 0.06 0.03 -0.06 -0.03 -0.01 -0.01 0.02 3 6 0.00 0.01 0.01 0.03 -0.05 -0.03 0.02 -0.03 -0.01 4 6 0.05 -0.03 -0.08 -0.01 0.00 0.09 -0.01 -0.03 0.03 (中略) 72 1 0.01 0.01 0.03 -0.07 -0.05 0.04 0.06 0.04 -0.05 73 1 -0.02 -0.01 0.13 -0.07 -0.10 0.06 0.05 0.05 -0.01
続いて、Forces (Hartrees/Bohr)です。これは、”***** Axes restored to original set *****” と検索すれば該当箇所が見つかります。この Forces (Hartrees/Bohr) を記述しておかないと Gauss View に認識してもらえません。
------------------------------------------------------------------- Center Atomic Forces (Hartrees/Bohr) Number Number X Y Z ------------------------------------------------------------------- 1 6 -0.000002784 -0.000006369 -0.000040343 2 6 -0.000031389 0.000015686 0.000017181 3 6 -0.000018626 -0.000025828 -0.000023635 4 6 0.000009298 -0.000003849 0.000003295 5 1 0.000020497 0.000015586 -0.000021494 6 1 0.000007990 0.000025326 -0.000030640 (以下略)
最後に振動計算の分子座標です。log ファイルの末尾に一行で書かれています。
\原子1, x 座標, y 座標, z 座標\ 原子2, x 座標, y 座標, z 座標\……
のように書き表されています。下の例で言うと
\Ni,1.53367,-0.36081,3.53557\P,0.70675,-2.3638,4.09998\
という部分に相当します。
Test job not archived. 1\1\G\OPT\R\C\S 5\0\\#P\ q\\t\\1,1\Ni,1.53367,-0.36081,3.53557\P,0.70675,-2.3638,4.09998\P,1.5 2849,0.461006,5.60677\C,2.01907,-0.829371,6.84903\H,2.15359,-0.353347, 7.83133\H,2.99917,-1.22872,6.54504\C,1.06624,-2.99547,5.80698\H,2.0835 8,-3.41301,5.73663\H,0.406098,-3.8503,6.0221\C,1.00734,-1.97312,6.9443 2\H,1.19068,-2.51432,7.88432\H,-0.0086774,-1.55726,7.0375\C,-1.10493,- 2.08083,4.07552\C,-1.62404,-1.47773,2.92211\C,-1.96946,-2.36518,5.1344 4\C,-2.97209,-1.15271,2.83668\H,-0.954728,-1.25987,2.08714\C,-3.32036, -2.0354,5.05131\H,-1.59544,-2.8432,6.03967\C,-3.82227,-1.42525,3.90683 (以下略)
GRRM の振動計算の log ファイルの読み方
GRRM の振動計算の log ファイルは、以下のようになっております。
************************************************************************* 1. Normal-Mode Analysis (FREQ) ************************************************************************* INPUT ORIENTATION Ni 1.533666668955 -0.360810433927 3.535573565612 P 0.706749991041 -2.363795114225 4.099979695742 P 1.528486377844 0.461005543655 5.606772701644 C 2.019068236966 -0.829370990574 6.849033745586 H 2.153587653477 -0.353347206822 7.831327900723 H 2.999169171415 -1.228720553456 6.545035777733 (中略) H 8.413312068503 -2.798341202305 -0.590379942791 Br 7.524312962905 -4.698661715174 -2.975150947622 ENERGY = -3576.910842407565 Spin(**2) = 0.000000000000 GRADIENT VECTOR 0.000000116780 0.000000005567 0.000000208948 -0.000000345253 0.000000030050 -0.000000911231 -0.000000487246 (中略) 0.000001633546 0.000001812082 0.000001855264 0.000000499825 HESSIAN MATRIX 0.089712444 0.016548407 0.189362469 -0.038398294 -0.008020183 0.206213247 -0.023090919 -0.008338985 0.003203989 0.288428870 -0.007608658 -0.041964634 0.000783867 -0.016059065 0.296973917 0.002905532 0.000875003 -0.021728225 0.001186394 0.007930761 -0.018187929 -0.000229555 -0.005309716 0.003391566 -0.004335371 (中略) 0.384017240 0.384051815 0.384946987 0.385384449 0.385652716 0.385951320 0.386083907 0.386138763 0.386632503 0.389838419 0.392462685 Normal Mode Analysis at Stationary Point ========================================================================= FREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQF Geometry (Origin = Center of Mass, Axes = Principal Axis), SYMMETRY = C1 Ni 0.914437218599 -0.376322197005 0.323068072649 P 0.647715915901 -2.598900939955 0.265652225350 P 2.622804426593 -0.405558068221 1.753418010255 C 2.392416293279 -1.707427716100 3.057646283926 (中略) C -6.973320750588 0.946853673939 0.378442563793 H -6.901442988827 1.921067693930 2.305134089713 Br -8.868733752528 0.776660885311 0.089559601163 0 1 2 Freq. : -232.22271622 8.05702727 14.28349012 Red. M : 9.87053291 8.29058976 6.93037416 Ni - x : -0.0271203698 -0.0033829990 -0.0054092453 Ni - y : 0.0499012983 0.0195834624 0.0253760842 Ni - z : -0.0282052927 0.0172022988 -0.0308288192 P - x : 0.0026880537 -0.0178134537 -0.0205830044 P - y : -0.0029656768 0.0184369403 0.0281588772 P - z : 0.0105252827 0.0296062478 -0.0380252693 (中略) Br - x : 0.0000153950 -0.0000007144 0.0000719413 Br - y : 0.0000018091 0.0000000420 0.0000163751 Br - z : 0.0000070027 -0.0000001978 0.0000316002 Thermochemistry at 298.150 K, 1.000 Atm E(el) = -3576.910842407565 ZPVE = 0.591235403944 Enthalpie(0K) = -3576.319607003621 E(tr) = 0.001416284860 E(rot) = 0.001416284860 (中略) S(rot) = 0.000064019245 (Symmetry number = 1) S(vib) = 0.000305387697 G-E(el) = 0.505519025493 Free Energy = -3576.405323382071 FREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQFREQF Normal termination of Global Reaction Route Mapping Programs (以下略)
まずは input ファイルの座標が書かれています。”INPUT ORIENTATION” で検索すると該当箇所が見つかります。座標は、Gaussian での振動計算の log ファイルと違い読みやすい形式で書かれています。
座標の下には GRADIENT VECTOR が書いてあります。これは、Gaussian での log ファイルの Forces (Hartrees/Bohr) に対応するものです。
GRRM では、非常に分かりづらい表記法が取られており、
原子1 の x 方向
原子1 の y 方向
原子1 の z 方向
原子2 の x 方向
原子2 の y 方向
原子2 の z 方向
……
のように縦一列で表記されています。
さらに下に行くと、HESSIAN MATRIXが書かれていますが、この部分は Gauss View で解析する際には不要な部分です。
HESSIAN MATRIX の下には、計算後の座標が書かれています。”Normal Mode Analysis at Stationary Point” で検索すると該当箇所が見つかります。
振動計算は一点計算であるため座標は変わらないのですが、コンピュータが計算しやすい値になるように分子を回転させています。
この座標情報の下に振動に関する計算結果が記されています。”Freq. :” で検索すると該当箇所が見つかります。これも非常に分かりにくく縦一列ずつ表記されています。
Gauss View で視覚化するには?
Gauss View で視覚化するためには、GRRM の log ファイルから必要な部分のみを抜き出し、適切な形式に並び変える必要があります。
大き分けて 3 つの操作が必要です。
- 振動に関する情報を抜き出し、x, y, z 座標を縦一列から横並びの形式へと変更する。
- Forces (Hartrees/Bohr) の情報を抜き出し、x, y, z 座標を縦一列から横並びの形式へと変更する。
- 座標情報を抜き出し、横一列表示へと変更する。
後は、一番最初に #p の行を入れるとか、最後に Normal termination を付け加えるなどの細かい作業があるだけです。
変換後の出力ファイルの末尾に組成式を書き込む必要があるのですが、ここは適当で大丈夫です。
プログラムで変換しよう!
このような操作は、手作業でやると膨大な時間がかかってしまうため、実際には簡単なプログラムを作って実行します。
こちらからダウンロードしてください。
今回のプログラムは、C++ で書いてあります。
g++ EQ_list.cpp
でコンパイルしてください。実際使用する際は、
./a.out sample_EQ_list.log
と入力して実行してください。
GaussView のファイルパスは、各自で書き換えておいて下さい。
プログラムの実行に当たってエラーなど出ましたら、コメント欄で質問してください。そうすれば他の人とも情報共有できます。
// Written by 計算化学.com 管理人 2/1/2017 // This program converts GRRM's Freq log file to Gaussian's format. // #include <iostream> #include <string> #include <fstream> #include <sstream> #include <iomanip> #include <cstring> #include <cstdio> using namespace std; int genso(string atom); //prototype declaration int main(int argc, char *argv[]){ // Input Error if(argc == 1){ cout << "\n *************************\n ERROR -- no argumanet --\n ./a.out [file_name]\n *************************\n\n Copyright (c) 2017\n 計算化学.com All Rights Reserved.\n" << endl; exit(1); } // Reading file name ifstream ifs(argv[1]); string line; int i(0), k(0), Geometry_startline(0), Initial_Structure(0), Gradient_Startline(0), Gradient_Endline(0); long Freq_startline[9999]; while (getline(ifs,line)){ i++; std::string::size_type Geometry_Find = line.find("Geometry"); std::string::size_type Freq_Find = line.find("Freq."); std::string::size_type Initial_Find = line.find("INPUT ORIENTATION"); std::string::size_type Gradient_Find = line.find("GRADIENT VECTOR"); std::string::size_type Hessian_Find = line.find("HESSIAN MATRIX"); if (Initial_Find != std::string::npos ){ Initial_Structure = i; }if (Gradient_Find != std::string::npos ){ Gradient_Startline = i; }if (Hessian_Find != std::string::npos ){ Gradient_Endline = i; }if (Geometry_Find != std::string::npos ){ Geometry_startline = i; }if (Freq_Find != std::string::npos ) { k++; Freq_startline[k] = i; } } int atom_number = Freq_startline[1] - Geometry_startline -3; //writing outputfile char filename_without_com[50]; int file_length (0); for(int i (0); argv[1][i] !='\0'; ++i) file_length +=1; for(int i = 0; i <= file_length-4; i++) filename_without_com[i] = argv[1][i]; filename_without_com[file_length-4] = '\0'; stringstream outputfile_name; outputfile_name << filename_without_com << "_convert.log"; ofstream ofs(outputfile_name.str()); ofs << "----------------------------------------------------------------------\n#P\n----------------------------------------------------------------------\n\n" << " and normal coordinates:" << endl; ifs.close(); ifs.open(argv[1]); i=0, k=1; int l(0); double vibration[99999][3]={}; string element[99999]; while (getline(ifs,line)){ i++; if(i == Freq_startline[k]-1){ int AA(0), BB(0), CC(0); sscanf(line.data(),"%d %d %d", &AA, &BB, &CC); ofs << " " << AA +1 << " " << BB +1 << " " << CC +1 <<endl; ofs << " A A A" << endl; }if(i == Freq_startline[k]){ double Frequencies[3]={}; sscanf(line.substr(10,47).data(),"%lf %lf %lf", &Frequencies[0], &Frequencies[1], &Frequencies[2]); ofs << "Frequencies -- " << Frequencies[0] << " " << Frequencies[1] << " " << Frequencies[2] << endl; ofs << " Atom AN X Y Z X Y Z X Y Z"<<endl; }if(i > Freq_startline[k] + 1 && i < (Freq_startline[k] + 1 + atom_number * 3) ){ l++; double a(0), b(0), c(0); sscanf(line.substr(10,47).data(),"%lf %lf %lf", &a, &b, &c); vibration[l][0] = a; vibration[l][1] = b; vibration[l][2] = c; element[l] = line.substr(0,2); }if(i == (Freq_startline[k] + 1 + atom_number * 3)){ for(int j=0; j<=atom_number-1; j++){ //writing element symbol if(j+1<10) ofs<< " " << j+1; else ofs<< " " << j+1; if (genso(element[j*3+1])<10) ofs << " " << genso(element[j*3+1]) << " "; else ofs << " " << genso(element[j*3+1]) << " "; //1つ目 if(vibration[j*3+1][0]>=0) ofs << " "; ofs << fixed << setprecision(2) << vibration[j*3+1][0] << " "; if(vibration[j*3+2][0]>=0) ofs << " "; ofs << fixed << setprecision(2) << vibration[j*3+2][0] << " "; if(vibration[j*3+3][0]>0) ofs << " "; ofs << fixed << setprecision(2) << vibration[j*3+3][0] << "\t"; //2つ目 if(vibration[j*3+1][1]>=0) ofs << " "; ofs << fixed << setprecision(2) << vibration[j*3+1][1] << " "; if(vibration[j*3+2][1]>=0) ofs << " "; ofs << fixed << setprecision(2) << vibration[j*3+2][1] << " "; if(vibration[j*3+3][1]>0) ofs << " "; ofs << fixed << setprecision(2) << vibration[j*3+3][1] << "\t"; //3つ目 if(vibration[j*3+1][2]>=0) ofs << " "; ofs << fixed << setprecision(2) << vibration[j*3+1][2] << " "; if(vibration[j*3+2][2]>=0) ofs << " "; ofs << fixed << setprecision(2) << vibration[j*3+2][2] << " "; if(vibration[j*3+3][2]>0) ofs << " "; ofs << fixed << setprecision(2) << vibration[j*3+3][2] << endl; } k++; } } ofs << "\n***** Axes restored to original set *****\n" << " -------------------------------------------------------------------\n" << " Center Atomic Forces (Hartrees/Bohr)\n" << " Number Number X Y Z\n" << " -------------------------------------------------------------------\n"; // writing gradient ifs.close(); ifs.open(argv[1]); double gradient[99999]; i=0, k=0; while (getline(ifs,line)){ i++; if(i > Gradient_Startline && i < Gradient_Endline){ k++; double a(0); sscanf(line.data(),"%lf", &a); gradient[k] = a; } } for(int j=0; j<=atom_number-1; j++){ //writing element symbol if(j+1<10) ofs<< " " << j+1; else ofs<< " " << j+1; //atomic number if (genso(element[j*3+1])<10) ofs << " " << genso(element[j*3+1]) << " "; else ofs << " " << genso(element[j*3+1]) << " "; //1つ目 if(gradient[j*3+1]>=0) ofs << " "; ofs << fixed << setprecision(10) << gradient[j*3+1] << " "; if(gradient[j*3+2]>=0) ofs << " "; ofs << fixed << setprecision(10) << gradient[j*3+2] << " "; if(gradient[j*3+3]>0) ofs << " "; ofs << fixed << setprecision(10) << gradient[j*3+3] << endl; } //initial coordinate ifs.close(); ifs.open(argv[1]); i=0; stringstream Coordinates; while (getline(ifs,line)){ i++; if(i == Initial_Structure){ ofs<<" Test job not archived.\n"<<" 1\\1\\G\\OPT\\R\\C\\S\n"<<" 5\\0\\\\#P\\"<<endl; Coordinates<< " q\\\\t\\\\1,1"; }if(i == Initial_Structure+1){ double a(0), b(0), c(0); sscanf(line.substr(3,54).data(),"%lf\t%lf\t%lf",&a,&b,&c); if(line.substr(1,1) == " ") Coordinates<<"\\"<<line.substr(0,1)<<","<<a<<","<<b<<","<<c; else Coordinates<<"\\"<<line.substr(0,2)<<","<<a<<","<<b<<","<<c; }if(i > Initial_Structure+1 && i < Initial_Structure + atom_number){ double a(0), b(0), c(0); sscanf(line.substr(3,54).data(),"%lf\t%lf\t%lf",&a,&b,&c); if(line.substr(1,1) == " ") Coordinates<<"\\"<<line.substr(0,1)<<","<<a<<","<<b<<","<<c; else Coordinates<<"\\"<<line.substr(0,2)<<","<<a<<","<<b<<","<<c; }if(i == Initial_Structure + atom_number){ double a(0), b(0), c(0); sscanf(line.substr(3,54).data(),"%lf\t%lf\t%lf",&a,&b,&c); if(line.substr(1,1) == " ") Coordinates<<"\\"<<line.substr(0,1)<<","<<a<<","<<b<<","<<c; else Coordinates<<"\\"<<line.substr(0,2)<<","<<a<<","<<b<<","<<c; } } Coordinates<<"\\\\Version=Fujitsu-XTC-G16RevB.01\\"; for(int j=0; j<=(Coordinates.str().length()/70) ; j++){ ofs << " " << Coordinates.str().substr(70*j,70)<<endl; } ofs<<"[X(C20H37O2)]\\\\\\@"<<endl<<endl; ofs<<"Normal termination of Gaussian 16"<<endl; stringstream buff; buff << "open /Applications/gv/gview.app " << outputfile_name.str() << " &\n"; //cout << buff.str().c_str() << endl; system(buff.str().c_str()); ofs.close(); } //End of Main Functional //////////////////////////////////////////////////////////////////////////////////////////////////////// // Function converting element name to atomic nuimber int genso (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","Ds","Rg","Cp","Uut","Uuq","Uup","Uuh","Uus","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; }
GRRM での振動計算はやめた方が良い?
GRRM は、裏で gaussian を走らせて計算しています。すなわち、GRRM で計算した結果を Gaussian の log 形式へと変換する作業は二度手間です。そもそも、GRRM と Gaussian のどちらで振動計算をしても結果は変わりません。Gaussian で計算した方が早いため、 GRRM での振動計算はお勧めしません。
【2017.03.28 追記】
GRRM と Gaussian の振動計算の結果は異なります。GRRM は独自に振動計算をしているそうです。しかし、なぜ結果が食い違うのかは謎です。。。どちらかが間違っているということもないと思います。でも、計算速度の観点から Gaussian を使ったほうが良いのではと思います。
【2019.03.14 追記】
管理人はプログラムの専門家ではないので、より良い方法がありましたらコメントまたはメールにて教えてください。よろしくお願い致します。
関連する記事
- AFIR の計算結果を素早く可視化するプログラム
- クリップボード上の分子を GaussView で表示【Python】
- 論文ってメールで送っても良いの?【著作権】
- 論文 PDF のダウンロードを効率化!【Webスクレイピング 第2回】
- Webスクレイピング で文献検索を効率化!第1回【Python】
- 【PTSB】Progdyn の使い方【遷移状態後枝分かれ】
- GRRMのinputファイルの作り方
- GRRMのlogファイル解析法 PART1
- GRRMのlogファイル解析法 PART2
- 読み取りにくい SI の座標を変換するプログラム
- ファイル名と原子数を読み取るプログラム(複数ファイル)
- 【河波保雄】はじめての量子化学計算―基礎と可視化 (I・O BOOKS)
- 分子を回転させよう!