H2 tutorial

From Ecal

Jump to: navigation, search

H2が箱に二個入った例。これはたぶんなかなかMTOのみでいい結果を得るのはしんどいです.

MTO+APW(PMT)法のrobustnessにより、ctrl fileが簡単につくることが可能になりました。 まずは、 SrTiO3の例ぜひやってほしいです。 説明が重複しないように、こちらの例は簡単に記述します。

Contents

ctrl fileの作り方

ctrl file作成のために最初に"ctrls.*"にcrystal structureを記入する必要があります。

ctrlでなくctrlsです。その後それをスクリプトctrlgen.pyで加工してctrl.*をつくります。

ここにctrlgen.py|バージョンです(2010Mar1stのパッケージはほんのすこし古い;以後 のecalパッケージに入れます)。

結晶構造のファイル"ctrls.{name}"の用意

結晶構造だけのファイルを作ります。 ctrl file作成のために最初に"ctrls.*"にcrystal structureを記入する必要があります。

ctrlでなくctrlsです。
%const a=5.0 c=15.0 d=1.398
STRUC   ALAT={a}
        PLAT= 1.0 0.0 0.0
               0.0 1.0 0.0 
               0.0 0.0 {c}/{a}
SITE    ATOM=H POS= 0.5*{d}/{a} 0.0 0.0
        ATOM=H POS=-0.5*{d}/{a} 0.0 0.0
        ATOM=H POS= 0.5*{d}/{a} 0.0 {c}/{a}/2
        ATOM=H POS=-0.5*{d}/{a} 0.0 {c}/{a}/2

これをctrls.hとしてセーブします。

  • SrTiO3の例ぜひ見てください.説明を書いています。原子の計算などでMT半径が計算できない場合があります。そのときはマニュアル設定をrmt.tmpに書かないといけない.この場合は分子なのでtouchingでMT半径が選ばれます。

ctrlgen.py scriptの実行

ctrlgen.py h
あるいは
~/ecal/lm7K/ctrlgen.py h

と実行します。これは内部で他のprogramを実行します。もしエラーが起きたらPATHが正しく設定されているのかをもう一度調べてください。lmfaとかにパスが通ってる必要があります。ctrls.hは変更されません。


ctrlgen.pyで生成されたctrl file

ctrls.hからctrlgen.pyにより作られたctrl.hは以下になっています。 このファイルはlmfa,lmfを実行するためのすべてのパラメタをすでに含んでいます。

### Do lmf --input to see all effective category and token ###
### It will be not so difficult to edit ctrlge.py for your purpose ###
VERS    LM=7 FP=7
             # version check. Fixed.
IO      SHOW=T VERBOS=35
             # SHOW=T shows readin data (and default setting at the begining of console output)
	     # It is useful to check ctrl is read in correctly or not (equivalent with --show option).
	     #
	     # lerger VERBOSE gives more detailed console output.
SYMGRP find  # 'find' evaluate space-group symmetry automatically.
             # Usually 'find is OK', but lmf may use lower symmetry
	     # because of numerical problem.
             # Do lmchk to check how it is evaluated.
%const kmxa=7  # kmxa=7 is good for pwemax=5 or lower.
               # larger kmxa is better but time-consuming (maybe not the critical part for large systems).
%const a=5.0 c=15.0 d=1.398
STRUC   ALAT={a}
        PLAT= 1.0 0.0 0.0
               0.0 1.0 0.0 
               0.0 0.0 {c}/{a}
  NBAS= 4  NSPEC=1
SITE    ATOM=H POS= 0.5*{d}/{a} 0.0 0.0
        ATOM=H POS=-0.5*{d}/{a} 0.0 0.0
        ATOM=H POS= 0.5*{d}/{a} 0.0 {c}/{a}/2
        ATOM=H POS=-0.5*{d}/{a} 0.0 {c}/{a}/2
SPEC
 ATOM= H Z= 1 R=0.699000
      RSMH= 0.500  EH= -0.108
      KMXA={kmxa} LMXA=4
      MMOM=0 0 0


% const pwemax=3 nk=2
BZ    NKABC={nk} {nk} {nk}  # division of BZ for q points.
      METAL=3   # METAL=3 is safe setting. For insulator, METAL=0 is good enough.
		# When you plot dos, set SAVDOS=T and METAL=3, and with DOS=-1 1 (range) NPTS=2001 (division) even for insulator.
		#   (SAVDOS, DOS, NPTS gives no side effect for self-consitency calculaiton).
                # 
                #BUG: For a hydrogen in a large cell, I(takao) found that METAL=0 for
                #(NSPIN=2 MMOM=1 0 0) results in non-magnetic solution. Use METAL=3 for a while in this case.
                # 
      BZJOB=0	# BZJOB=0 (including Gamma point) or =1 (not including Gamma point).
		#  In cases , BZJOB=1 makes calculation efficient.

ITER  CONV=1e-6 CONVC=1e-6 NIT=30
                # An other choice is
                # ITER MIX=A2,b=.5,n=3 CONV=1e-6 CONVC=1e-6 NIT=20
                # Practically results are independent from mixing procedure.
		
HAM   NSPIN=1   # Set NSPIN=2 for spin-polarize case; then set SPEC_MMOM (initial guess of magnetic polarization).
      FORCES=0  # 0: no force calculation, 1: forces calculaiton 
      GMAX=9    # this is for real space mesh. See GetStarted.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for VWN; GGA is not yet.
                # XCFUN=2 shows a bug for Hydrogen atom. 
		# (subs/evxc.F works only for XCFUN=1 if rho(up)=0 or rho(down)=0).

      PWMODE=11 # 10: MTO basis only (LMTO) PW basis is not used.
                # 11: APW+MTO        (PMT)
                # 12: APW basis only (LAPW) MTO basis is not used.

      PWEMAX={pwemax} # (in Ry). When you use larger pwemax more than 5, be careful
                      # about overcompleteness. See GetStarted.
      ELIND=-1  # this is only for accelaration of convergence. Not need to change.
      
OPTIONS PFLOAT=1 # Q=band (this is quit switch if you like to add)
                 # 

lmfa,lmfの実行

lmfを実行する準備が整います。

lmfa h
lmf h

を実行してください。とにかく結果の正しさは別にして「安全性の高い」計算ができる設定なので まずこのままやってください。

原子計算の場合は、real-space mesh (HAM_FTMESH or HAM_GMAX)がセンシティブです。 lmfでは、分子計算、原子計算がまだあんまりやられてないので情報をためていくようお願いしたいです.

ctrlを書き換えて収束チェックや格子定数のスクリプト計算

ctrl.hを以下のように書き換えてみました。以下のファイルでcsh ctrl.hとすると、計算が一気に連続して行えます。 最初の部分がcshでもbashでも使えるのでわかりやすいと思います。

#!/bin/csh
set kmxa=0
set pwemax=0
lmfa h > out_lmfa

# rule of thumb. kmxa > pwemax
foreach gmax   (9 12 15 20)
foreach kmxa   (4 8 12 16 20)
foreach pwemax (4 8 12 16 20) 
 lmf h -vkmxa=$kmxa -vpwemax=$pwemax -vkmxa=$kmxa -vgmax=$gmax > out_pwemax${pwemax}_kmxa${kmxa}_gmax${gmax}
 echo out_pwemax${pwemax}_kmxa${kmxa} $status
end
end
end
exit

ctrlstart ==============================
### Do lmf --input to see all effective category and token ###
### It will be not so difficult to edit ctrlge.py for your purpose ###
VERS    LM=7 FP=7
             # version check. Fixed.
IO      SHOW=T VERBOS=35
             # SHOW=T shows readin data (and default setting at the begining of console output)
	     # It is useful to check ctrl is read in correctly or not (equivalent with --show option).
	     #
	     # lerger VERBOSE gives more detailed console output.
SYMGRP find  # 'find' evaluate space-group symmetry automatically.
             # Usually 'find is OK', but lmf may use lower symmetry
	     # because of numerical problem.
             # Do lmchk to check how it is evaluated.
%const kmxa=20 # kmxa=7 is good for pwemax=5 or lower.
               # larger kmxa is better but time-consuming (maybe not the critical part for large systems).
%const a=5.0 c=15.0 d=1.398
STRUC   ALAT={a}
        PLAT= 1.0 0.0 0.0
               0.0 1.0 0.0 
               0.0 0.0 {c}/{a}
  NBAS= 4  NSPEC=1
SITE    ATOM=H POS= 0.5*{d}/{a} 0.0 0.0
        ATOM=H POS=-0.5*{d}/{a} 0.0 0.0
        ATOM=H POS= 0.5*{d}/{a} 0.0 {c}/{a}/2
        ATOM=H POS=-0.5*{d}/{a} 0.0 {c}/{a}/2
SPEC
 ATOM= H Z= 1 R=0.699000
      RSMH= 0.500  EH= -0.108
      KMXA={kmxa} LMXA=4
      MMOM=0 0 0

% const pwemax=2 nk=2 gmax=9
BZ    NKABC={nk} {nk} {nk}  # division of BZ for q points.
      METAL=0   # METAL=3 is safe setting. For insulator, METAL=0 is good enough.
		# When you plot dos, set SAVDOS=T and METAL=3, and with DOS=-1 1 (range) NPTS=2001 (division) even for insulator.
		#   (SAVDOS, DOS, NPTS gives no side effect for self-consitency calculaiton).
                # 
                #BUG: For a hydrogen in a large cell, I(takao) found that METAL=0 for
                #(NSPIN=2 MMOM=1 0 0) results in non-magnetic solution. Use METAL=3 for a while in this case.
                # 
      BZJOB=0	# BZJOB=0 (including Gamma point) or =1 (not including Gamma point).
		#  In cases , BZJOB=1 makes calculation efficient.

ITER  CONV=1e-6 CONVC=1e-6 NIT=30
                # An other choice is
                # ITER MIX=A2,b=.5,n=3 CONV=1e-6 CONVC=1e-6 NIT=20
                # Practically results are independent from mixing procedure.
		
HAM   NSPIN=1   # Set NSPIN=2 for spin-polarize case; then set SPEC_MMOM (initial guess of magnetic polarization).
      FORCES=0  # 0: no force calculation, 1: forces calculaiton 
      GMAX={gmax}    # this is for real space mesh. See GetStarted.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for VWN; GGA is not yet.
                # XCFUN=2 shows a bug for Hydrogen atom. 
		# (subs/evxc.F works only for XCFUN=1 if rho(up)=0 or rho(down)=0).

      PWMODE=11 # 10: MTO basis only (LMTO) PW basis is not used.
                # 11: APW+MTO        (PMT)
                # 12: APW basis only (LAPW) MTO basis is not used.

      PWEMAX={pwemax} # (in Ry). When you use larger pwemax more than 5, be careful
                      # about overcompleteness. See GetStarted.
      ELIND=-1  # this is only for accelaration of convergence. Not need to change.
      
OPTIONS PFLOAT=1 # Q=band (this is quit switch if you like to add)
                 # 

これを実行すると(まだ途中までなんですが)以下のような感じの結果です。つかったバージョンは、

commit e4abe8dc15ec57b6d7900891d3cd00df1dbabdb1
Author: Takao Kotani <takaokotani@yahoo.co.jp>
Date:   Mon Mar 1 14:28:25 2010 +0900

    ropbes.F is replaced.

です。以下が結果。固有値の収束チェックとかもお願いしたいです.1-shotGWならPMTでもできます(のはず。できなければ すぐ直します)。

takao@mar:~/ecal/lm7K/TESTsamples/H_atoms_miyake3$ grep 'c ' save.h 
c ehf=-4.483031 ehk=-4.4830943
c kmxa=4 pwemax=4 gmax=9 ehf=-4.4933326 ehk=-4.4933791
c kmxa=4 pwemax=8 gmax=9 ehf=-4.5100143 ehk=-4.5100441
c kmxa=4 pwemax=12 gmax=9 ehf=-4.5137302 ehk=-4.5137546
c kmxa=4 pwemax=16 gmax=9 ehf=-4.5151555 ehk=-4.5151777
c kmxa=4 pwemax=20 gmax=9 ehf=-4.5157507 ehk=-4.5157726
c kmxa=8 pwemax=4 gmax=9 ehf=-4.4933306 ehk=-4.4933789
c kmxa=8 pwemax=8 gmax=9 ehf=-4.5100146 ehk=-4.5100444
c kmxa=8 pwemax=12 gmax=9 ehf=-4.5137306 ehk=-4.5137549
c kmxa=8 pwemax=16 gmax=9 ehf=-4.5151559 ehk=-4.5151781
c kmxa=8 pwemax=20 gmax=9 ehf=-4.515751 ehk=-4.5157729
c kmxa=12 pwemax=4 gmax=9 ehf=-4.4933306 ehk=-4.4933789
c kmxa=12 pwemax=8 gmax=9 ehf=-4.5100146 ehk=-4.5100444
c kmxa=12 pwemax=12 gmax=9 ehf=-4.5137306 ehk=-4.5137549
c kmxa=12 pwemax=16 gmax=9 ehf=-4.5151559 ehk=-4.5151781
c kmxa=12 pwemax=20 gmax=9 ehf=-4.515751 ehk=-4.5157729
c kmxa=16 pwemax=4 gmax=9 ehf=-4.4933306 ehk=-4.4933789
c kmxa=16 pwemax=8 gmax=9 ehf=-4.5100146 ehk=-4.5100444
c kmxa=16 pwemax=12 gmax=9 ehf=-4.5137306 ehk=-4.5137549
c kmxa=16 pwemax=16 gmax=9 ehf=-4.5151559 ehk=-4.5151781
c kmxa=16 pwemax=20 gmax=9 ehf=-4.515751 ehk=-4.5157729
c kmxa=20 pwemax=4 gmax=9 ehf=-4.4933306 ehk=-4.4933789
c kmxa=20 pwemax=8 gmax=9 ehf=-4.5100146 ehk=-4.5100444
c kmxa=20 pwemax=12 gmax=9 ehf=-4.5137306 ehk=-4.5137549
c kmxa=20 pwemax=16 gmax=9 ehf=-4.5151559 ehk=-4.5151781
c kmxa=20 pwemax=20 gmax=9 ehf=-4.515751 ehk=-4.5157729
c kmxa=4 pwemax=4 gmax=12 ehf=-4.4949894 ehk=-4.4949896
c kmxa=4 pwemax=8 gmax=12 ehf=-4.5113217 ehk=-4.5113219
c kmxa=4 pwemax=12 gmax=12 ehf=-4.5149105 ehk=-4.5149105
c kmxa=4 pwemax=16 gmax=12 ehf=-4.5162939 ehk=-4.516294
c kmxa=4 pwemax=20 gmax=12 ehf=-4.5168912 ehk=-4.5168912
c kmxa=8 pwemax=4 gmax=12 ehf=-4.4949892 ehk=-4.4949894
c kmxa=8 pwemax=8 gmax=12 ehf=-4.511322 ehk=-4.5113221
c kmxa=8 pwemax=12 gmax=12 ehf=-4.5149108 ehk=-4.5149109
c kmxa=8 pwemax=16 gmax=12 ehf=-4.5162942 ehk=-4.5162943
c kmxa=8 pwemax=20 gmax=12 ehf=-4.5168915 ehk=-4.5168915
c kmxa=12 pwemax=4 gmax=12 ehf=-4.4949892 ehk=-4.4949894
c kmxa=12 pwemax=8 gmax=12 ehf=-4.511322 ehk=-4.5113221
c kmxa=12 pwemax=12 gmax=12 ehf=-4.5149108 ehk=-4.5149109
c kmxa=12 pwemax=16 gmax=12 ehf=-4.5162942 ehk=-4.5162943
c kmxa=12 pwemax=20 gmax=12 ehf=-4.5168915 ehk=-4.5168915
c kmxa=16 pwemax=4 gmax=12 ehf=-4.4949892 ehk=-4.4949894
c kmxa=16 pwemax=8 gmax=12 ehf=-4.511322 ehk=-4.5113221
c kmxa=16 pwemax=12 gmax=12 ehf=-4.5149108 ehk=-4.5149109
c kmxa=16 pwemax=16 gmax=12 ehf=-4.5162942 ehk=-4.5162943
c kmxa=16 pwemax=20 gmax=12 ehf=-4.5168915 ehk=-4.5168915
c kmxa=20 pwemax=4 gmax=12 ehf=-4.4949892 ehk=-4.4949894
c kmxa=20 pwemax=8 gmax=12 ehf=-4.511322 ehk=-4.5113221
c kmxa=20 pwemax=12 gmax=12 ehf=-4.5149108 ehk=-4.5149109
c kmxa=20 pwemax=16 gmax=12 ehf=-4.5162942 ehk=-4.5162943
c kmxa=20 pwemax=20 gmax=12 ehf=-4.5168915 ehk=-4.5168915
c kmxa=4 pwemax=4 gmax=15 ehf=-4.4950206 ehk=-4.4950207
c kmxa=4 pwemax=8 gmax=15 ehf=-4.5113457 ehk=-4.5113458
c kmxa=4 pwemax=12 gmax=15 ehf=-4.5149318 ehk=-4.5149318
c kmxa=4 pwemax=16 gmax=15 ehf=-4.5163145 ehk=-4.5163145
c kmxa=4 pwemax=20 gmax=15 ehf=-4.516912 ehk=-4.516912
c kmxa=8 pwemax=4 gmax=15 ehf=-4.4950203 ehk=-4.4950204
c kmxa=8 pwemax=8 gmax=15 ehf=-4.511346 ehk=-4.5113461
c kmxa=8 pwemax=12 gmax=15 ehf=-4.5149321 ehk=-4.5149321
c kmxa=8 pwemax=16 gmax=15 ehf=-4.5163148 ehk=-4.5163148
c kmxa=8 pwemax=20 gmax=15 ehf=-4.5169123 ehk=-4.5169123
c kmxa=12 pwemax=4 gmax=15 ehf=-4.4950203 ehk=-4.4950204
c kmxa=12 pwemax=8 gmax=15 ehf=-4.511346 ehk=-4.5113461
c kmxa=12 pwemax=12 gmax=15 ehf=-4.5149321 ehk=-4.5149321
c kmxa=12 pwemax=16 gmax=15 ehf=-4.5163148 ehk=-4.5163148
c kmxa=12 pwemax=20 gmax=15 ehf=-4.5169123 ehk=-4.5169123
c kmxa=16 pwemax=4 gmax=15 ehf=-4.4950203 ehk=-4.4950204
c kmxa=16 pwemax=8 gmax=15 ehf=-4.511346 ehk=-4.5113461
c kmxa=16 pwemax=12 gmax=15 ehf=-4.5149321 ehk=-4.5149321
c kmxa=16 pwemax=16 gmax=15 ehf=-4.5163148 ehk=-4.5163148
c kmxa=16 pwemax=20 gmax=15 ehf=-4.5169123 ehk=-4.5169123
c kmxa=20 pwemax=4 gmax=15 ehf=-4.4950203 ehk=-4.4950204
c kmxa=20 pwemax=8 gmax=15 ehf=-4.511346 ehk=-4.5113461
c kmxa=20 pwemax=12 gmax=15 ehf=-4.5149321 ehk=-4.5149321
c kmxa=20 pwemax=16 gmax=15 ehf=-4.5163148 ehk=-4.5163148
c kmxa=20 pwemax=20 gmax=15 ehf=-4.5169123 ehk=-4.5169123
c kmxa=4 pwemax=4 gmax=20 ehf=-4.4950241 ehk=-4.4950241
c kmxa=4 pwemax=8 gmax=20 ehf=-4.5113463 ehk=-4.5113464
c kmxa=4 pwemax=12 gmax=20 ehf=-4.5149322 ehk=-4.5149322
c kmxa=4 pwemax=16 gmax=20 ehf=-4.5163148 ehk=-4.5163148
c kmxa=4 pwemax=20 gmax=20 ehf=-4.5169124 ehk=-4.5169124
c kmxa=8 pwemax=4 gmax=20 ehf=-4.4950238 ehk=-4.4950239
c kmxa=8 pwemax=8 gmax=20 ehf=-4.5113466 ehk=-4.5113466
c kmxa=8 pwemax=12 gmax=20 ehf=-4.5149325 ehk=-4.5149325
c kmxa=8 pwemax=16 gmax=20 ehf=-4.5163152 ehk=-4.5163152
c kmxa=8 pwemax=20 gmax=20 ehf=-4.5169127 ehk=-4.5169127
c kmxa=12 pwemax=4 gmax=20 ehf=-4.4950238 ehk=-4.4950239
c kmxa=12 pwemax=8 gmax=20 ehf=-4.5113466 ehk=-4.5113466
c kmxa=12 pwemax=12 gmax=20 ehf=-4.5149325 ehk=-4.5149325
c kmxa=12 pwemax=16 gmax=20 ehf=-4.5163152 ehk=-4.5163152
c kmxa=12 pwemax=20 gmax=20 ehf=-4.5169127 ehk=-4.5169127
c kmxa=16 pwemax=4 gmax=20 ehf=-4.4950238 ehk=-4.4950239
c kmxa=16 pwemax=8 gmax=20 ehf=-4.5113466 ehk=-4.5113466
c kmxa=16 pwemax=12 gmax=20 ehf=-4.5149325 ehk=-4.5149325
c kmxa=16 pwemax=16 gmax=20 ehf=-4.5163152 ehk=-4.5163152
c kmxa=16 pwemax=20 gmax=20 ehf=-4.5169127 ehk=-4.5169127
c kmxa=20 pwemax=4 gmax=20 ehf=-4.4950238 ehk=-4.4950239
c kmxa=20 pwemax=8 gmax=20 ehf=-4.5113466 ehk=-4.5113466
c kmxa=20 pwemax=12 gmax=20 ehf=-4.5149325 ehk=-4.5149325
c kmxa=20 pwemax=16 gmax=20 ehf=-4.5163152 ehk=-4.5163152
c kmxa=20 pwemax=20 gmax=20 ehf=-4.5169127 ehk=-4.5169127
Personal tools