QandA

From Ecal

Jump to: navigation, search

Contents

SrTiO3---bulkで結果が変(解決済み)

ctrls.svo

%const  au=0.529177                       #<--   These % sections are to define variables.   
%const  a1= 20.0  da=0
% show vars                               
HEADER  TAstack
STRUC   
        ALAT=7.2644 
        PLAT=1 0 0 
             0 1 0 
             0 0 1 

SITE    ATOM=V  POS=0.0 0.0 0.0 
        ATOM=Sr POS=0.5 0.5 0.5 
        ATOM=O  POS=0.5 0.0 0.0 
        ATOM=O  POS=0.0 0.5 0.0 
        ATOM=O  POS=0.0 0.0 0.5 

syml.svo

21 .5 .5 .5    0 0 0 
21  0 0 0   .5 0 0 
21 .5 0 0    .5 .5 0
21 .5 .5 0   0 0 0
0 0 0 0   0 0 0 

フェルミレベル付近のdispersionが変(がたがた)になる。pwemaxを変更、overlapが小さくなるSrのs,p channelを除いてみるなどしましたが相変わらず変です。


REPLY: ぼくがこっちでやると正常にできました.

c ehf=-8696.073227 ehk=-8696.073238

job_band

#!/bin/bash
ddd=svo
lmdir=../..
rm -rf wkp.$ddd
$lmdir/lmf $ddd --quit=band    |tee llmf_ef
$lmdir/lmf --band:fn=syml $ddd 

echo -15,15,15,15|$lmdir/plbnds -fplot -ef=0 -scl=13.605 $ddd
$lmdir/fplot -f plot.plbnds
echo 'Check ps.dat(postscript)!'

実行すると、バンドはFile:Svosample.pdf。ctrl.svoはCtrlsvoSvo3_caseにもあります。

ーーー>これは解決しました。ファイルに一部古いものが混ざっていたせいでした。--Kino 16:50, 19 October 2009 (JST)





ehfとehkが一致しない(なおしちゅう)

ひとつには、mixrhoのinterpolationの際(bndfp.Fがよぶ)、

mixrho: (warning) scr. and lin-mixed densities had 5520 and 5520 negative point

が表示されるとそういうことがおこるようです。bndfpのiterationの最後で

       call mixrho(ssite,sspec,slat,nsp,iter,sstrn(i1:i2),qval-qbg,
    .    elind,orhat1,w(oorhat),k1,k2,k3,dmxp,w(osrout),w(osmrho))

がコールされてますが、ここでw(osrout)があたらしく作られた電子密度のsmooth part(実空間メッシュの上での密度)、w(osmrho)が、mixの後のsmooth partの電子密度です。mixingのためには全回のiterationの電子密度な どもいりますがそれはmixm.*ファイルにあり、それをこのルーチンがアクセスしてます。で、上のメッセージはこのmixrhoがだしています。調べると「w(osrout)がただしく正定値なときにでも、mixした結果が負になりそれからむりに負の部分をカットしたものがw(osmrho)になる」場合に出ます.たちがわるいのは、それがおこるとehkとeksがずれたまま収束するということがおこります。

bug in lmchk

       if (alat .eq. 1) write(stdo,345) ib,ntab(ib+1)-ntab(ib),rcut
 345 format(' pairc, ib=',i3,':',i4,' neighbors in range',2f7.3)     <--- not f7.3 but 2f7.3

When did I include this bug? It was working... --Tkotani2 03:06, 26 September 2009 (JST)

input statement requires too much data, file=...rst.*

Q. intel fortanで実行すると

forrtl: severe (67): input statement requires too much data, unit 12, file  /home/kino/.../rst.svo

というrstファイルの読み込みに関するエラーが出る場合があります。

A. intel fortranのバグの様です。lapack,blasともにgfortranなど別のコンパイラでコンパイルし直してください。


bug for AnyQ mode

search anyq in fpgw/gwsrc/mkqg.F At L387 or around, there is a block "if(anyq) then...". At the end of this if block, you need to release qsave as

       allocate(q0i(3,nq0i))
       q0i=qsave
       deallocate(qsave) <--- this shuold be inserted.
       close(ifqpnt)
     endif

--Tkotani2 13:50, 25 September 2009 (JST)


zgemm (called from matm or so)

Be careful to use acml, atlas, or intel math libraries. zgemm can give wrong results for large size of array (when size is more than several hundred or so). --Tkotani2 12:22, 25 September 2009 (JST)

removed gwd/sugwin.F

We removed gwd/sugwin.F. It is not nessesary. --Tkotani2 12:22, 25 September 2009 (JST)



酸素原子:sigletとtriplet(解決)

(ここで言うsingletはスピンがない場合(縮退の場合)であり、Antiferroにオーダーしたものでない。あくまでバーチャルなもの)



酸素原子で2pがスピンが縮退した場合の解、その一般化について(物理的描像)

スピンが縮退した場合(up downを分けて考えない場合)、 2px,2py,2pzに4個の電子を詰めることを考えればよい。スピンがない場合には、 2px,2pyに二個づつ詰め、2pzは空にする、という詰めかたが考えられる。しかし、 このとき、クーロン反発力のため、HOMOである2px,2pyは2pz軌道より高いエネルギー準位となる。それで、HOMOより低い2pzに穴をあけた状態でないと解が求まらないことになる。

別の考え方として,fractionalな占有数を許す方法(有限温度で解くなど)がある。この時には2px,2py,2pzに4/3個ずつ積めていくやり方が考えられる。しかしこの解は、電子移動に対して全エネルギーは不安定である。2px,2pyから2pzに電子を移したほうが全エネルギーが下がるのである。ただし、このとき2pz軌道のほうがエネルギー準位はさがってしまう(結局HOMO軌道の下に穴があく)。 この場合、有限温度で解く方法ではある温度以下では「解なし」となる。

一般的に言って、運動エネルギー+静電エネルギーを最小にする問題を微分=0で解こうとする場合、この問題が発生しうる。すなわち運動エネルギーのゲインを犠牲にしてでも(HOMOの下に穴が開いた状態ででも)、クーロンエネルギーを下げるような電子配置を選んだほうが全エネルギーが得になるような状況である。フレンケルエキシトンと磁気的励起が、同一のものになる。逆に言えば、こういうものをコアにすれば、磁気ー光ー電気伝導の強い結合をつくることができる。(調べること)。

  • 「HOMOの下に穴があいたような状況」は状況によっては、たとえ基底状態であっても存在しうるし、格子変形によって実現されたりしうる。
  • 誘電応答はどうなるか?線型応答のレベルでstaticには正定値であるはず。RPAは有効なのか?
  • また一方で、2原子分子などでは、HOMOがfractionalな占有数を持つ状態が電子移動に関して安定であったりもするようだ。
  • 占有数をコントロールする。ある極小解Aがあるとする。占有数のconstraintを連続的に変化させていけば、連続的に変形する解が得られるだろう。そして、他の極小解Bへ行き着くだろう。ただ、占有数constraintの変形経路が違うと別の極小解B'に達するかもしれない。

Gaussianの結果

BE= binding energy of O2
dE=E(triplet)-E(singlet)

O pbepbe/6-311+G(d,p)
2Sz+1  E_tot
1 -74.8975105519 au
3 -75.0033609885 au   dE= .1074765926 au = 2.9246 eV

O2 pbepbe/6-311+G(d,p)
2Sz+1  E_tot
1 -150.169318957  re=0.610476*2 ang
3 -150.229513783  re=0.609845*2 ang  BE=.2227918060  = 6.0625 eV

O LSDA/6-311+G(d,p) 
2Sz+1  E_tot
1 -74.5704620153
3 -74.6737197564  dE=.1032577411 au =  2.8098 eV

O2 LSDA/6-311+G(d,p) 
2Sz+1  E_tot
1  SCF NG
3 -149.624418359 re=0.601990*2 ang BE=.2769788462 = 7.5370 eV

1au=27.2116eVとして変換。singletの計算はSCF(tight)では標準settingでは収束していないだけ。

lmfの結果

  • 球対称化した電子密度でのO原子の計算。

singlet:-149.0574 triplet:-149.1664 差は、0.11Ry程度。

  • triplet:以下のctrlで、NSPIN=2とすればtripletは計算できる。
c mmom=2 ehf=-149.1698502 ehk=-149.1698584

バンドギャップは0.15eVと小さい。METAL=0(insulator)だと収束が安定しない.METAL=3だと よい(why?;たぶん最初のiterationで縮退がただしくほどけるから)。 が、収束安定性はもうすこし上げたい感じ(エネルギは1e-5Ryまで収束するがデフォルト収束条件1e-6Ryに達しないで振動する)だが基本的には問題がないと思われる.

  • singlet:穴をあけて収束させる方法でないと収束しない.

「酸素singletはスレーターの遷移状態計算のようなことを考えないとできない。穴を開けて計算するってこと です。」 というものです。lmfは基本的には(数値的にもっと安定させたいし、小さいMT半径の問題や 他の分子テストはこれからにせよ)、それなりにちゃんと動いてるように思えます。

(1)まず、r4zで収束させる。ただし、たとえばzhev.Fで、
   3つの(P状態の)準位のうち、いちばん深いものを空にするようにする。
   たとえば、以下の文をzhev_tkでreturnする直前にいれました.
   (insulatorだがmetal=3で回す;unoccupied(4番め)を計算させるため)。

      if(abs(e(3)-e(4))<0.001) then
         eee= e(2)
         e(2)= e(3)
         e(3)= e(4)
         e(4)= eee+1d0
         allocate(zzz(n))
         zzz   = z(:,2)
         z(:,2)= z(:,3)
         z(:,3)= z(:,4)
         z(:,4)= zzz
         deallocate(zzz)
      endif

   (エネルギー分散はあまりないので、nk=2でもe(3)-e(4)で判定できた)。

(2)そのあと、SYMOPS Eにする。
 そうすると、c ehf=-148.9604801 ehk=-148.960482を得ます(nk=2)
            c ehf=-148.9603921 ehk=-148.960394(nk=1)
 となる。
 固有値は -1.6051 -0.6277 -0.4892 -0.4889 ...(nk=1)となる。
  占有数は  2        0      2       2     です。
  すなわち-0.6277の波動関数を空にしてます。
  (-0.4892と-0.4889はおそらく本来は縮退してるだろうとおもわれる状態ですがわずか
 に数値誤差があるようです)。

(3)で、そのような占有数でのtripletでの計算は、
  c mmom=2 ehf=-149.1700881 ehk=-149.1700914 nk=2
  c mmom=2 ehf=-149.1698502 ehk=-149.1698584 nk=1
 (すこしだけ数値不安定で、CONV=1e-5にしてやったが、.17009まで収束してる)
 となります。これで、singlet,tripletの差は2.85eVとなります。

用いたctrl.o

### 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=42
             # 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 E #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  au=0.529177                       #<--   These % sections are to define variables.   
%const  a1= 9.0/au  da=0
%const distance=1.2
STRUC   ALAT=1
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}-4  # PLAT: primitive vector in ALAT. ALAT  (in  a.u.)  
                                              # DALAT gives "shift" of ALAT  (in  a.u.)  
  NBAS= 1  NSPEC=1
SITE						
  ATOM=O POS=0 0 0
SPEC
 ATOM= O Z= 8 R=2
      RSMH= 0.916 0.930  EH= -1.367 -0.434
      KMXA={kmxa} LMXA=4
      MMOM=0 2 0

% const pwemax=2 nk=1
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.
      TETRA=0 
      N=-1
      W=0.001
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)
                 # 


  • 注意:smearingだと(BZ TETRA=0だと)Wの値はrstファイルから読まれるーー>新バージョンではctrlから読まれるようにすること. lmfp.F L506あたり。


frozen core近似のせいかもしれません。--Kino 19:56, 13 September 2010 (JST)

計算の限界は?

  1. intel corei7(nehalem). 8CPU, 32core上でlmf-MPI 30 MPIを使って20.0Ang^3のsupercellでの30原子(thymime-adenosine,パラメタはctrlgen.pyで作成)を含む系が24時間で一度もiterationが回りません。なおこの場合の必要メモリーはwksize=800 000 000です。)non MPIのlmfなら3回iteration程度まわります。


ctrlgenは大文字のファイル名を受け付けるがlmfaは小文字だけ

整合性がないです。lmfは基本はファイル名などはすべて小文字にするように なってますがときどき(小谷の加えたものなど)大文字の場合もあります。 とにかく、ctrlファイルは小文字オンリーであることを明確にすべきですね(将来なおす)。


BZ_ZVAL Number of electrons to accumulate in BZ integration のdefault=0となってますが、正しいのでしょうか。

たぶんこのZVALはlmfモードではつかわれてないとようです(要チェック)。これを消すには m_rdctrl.FのLMF switchなどで、call tkadd(" BZ_ZVAL~")をコメントアウトすればいいだけ。 ですが、当面ほっときます(確認後消すこと)。

おそらくBZ_ZBAK(back ground charge)のみが生きてます。

fp/lmfp.F:233:      qbg    = dgets('ctrl zbak',sctrl)

です。要チェック。

(以下メモ:要チェック:lmf --inputで必要な入力とオプショナル入力が示されてます。m_rdctrl.F内でgtvを読んでるところで読み込まれています。lmfなどに関連した、tkadd("BZ_ZVAL~",...)などとなってるのでオプショナルでデフォルトはたしかにゼロに設定されてます。この値は、rdctrl2.F内でbz_zvalとして引用されbzのzvalにパックされる。ですがbzwtsf (bndfp.Fからコールされる)がフェルミエネルギーを決めるがこのzvalはzの和とBack ground chargeで決まるんだと思う。)

補足:BZ_ZBAK=1が(system) + 1、BZ_ZBAK=-1が(system) − 1に対応。


SrTiO3_tutorialのサンプルにかんして.バンド図から見てO2s,Sr4pを越えない内核準位は凍結してあるようですがこの凍結範囲はどこで制御するのか?例えばコアレベルシフトの計算をどのようにするのか?

SrTiO3_tutorialのサンプルでは、ctrls.srtio3から生成されたctrl.srtio3でのSrの指定は

SPEC
ATOM= Sr Z= 38 R=3.616323
     RSMH= 2.411 2.411 2.094  EH= -0.100 -0.100 -0.100  PZ=0 14.9 P=0 5.3
     KMXA={kmxa} LMXA=4
    MMOM=0 0 0 0

となってます。このctrlファイルをつかって計算は進行するわけです。まず、 lmfa srtio3 |grep conf してもらうと、

conf:SPEC_ATOM= Sr : --- Table for atomic configuration ---
conf int(P)z = int(P) where P is replaced by PZ if it is semicore
conf:   isp l    int(P) int(P)z    Qval    Qcore   CoreConf
conf:    1  0       5   5        2.000    8.000 => 1,2,3,4,
conf:    1  1       5   4        6.000   12.000 => 2,3,
conf:    1  2       4   4        0.000   10.000 => 3,
conf:    1  3       4   4        0.000    0.000 =>
conf:    1  4       5   5        0.000    0.000 =>
conf:-----------------------------------------------------
           (以下略)

となります。この表が電子状態計算でつかわれる原子のconfigurationです。 int(P)がvalenceの主量子数で、5s,5p,4d,4f,5gとなっています。 Qval、Qcoreが初期の原子計算での各チャンネルでの電子数です。総和は38になってます。 (この初期計算での電子数も変えれますがそれにはQ=というオプションをつかいます。  lmf --inputで、ctrlファイルで有効なオプションがみれます)。

pチャンネル (l=1)のint(P)zだけint(P)とずれていますが、 これはlocal orbitalをそのチャンネルで使ってることを意味してます。 要するに4pに対応するMTO基底をいれています (これがPZ=0 14.9の意味です。最初のゼロはsチャンネルではつかわないという意味)。

もしPZ=0 4.9とするとMuffin-Tin内で閉じた基底(これが通常の意味でのlocal orbital) を入れれます。これでも、場合によるけど、そんなには悪くはないです。  (単純には最初からこの状態で、lmfa lmfをやりなおすと  全エネルギーがgrep 'c ' save.srtio3すると、

c ehf=-8505.3551189 ehk=-8505.3551211 <--- PZ=0 14.9 のとき
c ehf=-8505.3544806 ehk=-8505.3544831  <--- PZ=0 4.9  のとき

 でわずかにPZ=14.9のほうが良さげです。ただ、pwemax=3をもっとおおきくすると  差は小さくなるでしょう。)

で、それらのPZに対応するチャンネルは対角化にいれてるので、valenceの電子と 同等に計算してるわけです。そういうやりかたが時間はかかるけど浅い側のsemi coreレベル を正確に解くにはいい方法です。lmfでは、各lごとにPZに対応する軌道を一個ずつ入れれます。

たとえば、Sr 4sもvalenceにするには、PZ=14.9 14.9 P=5.3 5.3とでも してください。Pは"continuously variable" principal quantum Pと呼んでいます。  http://titus.phy.qub.ac.uk/Programs/LMTO/lmto.html#section2 これは原子のサイトでのradial schoredinger eqを解くときの対数微分値 (どのエネルギーでradial関数を解くか?といってもおなじこと)を決めています。 通常の計算では(IDMOD=というのを指定しなければ)、 「深い方の軌道のPはself-consistentに自動決定される。また浅い方は固定」ということになってます。 ですので、PZ=14.9 P=5.3だと、「4sに対応するMTOのPは自動決定(これがPZと呼ばれてます)。 Pの方は固定」ということになってます。この場合radial関数としては、 「PZの値で解いたもの一個+Pの値でといたもの+Pの値で解いたもののエネルギー微分」 の3つのものがつかわれます。P=5.3をP=5.4にしたら多少は最終結果に影響がでます( これはいまのコードでは避けれないですがわずかです)。PZのほうは初期値なので影響しませんが。

えっと、説明になってるでしょうか?もっときちんと説明すべきなのですが (それと改善の余地もある部分ですし)、むしろ下の方がポイントなのかもしれない。

とにかくこのコンソール出力の説明をきっちりつけないといけないのですが、 現状では断片的なドキュメントしかないです。。 http://pmt.sakura.ne.jp/wiki/index.php?title=EcalJ にすこしずつかいていこうとしてますが、 ここに引用してる基礎文献をざっとみてもらうのがまずはいいかとおもいます。

4

もっと深いコアのレベルはどうか? とにかく「標準(デフォルト)」では深いコアはfrozen coreとしてますので、 計算していません。それをするには、 http://titus.phy.qub.ac.uk/Programs/LMTO/fp.html にかいてあるLFOCA=をつかいます。LFOCA=0とするとcoreの波動関数を MTの端でphi=phidot=0となるようにして計算します。各原子のSPECごとに指定できます。 たとえば、

ATOM= Sr Z= 38 R=3.616323
     RSMH= 2.411 2.411 2.094  EH= -0.100 -0.100 -0.100  PZ=0 14.9 P=0 5.3
     KMXA={kmxa} LMXA=4 
    MMOM=0 0 0 0 LFOCA=0

とします。(最初に表示されるSPEC_LFOCAをみてきちんとよみこまれたか確認してください)

そうするとsc計算内でcoreも計算されます。これはlmv7-lmfp-bndfp-mkrout-getcor.F L250の

C   --- Make new core density and core eigenvalue sum ---
       if (lfoc .eq. 0) then
         call pshpr(ipr-11)
         call getcor(0,z,a,pnu,pnz,nr,lmxa,w(orofi),w(ov01),kcor,lcor,
    .      qcor,smec,smtc,w(orhoat(3,ib)),ncore,0d0,0d0)

あたりで計算されてますが、(lfocがLFOCA)、このecoreをみるには

IO      SHOW=T VERBOS=60

とVERBOSEをあげてみてください。すると、

state  chg          ecor0          ecore          tcore    nre  rho(rmax)
1s    2.00   -1168.670779   -1168.670779    1487.767467    391  0.00000
2s    2.00    -157.828824    -157.828824     301.089848    449  0.00000
3s    2.00     -23.949615     -23.949615      75.150042    451  0.00000
4s    2.00      -2.523191      -2.521687      14.159545    451  0.00004
2p    6.00    -140.279739    -140.279739     287.781057    445  0.00000
3p    6.00     -18.292372     -18.292372      68.362332    451  0.00000
3d   10.00      -8.926275      -8.926275      57.768402    451  0.00000

という表がでてきます。このecoreがコアレベルです。ecor0とecorがもしずれてるなら それはcoreが十分にMT内に局在してないということです。単位はRyでvalenceのレベルと同じくRyです。 これは、getcorから呼ばれるrhocorで書かせてます(atomsr.F内)。

またcore-level spectroscopyというオプションもあるのですが(EELS, --cls)つかったことないので、 よくわかってません。fp.html によると、 suclst.Fをみろとかいうことですが、make checkのサンプルでもfe,crnが対応してるようです。 fp/testで./test.fp crnすると実行されるのがみれます。ただこのサンプルはあんまりきれいじゃなくて みてたら何がなんだかわからなくていやになるかもしれないです。時間もかかる。 いるなら、ぼくが先にクリーンアップしてわかりやすくします。いずれやらないといけないし。


5

とにかくコアレベル表示させるひとつの方法はLFOCA=0です。mkrout.F L248を

call pshpr(ipr+11)

とプラスに直してもらうと、verbose=35でもcoreレベルが表示されます。

実際、まちがいなのかもしれない(これはpush popでgetprでとれるiprintの値の上げ下げを してるんです。それをみてどれだけ出力するかきめてます)。どっちにせよ、表示がかわるだけです。

また、収束計算させた後、LFOCA=0をくわえて、non-self-consitentに計算 させてもいいかもしれないです。そのためにはOPTIONS HF=tとするかNIT=1として (これでiterationは一回で終わる)lmf --rs=1,0 srtio3で実行します。 これで、電子密度のrstファイルは上書きされないです。lmf --help参照。

Personal tools