自由エネルギーの計算

ここでは最適化された構造を用いthermochemistryモジュールを用いて自由エネルギーを計算する方法を記述します。

  • Pythonスクリプト (ase-qe-opt+vib.py)
    from ase.io import read
    from ase.io import write
    from ase.calculators.espresso import Espresso
    from ase.optimize import BFGS
    from ase.vibrations import Vibrations
    from ase.thermochemistry import IdealGasThermo
    
    input_data = {'control': {'outdir' : './tmp/',
                              'prefix' : 'h2'
                            },
                             'system': {'ecutwfc': 80,
                             'ecutrho': 800},
                             'electrons': {'diagonalization' : 'cg',
                             'conv_thr'        : 1.0e-12,
                             'mixing_beta'     : 0.3}
                  }
    
    pseudopotentials = {'H':'H.pbe-kjpaw_psl.1.0.0.UPF'}
     
    calc = Espresso(pseudopotentials=pseudopotentials,input_data=input_data,
                               tprnfor=True)
    
    h2=read('scf.in',format='espresso-in')
    h2.calc = calc
     
    e_h2 = h2.get_potential_energy()
     
    vib = Vibrations(h2)
    vib.run()
    vib.summary(log='H2_vib_summary.txt')
    vib.write_mode(-1)
    vib_energies = vib.get_energies()
     
    thermo = IdealGasThermo(vib_energies=vib_energies,
                                                potentialenergy=e_h2,
                                                atoms=h2,
                                                geometry='linear',
                                                symmetrynumber=2,spin=0)
    
    G = thermo.get_gibbs_energy(temperature=298.15, pressure=101325.)
    ここで計算に使用する構造は以下のようにQEフォーマットのscf.inというファイルに記述されています。
    h2=read('scf.in',format='espresso-in')
    先ず全エネルギーのSCF計算を行っています。
    e_h2 = h2.get_potential_energy()
    その後、自由エネルギーの計算に必要な(調和)振動数の計算を実行しています。
    vib = Vibrations(h2)
    vib.run()
    vib.summary(log='H2_vib_summary.txt')
    vib.write_mode(-1)
    ここで得られた振動数と零点エネルギーは以下のようにH2_vib_summary.txtに出力されます。
    ---------------------
      #    meV     cm^-1
    ---------------------
      0    0.0i      0.0i
      1    0.0       0.0
      2    0.0       0.0
      3    4.5      35.9
      4    4.5      35.9
      5  535.2    4316.3
    ---------------------
    Zero-point energy: 0.272 eV
    そして自由エネルギーに必要な計算を実行しています。
    thermo = IdealGasThermo(vib_energies=vib_energies,
                                                potentialenergy=e_h2,
                                                atoms=h2,
                                                geometry='linear',
                                                symmetrynumber=2,spin=0)
    symmetrynumberはEssentials of Computational Chemistryを参照して選ぶと良いでしょう。
    最終的に298.15 K、101325 Pa (101.325 kPa (~0.1 MPa))における自由エネルギーを標準出力に以下のように出力しています。
    G = thermo.get_gibbs_energy(temperature=298.15, pressure=101325.)
    Enthalpy components at T = 298.15 K:
    ===============================
    E_pot                -31.747 eV
    E_ZPE                  0.268 eV
    Cv_trans (0->T)        0.039 eV
    Cv_rot (0->T)          0.026 eV
    Cv_vib (0->T)          0.000 eV
    (C_v -> C_p)           0.026 eV
    -------------------------------
    H                    -31.389 eV
    ===============================
    
    Entropy components at T = 298.15 K and P = 101325.0 Pa:
    =================================================
                               S               T*S
    S_trans (1 bar)    0.0012188 eV/K        0.363 eV
    S_rot              0.0001341 eV/K        0.040 eV
    S_elec             0.0000000 eV/K        0.000 eV
    S_vib              0.0000000 eV/K        0.000 eV
    S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV
    -------------------------------------------------
    S                  0.0013518 eV/K        0.403 eV
    =================================================
    
    Free energy components at T = 298.15 K and P = 101325.0 Pa:
    =======================
        H        -31.389 eV
     -T*S         -0.403 eV
    -----------------------
        G        -31.792 eV
    =======================
トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS