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

構造最適化と振動数の計算 †

ここではQuick Startに引き続き水素分子の構造最適化と振動数の計算方法を示します。

• 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
  
  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('relax.in',format='espresso-in')
  h2.calc = calc
  BFGS(h2).run(fmax=0.01)
  write('H2_opt.xyz',h2)
  etot = h2.get_potential_energy()
  print('Total energy: %12.8f eV' % etot)
  
  vib = Vibrations(h2)
  vib.run()
  vib.summary(log='H2_vib_summary.txt')
  vib.write_mode(-1)

  ここで計算に使用する構造は以下のようにQEフォーマットのrelax.inというファイルに記述されています。

  h2=read('relax.in',format='espresso-in')

  このスクリプトでは以下で構造最適化を実行しています。

  BFGS(h2).run(fmax=0.01)

  以下で最終的な全エネルギーを標準出力に出力しています。

  etot = h2.get_potential_energy()
  print('Total energy: %12.8f eV' % etot)

  その後、最適化された構造を使って以下のように（調和）振動数の計算を実行しています。

  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

参考 †

• Vibrational analysis