第一原理分子動力学プログラム STATE Senri Wiki

自由エネルギーの計算 †

ここでは最適化された構造を用いthermochemistryモジュールのIdealGasThermoを用いて分子（理想気体）の自由エネルギーを計算する方法を記述します。

• Pythonスクリプト (ase-vib+thermo.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
  =======================

参考 †

• Thermochemistry