SrTiO3 tutorial

From Ecal

Jump to: navigation, search

MTO+APW(PMT)法のrobustnessにより、ctrl fileが簡単につくることが可能になりました。 ここではMTO+APWのctrl file作成の最小限の説明を行います。 MarksOriginalDoc/fp.htmlも参考になります。MTO+APWではfplmtoにあったfloating orbitalやorbital optimzationは不要です。

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

ctrlでなくctrlsです。

そしてctrlgen.pyを"ctrlgen.py srtio3"として実行し,ctrls.srtio3をctrl.srtio3に変換させます。ctrl.srtio3が以降の計算で用いられます。



注意:ctrlgen.pyの新バージョン(2010、7月)は以下の記述とすこしちがうctrlを生成します。 たとえばSrTiO3では、ctrlgen.pyによってMTO基底が

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

となります。計算するとき、これを

    RSMH=  0 0 1.808 1.808 EH=  -0.5 -0.5 -0.5 -0.5
    RSMH2= 0 0 1.808 1.808 EH2= -2 -2 -2 -2

になおさないと計算できません。正当な理由があるのですが自動化できてません。 (単純には、Srのs軌道、p軌道は、平面波で表されるのでPMT=MTO+APWの計算法では線形独立性が 損なわれるのでそういう軌道をいれないようにする必要があります。上の置き換えはそれを意味します)。



Contents

ctrl fileの作り方

この節のすべてのファイルはTESTsamples/SrTiO3_tutorialにあります。

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

結晶構造だけのファイルを作ります。

%const da=0 au=0.529177                       #<--   These % sections are to define variables.   
%const d0=1.95/au a0=2*d0 v0=a0^3 vfrac=1 v=v0*vfrac a1=v^(1/3) 
% show vars                               
HEADER  SrTiO3 cubic                          #<--   Any header    				  
STRUC   ALAT={a1}	      
        DALAT={da} PLAT=1 0 0  0 1 0  0 0 1   # PLAT: primitive vector in ALAT. ALAT  (in  a.u.)  
                                              # DALAT gives "shift" of ALAT  (in  a.u.)  
SITE						
  ATOM=Sr POS=1/2 1/2 1/2	# Atom SPEC and position in cartesian coordinate    		
  ATOM=Ti POS= 0   0   0+0			
  ATOM=O  POS=1/2  0   0			                                                                                    
  ATOM=O  POS= 0  1/2  0		
  ATOM=O  POS= 0   0  1/2
SPEC
  ATOM=Sr Z=38 
  ATOM=Ti Z=22 
  ATOM=O  Z=8  

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

NOTE.
1 "% const" sectionは変数を宣言します。"% show vars"はprogram (lmf, lmfa, lmchk)に開始時に変数を表示するように要求します。この例が示すように、定義された変数は{a1}のように{something}として使用されます。
lmfを呼び出すときのcommand line argumentにより、lmfなどの実行時に"--v{varilable}={value}"と% const sectionで定義された変数の値を置き換えることができます。たとえば、"lmchk --va0=2.5 srtio3"は% const sectionで定義されたa0=2*d0の代わりにa0=2.5を用います。(これを使うときは、lmfにおいてはこれはsave.* に記録されます。/Li_Latticeをみて調べてください。lattice constantは save.*の中に再度記録されます。)
2. lattice constantはALAT+DALATで与えられます。lmf内で使われる計算のcutoffはALATだけで決められます。lattice constantのなめらかな関数としてtotal energyを得るためには、格子定数を変える時にDLATを使う必要があります。一般的に言って、変形しない結晶構造に対してそのようなcutoffを止めて cellを変形させる方法があります。lmf -inputをしてください。
cagetoryに対するtagを覚えてください。たとえば、SPECは行頭から始まらないといけません。tokenに対するtagは、たとえばATOM=, Z=は行頭から始まってはいけません。
3 STRUC_NBAS=とSTRUC_NSPEC=は必要ありません。更に、もしあなたが原子に対して"standard SPEC"で定義されている標準の名前を使う場合はSPECセクションを書かないでください。標準の名前はctrlgen.pyをctrlgen.py --helpatomnameで実行すると表示されます。
4 ctrlgen.py --helpしてみてください。

ctrlgen.py scriptの実行

ctrlgen.pyはまだあまりテストされていません。おそらくうまく動きますが、将来さらに改善する予定です。 何か妙な動きがあったら教えてください。これがctrlgen.pyバージョンです。 (2010Mar1stのパッケージにいれてない;以後は入れます)。

ctrlgen.py srtio3

と実行します。これは内部で他のprogramを実行します。もしエラーが起きたらPATHが正しく設定されているのかをもう一度調べてください。

これはpython scriptで、ctrl.srtio3を作ります。ctrls.srtio3は変更されません。 ここではctrlgen.pyの内部動作を説明します。(ここはskipしても構いません。)

ctrlgen.py

MT半径の生成

内部的には(L733: ctrlgen.py)"lmchk --getwsr tmp > llmchk_getwsr"を呼びます。(ctrl.tmpはctrls.srtio3 + default sectionです。)これはrmax.tmp に保存された推奨されるMT半径を生成します。 どうMT半径を決めるかはfp.htmlとlmto.htmlの"section8を見てください。("lmchk srtio3"は結晶構造のgeometricalな情報を示します。

RSMH, EH, KMXA, PZ, Pのdefault値

ctrlgen.pyはctrl.tmp2として名付けられたctrl fileにRを加えます。そして"lmfa srtio3 > llmfa.tmp2”をします。(L782: ctrlgen.py) それはRSMH,EH,KMXA,pZ,Pの推奨値を含む"mtopara.tmp2"を生成します。

Sr@ RSMH= 2.411 2.411 2.094  EH= -0.100 -0.100 -0.100 KMXA=7   PZ=0,14.9  P=0,5.3   
Ti@ RSMH= 1.393 1.393 1.055  EH= -0.100 -0.100 -0.100 KMXA=7    
O@ RSMH= 0.900 0.900  EH= -1.360 -0.512 KMXA=7    

これはctrl fileの一部として使われます。

Note: 出力(llmfa.tmp2)は"conf: sections"を含みます。それはそれぞれのSPEC_ATOMに対して原子の配置を示します。srtio3の場合はSrに対してそれは

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     5        0.000   18.000 => 2,3,4,   
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:-----------------------------------------------------   

です。 "grep conf llmfa.tmp2"でこれを見ることがでいます。ここではint(P)はvalence orbitalの主量子数を意味します。P自体はinteger部分が普通の主量子数に対応した連続した主量子数です。 MarksOriginalDoc/lmto.htmlをみてください。

RSMHとEHの組はenvelope 関数(smooth Hankel function)を明記します。fp.htmlをみてください。これらのPZ,Pはlocal orbitalがP=4.9でvalence orbialがp channelの3.5をもつことを意味します。(最初の0はs channelにはlocal orbitalがないことを意味します。)PZの2桁目の1はlocal orbitalが拡張されたlocal orbitalであることを意味します。(論理的にはそれはlocal orbialではなくMTOです。なぜならそれはextended(拡張された?広がった?)でaugumentされているからです。)

このmtopara.srtio3のdefaultの設定はsubs/writebasis.Fのwritebasis subroutineで与えられます。それはdefaultの設定を用います。

 ===Possible semicore is (Takao thinks) =========
 ----- these are important ---
  Zn,Ga,Ge    : 3d  PZ=0,0,13.9    P=0,0,4.2       
  Cd,In,Sn,Te : 4d  PZ=0,0,14.9    P=0,0,5.2       
  Hg,Tl,Pb,Bi : 5d  PZ=0,0,15.9    P=0,0,6.2       
  Lu,Hf,Ta,W  : 4f  PZ=0,0,0,14.9  P=0,0,0,5.15       
 ---- somehow meaningful, but can be negligible ---    
  Na,Mg   : 2p  PZ=0,12.9      P=0,0,3.2       
  K,Ca    : 3p  PZ=0,13.9      P=0,0,4.2       
  Rb,Sr   : 4p  PZ=0,14.9      P=0,0,5.2          
  Cs,Ba   : 5p  PZ=0,15.9      P=0,0,6.2          
 --- less important ---    
    Li : 1s         PZ=11.9    P=2.6   Q=0,1  (this is commented as #PZ.    
     In this case Q is needed. See explanation shown by lmfa without Q,   
     and SPEC_ATOM_Q part shown by "lmfa --input")  
==================================================

最後に、ctrlgen.pyはdefaultのcategoryをctrl.srtio3に加えます。これでlmfa,lmfを実行できる状態になりますが、後はctrl.srtio3を手で編集できます。

必要ならばctrl fileの編集

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

# Do lmf --input to see all effective category and token #
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 da=0 au=0.529177
%const d0=1.95/au a0=2*d0 v0=a0^3 vfrac=1 v=v0*vfrac a1=v^(1/3)
% show vars
HEADER  SrTiO3 cubic 
STRUC   ALAT={a1} DALAT={da} 
        PLAT=1 0 0  0 1 0  0 0 1
  NBAS= 5  NSPEC=3
SITE
  ATOM=Sr POS=1/2 1/2 1/2
  ATOM=Ti POS= 0   0   0+0
  ATOM=O  POS=1/2  0   0
  ATOM=O  POS= 0  1/2  0
  ATOM=O  POS= 0   0  1/2
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

 ATOM= Ti Z= 22 R=2.089960
      RSMH= 1.393 1.393 1.055  EH= -0.100 -0.100 -0.100
      KMXA={kmxa} LMXA=4
     MMOM=0 0 0 0

 ATOM= O Z= 8 R=1.595007
      RSMH= 0.900 0.900  EH= -1.360 -0.512
      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).
      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)

NOTE 1.(again): category, e.g.,SPEC は行頭にないといけません。 token, e.g, ATOM=, Z= は行頭にあってはいけません。 NOTE 2. #以降はコメントです。

ctrl fileを変更する必要があるかもしれません。 他のcategory-tokenはたとえばlmf --inputでみれます。 ctrlgen.pyは最小限のパラメタしか書きません。

NOTE: パラメタが正しいかどうかはlmf --show srtio3で確認できます。outputの最初の方にdefaultなどを補完したinput parameterが表示されます。

  • metalに対しては大きなNKABC=n1 n2 n3が必要です。
  • テストした限りではKMXA=7 はPWEMAX=5 (PW cutoff by 5Ry)に対して十分安全な(収束した?)値です。 計算時間を短くするにはもっと小さな値を使えます。
  • Mixing parametes ITER_MIX はやや複雑なパラメタです lmf はmixing scheme を自動的に選択します。 しかし、自分で 例えば"ITER MIX=A2,b=.5,n=3"と設定することもできます。
  • Feのようなmagneticな系はHAM_NSPIN=2とします。同時にMMON=0,0,2などとmagnetic momentの初期値を与えます。(この場合はd channelだけ2Sz=2) もしPZを変更したらlmfaを再度実行してください。repeat lmfa for safe; when PZ is added, need to generate density PZ as valence。
  • すべての固有値を見るにはVERBOSE=35 or 40 で十分です。
  • 計算のはじめではHAM_PWMODE=11 HAM_PWEMAX=3を進めます。 収束したら、PWEMAX to 4, 5 or soとして再度計算できます。 大きなPWEMAX には大きなKMXAが必要です。 KMXA.ge.6 for SrTiO3 for PWEMAX=5です。この場合、KMXA=20(KMAXA for PWEMAX=5の収束のテストのためです) としたときと差が少しあるのがわかります。KMXA <7を用いた場合計算がiterationの途中で止まることがあります。it gives very unstable divergent behevior (out of convergence path; then the Hamiltonan becomes very strange; we will have to improve KMAX-related pat ---> "Gaussian Polynomial" expansion of smooth Hankel function at atomic site).
  • FTMESH はcharge densityのreal-space meshを意味します。 あとの節を見てください。FTMESHの代わりにGMAXもreal space meshの指定に使えます。GMAX=9は普通は収束した値を与えます。あとのHAM_FTMESH節を見てください。

lmfの実行

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

lmfa srtio3

を実行してください。

directoryには前回行った結晶の結果が残っているかもしれません。 atomの結果から計算を行いたい場合は mix.*は前回の計算でのmixing data、 moms.*はmomentum data、 rst.*はcharge/potentailデータ? ですので、それらを消去してください。 もしそれらがあると、atom dataより優先して読まれます。 lmf --rs(something)とすると強制的にlmfaの結果であるatomic dataを読みます。(also see "lmf --help" and MarksOriginalDoc/fp.html)

lmf srtio3   

OPTIONS_PFLOAT

OPTIONS_PFLOAT=1 はold versionのbug fixを行いますので用いないといけません。 これがないとradial functionが占有状態の重心で計算されません。 (しかし、前の版ではどういうわけかうまく(正しく?安定して?)動いていました。) OPTIONS_PFLOAT=1は普通安定で、古い版より速く収束します。

l-cutoff? (STRUC_NL; SPEC_ATOM_LMXA,LMX,LMXL)の決め方

いくつかのl-cutoffが存在します。"lmf --input"でわかるとおり

SPEC_ATOM_LMXA: l-cutoff for augmentation. LMXA=NL+1 for defaults.
SPEC_ATOM_LMXL: lmax for which to accumulate rho,V in sphere
SPEC_ATOM_LMX : l-cutoff for basis (head part specfied. this corresponds to # of parameters in EH= (the same that for RSMH=).
STRUC_NL : global default lmax+1 for basis and augmentation

Specification by SPEC_ATOM is stronger than the global default by STRUC_NL. This default is not effective if SPEC_ATOM_LMXA are spefcified for all the SPECs.

普通はLMXとLMXLは指定する必要はありません、なぜなら、それらは自動的にdefaultに選択されるからです。 LMX= はEH=,RSMH= sectionsの数から. LMXL はLMXAから決められます。

I(takao) think a possible (automatic) choice of LMXA is to set LMXA= {# of EH channel} + 2 (or +1 ?) in MTO+APW mode. This is given by ctrlgen.py. Need to be tested.

それらが意味するところはMarkOriginalDoc/fp.html を見てください。

real-space mesh (HAM_FTMESH or HAM_GMAX)の決め方

electron densityのsmooth partはreal space meshであらわされます。 "HAM_FTMESH {3 integers}"は格子ベクトルの分割数です。もしくはwave vectorの最大値であるHAM_GMAXを使います。(GMAXの方がFTMESHより優先されます。) MarksOriginalDoc/fp.htmlでこれらの単語を探してください。 FTMESHはreal spece meshを明示的に指定するのに役立ちます。(lattice constantを変えながらtotal energyを計算するときなど)


FTMESHが小さすぎるとtotal energyが収束しません。FTMESHが大きすぎると、現在のversionでは数値的な問題が起きます。将来これは直さないといけません。 そのため、FTMESHの関数としてtotal energyの変化を観測する必要があるかもしれません。 しかしFTMESHに関する収束は速いことが分かっています。

There is another experimental observation to choose FTMESH in a reasonable manner; lmfを実行すると最初の方に以下のような出力が現れます。

sugcut:  make orbital-dependent reciprocal vector cutoffs for tol= 1.00E-05    
spec      l    rsm    eh     gmax    last term    cutoff    
 La       0    2.00  -0.10   3.393    1.16E-05    1055     
 La       1    1.25  -1.14   5.825    1.00E-05    5343     
 La       2    1.68  -0.10   4.551    1.01E-05    2509     
 La       3    1.21  -0.10   6.862    1.00E-05    8617     
 Ga       0    1.90  -0.55   3.574    1.01E-05    1237     
 Ga       1    1.86  -0.10   3.854    1.00E-05    1511     
 Ga       2    0.73  -0.10  11.052    1.01E-05   36243     
 O        0    0.88  -1.32   7.747    1.00E-05   12461     
 O        1    0.83  -0.38   8.930    1.00E-05   19109    

"last term"が~1e-5になるようにGMAX or FTMESH を調整する必要があります。 comment: tol=HAM_TOL(ctrlgen.pyのdefaultは1e0-6です)。last termはHAM_GMAXを大きくすると小さくなります。

"last term"is ~1d-5 かそれ以下ならreal space meshが十分細かく、計算がmeshに関して収束していることを示します。 

Overcompleteness

大きなPWEMAX (>~5Ry or so)を用いるときは注意深くやらねばなりません。 それはbasis setのlinear-dependencyが悪くなる問題、overlap matrixがとても小さな固有値をもつこと、 を発生させます。 この場合overlap matrixの固有値はとても奇妙になりえます。MTOとAPWが強くcancelするlinear combinationです。

[overlap matrix の固有値はKohn-Sham equationsの固有値と違います。]

その場合、空間を制限する(つまりある基底を除く)ことでHamiltonianの次元=MTOs+APWsで張られたHilbert空間を減らす必要があります。その操作がないと、計算が途中で異常終了するか、最終結果がとても奇妙な固有関数またはtotal energyになります。

(A)と(B)二つのlinear dependency problem を克服する方法を説明します。 (さらにOVNCUTを用いる方法がありますがここでは説明しません。)

HAM_OVEPS

一つの方法はHAM_OVEPS=1d-6 (or so)とする方法です。 [OVEPS=#. もし#が正の数だったら, lmfはoverlap matrixを対角化して, eigenvaluesが#より小さいsubspaceを除きます。] Note: lm-7.0betaKではMTOのsubspaceは保存されます、and within the space spanned by APWs orthogonalized to the space of MTO, overlap matrixの固有値を計算します。 (See subroutine zhev_tk in slatsm/zhev.F). この手続きはlm-7.0betaとは違います。 lm-7.0betaはMTOのsubspaceを保存しません。 fp.htmlも見てください。

lmfの出力で, "zhev_tk: ovlmat=" がoverlap matrixの固有値を示します。 例えば、 Li/llmf, これは"lmf li"の出力です, は

zhev_tk: ovlmat=   
   1  0.66D-05    2  0.14D-01    3  0.14D-01    4  0.14D-01    5  0.33D+00    
   6  0.10D+01    7  0.10D+01    8  0.10D+01    9  0.10D+01   10  0.10D+01    
  11  0.10D+01   12  0.10D+01   13  0.10D+01   14  0.92D+04   15  0.10D+05    
  16  0.10D+05   17  0.10D+05   18  0.11D+05    

となります。 HAM_OVEPS=1d-6とすると、 overlap matrixが0.66D-05の固有値に対応する固有関数がHilbert spaceから除かれます。いくつかの固有値はほぼ1です。これらが線形独立なAPWであることを意味します。 他の 1D+5をもつ基底はMTOsです。

zhev_tkでoverlap matrixを計算するとき、overlap matrixのdiagonal part の規格化を<APW(q+G)|APW(q+G)>=1、<MTO(q)_i|MTO(q)_i>=1d+5としています。ここで{APW(q+G)} はAPW setで、{MTO_i(q)}はMTOのsetです。

このweighted procedureでMTOのHilbert spaceを容易にpreserveすることができます。

あるchannelのMTOの除き方

MTOs のいくつかを明示的に除くことでHilbert spaceを減らすことができます。 この方法は実際上(A)よりいい方法です。なぜなら、(A)はlattice constantを変えたときにtotal energyの突然の変化を引き起こすからです。

"Li 2s"を除くにはもともとの"RSMH=1.6 1.6 EH=-0.1 -0.1"でなく"RSMH=0 1.6 EH=0 -0.1"とします。 すると、最も小さい0.66D-05の固有値はなくなります。これは0.66D-05の固有値がMTO のLi 2sから来ていて すでにAPWの重ね合わせによりうまく表現できている、という事実に由来します。 小さな固有値があったらそれがどの成分か分かる仕組みが欲しい。--Kino 15:32, 14 October 2009 (JST)

[ あるチャンネルのMTOを除くためには、 RSMH と EHのその成分を上で説明した通り0にします。]

Note: ある原子の全てのMTOを除くにはRSMH= 0 0 0 0 0 EH= 0 0 0 0 0と多くの0を書く必要があります。これはバグです。

収束のチェック

少なくともGamma pointですべてのKohn-Sham固有値を見る方がいいと思います。出力ではeigenvalue=のあとに表示されます。 linear-dependencyが悪くなると上のいくつかの固有値は高すぎる値になります。 しかしながらlmfの結果は収束している限りそれほど悪くなりません。というのはLDAのself-consistencyはFermi energyの下だけに依っているからです。

一般的に言って、minimum MTOとPWEMAX=5はtotal energyに関してそんなに悪い値を与えません。 (total energyはlattice constantを決めるときなどに用いられます。 しかし、それぞれの系の数値的な収束は系によります。) PWEMAX=10かそれ以上を使うときは(A)か(B)の手続きを用いる必要があります。 そのときは大きなKMXAを用いることを忘れないでください。(KMXA=15程度はPWEMAX=10には必要です。)


(note: 再計算するには, mix* moms* rst*を消して計算を繰り返してください。)

原子位置の緩和

原子位置を緩和させるにはLaGaO_323の例を見てください。 これはunit cellに20個の原子を持ちます。 以下を加えることで原子をrelaxさせられます。

HAM  FORCES=1   
DYN  MSTAT[MODE=5 HESS=t XTOL=.001 GTOL=0 STEP=.015]  NIT=30   

"grep ATOM log.lagao" でどう原子が動いたのかが表示されます。 よりよい力を計算するにはHAM_FRZWF=1を加える必要があります。

Personal tools