Ecalquestion

From Ecal

Jump to: navigation, search

Contents

ctrlgen.pyが生成するctrlにおいてsemi-coreにMTOを入れるかどうかどこで決めてるか?

subs/writebasis.FでPZセッティングを決め、mtopara.*に書き出させてます。 subs/writebasis.Fの

c d channel local orbital
      if(29.5d0<zz .and. zz<32.5d0) then !3d
        szz='  PZ=0,0,13.9 P=0,0,4.2'

のところです。writebasisはsubs/freeat.Fから呼んでます。 これをctrlgen.pyで読み込んでそのまま書かせてます。

lmf-MPIでのモーメントファイル

mpirun -n 2 lmf srtio3 を走らせた後、momentファイルができないのでlmdosでDOS書けない。現状ではlmfを再度走らせないといけない。

subs/swich.Fのfullmesh=T

空間メッシュを適当にカットするのを止める方法ーーー>正定値がきちんと確保できるはず。GGAの精度などに有利? まだなにかおかしい場合がある。たぶんgvlist.Fのshortbzmはkotaniのshortn3のアルゴリズムにおきかえないといけない。


原子の計算が収束しない

  • lmfa計算では、たとえばLa原子で、以下のような

設定で計算できる。(MOMとQに注目。このときには、fに3個の電子をfullに分極した状態を仮定して球対称で原子計算させている)計算できる。MOMがモーメントで、Qがcharge.lmfaのコンソール出力で、Pl,Qlを検索。NSPIN=2の場合でないとMOMには意味がない。(この場合、ちょっとワーニングがでておかしいが)。

SPEC
 ATOM= La Z= 57 R=2.5
      RSMH= 1.667 1.667 1.552  EH= -0.100 -0.100 -0.100
      KMXA=12 LMXA=4
      MMOM=0 0 0 1 0
      Q=   2 0 0 1 0 
  • 原子計算は、かならずしもinsulatorで収束しない。ある状態が占有されるとクーロン相互作用のせいでその状態のレベルが上がる。それで、次のiterationでは、空になってしまうということがおこる。それで

収束サイクルにおいてシーソー現象がおきる。酸素原子でもsiglet状態だとそういうことがおきた。 (どういう場合におこるか調べ上げる必要あり;文献さがす)。落ち着きどころは、ちょうどバランスした非整数個の占有状態ということになる。これは、励起状態などでもおなじ。と、いうわけで固体のコードではちょっとした工夫がないと異方性をいれた原子計算がきっちりできない。

  • MT半径が小さすぎる(なんとかコードを改善したい。とくにノードカウントで落ちる部分)。たとえ収束しても平面波数がたくさん必要だがなんとかならないか?
  • 初期条件、QとMOMをいれてからlmfaを動かしてセットアップ
  • Metalスイッチ。smearingが有効かも(BZ_METAL=3,BZ_TETRA=Fにする).smearingの幅をW=0.001とかW=.005
  • MTが小さいとき、EFMAX、NEVMXを大きくとらないと”それらをおおきくせよ”というメッセージを出して終了する。
  • 対称性。歪んだ解を求めるには,"SYMOPS E".
  • HAM_NSPIN=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.
      TETRA=F
      EFMAX=6
      NEVMX=15
      W=0.001

SrVO3の場合

svo3 case とくにSrのRSMHを小さめにとるべき。2/3*Rmtルール.

RSMHの設定

ctrlgen.pyでは内部で、lmfa --getwsrをよび(ctrlgen.pyの内部を見てください)、 RSMHなどを生成している。この過程で、RSMHに<2/3*RMTという拘束をかけるようにしている がこれでいいのかどうか。。。(さらには2/3*RMTが1より小さいときには、RSMH<1となるようにしてる:9月26日にから試してる;これでよいのかどうか疑問。いろいろためしていかないと。)。 (freeat.F L850あたり.しかしあまり論拠がないのと、L757でのrlim1の拘束との矛盾(マークのオリジナルでは rlim1>0.3.いろいろためしてるが現状ではマークのオリジナル(ただしdではrmt*2/3に戻してる)。プログラムの一部は

C   ... Possibly constrain rsm
ctakao
        rmtmx= rmt*2d0/3d0
        if(rmtmx<.5d0) rmtmx=.5d0
        if(l>=2) rmtmx=rmt

        if (mod(isw,10) .eq. 1 .and. rsm .gt. rmtmx) then
        if (ipr .ge. 20)
     .  write(stdo,
     .  '('' ...rsm exceeded rmtmx.. repeat with rsm= rmtmx ='',f8.5)') rmtmx
        rsm = rmtmx

--- 追記:現在(Aug2010以後)は、上述と違う考え方の計算法を研究している。あたらしいctrlgen.pyを見てもらうと分かるが(まだ試験中だが)、RSMH,RSMH2と二つの基底を各lに入れ、RSMH=1/2*RMTととっている。これにより計算が安定する。このときは、MT外はsmooth HankelはほぼHankelと一致している。それで、smooth Hankelの役割は概念的に以前のものと違ってきている(いまや単にsmoothなcharge densityをつくるためのものという役割しかない)。

場合によっていろんな収束があるようにみえる

ひとつの原因にwarningで、

(warning) failed to find proper node count for l=2 ebar=-0.3216: pnu not

などとなる場合。要注意です.あたえられたenuでradial関数を計算したときに ノードの数が一致しなくなっている。


追記:(以下の記述は正しいかどうか調べないといけない)。 おそらくノード数にはこだわらなくていいと思うから、このwarningは気にしなくても いいとおもう。ただ、PZがはいるとき、ちょっと問題がありそう。これがその原因かもしれない。 『semi-coreにlocal orbitalを用いるとき:ebarが深すぎるようにみえる』 (ebar=hbyl/qbylでpnunew.Fで計算されているが、local orbitalのDOSとvalence channelのDOSの分離がうまくてきていない)。

iodine case(長柄)

%const da=0 au=0.529177                                                  
%const d0=2.11/au a0=d0 v0=a0^3 vfrac=1 v=v0*vfrac a1=v^(1/3)            
% show vars                                                              
HEADER  I Cmca                                                           
STRUC   ALAT=13.485085627954255                                          
        PLAT= 1.0000000 0. 0.   0. 0.6566704 0.    0. 0. 1.3710762       
        DALAT=0                                                          
                                                                         
SITE    ATOM=I POS=       0.0000000      0.1026299      0.1605894        
        ATOM=I POS=       0.0000000     -0.1026299     -0.1605894        
        ATOM=I POS=       0.0000000      0.4309651      0.5249487        
        ATOM=I POS=       0.0000000      0.2257053      0.8461276

これからctrlgen.py(Aug2010 version)でctrlを作って計算できるが、PZ=15.9を入れて計算すると、6sのebarが-0.718となる。これはほとんど5s軌道のebar=-0.80と同じ。なにか変。 (一方で、PZなしだと、6sのebarが-0.22程度だからそれぐらいになるべき…)。

PZをもちいるときの改良予定

PZももちいるときには、P(pnu)の値を現状では固定してますが、このあたり改良の余地あり。 上述のポイントと関係してます。

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がずれたまま収束するということがおこります。


memo

TETRA=0のとき、2.5Ryより上の固有値計算してない?****になってる(たぶん問題ない)。


RSEQ: bad nodesで終了

RSEQ : nit gt 80 and bad nodes for l=1.  Sought 0 but found 1.  e=-1.3247D+00
Exit -1 RSEQ: bad nodes

で終了。

### 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 op
tion).
             #
             # 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  au=0.529177                       #<--   These % sections are to define variables.
%const  a1= 15.0  da=0
% show vars
HEADER  H2O
STRUC   ALAT=1/{au}
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}   # PLAT: primitive vector in ALAT. ALAT  (in
 a.u.)
                                              # DALAT gives "shift" of ALAT  (in  a.u.)
  NBAS= 13  NSPEC=3
SITE
 ATOM= C POS= -.330862 -1.356672 -.474487
 ATOM= C POS= -1.349851 -.399686 -.445168
 ATOM= C POS= -1.030612 .961458 -.440462
 ATOM= C POS= .307703 1.365816 -.464868
 ATOM= C POS= 1.326896 .409021 -.494189
 ATOM= C POS= 1.007555 -.952221 -.498452
 ATOM= H POS= -.581419 -2.423371 -.505843
 ATOM= H POS= -2.399256 -.716233 -.454166
 ATOM= H POS= -1.829831 1.711647 -.445891
 ATOM= H POS= .557003 2.432982 -.488946
 ATOM= H POS= 2.375146 .726206 -.540558
 ATOM= H POS= 1.805705 -1.701941 -.548427
 ATOM= Na POS= .044334 -.017780 1.808144
SPEC
 ATOM= C Z= 6 R=1.312815
      RSMH= 0.875 0.875  EH= -0.395 -0.114
      KMXA={kmxa} LMXA=3
     MMOM=0 0 0 0

 ATOM= H Z= 1 R=0.758599
      RSMH= 0.506  EH= -0.111
      KMXA={kmxa} LMXA=2
     MMOM=0 0 0 0

 ATOM= Na Z= 11 R=3.730086
      RSMH= 2.487 2.487  EH= -0.100 -0.100  PZ=0 12.9 P=0 3.3
      KMXA={kmxa} LMXA=4
      MMOM=0 0 0


% const pwemax=3 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=20
01 (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.
      ZBAK=1

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 magne
tic polarization).
      FORCES=1  # 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 Vosko-Ceperly-Alder; 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)
                 #
DYN     MSTAT[MODE=5 HESS=t XTOL=.001 GTOL=0 STEP=.015]  NIT=30
### 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  au=0.529177                       #<--   These % sections are to define variables.   
%const  a1= 10.0  da=0
% show vars                               
HEADER  O
STRUC   ALAT=1/{au}
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}   # PLAT: primitive vector in ALAT. ALAT  (in  a.u.)  
                                              # DALAT gives "shift" of ALAT  (in  a.u.)  
  NBAS= 1  NSPEC=1
SITE                                            
  ATOM=Fe POS=0 0 0
SPEC
 ATOM= Fe Z= 26 R=3.0
      RSMH=  2.000   EH=  -0.100 
      KMXA={kmxa} LMXA=4
      MMOM=0 0 4


% 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=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=9    # this is for real space mesh. See GetStarted.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for Vosko-Ceperly-Alder; 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)

計算の限界は?

  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程度まわります。

RSMH,EHの値が全く同じ

ctrlgen.pyで

 RSMH= 0.667 0.667  EH= -2.500 -2.500

と同じ値になるがいいんでしょうか。


ーー>これになる場合の入力をpreではりつけてください。見てみます。

たぶん、EH=-2.5、RSMH= 0.667はプログラムでの限度値になってます。あまりよろしくない。

### 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 ou
tput)
             # It is useful to check ctrl is read in correctly or not (equivalent with --s
how 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  au=0.529177                       #<--   These % sections are to define variables.

%const  a1= 10.0  da=0
% show vars
HEADER  O
STRUC   ALAT=1/{au}
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}   # PLAT: primitive vector in ALAT. ALAT
  (in  a.u.)
                                              # DALAT gives "shift" of ALAT  (in  a.u.)
  NBAS= 1  NSPEC=1
SITE
  ATOM=Al POS=0 0 0
SPEC
 ATOM= Al Z= 13 R=1.0
# RSMH= 0.667 0.667  EH= -2.500 -2.500
      RSMH= 0.667  EH= -2.500
      RSMH= 0.667  EH= -2.500
      KMXA={kmxa} LMXA=4
      MMOM=0 1 0


% const pwemax=3 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) N
PTS=2001 (division) even for insulator.
                #   (SAVDOS, DOS, NPTS gives no side effect for self-consitency calculaito
n).
                #
                #BUG: For a hydrogen in a large cell, I(takao) found that METAL=0 for
                #(NSPIN=2 MMOM=0 1 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.

     EFMAX=5

ITER  CONV=1e-6 CONVC=1e-6 NIT=60 MIX=A2,b=.2,n=3
                # 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=9    # this is for real space mesh. See GetStarted.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for Vosko-Ceperly-Alder; 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)

他にもRSMH, EHが同じになってしまう例

### 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 ou
tput)
             # It is useful to check ctrl is read in correctly or not (equivalent with --s
how 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  au=0.529177                       #<--   These % sections are to define variables.

%const  a1= 12.0  da=0
% show vars
HEADER  O
STRUC   ALAT=1/{au}
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}   # PLAT: primitive vector in ALAT. ALAT
  (in  a.u.)
                                              # DALAT gives "shift" of ALAT  (in  a.u.)
  NBAS= 1  NSPEC=1
SITE
  ATOM=Fe POS=0 0 0
SPEC
 ATOM= Fe Z= 26 R=2.0
      RSMH= 1.333 1.333 0.900  EH= -0.100 -0.100 -0.179
      KMXA={kmxa} LMXA=4
      MMOM=0 0 4


% const pwemax=3 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) N
PTS=2001 (division) even for insulator.
                #   (SAVDOS, DOS, NPTS gives no side effect for self-consitency calculaito
n).
                #
                #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 MIX=A2,b=.5,n=3

                # 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=9    # this is for real space mesh. See GetStarted.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for Vosko-Ceperly-Alder; 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)


zhev_tk: ovlmat= の値がNaNになる

少し前の

RSEQ (warning) eval for l=? node=? did not converge: de=??? e= ???

という部分が収束していません。ではRSEQをなんとかするには?

answer:まだなし。 ーーー>これはバグです。これも入力サンプルがほしいです。

### 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 ou
tput)
             # It is useful to check ctrl is read in correctly or not (equivalent with --s
how 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  au=0.529177                       #<--   These % sections are to define variables.

%const  a1= 12.0  da=0
% show vars
HEADER  O
STRUC   ALAT=1/{au}
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}   # PLAT: primitive vector in ALAT. ALAT
  (in  a.u.)
                                              # DALAT gives "shift" of ALAT  (in  a.u.)
  NBAS= 1  NSPEC=1
SITE
  ATOM=La POS=0 0 0
SPEC
 ATOM= La Z= 57 R=1.5
      RSMH= 1.000  EH= -0.100
      KMXA={kmxa} LMXA=4
      MMOM=0 0 2


% const pwemax=3 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) N
PTS=2001 (division) even for insulator.
                #   (SAVDOS, DOS, NPTS gives no side effect for self-consitency calculaito
n).
                #
                #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=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=9    # this is for real space mesh. See GetStarted.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for Vosko-Ceperly-Alder; 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)
                 #

収束しにくい場合のおすすめのパラメタは

例えば、MAX=A2,b=.5,n=8

収束しない例

### 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 ou
tput)
             # It is useful to check ctrl is read in correctly or not (equivalent with --s
how 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  au=0.529177                       #<--   These % sections are to define variables.

%const  a1= 12.0  da=0
% show vars
HEADER  O
STRUC   ALAT=1/{au}
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}   # PLAT: primitive vector in ALAT. ALAT
  (in  a.u.)
                                              # DALAT gives "shift" of ALAT  (in  a.u.)
  NBAS= 1  NSPEC=1
SITE
  ATOM=Fe POS=0 0 0
SPEC
 ATOM= Fe Z= 26 R=2.5
      RSMH= 1.667 1.667 1.023  EH= -0.150 -0.100 -0.233
      KMXA={kmxa} LMXA=4
      MMOM=0 0 4


% const pwemax=3 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) N
PTS=2001 (division) even for insulator.
                #   (SAVDOS, DOS, NPTS gives no side effect for self-consitency calculaito
n).
                #
                #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 MIX=A2,b=.5,n=3

                # 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=9    # this is for real space mesh. See GetStarted.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for Vosko-Ceperly-Alder; 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)

H2が収束しない

ctrlgen.pyの値では収束しません。

     RSMH= 0.482 EH= -0.100

と例えばEHが大きい方だけをとってください。 (よくわからん。もうすこし調べる)。

追記 遷移金属原子などでH2と同じく収束しない場合があります。RSMH,EHの数を減らして、必ずしも収束しやすくなるわけではありません。

lmfaがfailed to match log der=??? to envelope, で落ちる

mtchr2 mode 1  l = 1  r =  1.000000
l  it  ir      Rsm         Eh        phi'/phi     target
1   1  -1    1.000000  -10.000000   -6.341E-01  -8.428E-01
1   2  -4    1.000000   -0.020000   -6.730E-02  -8.428E-01
1   3  -5    1.000000   -5.010000   -5.042E-01  -8.428E-01
exit mtchr2 info =-1
Exit -1 : ftfalo : failed to match log der=-0.843 to envelope, l=1

原子、分子でRがその原子種に対して小さいときに起きます。

answer: まだなし

### 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 ou
tput)
             # It is useful to check ctrl is read in correctly or not (equivalent with --s
how 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  au=0.529177                       #<--   These % sections are to define variables.

%const  a1= 12.0  da=0
% show vars
HEADER  O
STRUC   ALAT=1/{au}
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}   # PLAT: primitive vector in ALAT. ALAT
  (in  a.u.)
                                              # DALAT gives "shift" of ALAT  (in  a.u.)
  NBAS= 1  NSPEC=1
SITE
  ATOM=Ca POS=0 0 0
SPEC
 ATOM= Ca Z= 20 R=1.0
      RSMH= 0.667  EH= -0.100  PZ=0 13.9 P=0 4.3
      KMXA={kmxa} LMXA=4
      MMOM=0 0 0


% const pwemax=3 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) N
PTS=2001 (division) even for insulator.
                #   (SAVDOS, DOS, NPTS gives no side effect for self-consitency calculaito
n).
                #
                #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 Vosko-Ceperly-Alder; 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)

match phi'/phi=??? to envelopeで落ちる

elocp:
Fit local orbitals to sm hankels, species K, rmt=1.5
 l  type    Pnu      Eval        K.E.       Rsm       Eh      Q(r>rmt)    Fit K.E.
Exit -1 : mtchre : failed to match phi'/phi=-2.052 to envelope, l=1

原子、分子でRがその原子種に対して小さいときに起きます。

answer: まだなし。

## 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 ou
tput)
             # It is useful to check ctrl is read in correctly or not (equivalent with --s
how 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  au=0.529177                       #<--   These % sections are to define variables.

%const  a1= 10.0  da=0
% show vars
HEADER  O
STRUC   ALAT=1/{au}
        DALAT=0 PLAT={a1} 0 0  0 {a1} 0  0 0 {a1}   # PLAT: primitive vector in ALAT. ALAT
  (in  a.u.)
                                              # DALAT gives "shift" of ALAT  (in  a.u.)
  NBAS= 1  NSPEC=1
SITE
  ATOM=K POS=0 0 0
SPEC
 ATOM= K Z= 19 R=1.5
      RSMH= 1.000 1.000  EH= -0.100 -0.100  PZ=0 13.9 P=0 4.3
      KMXA={kmxa} LMXA=4
      MMOM=1 0 0


% const pwemax=3 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) N
PTS=2001 (division) even for insulator.
                #   (SAVDOS, DOS, NPTS gives no side effect for self-consitency calculaito
n).
                #
                #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=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=9    # this is for real space mesh. See GetStarted.
      REL=T     # T:Scaler relativistic, F:non rela.

      XCFUN=1   # =1 for Vosko-Ceperly-Alder; 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)
                 #

efmax=2<...でおちる

BZWTS : --- Non-metal sampling ---
Fermi energy: 3.347016;  12 electrons;  Sum occ. bands: 27.249626
VBmax = 3.347016  CBmin = 3.588149  gap = 0.241133 Ry = 3.27941 eV
Exit -1 BZWTS: efmax=2 < 3.347016
wkinfo:  used  9954 K  workspace of 80000 K   in  15 K calls

これはBZのEFMAX=5などとする。フェルミエネルギーをMT半径の端から計算してるのでMT半径が小さいと見掛け上フェルミエネルギーが上がったりする。

[上の例だと efmax=3.347016 にすればいいのでその値にプログラム内部ですればいいんですよね。] lmf --inputとしてBZ_EFMAXを探してください.これの値を上げるにはctrlファイルに書くとよい. ですが、m_rctrl.Fで

     xx = 2d0
     if (ltbe) xx = 5d0
     nm='BZ_EFMAX'; call gtv(trim(nm),tksw(prgn,nm),bz_efmax,
    .  def_r8=xx,note='Find evecs up to efmax')

のxxを上げてもいいかもしれない.これで固有ベクトルをEFMAXまでしか求めないことに なってますが(LDA計算の時)、そうなってるかどうかもういちどみないとわかりません。 そもそもLDA計算ならsc計算においては占有バンドのみ求めればいいわけですが(課題にします)。

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に対応。

H2の格子定数が小さいときに止まる

いまのバージョンのちょっとしたバグ。subs/hansr.FのL240を

c      parameter (nmax=40,nm1=14,nm2=20)
      parameter (nmax=60,nm1=14,nm2=20)

としてください。makeしなおして、自分のbinにlmfa,lmfをインストールしてください。 あとは原子間距離0.6でできます。distance=0.5angstromeだとEfermi>2というメッセージで止まります。 (たぶん安全弁としてつけてるだけだとおもいますが)。生成されたctrl.h2でnk=1でもいいのかもしれないです。 また、僕の場合は、ctrl.h2において

#!/bin/tcsh
foreach distance (0.6 0.65 0.7 0.75 0.80 0.85 0.90)
  rm -f rst.h2 mixm.h2 moms.h2
  lmfa -vdistance=$distance h2 >llmfa.${distance}
  lmf  -vdistance=$distance h2 >llmf.${distance}
end
exit

ctrlstart ==============================
%const distance=0.60 
 ーーー中略 ーーー
SITE						
 ATOM=H POS=0 0 0
 ATOM=H POS=0 0 {distance}/{au}

などとしたのち、ctrl.h2をシェルスクリプトとして実行してます。 また、lmchkがrmax=0を示すのがなにか変です。次に直します。

小谷


また、この修正をほどこす前に「手でctrlを細工して得た結果」は grep 'c ' save.h2で

c distance=.6 ehf=-2.2026354 ehk=-2.2034157
c distance=.65 ehf=-2.2279434 ehk=-2.2286852
c distance=.7 ehf=-2.2413074 ehk=-2.2419823
c distance=.75 ehf=-2.2461606 ehk=-2.246794
c distance=.8 ehf=-2.2450245 ehk=-2.245609
c distance=.85 ehf=-2.2397415 ehk=-2.240255
c distance=.9 ehf=-2.2314129 ehk=-2.2318795

でした。ehfとehfがすこしずれててきもちわるいですが。。。要チェック。

Personal tools