O2tutorial

From Ecal

Jump to: navigation, search

←戻る

Contents

O2分子の計算の仕方

まず、インストールし、make checkをとおしておくこと。「make install」で必要なものが自分のbin(~/bin)にコピーされる。

構造ファイルを書く

構造ファイルctrls.o2は、たとえばこんな感じで書く。

% const disexp=2.282 a=25.0 dis=0.0  dd=(disexp+dis)/a/sqrt(3.0)       
% const nk=1 da=0                                             
STRUC   ALAT={a} DALAT={da} PLAT= 0 .5 .5   .5 0 .5    .5 .5 0         
SITE    ATOM=O POS=-.5*{dd} -.5*{dd} -.5*{dd}                          
        ATOM=O POS=.5*{dd} .5*{dd} .5*{dd}                    

ここで、disexpなどは定数で、それを引用する形で、構造と酸素Oの位置が指定される。 (このdisを変化させながら計算を繰りかえすることができる。以下で説明する)。 ATOM=OのOが何であるかは、デフォルトで以下のctrlgen.pyに よって認識される。(ctrlgen --helpを参照。ctrlgen --helpatomnameを打ってみたらいいといわれる)。 この構造はfccでcell volumeは3906.25au^3で、一分子あたり、8angstrom^3程度である。

  •  %で始まる行によって変数を定義する。これらの変数は、プログラム起動時にプログラムの引数によって置き換えることが可能である。
  • ALATの単位はa.u.である。上の例にあるように、簡単な代数計算をconstのなかで行わせることができる。これが基本ユニットとなる。これを単位として構造を指定する。
  • PLATはprimitive unit cell. 3つずつに分かれている。ALATを単位とする。
  • ATOM=OのOは原子の種類を意味する。なにも指定しない時は、ctrlgen --helpatomnameで指定されるZをもつとされる。たとえばfractionalな値を入れるとかの場合、ATOM=MYATOM1などと指定できる。だがそのときには

 SPEC ATOM=MYATOM1 Z=8

      ATOM=MYATOM2 Z=8

などとしてSPEC(原子種)指定する必要あり。等価位置になくても区別する必要はない(ただ初期値の磁気モーメントが違うときなどはFe1 Fe2などとしてわける必要がある。その場合もいったん、以下のctrlgen.pyで生成してからそれをeditすることも考えられる。)。

ctrlgen.py o2を実行して、ctrl.o2を生成する。

ctrlgen.py o2

を実行する。するとMT半径を自動設定し(このばあいはtouching MTになる)、 ctrl.o2ができる。これがテンプレートになる。あるいは手でMT半径を 設定したい場合は、rmt.tmpファイルをeditして、

O 1.0

として、そのあと

ctrlgen.py o2 --readrmt

とする。この場合はrmt.tmpであたえたMT半径を採用したctrl.o2ができる。これは、

### This is generated by ctrlgen.py from ctrls ; ver.2.0 tkotani 13July2010
### 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=5  # kmxa=5 is good for pwemax=3 or lower.
               # larger kmxa is better but time-consuming. A rule of thumb: kmxa>pwemax in Ry.
% const disexp=2.282 a=25.0 dis=0.0  dd=(disexp+dis)/a/sqrt(3.0)
% const nk=1 da=0
STRUC   ALAT={a} DALAT={da} PLAT= 0 .5 .5   .5 0 .5    .5 .5 0
  NBAS= 2  NSPEC=1
SITE    ATOM=O POS=-.5*{dd} -.5*{dd} -.5*{dd} 
        ATOM=O POS=.5*{dd} .5*{dd} .5*{dd} 
SPEC
 ATOM= O Z= 8 R=1.0
     RSMH=  0.500 0.500 0.500 0.500 EH=  -0.5 -0.5 -0.5 -0.5
     RSMH2=  0.500 0.500 0.500 0.500 EH2= -2 -2 -2 -2
     KMXA={kmxa}  LMX=3 LMXA=4
     MMOM=0 0 0 0

% const pwemax=2 nk=2 nit=30 gmax=12
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.

      #Setting for molecules. No tetrahedron integration. (Smearing))
      #TETRA=0 
      #N=-1
      #W=0.001

ITER  CONV=1e-6 CONVC=1e-6 NIT={nit}
                # 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. (Real spece mesh for charge density).
                # Instead of GMAX, we can use FTMESH.
                # You need to use large enough GMAX to reproduce smooth density well.
                # Look into sugcut: shown at the top of console output. 
                # It shows required gmax for given tolelance HAM_TOL.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for VWN.
                # =2 Birth-Hedin (if this variable is not set).
		#    (subs/evxc.F had a problem when =2 if rho(up)=0 or rho(down)=0).
                # =103 PBE-GGA

      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.
      #STABILIZE=1e-10  # Default is 1e-8 for robust convergence. But 1e-10 can give slightly lower total energy.

      #ELIND=-1  # (unused now) this for accelaration of convergence. maybe not efficient. See mixrealsmooth() in switch.F
      
OPTIONS PFLOAT=1 # Q=band (this is quit switch if you like to add)
                 # 

となる。これはあくまで計算の入力ファイルのテンプレートである。原則的に、ここに書かれているパラメーターの意味は知っておく必要がある。これ以外にもオプションのパラメーターはあり、lmf --inputで簡単な説明が表示される。 従来、難しかったのはSPECセクションであったが、小谷らの提案する新方式では(論文はこれから)、

    RSMH=  0.500 0.500 0.500 0.500 EH=  -0.5 -0.5 -0.5 -0.5
    RSMH2=  0.500 0.500 0.500 0.500 EH2= -2 -2 -2 -2

は基本的に固定でかまわない。ちょうどRSMH=R/2に固定されている。(これ等の数字はs,p,d,f軌道についてならんでいる)。 それでs+p+d+fで16基底を二セットで計32基底取っていることに対応している。

ctrl.o2を修正する

O2の場合は、

  1. 磁性をもつのでNSPIN=2とする。
  2. k点の数は一個でいいので、nk=1とする。
  3. MMOM=0 0 0 0は原子計算で初期値をつくるときの磁気モーメント指定であり、0 2 0 0とする。(あくまで原子の計算に関する初期値設定、s,p,d,fの順)。これで初期値の酸素原子の電子数がp chanelで分極したものになる。lmfa o2 |grep confと打ってみればよい。
  4. 分子計算なのでinsulator計算でよいはず。ただ、レベルが計算途中でいれかわったりすることを考えるとsmearingをもちいたmetal計算の方がいい。このためにはTETRA=0以下の3行のコメントをはずす。すなわちBZセクションの
     #Setting for molecules. No tetrahedron integration. (Smearing))
     TETRA=0  
     N=-1
     W=0.001

lmfa,lmfを実行

順次実行する。

lmfa o
lmf o

なんらかの結果を得る。lmfaは原子計算でありすぐにおわる。これで計算した電子(スピン)密度を重ね合わせてlmf(バンド計算のメイン部分)の初期値とする。全エネルギーはsave.o2に書かれる。


原子間距離を変えながら実行

これはctrlファイルの冒頭に種得るスクリプトを書き込む方法で可能になる。 このためには、

#!/bin/tcsh                                                                                      
set nk = 1
set dis=0
lmfa o2 -vnk=$nk -vdis=$dis > out_lmfa 
foreach dis (-0.25 -0.20 -0.15 -0.1 -0.05 -0.025 0 0.025 0.05 0.1 0.15 0.2 0.25 0.30 0.35 )
lmf --rs=1,1,1,0,0 o2 -vnk=$nk -vdis=$dis >out_dis${dis}_nk${nk}
end            
exit 
     
ctrlstart ==============================                    

をctrlファイルの最初に書き込む。そうすると、このctrlファイル自身をシェルスクリプトとみなせるから, たとえばtcsh ctrl.genで実行ができる。あるいは、これをchomod +xとして実行形式にしてしまってもいい。 実際のctrlファイルの中身はctrlstart 以下のものである。ctrlファイルで定義された dis変数が置き換えられながら計算を順次実行していくことになる。

それで、結局、最終的に得られるファイルは以下のものであり、これをtcsh ctrl.o2で実行することになる。

#!/bin/tcsh
set nk = 1
set dis=0
lmfa o2 -vnk=$nk -vdis=$dis > out_lmfa
foreach dis (-0.25 -0.20 -0.15 -0.1 -0.05 -0.025 0 0.025 0.05 0.1 0.15 0.2 0.25 0.30 0.35 )
lmf --rs=1,1,1,0,0 o2 -vnk=$nk -vdis=$dis >out_dis${dis}_nk${nk}
end
exit

ctrlstart ==============================
### This is generated by ctrlgen.py from ctrls ; ver.2.0 tkotani 13July2010
### 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=5  # kmxa=5 is good for pwemax=3 or lower.
               # larger kmxa is better but time-consuming. A rule of thumb: kmxa>pwemax in Ry.
% const disexp=2.282 a=25.0 dis=0.0  dd=(disexp+dis)/a/sqrt(3.0)
% const nk=1 da=0
STRUC   ALAT={a} DALAT={da} PLAT= 0 .5 .5   .5 0 .5    .5 .5 0
  NBAS= 2  NSPEC=1
SITE    ATOM=O POS=-.5*{dd} -.5*{dd} -.5*{dd} 
        ATOM=O POS=.5*{dd} .5*{dd} .5*{dd} 
SPEC
 ATOM= O Z= 8 R=1.0
     RSMH=  0.500 0.500 0.500 0.500 EH=  -0.5 -0.5 -0.5 -0.5
     RSMH2=  0.500 0.500 0.500 0.500 EH2= -2 -2 -2 -2
     KMXA={kmxa}  LMX=3 LMXA=4
     MMOM=0 2 0 0



% const pwemax=2 nk=1 nit=30 gmax=12
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.

      #Setting for molecules. No tetrahedron integration. (Smearing))
      TETRA=0 
      N=-1
      W=0.001

ITER  CONV=1e-6 CONVC=1e-6 NIT={nit}
                # 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=2   # 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. (Real spece mesh for charge density).
      #FTMESH=55 # Instead of GMAX, we can use FTMESH.
                # You need to use large enough GMAX to reproduce smooth density well.
                # Look into sugcut: shown at the top of console output. 
                # It shows required gmax for given tolelance HAM_TOL.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for VWN.
                # =2 Birth-Hedin (if this variable is not set).
		#    (subs/evxc.F had a problem when =2 if rho(up)=0 or rho(down)=0).
                # =103 PBE-GGA

      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  # (unused now) this for accelaration of convergence. maybe not efficient. See mixrealsmooth() in switch.F
      
OPTIONS PFLOAT=1 # Q=band (this is quit switch if you like to add)
                 # 

収束チェック

この結果を基本として、収束チェックをしていく。

  • pwemax 平面波基底の数。原則的には、これを増やすことで収束チェックができる。kmxaは平面波をMTサイトで展開する時の多項式の次元だが、経験則としてはkmxa>pwemax+2ぐらい以上になっていれば問題ない。pwemaxを大きくしていくと、基底関数の線形独立性のなさが問題になる。

デフォルト設定で異常終了するときなどには小谷までご一報願いたいです。対応します。

  • コンソール出力でovlmatという文字列以下に基底関数のかさなり積分の固有値が表示される。これが非常に小さくなってくると、固有値が安定せず異常に小さなeigenvalueとかが表示され計算が破綻することがありえる。現在テスト中の手法は、新規開発中のautomatic stabilizerの手法である。で、もし収束が安定しない場合。*R. これはMT半径。Rを小さくしてみる。RSMHはR/2になるようにとる。ただし、あまりちいさすぎるのはこまるのでRを小さくした場合でも0.5を最低値とする。2原子分子の結合半径よりも小さく取ったりもできる。PMTの方法論ではMT半径が多少重なってても気にしなくてよい(現状のGWではまだその点で改良の余地あり)。
  • そもそも、pwemax=2でもそんなに悪い計算ではない。pwemax=4と比較してそんなにかわらないならよいとする。pwemax=6が計算が安定しない場合は、

HAM_STABILIZE=1e-6ぐらいにとるのもありえるかもしれない。

結果

上のctrlファイルでの基本結果は(pwemax=2)、commit: 17b7b17で

c nk=1 dis=-.25 mmom=2 ehf=-298.8136525 ehk=-298.8136524                                        
c nk=1 dis=-.2 mmom=2 ehf=-298.8404898 ehk=-298.8404897                                         
c nk=1 dis=-.15 mmom=2 ehf=-298.8593074 ehk=-298.8593074                                        
c nk=1 dis=-.1 mmom=2 ehf=-298.8714217 ehk=-298.8714217                                         
c nk=1 dis=-.05 mmom=2 ehf=-298.8779182 ehk=-298.8779182                                        
c nk=1 dis=-.025 mmom=2 ehf=-298.8793468 ehk=-298.8793468
c nk=1 dis=0 mmom=2 ehf=-298.8796959 ehk=-298.8796959
c nk=1 dis=.025 mmom=2 ehf=-298.8790559 ehk=-298.8790559
c nk=1 dis=.05 mmom=2 ehf=-298.8775109 ehk=-298.8775109                                         
c nk=1 dis=.1 mmom=2 ehf=-298.8720131 ehk=-298.872013                                           
c nk=1 dis=.15 mmom=2 ehf=-298.8637733 ehk=-298.8637732                                         
c nk=1 dis=.2 mmom=2 ehf=-298.853298 ehk=-298.8532977                                           
c nk=1 dis=.25 mmom=2 ehf=-298.8410331 ehk=-298.8410329                                         
c nk=1 dis=.3 mmom=2 ehf=-298.8273627 ehk=-298.8273626 
c nk=1 dis=.35 mmom=2 ehf=-298.8126029 ehk=-298.8126029

となる。対応する酸素原子の計算と比べると結合エネルギーは7.58eV =(-298.8796959+2*149.16113)*13.605程度である。 これはたまたま見つけた545030.pdf(グーグルする)のp.6の右側に、O2のbinding energyが、  LDA 7.53eV  GGA 6.31eV  expr. 5.2eV と書いてあるLDAの結果とよく一致している。GGAでも同等の一致をしめす(XCFUN=103とする)。

pwemax=4とすると、

c nk=1 dis=-.25 mmom=2 ehf=-298.8357858 ehk=-298.8357858                                        
c nk=1 dis=-.2 mmom=2 ehf=-298.8625796 ehk=-298.8625796                                         
c nk=1 dis=-.15 mmom=2 ehf=-298.8813641 ehk=-298.8813641                                        
c nk=1 dis=-.1 mmom=2 ehf=-298.8934565 ehk=-298.8934564                                         
c nk=1 dis=-.05 mmom=2 ehf=-298.8999465 ehk=-298.8999465                                        
c nk=1 dis=-.025 mmom=2 ehf=-298.9013794 ehk=-298.9013792                                       
c nk=1 dis=0 mmom=2 ehf=-298.9017383 ehk=-298.901738                                            
c nk=1 dis=.025 mmom=2 ehf=-298.9011135 ehk=-298.9011132                                        
c nk=1 dis=.05 mmom=2 ehf=-298.8995887 ehk=-298.8995887                                         
c nk=1 dis=.1 mmom=2 ehf=-298.8941449 ehk=-298.8941449                                          
c nk=1 dis=.15 mmom=2 ehf=-298.8859699 ehk=-298.8859699                                         
c nk=1 dis=.2 mmom=2 ehf=-298.8755592 ehk=-298.8755592                                          
c nk=1 dis=.25 mmom=2 ehf=-298.8633478 ehk=-298.8633477                                         
c nk=1 dis=.3 mmom=2 ehf=-298.84971 ehk=-298.8497096                                            
c nk=1 dis=.35 mmom=2 ehf=-298.8349587 ehk=-298.834958                

を得た。このときにも結合エネルギーは7.58eV程度 =(-298.9017+2*149.1723127)*13.605である。(絶対値はだいぶと違うが)

GGA計算

このときにはXCFUN=103とする。lmfaからスタートしないといけない。lmfaがコアの初期値を作るのだがそれがGGA計算されていないといけない。

注意事項

lmfはデフォルトでは、そのディレクトリに存在するrstファイルに書かれている原子位置posを読み込んでスタートする。 (lmf --help).もし位置をctrlファイルから読み込みたいなら、たとえば--rs=1,1,1,0,0などとする必要がある。

Personal tools