gaussian scan の結果を解析するプログラム

遷移状態構造の探索などを行う際に gaussian の modredundant の scan コマンドを使うことが多いと多います。
通常は、log ファイルを gauss view で開いてエネルギーの山を探し、gauss view 上でその構造を save し、次の計算のインプットファイル作っていくと思います。

しかし、scan コマンドでは複数の構造を探索するために log ファイルが結構な大きさになってしまうことがあります。

手元にある計算機から log ファイルをダウンロードする際は転送はあっという間です。しかし、出張先などで wi-fi 環境がよくない際には数百 MB のファイルのダウンロードには時間がかかってしまうため、やる気が失せてしまいます。。。

ダウンロードなしで、terminal 上で計算の進行状況を確認し、かつ目的の構造を抜き出すプログラムがあれば良いのになぁということで、簡単なプログラムを書いてみました。

環境

Mac
gcc 7.3.0
Gauss View 6
Gaussian16

Linux (Ubuntu server 16.04)
gcc 5.4.0
Gauss View 6
Gaussian16

プログラムは C++ で書きました。

必要な機能

  1. terminal 上で エネルギーダイヤグラムを表示
  2. local maxima または local minima の座標を保存

local maxima/minima の座標は com 形式のファイルで保存。ファイルの大きさは数 KB なので、ダウンロードするにしても一瞬で終わります。

実装

以下のように実装しました。

#include <iostream>
#include <fstream>

using namespace std;

int main(int argc, char *argv[]){
    
    if(argc == 1){
        cout << "\nERROR -- no argumanet --\nirc [file_name]\n" << endl;
        exit(1);
    }
    
    ifstream ifs(argv[1]);
    string str;
    long l;
    int k(0);
    double energy[100];
    double value;
    
    while(getline(ifs,str)){
        l++;
        std::string::size_type SCF = str.find("SCF Done");
        std::string::size_type OPT = str.find("Optimized Parameters");
        if (SCF != std::string::npos){
            sscanf(str.substr(24,17).data(),"%lf",&value);
            energy[k] = value;
        }if (OPT != std::string::npos){
            k++;
        }
    }
    
    double min_energy = energy[0];
    double max_energy = energy[0];
    int min_num(0), max_num(0);
    
    for(int i(0); i<k;i++){
        if(min_energy > energy[i]){
            min_energy = energy[i];
            min_num = i;
        }if(max_energy < energy[i]){
            max_energy = energy[i];
            max_num = i;
        }
    }
    
    double step = (energy[max_num] - energy[min_num])*627.51 / 19;
    
    int **graph = new int *[20];
    for(int i(0); i<20; i++) graph[i] = new int[k];
    
    // preparing the output
    for(int i(0);i<20;i++){
        for(int j(0);j<k;j++){
            if((energy[j]-energy[min_num])*627.51 >= step * i && (energy[j]-energy[min_num])*627.51 < step * (i+1)){
                graph[i][j] = 1;
                //printf("i = %d, j = %d\n",i,j);
            }else graph[i][j] = 0;
        }
    }
    
    //Output
    if(k>10){
        for(int j(1);j<10;j++) printf("  %d ",j);
        for(int j(10);j<=k;j++) printf(" %d ",j);
    }else{
        for(int j(1);j<=k;j++) printf("  %d ",j);
    }
    
    cout << endl;
    
    for(int i(19);i>=0;i--){
        for(int j(0);j<k;j++){
            if(graph[i][j] == 1) printf("-**-");
            else printf("----");
            
            if(j == k-1) cout << endl;
        }
    }
    
    //free memory
    for(int i(0); i<20; i++){
        delete[] graph[i];
        graph[i] = 0;
    }
    delete[] graph;
    graph = 0;
    
    if(k>10){
        for(int j(1);j<10;j++) printf("  %d ",j);
        for(int j(10);j<=k;j++) printf(" %d ",j);
    }else{
        for(int j(1);j<=k;j++) printf("  %d ",j);
    }
    cout << "\n\n";
    
    printf("min(%d)\t= %10lf\nmax(%d)\t= %10lf\n",min_num+1,min_energy,max_num+1,max_energy);

    cout << "save 1.local maxima structure, 2.local minima structure?\n Type 1 or 2 (Type 0 to escape)\n";
    int option(0);
    cin >> option;
    
    int start_line(0), end_line(0);
    int a, b, c;
    double x, y, z;
    
    if(option == 1){
        
        ofstream ofs("max_structure.com");
        ofs << "%mem=8GB\n%nprocshared=8\n#p B3LYP/6-31+G** opt pop=none freq=noraman\n\ntitle\n\n0 1\n";
        l=0; k=0;
        ifs.clear();
        ifs.seekg(0, ios_base::beg);
        while(getline(ifs,str)){
            l++;
            std::string::size_type Number = str.find(" Number     Number       Type");
            std::string::size_type OPT = str.find("Optimized Parameters");
            std::string::size_type Rotation = str.find("Rotational constants");
            
            if(Number != std::string::npos)start_line = l+1;
            if(Rotation != std::string::npos)end_line = l-1;
            if(OPT != std::string::npos){
                k++;
                if(k == max_num+1) break;
            }
        }
        
        l=0;
        ifs.clear();
        ifs.seekg(0, ios_base::beg);
        while(getline(ifs,str)){
            l++;
            if(l > start_line && l < end_line){
                a = 0;
                sscanf(str.data(),"%d %d %d %lf %lf %lf",&a, &b, &c, &x, &y, &z);
                if (a != 0) ofs << b << "\t" << str.substr(35,35) << endl;
            }
        }
        ofs << "\n\n\n\n";
        ofs.close();
    }if(option == 2){
        ofstream ofs("min_structure.com");
        ofs << "%mem=8GB\n%nprocshared=8\n#p B3LYP/6-31+G** opt=(calcfc,ts,noeig,maxcyc=500,tight) iop(1/8=10) pop=none freq=noraman\n\ntitle\n\n0 1\n";
        l=0; k=0;
        ifs.clear();
        ifs.seekg(0, ios_base::beg);
        while(getline(ifs,str)){
            l++;
            std::string::size_type Number = str.find(" Number     Number       Type");
            std::string::size_type OPT = str.find("Optimized Parameters");
            std::string::size_type Rotation = str.find("Rotational constants");
            
            if(Number != std::string::npos)start_line = l+1;
            if(Rotation != std::string::npos)end_line = l-1;
            if(OPT != std::string::npos){
                k++;
                if(k == min_num+1) break;
            }
        }
        
        l=0;
        ifs.clear();
        ifs.seekg(0, ios_base::beg);
        while(getline(ifs,str)){
            l++;
            if(l > start_line && l < end_line){
                a = 0;
                sscanf(str.data(),"%d %d %d %lf %lf %lf",&a, &b, &c, &x, &y, &z);
                if (a != 0) ofs << b << "\t" << str.substr(35,35) << endl;
            }
        }
        
        ofs << "\n\n\n\n";
        ofs.close();
    }if(option == 0){
        cout<<"\n\nCopyright (c) 2018\n計算化学.com  All Rights Reserved.\n\n";
        exit(0);
    }
    
    ifs.close();
    cout<<"\n\nCopyright (c) 2018\n計算化学.com  All Rights Reserved.\n\n";
    
}

実行方法と実行結果

まずは、プログラムをコンパイルします。

g++ modredundant_analysis.cpp -omodredundant_analysis

これを script ディレクトリに置き、path を通して置きます。

modredundant_analysis test.log

のようにして、実行します。実行しやすいようにもっと短い名前にした方が良いと思いますが、ここれは説明のために、長い名前にしています。

以下のように terminal 上に結果が表示されます。

計算化学.com[~/]$modredundant_analysis test.log 
  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
-----------------**-----------------------------------------------------------------
-------------**------**--**---------------------------------------------------------
---------**------------------**-----------------------------------------------------
---------------------------------**-------------------------------------------------
-----**-----------------------------------------------------------------------------
-------------------------------------**---------------------------------------------
-**---------------------------------------------------------------------------------
-----------------------------------------**-----------------------------------------
------------------------------------------------------------------------------------
---------------------------------------------**-------------------------------------
------------------------------------------------------------------------------------
-------------------------------------------------**---------------------------------
------------------------------------------------------------------------------------
------------------------------------------------------------------------------------
-----------------------------------------------------**-----------------------------
------------------------------------------------------------------------------------
---------------------------------------------------------**-------------------------
-------------------------------------------------------------**------------------**-
-----------------------------------------------------------------------------**-----
-----------------------------------------------------------------**--**--**---------
  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 

min(18)	= -977.108234
max(5)	= -977.105380
save 1.local maxima structure, 2.local mimma structure?
 Type 1 or 2 (Type 0 to escape)

Gauss View 6 では、二点間距離を横軸にしてエネルギーダイヤグラムを表示するため、左右反転している場合があります。おそらく、二点間距離を近づいけていく場合に左右反転しているのだと思います。scan のプロット数が 30 を越えるとラップトップ画面上では見にくくなると思います。

save 1.maximum structure, 2.mimmum structure? と聞かれますので、1 または 2 を入力して Enter を押します。
1 を入力した場合は、local maxima が max_structure.com という名前で保存されます。
2 を入力した場合は、local minima が min_structure.com という名前で保存されます。

ちなみに、output ファイルは、あまり丁寧に形を整えていませんが、一応 gauss view 6 で開くことが出来ます。

記事中に間違い等ある場合は、コメント欄、twitter またはメールにてお知らせいただけると幸いです。

関連する記事

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