GAMESSを用いてpythonで実装する変分型遷移状態理論探索プログラム

初めに

GAMESSを使って2次元のポテンシャルエネルギー曲面を描く必要が出てきたのでここに備忘録として記すこととする。ちなみにこのスクリプトのテスト環境はUbuntu。

下準備

スクリプトを動かす前に下準備を整えておく。
まずはGAMESSのスクリプトファイル"rungms"のパスを通しておく。
設定ファイル.bashrcにパスを通しておこう。

~/.bashrc

export GMSPATH='$HOME/gamess'
export GMSSCR='$HOME/scr'

GAMESS計算を自動化するモジュール

vtst.py

import os
import subprocess
import numpy as np

# A-B -> A + B 
#
# number:A = na  Fixed
# number:B = nb  Free
# 
# distance: A-B = r, rab
#

def func_split(s):
    """Separate input data into $DATA group and others."""

    list_line = s.splitlines()
    flag_data_found = False

    for i, line in enumerate(list_line):
        if "$DATA" in line.strip():
            n_line_data = i
            flag_data_found = True
        if flag_data_found and (line.strip() == "$END"):
            n_line_end = i
            break        
    else:
        print("No $Data Group!!!!!")
        raise SyntaxError()

    s_data = "\n".join(list_line[n_line_data:n_line_end+1])
    
    if n_line_end == len(list_line)-1:
        s_others = "\n".join(list_line[:n_line_data])
    elif n_line_data == 0:
        s_others = "\n".join(list_line[n_line_end+1:])
    else:
        s_others = "\n".join(list_line[:n_line_data]+list_line[n_line_end+1:])

    return s_data, s_others


def func_calc_r(s, na, nb):
    """Calculate distances between nuclei."""
    list_data = s.splitlines()[3:]
    xa = np.array([float(x) for x in list_data[na-1].split()[2:]])
    xb = np.array([float(x) for x in list_data[nb-1].split()[2:]])
   
    r = np.linalg.norm(xa-xb)
    
    return r 


def func_sub(s, cof, na, nb, dict_group):
    """Change the coordinates of atoms B."""
    list_data = s.splitlines()[3:]
    xa = np.array([float(x) for x in list_data[na-1].split()[2:]])
    xb = np.array([float(x) for x in list_data[nb-1].split()[2:]])

    for nb_i in dict_group["b"]:    
        data_b = list_data[nb_i-1].split()[:2]
        xb_i = np.array([float(x) for x in list_data[nb_i-1].split()[2:]])   
        xb_new = [str(x) for x in list(xb_i + cof*(xb-xa))]    
        list_data[nb_i-1] = "  ".join(data_b + xb_new)

 
    return "\n".join([" $DATA", "Substituted coordinates", "C1"]+list_data)
    

def find_temp(s):
    list_block = s.split("$END")
    for block in list_block:
        if ("$FORCE" in block.upper()) and ("TEMP(1)=" in block.upper()):
            block_force = block
            break
    else:
        print("No $FORCE group data.")
        raise SyntaxError
    
    list_word = block_force.split("TEMP(1)=")[1].split()
    list_temp = []
    for word in list_word:
        if (word!=word.upper()) or (word!=word.lower()):
            break
        else:
            list_temp.append(float(word))
        
    return np.array(list_temp)


def find_pot_energy(s):
    list_line = reversed(s.splitlines())
    for line in list_line:
        if "TOTAL ENERGY =" in line:
            energy_pot = float(line.split()[3])*2625.50
            return energy_pot
            break
    else: 
        print("No Potential Energy Data!!!")
        raise SyntaxError()


def find_free_energy(s):
    list_data = s.split("THERMOCHEMISTRY AT T=")[1:]
    list_energy = []
    list_temp = []
    for s_data in list_data:
        list_line = s_data.splitlines()
        temp = float(list_line[0].split()[0])
        list_temp.append(temp)
        for i, line in enumerate(list_line):
            if "KJ/MOL    KJ/MOL    KJ/MOL   J/MOL-K   J/MOL-K   J/MOL-K" in line:
                index_unit = i
                break
        else:
            print("No Thermochemistry Data!!!")
            raise SyntaxError()

        energy_g = float(list_line[index_unit+5].split()[3])                
        list_energy.append(energy_g)
 
    return np.array(list_energy), np.array(list_temp)
            


def main_energy_diss(input_file, na, nb, dict_group, ncore=1):
    """main function"""
    with open(input_file) as finput:
        s = finput.read()

    #Separate input data into $DATA group and others.
    s_data, s_others = func_split(s)

    #check for C1 symmetry
    if s_data.splitlines()[2].strip().upper() != "C1":
        raise SyntaxError()

    #calculate distances between nuclei.
    rab = func_calc_r(s_data, na, nb)

    path_gms = os.environ["GMSPATH"]
    path_gms_scr = os.environ["GMSSCR"]


    #Define coefficients about distances between nuclei.
    array_cof_ab = np.arange(0.0, 10.0, 0.05)/rab
    
    num_temp = len(find_temp(s))
    array_energy = np.zeros([num_temp, len(array_cof_ab)])

    #Run Quntum Chemistry Calculation on each grid.
    for i, cof_ab in enumerate(array_cof_ab):
        print(f"i={i} Rab={rab*(1.0+cof_ab)} ANGS Run!!!")

        #Change the coordinates of atoms B
        s_data_new = func_sub(s_data, cof_ab, na, nb, dict_group)
            
        newfile = f"{input_file[:-4]}_{i:04d}.inp"
        with open(newfile, "w") as fnew:
            fnew.write("\n".join([s_data_new, s_others]))
        output = subprocess.run(f"{path_gms}/rungms {newfile} 00 {ncore}".split(), capture_output=True, text=True).stdout
        with open(f"./dir_log/{input_file[:-4]}_{i:04d}.log", "w") as flog:
            flog.write(output)

        os.system(f"mv $GMSSCR/{newfile.split('.')[0]}.dat ./dir_log/")
        os.system(f"rm -f $GMSSCR/{newfile.split('.')[0]}.*")
        os.system(f"rm -f ./{newfile}")
        energy_pot = find_pot_energy(output)
        array_energy_g, array_temp = find_free_energy(output)
        array_temp_energy = array_energy_g + energy_pot
        print(f"energy = {array_temp_energy}")
        for j, energy in enumerate(array_temp_energy):
            array_energy[j][i] = energy
            
 
    array_rab = (1.0 + array_cof_ab)*rab
    
    np.savetxt("./data_Gibbs_energy.txt", array_energy)
    np.savetxt("./data_diss_rab.txt", array_rab)
    np.savetxt("./data_diss_temp.txt", array_temp)

    return 0

計算を実行する

今回は以下の反応の2次元ポテンシャル曲面を描く。
GAMESSのインプットファイルと実行スクリプトは以下の通り。

opt.inp

 $BASIS DIFFS=.T. DIFFSP=.T. GBASIS=N31 NGAUSS=6 $END
 $CONTRL MAXIT=200 MULT=1 RUNTYP=OPTIMIZE SCFTYP=RHF $END
 $DATA 
Title
C1
C     6.0     0.09348     0.08911     0.03863
H     1.0    -0.97742    -0.11856     0.03364
H     1.0     0.88571    -0.65057    -0.04802
O     8.0     0.46597     1.26945     0.15245
O     8.0    -0.45055     2.23124     0.26365
 $END
 $DFT DC=.T. IDCVER=4 $END
 $FORCE METHOD=ANALYTIC
  TEMP(1)=0.001 200 300 400 500 600 700 800 900
 $END
 $SCF DAMP=.T. DIRSCF=.T. SOSCF=.T. $END
 $STATPT HESS=CALC HSSEND=.T. IHREP=10 NSTEP=500 OPTTOL= 1.000000E-04 
  IFREEZ(1)=10, 11, 12, 13, 14, 15
 $END
 $SYSTEM MWORDS=5 $END

run.py

import vtst as diss

def main():
    dict_group= {"a":[1,2,3,4], "b":[5]}
    diss.main_energy_diss("opt.inp", na=4, nb=5, dict_group=dict_group, ncore=2)


if __name__ == "__main__":
    main()

得られた自由エネルギー曲線のデータをプロットし、最大値を算出する。

plot.py

import numpy as np
import matplotlib.pyplot as plt

x = np.loadtxt("data_diss_rab.txt")
x = x.flatten()
y = np.loadtxt("data_Gibbs_energy.txt")
temp = np.loadtxt("data_diss_temp.txt")
print(temp)
print(x)
print(y)

print(f"G_MAX={np.max(y)} kJ/mol  R(A-B)={x[y==np.max(y)]} ANGS")

plt.xlabel('R(O-O) ANGS', fontsize=24)
plt.ylabel("Gibbs Energy kJ/mol", fontsize=24)
plt.title("CH2OO -> CH2O + O(1D)", fontsize=20)
plt.plot(x, y)
#plt.legend(loc="best")
plt.show()

#plt.savefig("pot_colormap.png", dpi=150)

ご利用は自己責任でお願いします。

opt.inp

 $BASIS DIFFS=.T. DIFFSP=.T. GBASIS=N311 NDFUNC=1 NGAUSS=6 NPFUNC=1 $END
 $CONTRL DFTTYP=M06-2X MAXIT=200 MULT=1 RUNTYP=OPTIMIZE SCFTYP=UHF $END
 $DATA 
Title
C1
S    16.0    -0.87134    -0.39017    -0.07683
O     8.0    -1.25733     0.96906     0.05849
O     8.0     0.73125    -0.51869     0.01231
O     8.0    -1.00351    -1.02824     1.46306
O     8.0     0.20743    -1.54636     1.87957
C     6.0     1.14390    -0.67642     1.34613
H     1.0     1.13126     0.27582     1.85903
H     1.0     2.10729    -1.15635     1.36504
 $END
 $FORCE METHOD=SEMINUM $END
 $SCF DAMP=.T. DIRSCF=.T. SOSCF=.T. $END
 $STATPT HESS=CALC HSSEND=.T. IHREP=10 NSTEP=500 OPTTOL= 1.000000E-04 $END
 $SYSTEM MWORDS=5 $END