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

Quick start †

水素分子を例にとりASEを用いてQEを実行する方法を説明します。

従来のQEの実行方法 †

先ず最初に従来の方法でQEを実行しましょう。 このチュートリアルではxeシステム (xe1) 上で8コア、あるいはxe2で12コア使用することを想定しています。

• 入力ファイル (h2_scf.in)

  &control
     calculation   = 'scf'
     restart_mode  = 'from_scratch'
     pseudo_dir    = '/home/ikutaro/QE/pseudo/'
     outdir        = './tmp/'
     prefix        = 'h2'
     tprnfor       = .true.
     forc_conv_thr = 1.d-4
  /
  &system
     ibrav   = 1
     A       = 10.d0
     nat     = 2
     ntyp    = 1
     nbnd    = 4
     ecutwfc = 30.0
     ecutrho = 120.0
  /
  &electrons
     diagonalization = 'david'
     conv_thr        = 1.0e-12
     mixing_beta     = 0.1
  /
  ATOMIC_SPECIES
   H 0.0000  H_ONCV_PBE-1.0.upf
  ATOMIC_POSITIONS (angstrom)
  H      0.000000000000      0.000000000000      0.00000000000
  H      0.000000000000      0.000000000000      0.74000000000
  K_POINTS (gamma)

• ジョブスクリプト (qsub-qe.sh)

  #$ -S /bin/bash
  #$ -cwd
  #$ -q xe1.q
  #$ -pe x8 8
  #$ -N H2
  #
  module load intel/2020.2.254
  module load intelmpi/2020.2.254
   
  INPUT_FILE='h2_scf.in'
  OUTPUT_FILE='h2_scf.out'
   
  PW_DIR=/home/ikutaro/QE/src/qe-6.6/bin
   
  MPI_COMMAND=mpirun
  PW_COMMAND=pw.x
  PW=$PW_DIR/$PW_COMMAND
   
  I_MPI_PIN=1
  I_MPI_FABRICS=shm:ofa
  OMP_NUM_THREADS=1
   
  if [ $OMP_NUM_THREADS -gt 1 ];
  then
  cat $PE_HOSTFILE | awk '{ print $1":"$2/ENVIRON["OMP_NUM_THREADS"] }' > hostfile.$JOB_ID
  fi
   
  $MPI_COMMAND $PW < $INPUT_FILE > $OUTPUT_FILE

• ジョブのサブミット

  $ qsub qsub-qe.sh

  ジョブが正常に終わるとh2_scf.outというファイルが出力されます。

ASEを用いたQEの実行方法 †

Single point (SCF) †

入力ファイルの指定の仕方は一通りではありませんが、以下の例ではQEの入力ファイル (h2_scf.in) を用いて構造を入力しています（通常のQE計算の結果やリソースを可能な限り再利用する、という方針です）。ASEを使用してQEの入力ファイルに記述された構造をASEのフォーマットに変換することも可能でしょう。

• ASE用pythonスクリプト (ase-qe.py)

  from ase.io import read
  from ase.io import write
  from ase.calculators.espresso import Espresso
  qe_input_data = {'control': {'outdir' : './tmp/',
                               'prefix' : 'h2'
                              },
                   'system':  {'nbnd'   : 4,
                               'ecutwfc': 30,
                               'ecutrho': 120
                              },
                   'electrons': {'diagonalization' : 'david',
                                 'conv_thr'        : 1.0e-12,
                                 'mixing_beta'     : 0.1
                                }
                   }
  pseudopotentials = {'H':'H_ONCV_PBE-1.0.upf'}
  calc = Espresso(pseudopotentials=pseudopotentials,input_data=qe_input_data)
  h2=read('h2_scf.in',format='espresso-in')
  h2.calc = calc
  etot=h2.get_potential_energy()
  print('Total energy: %12.8f eV' % etot)

  このスクリプトでは以下でh2_scf.inを読み込んでいます。

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

  構造以外のパラメーターはqe_input_dataで指定しています。パラメータは必要に応じて追加します。 また

  print('Total energy: %12.8f eV' % etot)

  のようにprintを実行しないと全エネルギーは出力されないので注意が必要です。
• ジョブスクリプト

  #$ -S /bin/bash
  #$ -cwd
  #$ -q xe1.q
  #$ -pe x8 8
  #$ -N H2
  #$ -o H2_out
  #$ -e H2_err
   
  module load intel/2020.2.254
  module load intelmpi/2020.2.254
   
  export ASE_ESPRESSO_COMMAND="mpirun ${HOME}/QE/src/qe-6.6/bin/pw.x -in PREFIX.pwi > PREFIX.pwo"
  export ESPRESSO_PSEUDO="/home/ikutaro/QE/pseudo"
   
  python3 ase-qe.py

  ASE_ESPRESSO_COMMANDとESPRESSO_PSEUDOを適切に指定する必要があります。この例ではQE (pw.x) の通常の入力はh2_scf.inではなくPREFIX.pwi、出力はPREFIX.pwoとなっています。計算に問題が起こった時はASEにより生成されたこれらのファイルの内容を確認すると良いでしょう。
• ジョブのサブミット

  $ qsub qsub-ase.sh

  計算が正常に実行されると、H2_outに

  Total energy: -31.48786645 eV

  が出力されます。

構造最適化 †

次に構造最適化を行いましょう。ここでBFGS法を用いることにします。

• 入力ファイル (ase-qe-opt.py)

  from ase.io import read
  from ase.io import write
  from ase.calculators.espresso import Espresso
  from ase.optimize import BFGS
  input_data = {'control': {'outdir' : './tmp/',
                            'prefix' : 'h2'
                          },
                'system': {'nbnd'   : 4,
                           'ecutwfc': 30,
                           'ecutrho': 120},
                'electrons': {'diagonalization' : 'david',
                              'conv_thr'        : 1.0e-12,
                              'mixing_beta'     : 0.1}
                }
  pseudopotentials = {'H':'H_ONCV_PBE-1.0.upf'}
  calc = Espresso(pseudopotentials=pseudopotentials,input_data=input_data,
                  tprnfor=True)
  h2=read('h2_scf.in',format='espresso-in')
  h2.calc = calc
  opt = BFGS(h2)
  opt.run(fmax=0.02)
  write('H2_opt.xyz',h2)
  etot=h2.get_potential_energy()
  print('Total energy: %12.8f eV' % etot)

• ジョブスクリプト

  #$ -S /bin/bash
  #$ -cwd
  #$ -q xe1.q
  #$ -pe x8 8
  #$ -N H2
  #$ -o H2_out
  #$ -e H2_err
   
  module load intel/2020.2.254
  module load intelmpi/2020.2.254
   
  export ASE_ESPRESSO_COMMAND="mpirun ${HOME}/QE/src/qe-6.6/bin/pw.x -in PREFIX.pwi > PREFIX.pwo"
  export ESPRESSO_PSEUDO="/home/ikutaro/QE/pseudo"
   
  python3 ase-qe-opt.py

  計算が正常に終了するとH2_outに以下のような結果が出力されます。

  Total energy: -31.48786645 eV
        Step     Time          Energy         fmax
  BFGS:    0 15:43:16      -31.487866        0.7746
  BFGS:    1 15:43:18      -31.496494        0.0197
  Total energy: -31.49649423 eV

  通常通りQEで構造最適化してみてエネルギーや構造、最適化にかかる時間やステップ数なども比較してみましょう。

参考 †

• ASE web page (espresso)
• 原子分子シミュレーションPythonパッケージAtomic Simulation Environment (ASE)を使ってみる　その２：Quantum Espressoで第一原理計算