Scratchpad Wiki

PMT: V: SrTiO3 tutorial

トーク0 シェアする

Thanks to the robustness in the MTO+APW(PMT) method, it is possible to set up a ctrl file easily. Here are minimum instructions to set up a ctrl file for MTO+APW. (Read also MarksOriginalDoc/fp.html. The floating orbitals or orbital optimization are unnecessary for MTO+APW.)

To set up a ctrl file, you first have to write down crystal structure in the ctrls.* file, e.g. ctrls.srtio3. Then you can run, e.g., " srtio3". convert ctrls.srtio3 into ctrl.srtio3, which can be used by the following lmf calculaiton (do lmfa in advance). We explain these procedure below.

We expect "" should work well, and we will improve it in future. But it is not well tested yet. So let me know something strange happens.

How to make a ctrl file 編集

All the files in this section is in TESTsamples/SrTiO3_tutorial.

Prepare crystal structure file. "ctrls.{name}" 編集

Let's make a file which contains only a crystal structure.

%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}	      # NBAS: number of site, NSPEC: number of atomic species  
        DALAT={da} PLAT=1 0 0  0 1 0  0 0 1   # PLAT: primitive vector in ALAT. ALAT  (in  a.u.)  
  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
  ATOM=Sr Z=38 
  ATOM=Ti Z=22 
  ATOM=O  Z=8  

Save it as ctrls.srtio3

1. "% const" section defines variables, "% show vars" requests programs(lmf, lmfa, lmchk) to show defined variables at the begining of console output. As this example shows, defined variables are used as {something}, e.g, {a1}.
By command line arguments when you invoke lmf, it is possible to replace values of variables defined in the % const section by "--v{variable}={value}" argument when you run lmf or so. For example, " lmchk --va0=2.5 srtio3" uses a0=2.5 instead of a0=2*d0 defined in the % const section. (When you use this, this is recorded in save.* file in the case of lmf. See Li_Lattice/ sample and examine it; the lattice constant is reecorded in save.* file.)
2.The lattice constant is given by ALAT+DALAT. The computational cutoffs internally used within lmf is determined by ALAT only. In order to obtain total energy as a smooth function of lattice constant, it is necessary to to use DALAT when you change lattice constant. Generally speaking, there are ways to deform cells with keeping such cutoffs for crystal structure without deformation.Do lmf --input
Remember The tag for category, e.g, SPEC must start at the beginning of the line. The tags for tokens, e.g, ATOM=, Z= must not start at the beginning of the line.
3. STRUC_NBAS= and STRUC_NSPEC= are not necessary. Furthermore, you don't need to supply SPEC if you use standard name for atoms defined in "standard SPEC", which is shown when you run without arguments.

Run a script 編集

Then run srtio3

It calls other programs. If an error occurs, probably your PATH is wrong.

This is a python srcipt. It generates ctrl.srtio3 file (ctrls.srtio3 is unchanged). Here we explain how works (you can skip this explaination below).


Generate MT radius編集

Internally( it calls "lmchk --getwsr tmp >llmchk_getwsr". (ctrl.tmp is ctrls.srtio3 + some default section). This generates suggested MT radius which are stored in rmax.tmp file. See fp.html and lmto.html#section8 for discussion about how to determine R. ("lmchk srtio3" shows geometrical informations of its crystal structure).

Default setting of RSMH, EH, KMXA, PZ, P編集

The adds these obtained R to a ctrl file named as ctrl.tmp2. Then it does "lmfa srtio3 >llmfa.tmp2" ( It generates 'mtopara.tmp2', which contains suggested values for RSMH, EH, KMXA, PZ and P as

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    

This is used as a part of ctrl file.

Note: Console output (llmfa.tmp2) includes "conf: sections"; it shows atomic configurations for each SPEC_ATOM. In the case of srtio3, for Sr, it is

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 =>    

(you can see this by "grep conf llmfa.tmp2".) Here int(P) means the principle quantum number for valence orbitals. (P is the continuous principle quantum number; integer part corresponds to the usual principle quantum number. See MarksOriginalDoc/lmto.html.

A pair of RSMH and EH specify an envelope function (the smooth Hankel function). See fp.html. These PZ and P means that local orbital is with P=4.9, and the valence orbital is with 5.3 for p channel (first zero means no local orbital for s channel). 1 of 10th ditit in PZ means the local orbital is the extended local orbital (Logically speaking, it is not a local orbital but a MTO because it is extended and augmented).

This default setting of mtopara.srtio3 is given in a subroutine writebasis in subs/ioorpb.F. It uses a default setting setting:

 ===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,13.2       
  K,Ca    : 3p  PZ=0,13.9      P=0,0,14.2       
  Rb,Sr   : 4p  PZ=0,14.9      P=0,0,15.2          
  Cs,Ba   : 5p  PZ=0,15.9      P=0,0,16.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")  

Finally, add some default categories to ctrl.srtio3. Look into it. We have to edit it by hand as follows.

Edit the ctrl file, if neccesary 編集

This file contains all the parameters to run lmfa and lmf. The ctrl.srtio3 generated by from ctrls.srtio3 is given as

# Do lmf --input to see all effective category and token #
VERS    LM=7 FP=7
             # version check. Fixed.
             # 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
  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
 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): The categories,e.g.,SPEC must start at the beginning of the line. The tokens,e.g, ATOM=, Z= must not start at the beginning of the line. NOTE 2. In each line, columns after "#" is neglected (these are comments).

You may need to edit this ctrl file. There are other category-token to control lmf in detail (do lmf --input to see it). only inlcudes minimum things.

NOTE: To check that lmf reads parameters correctly, do lmf --show srtio3. At the beginning of console output, supplied input is shown.

  • For metal, we need large NKABC=n1 n2 n3
  • KMXA=7 is safe setting up to PWEMAX=5 (PW cutoff by 5Ry), as far as we tested. In cases you can use smaller number to reduce the computational time.
  • Mixing parametes ITER_MIX is a little difficult to understand. lmf chooses a mixing scheme automatically. But you can set it by yourself, e.g. as "ITER MIX=A2,b=.5,n=3".
  • For the magnetic case like Fe, you have to use HAM_NSPIN=2 (this expression means NSPIN=1 in HAM category). Further, add initial magnetic momemnt specification MMOM=0,0,2 (this means polarized in d channel). We now have complete ctrl file. You can run lmfa and lmf (repeat lmfa for safe; when PZ is added, need to generate density PZ as valence).
  • VERBOSE=35 or 40 might be better to observe all eigenvalues.
  • To start calculations, we recommend HAM_PWMODE=11 HAM_PWEMAX=3 After it is converged, you may enlarge PWEMAX to 4, 5 or so. Note that larger PWEMAX requires larger KMXA. For example, we have to use for SrTiO3 for PWEMAX=5. In this case, you will see only a little difference even when you use KMXA=20 (this is a convergence check about KMXA for PWEMAX=5). When we used KMXA <7, the calculation stoped suddenly during iteration cycle; 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 defines real-space mesh for charge density. See a section below for setting. Instead of FTMESH, you can use GMAX to specify the real-space mesh. Usually GMAX=9 gives almost converged results. See HAM_FTMESH section below.

Run lmf 編集

Now you are ready to run lmf. Do

lmfa srtio3

Then you need to remove mix.* (for mixing data in previous itteration), moms.*(momentum data in previous itteration), rst.* (this contains self-consistent results) if you want to start over. If they exist, they are read by lmf. (This is a defalut setting. To change this beheviror, use --rs{something} options; do "lmf --help" and read MarksOriginalDoc/fp.html).

lmf srtio3   


You have to use OPTIONS_PFLOAT=1 (bug fix of old versions). Without it, radial functions are not calculated at the center of gravity of the occupied states. (But the previous version somehow worked well). OPTIONS_PFLOAT=1 is usually stable and gives more rapid convergence than old versions.

How to determine l-cutoff? (STRUC_NL; SPEC_ATOM_LMXA,LMX,LMXL)編集

There are some kinds of l-cutoffs. As "lmf --input" shows, they are

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.

Usually you neither need to set LMX nor LMXL because they are automatically choosen in defaults (LMX= is given by the numbers of EH=,RSMH= sections. LMXL is chosed to be 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 Need to be tested.

Read MarkOriginalDoc/fp.html to know what they mean.

How to determine real-space mesh (HAM_FTMESH or HAM_GMAX). 編集

The smooth part of electron density is represented on a real-space mesh. "HAM_FTMESH {3 integers}" is the number of divisions along the three lattice vectors. Or you can use HAM_GMAX, which is the maximum size of wave vector. (if GMAX exists, FTMESH is neglected). Search these words in MarksOriginalDoc/fp.html. FTMESH is useful to set number of real-space mesh explicitly (when you require smoothness of total energy when you change lattice constant).

Too small FTMESH will not give converged total energy. If FTMESH is too large, however, current version causes a numerical problem. It will be fixed in future. Thus, you may need to observe the change of the total energy (and eigenvalues) as function of FTMESH. Look into examples. In anyway the convergence as a function of FTMESH is rapid.

There is another experimental observation to choose FTMESH in a reasonable manner; when you run lmf (run and stop it within seconds), console output near its begining shows that

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    

look into "last term", you may need to tune GMAX or FTMESH so that "last term" is ~1e-05. (When "last term"is ~1d-5 or below, it indicates that the real-space mesh is dense enough, that calculations are converged well as for the mesh.)

Overcompleteness 編集

When you use large PWEMAX (more than ~5Ry or so), be careful! It can cause a problem of poor linear-dependency of the basis sets; the overlap matrix of a basis set can have very small eigenvalues. The eigenfunction of the overlap matrix corresponding to such a small eigenvalue can be very strange; it is the linear combination of MTOs and APWs where there exist strong cancellations each other.

[Don't confuse eigenvalues of the overlap matrix with those of the secular equation (eigenvalues of the Kohn-Sham equations).]

In such a case, it is necessary to reduce the dimension of the Hamiltonian (the Hilbert space spanned by MTOs+APWs) by eliminating some subspace (or removing some of basis functions). Without such a procedure, calculation stops suddenly, or the final results show very strange eigenvalues and total energies.

There are two way (A) and (B) as explained below to overcome this linear dependency problem (OVNCUT is another way, but not explained here).


A way is to set HAM_OVEPS=1d-6 (or so). [OVEPS=#. If # is a positive number, lmf will diagonalize the overlap matrix, and eliminate the subspace which has eigenvalues smaller than #.] Note: In lm-7.0betaK, the space of MTO is preserved, and within the space spanned by APWs orthogonalized to the space of MTO, we calculate eigenvalues of overlap matrix. (See subroutine zhev_tk in slatsm/zhev.F). This procedure is a little different from lm-7.0beta (it does not preserve space of MTO). See fp.html.

In the console output of lmf, "zhev_tk: ovlmat=" shows the eigenvalues of the overlap matrix. E.g., Li/llmf, which is an output of "lmf li", shows

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    

If we set HAM_OVEPS=1d-6, an eigenfunction corresponding to the eigenvalue 0.66D-05 of the overlap matrix is removed from the Hilbert space. Some of eigenvalues are ~1, this means these corresponds to APW which are linearly-independent. Other basis with 1D+5 are MTOs (when we calculate overlap matrix in zhev_tk, we set normalization (diagonal part of overlap matrix) as <APW(q+G)|APW(q+G)>=1, and <MTO(q)_i|MTO(q)_i>=1d+5, where {APW(q+G)} denotes the set of the APW, and {MTO_i(q)} denotes that of the MTO. With this weighted procedure, we can easily preserve the Hilbert space of MTO).

Remove MTOs for the channel 編集

Removing some of MTOs explicitly can also reduce the Hilbert space. This method is better than (A) in practice since (A) can cause a sudden change of total energy when you change lattice constant(s).

To remove the "Li 2s", set "RSMH=0 1.6 EH=0 -0.1" instead of original setting as "RSMH=1.6 1.6 EH=-0.1 -0.1". Then the smallest eigenvalue 0.66D-05 disappears. This is because the very small eigenvalue 0.66D-05 originates from the fact that the MTO of Li 2s is already well represented by the superposition of APWs.

[To remove MTOs of a channel, set zero for RSMH and EH for the channel as we explained above.

Note: To remove all MTO from an atom use many zero (number of LMXA) as RSMH= 0 0 0 0 0 EH= 0 0 0 0 0 ; this is a minor bug. ]

Checking convergence編集

It is better to watch all the Kohn-Sham eigenvalues (shown after "eigenvalue=" in console output;) at least at the Gamma point or so; some of eigenvalues from the top becomes too high when the linear-dependency becomes poorer; however, the convergence of lmf can be not so bad as far as it converges since the self-consistency in LDA is only related to the eigenfunctions below the Fermi energy.

Generally speaking, PWEMAX=5 with minimum MTOs gives not so bad result as for total energy (the energy which can be used to determine lattice constant or so; but the numerical convergence for each systems are dependent on sysmtems). When you use PWEMAX=10 or more, you must need to use procedure (A) or (B). Don't forget to use large KMXA (KMXA=15 or so is needed for PWEMAX=10)

(note: to start over, remove mix* moms* rst*, and repeat it).

How to relax atomic positions 編集

For relaxation of atoms (atomic force is calculated and relaxed), see a sample, LaGaO_323. This is a perovskite with twenty atoms per cell. What you have to do is to add

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

"grep ATOM log.lagao" ---> it shows how the atomic positions are moved. You may need to add HAM_FRZWF=1 to calculate better force.