GAMESSとmatplotlibを使ってポテンシャルエネルギー曲面を描く

初めに

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

下準備

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

~/.bashrc

export GMSPATH='$HOME/gamess'

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

pot_dim2.py

import os
import subprocess
import numpy as np

# A + BC -> A-B-C -> AB + C 
#
# number:A = na  Free
# number:B = nb  Fixed
# number:C = nc  Free
#
# distance A-B rab
# distance B-C rbc
#

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, nc):
    """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:]])
    xc = np.array([float(x) for x in list_data[nc-1].split()[2:]])
    rab = np.linalg.norm(xa-xb)
    rbc = np.linalg.norm(xb-xc)
    return rab, rbc


def func_sub(s, cof_ab, cof_bc, na, nb, nc, dict_group):
    """Change the coordinates of atoms A and C."""
    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:]])
    xc = np.array([float(x) for x in list_data[nc-1].split()[2:]])

    for na_i in dict_group["a"]:    
        data_a = list_data[na-1].split()[:2]
        xa_i = np.array([float(x) for x in list_data[na_i-1].split()[2:]])   
        xa_new = [str(x) for x in list(xa_i + cof_ab*(xb-xa))]    
        list_data[na-1] = "  ".join(data_a + xa_new)

    for nc_i in dict_group["c"]:    
        data_c = list_data[nc-1].split()[:2]
        xc_i = np.array([float(x) for x in list_data[nc_i-1].split()[2:]])   
        xc_new = [str(x) for x in list(xc_i + cof_bc*(xc-xb))]    
        list_data[nc-1] = "  ".join(data_c + xc_new)

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

def find_energy(s):
    list_line = s.splitlines()
    for line in list_line:
        if "TOTAL ENERGY =" in line:
            energy = float(line.split()[3])
            break
    else:
        print("Not Found Energy.")
        print("Check your the latest logfile.")
        raise SyntaxError()
    return energy
            


def main_pot(input_file, na, nb, nc, 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, rbc = func_calc_r(s_data, na, nb, nc)
    

    #Define coefficients about distances between nuclei.
    n_grid_ab = 5
    n_grid_bc = 10
    array_cof_ab = np.linspace(0.0, 0.9, n_grid_ab)
    array_cof_bc = np.linspace(0.0, 0.9, n_grid_bc)

    path_gms = os.environ["GMSPATH"]

    array_energy = np.array([])
   
    #Run Quntum Chemistry Calculation on each grid.
    for i, cof_ab in enumerate(array_cof_ab):
        for j, cof_bc in enumerate(array_cof_bc):
            print(f"i={i} j={j} Rab={rab*cof_ab} ANGS Rbc={rbc*cof_bc} ANGS Run!!!")

            #Change the coordinates of atoms A and C
            s_data_new = func_sub(s_data, cof_ab, cof_bc, na, nb, nc, dict_group)
            
            newfile = f"{input_file[:-4]}_{i:02d}_{j:02d}.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:02d}_{j:02d}.log", "w") as flog:
                flog.write(output)

            os.system(f"rm -f $HOME/scr/{newfile[:-4]}.*")
            os.system(f"rm -f ./{newfile}")
            energy = find_energy(output)
            print(f"energy = {energy} HARTREE")
            array_energy = np.append(array_energy, energy)
            

    array_rab = array_cof_ab*rab
    array_rbc = array_cof_bc*rbc
    
    np.savetxt("./data_energy.txt", array_energy)
    np.savetxt("./data_rab.txt", array_rab)
    np.savetxt("./data_rbc.txt", array_rbc)

    return 0

計算を実行する

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

ea_hbr.inp

!   File created by MacMolPlt 7.7
 $CONTRL SCFTYP=ROHF RUNTYP=ENERGY MAXIT=200 MULT=2 $END
 $SYSTEM TIMLIM=525600 MEMORY=1000000 $END
 $BASIS GBASIS=STO NGAUSS=6 $END
 $SCF DIRSCF=.TRUE. $END
 $DATA 
Title
C1
Br    35.0     1.27681     0.01778     0.00000
H     1.0     1.27681     1.52778     0.00000
H     1.0     1.30098     5.04680     0.01943
 $END

run.py

import pot_dim2
# A + BC -> A-B-C -> AB + C 
#
# number:A = na  Free
# number:B = nb  Fixed
# number:C = nc  Free
#
# distance A-B rab
# distance B-C rbc

def main():
    dict_group= {"a":[3], "b":[2], "c":[1]}
    pot_dim2.main_pot("ea_hbr.inp", na=3, nb=2, nc=1, 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_rab.txt")
y = np.loadtxt("data_rbc.txt")


nx, ny = len(x), len(y)

X, Y = np.meshgrid(x, y)

Z = np.loadtxt("data_energy.txt").reshape(ny, nx)   # 等高線はZに対して作られる。
# 等高線図の生成。
cont=plt.contour(X,Y,Z, colors=['black'])
cont.clabel(fmt='%1.1f', fontsize=14)


plt.xlabel('R(A-B)', fontsize=24)
plt.ylabel('R(B-C)', fontsize=24)
plt.title("AB+C->A-B-C->A+BC")

plt.pcolormesh(X,Y,Z, cmap='cool') #カラー等高線図
pp=plt.colorbar (orientation="vertical") # カラーバーの表示 
pp.set_label("Potential Energy (Hartree)",  fontsize=24)

plt.show()

得られた結果は以下の通り。
ご利用は自己責任でお願いします。