The full-potential program lmf was originally adapted from a program nfp written by M. Methfessel and M van Schilfgaarde. The method is described in some detail in the following reference (click here to view postscript file):
The following features are unique to this method:
The full-potential program builds on the basic package (documented in file lmto.html) which contains an implementation of a tight-binding LMTO program in the Atomic Spheres Approximation (ASA), and shares most things in common with it, including a number of auxiliary programs useful to both ASA and this one. For example, both methods are linear augmented-wave methods, and the wave functions inside the augmentation spheres are equivalent in the two cases. The reader is strongly advised to review that documentation before proceeding with this one. There are some differences in input between the two programs. A description of tokens available to this method are documented in the file tokens.html; which specifies which tokens can and cannot be used with this method. There is also a description of tokens specific to this method located in Section 3. See Building_FP_input_file.html or FPtutorial.html for a tutorial that helps you build an input file, and explains the output, and ASAtutorial.html for a corresponding tutorial for the ASA programs.
One important difference between the ASA and FP methods is that the FP method has no neat parametrization of total density in terms of the ASA energy moments Q0,Q1,Q2, or the representation of the potential by a few potential parameters, as in the ASA (see ASA overview in lmto.html). However, the basis within the augmentation spheres is defined from the spherical average of the potential, just as in the ASA, and the linearization proceeds in the same way. Program lmf uses the same ``continuously variable principal quantum numbers'' P described to establish a mapping between the linearization energy and logarithmic derivative at the MT boundary, and to float the linearization energy to band center-of-gravity. Thus, the description of P and IDMOD equally apply to the ASA and this program.
A second important difference is that the basis set is more complicated, and in its current form, the user must choose parameters defining the basis. This complication is the most onerous part of the present method, and there are plans for redesign; but at present is necessary to treat the interstitial region reliably. The inputs are described below; see also FPoptbas.html for a tutorial about choosing the basis optimally. Note also with the incorporation of the APW basis, other options are available.
The FP package adds two executable programs to the basic one:
lmf in an implementation of the local-density approximation, and generates density-of states, and total energy as a part of the self-consistency cycle.
Other extensions are:
Energy bands. You can run lmf in band mode to generate energy bands along lines or planes for, generating, e.g. Fermi surfaces. k-point specifications and invocation is described in generating-energy-bands.html. Particularly useful are the color weights.
Partial DOS and Mulliken analysis. lmf can generate partial dos within augmentation spheres, and construct a Mulliken analysis. DOS can be resolved by site, by site and l, or by site and lm. These options are invoked through in command-line switches. For illustrations, invoke
fp/test/test.fp co 2 fp/test/test.fp fe 2
Charge density. lmf can generate the charge density (smooth or not, with or without cores), and the contribution to the density from a selected window of states. See --wden and --window in command-line switches
Core-level spectroscopy. lmf can generate EELS spectra, which involve matrix elements between core and valence electrons. The EELS option is invoked with a command-line switches For illustrations, invoke
fp/test/test.fp fe 2 fp/test/test.fp crn 2
LDA+U. The LDA+U functional was built into lmf in v6.15 and later, by Walter Lambrecht. LDA+U needs in addition to the LDA parameters U and J for selected orbitals, which are empirical. The LDA+U constructs an additional potential for a particular l subblock (m=-l..l) from the U supplied by the user, and the density-matrix, which is generated by lmf. Two modifications must be added to the input, which are described here. In a strictly LDA calculation, complete information is carried by the density, contained in the restart file, rst.ext. In the LDA+U case, complete information is carried by density and on-site density matrices, which are contained in file dmats.ext. An example that illustrates LDA+U method is ErAs, which you can run by
fp/test/test.fp erasErAs is an interesting case because LDA puts all 4 minority f electrons in a single extremely narrow band at the Fermi level. In LDA+U the minority f is spit into a 4- and 3-manifold; see PRB 67, 035104 (2003).
GW. lmf is designed to work in coordination with a GW package by T. Kotani (the GW package comes comes separately). lmf acts both as a driver for the GW package and also can be used in a self-consistent GW cycle. An extra driver lmfgw is compiled as part of this extensions. Use of this driver is described in the GW driver documentation. You need the extension package GW.version.tar.gz. Also you will need the GW package itself. For illustrations of the driver invoke
gw/test/test.gw si gw/test/test.gw gas
Spin-Orbit coupling. lmf solves the scalar Dirac hamiltonian. The dominant difference between the full Dirac hamiltonian and the scalar one is the spin orbit coupling, which can be added as a term λL·S to the scalar Dirac Hamiltonian. Starting with v6.15, lmf can add λL·S to the scalar Dirac hamiltonian (courtesy of A. Chantis). It is possible to include the full L·S or just the L_{z}·S_{z} part. Beginning with v7.9, L_{z}·S_{z} + (L·S−L_{z}·S_{z}) can be added where the last term is treated in an approximate manner. The approximation turns out to be rather good. See here for some description and analysis of all three approximations. L·S in all three forms can be combined with the self-energy read by a QSGW calculation. In the L_{z}·S_{z} and approximate L·S forms the effect of SO coupling can be passed through to the GW code, to include its effect on the self-energy.
Use HAM_SO to turn on SO coupling. For illustrations of all three kinds of approaches, try
fp/test/test.fp felz fp/test/test.fp gasls
Local orbitals. In v6.12 and later, local orbitals may be added to the basis set. These orbitals are important when energy bands over a very wide energy window are required, when high accuracy is needed for shallow (semi)core states, or for energies far above the Fermi energy. Examples of the former occur in oxides: bond lengths are small and cations with shallow p orbitals extend somewhat beyond the augmentation radius.
Local orbitals are presented in one of two flavors. The first,
conventional type of local orbital is constructed by solving the
radial Schrodinger at a different linearization energy than the usual
valence states, and then subtracting off a particular linear of the
valence wave function phi and energy derivative phidot such that the
local orbital's value and slope vanish at the augmentation
radius. Thus
fp/test/test.fp gas fp/test/test.fp cuThe Ga 3d semicore the high-lying As 5s state are included as local orbitals. In the Cu case, the high-lying Cu 4d is included, which is important in GW calculations.
The second, extended, kind of local orbital can only be used for semicore states. Instead of artificially subtracting off some linear of the phi and phidot to make the orbital vanish at the augmentation radius, a smooth Hankel tail is attached to the orbital. The smoothing of tail is constructed to match as well as possible the kinetic energy of the semicore state. This type of orbital has the advantage that the valence envelope function need not `carry' the tail of the semicore state. Its drawback is that more things can "go wrong," namely it may fail to do a good job of fitting the kinetic energy. an example that illustrates the second kind of local orbital is
fp/test/test.fp srtio3The Sr 4p and Ti 3p semicore states are included as local orbitals. In the first part of the test, they are included as local orbitals of the first type; then the last step is recalculated using local orbitals of the second type.
Floating Orbitals. In v6.15 and later, floating orbitals may be added to the basis set. These orbitals can be important when very accurate calculations are needed in open systems, e.g. when reliable energy bands are needed for a wide energy window. These orbitals differ from the usual smooth Hankels in that they are not centered at an atom. They are augmented just as the other orbitals, but there is no augmentation radius and thus no "head" sphere.
The following illustrate the inclusion of floating orbitals in the basis:
fp/test/test.fp te fp/test/test.fp gaslc
Augmented Plane Waves. Starting with v6.17, Augmented Plane Waves can also be included in the basis. They play the same role as floating orbitals do, but APWs are superior because they are simpler to use (there are no parameters and they do not need to be located at any particular site), and the control over convergence is more systematic.
The following illustrate the inclusion of APWs in the basis:
fp/test/test.fp te fp/test/test.fp srtio3 fp/test/test.fp felz 4
Parallel Implementation. Two separate parallel versions of lmf have been made (courtesy of A. T. Paxton). One parallelizes over k-points, and is the most efficient for scaling; the other parallelizes many points at a lower level. Installation is not automatic, however. See instructions for installation.
The full-potential reads from the same input file as the basic package, documented in lmto.html. File input-file-style.html describes in generic terms how the input file is grouped into categories and tokens; tokens.html documents tokens shared by both the basic package and this one. Building_FP_input_file.html documents a (nearly) automatic input file generator. This section documents tokens specific to lmf.
Floating orbitals. To specify a floating orbital, add a new species with atomic number 0 and augmentation radius 0. No augmentation parameters (LMXA,LMXL, etc) are used, bu you must also specify the envelope function. In the following example
SPEC .. ATOM=X Z=0 R=0 LMX=1 RSMH=1.5,1.5,1.5,1.5 EH=-.3,-.3,-.3,-.3species X is specified as a floating orbital. Although parameters RSMH and EH are supplied for spdf orbitals, the LMX=1 token specifies that just s and p functions are to used. In the SITE category, you specify where to place the orbitals, e.g.
SITE .. ATOM=X POS=1/2 1/2 1/2 ATOM=X POS=3/4 3/4 3/4You can in principle put the floating orbitals anywhere at all. In practice, it makes little sense to put them too close to the regular orbitals or too close to each other: numerical instabilities can arise if you do so. For an example that illustrates the floating orbitals, try elemental Te, which has a very open structure:
fp/test/test.fp teNine sites with floating sp orbitals are added to the three actual atoms, increasing the basis from from 39 to 75 orbitals. The test shows that addition of the floating orbitals lowers the total energy by 4 mRy and minimally affects the forces. (Note that this test also includes APWs, which in this case offers a more convenient and efficient method to make the basis complete.)
Local orbitals. The local orbital is specified with the PZ= token, which is its analog of the continuously variable principal quantum number P. Conventional local orbitals may be used either for semicore or high-lying states. In this case, specify the integer part of PZ must either as one less than the integer part of P for a semicore state or one greater for a high-lying state. As for the fractional part, it is recommended that you specify 0.9 or so for the semicore case; while for high-lying states 0.3 or so is reasonable (a larger number specifies a higher linearization energy).
Extended local orbitals:. Extended local orbitals have a smooth Hankel envelope attached to the radial wave function which spills out into the interstitial. The tail is matched value and slope to the radial function, but lmf will try to vary both RSMH and EH to match the kinetic energy as well. To specify that a semicore local orbital is of the extended type, simply add 10 to token PZ=.
Hazard for conventional local orbitals:. Conventional local orbitals applied to semicore states must rely on the valence envelope functions for representation of these deep states in the interstitial. Therefore be advised to choose the smoothing radius for at least one envelope function rather small (typically around 1 or so), with a corresponding energy rather negative (typically around -1 or so).
Hazard for extended local orbitals:. It may happen that lmf cannot find any reasonable combination of RSMH and EH that matches the envelope's value and slope to the radial wave functions. Since this matching condition is required, lmf stops with a message like
Fit local orbitals to sm hankels, species Sr, rmt=3.1 l type Pnu Eval K.E. Rsm Eh Q(r>rmt) Fit K.E. Exit -1 : mtchre : failed to match phi'/phi=-5.127 to envelope, l=1This usually happens when the fractional part of PZ is too large (it was set to 0.98 in this case). The solution is to change the fractional part of PZ.
lmf will try to vary both RSMH and EH to match (value,slope,kinetic energy), and it may happen that lmf can match the first two but not the kinetic energy of the envelope to the radial wave function. This doesn't cause an error but it can reduce the accuracy of the orbital. It is seen in the following output:
l type Pnu Eval K.E. Rsm Eh Q(r>rmt) Fit K.E. 1 low 4.900 -1.708558 -0.923699 1.00000 -0.30708 0.07250 -0.289954*The true and fit K.E. differ. The solution is again to alter PZ (4.90 in this case). It is not essential, however. If PZ floats in the course of a self-consistent cycle, lmf usually picks a good value for PZ.
At present, the P with the higher principal quantum number is automatically frozen to its input value (that is, P is frozen when PZ is a semicore orbital, and PZ is frozen when it is a high-lying orbital). Thus, only the P corresponding to the deepest state is allowed to float (see IDMOD description).
APWs. It is very simple to add augmented plane waves to the basis. APW's and smoothed LMTO's can be naturally integrated together; a basis with both kinds of orbitals is called the "Plane Muffin-Tin" or PMT basis.
To include APWs add these tokens to category HAM:
PWMODE=1 PWEMIN=# PWEMAX=#PWMODE=1 tells lmf to include APWs in the basis (The default is PWMODE=0 which means no APWs will be included). PWEMAX is an energy cutoff (in Rydbergs): G-vectors whose energy fall below PWEMAX will be included. Because MTOs can also be included, and they already rather accurately describe low-energy states, it makes sense to include a minimum-energy cutoff as well, i.e. include APWs that satisfy PWEMIN < E < PWEMAX .
There are two further considerations the user should take note of. First, when many APWs are used, the overlap matrix can become nearly singular. When this happens the secular matrix becomes unstable, which can lead to somewhat wrong or even nonsensical results. To protect against this, add to HAM the following:
OVEPS=#If # is a positive number, lmf will diagonalize the overlap matrix, and eliminate the subspace which has eigenvalues smaller than #. This procedure will largely eliminate these instabilities. A reasonable choice is OVEPS=1E-7.
Second, high-energy APWs oscillate much more rapidly than MTOs do. That means polynomials of higher order may be needed to augment the envelope functions. The order of polynomial cutoff is set by (species-specific) token KMXA in the SPEC category. When using PWEMAX larger than 3 or 4, take care that the results are converged wrt KMXA. As a rough guide, the dimensionless integer KMXA should be approximately as large as PWEMAX is (in Rydbergs).
You can see the interplay between MTOs and APWs in the Figure, which shows the total energy E in SrTiO_{3} as a function of the number of APWs included in the basis. (The upper cutoff corresponds to PWEMAX=10.)
Curve "1" includes just s and p MT orbitals on the O sites, sufficient for a crude representation of the valence bands (no orbitals are available for conduction bands). Also included are local orbitals to represent Sr 4p and Ti 3p semicore states, as they are too extended to be appoximated as core states. A large number of APWs is need to get a good total energy: something like 150 APWs are needed to converge E to within about 50 mRy of the converged result. (Had the O s and p not been present, many more APWs would have been required.)
Curve "4" corresponds to an extreme tight-binding basis, consisting of sp orbitals on Sr, spd orbitals on Ti (the conduction band is mainly Ti d), and sp orbitals on O. The total energy of the MTO basis alone (no APWs) is rather crude --- more than 200 mRy underbound. However, only 25 orbitals (plus 6 for the local orbitals) are included in this basis. The energy drops rapidly as low-energy APWs are included: adding about 40 APWs is sufficient to converge E to about 50 mRy. As more APWs are added, as more APWs are added, the gain in energy becomes more gradual; indeed convergence is very slow for large E.
Curve "5" differs from curve "4" only in that a Sr d orbital was added. With the addition of these 5 orbitals, the MTO-only basis is already rather reasonable. This would be the smallest acceptable MTO-basis. As in the Curve "4" basis, there is initially a rapid gain in energy as the first few APWs are added, followed by a progressively slower gain in energy as more APWs are added. The blue curve close to Curve "5" uses the same MTO basis, but PWEMIN=1. Interestingly, the blue and black curves almost superimpose on each other. It means the gain in energy seems to depend on the number of APWs, but not details of their shape.
Curve "6" is a standard LMTO minimum basis: spd orbitals on all atoms. Comparing Curve "5" or Curve "6" to Curve "1" shows that the MTO basis is vastly more efficient than the APW basis in converging the total energy. This is true until a minimum basis is reached. Beyond this point, the gain APWs and more MTOs improve the total energy with approximately the same efficiency, as the next tests show.
Curve "8" is a standard LMTO larger basis: spdspd orbitals on Sr and Ti, and spdsp on O. Comparing curves "6" and '8" shows that the efficacy of any one orbital added to to the standard MTO minimum basis is rather similar in the APW and MTO cases. Thus, increasing the MTO basis from 51 to 81 orbitals in the MTO basis lowers the energy by 33 mRy; adding 33 APWs to the minimum basis (curve "6") lowers the energy by 36 mRy. APWs, however, have the advantage that they are simpler to include.
Curve "11" enlarges the MTO basis still more, with Sr: spdspd, Ti: spdspd, O: spdspd. Also local orbitals are used to represent the high-lying Ti 4d and O 3s and 3p states. For occupied states, these orbitals have little effect, but they are important for unoccupied states higher than about 2 Ry. Comparing the last curves, it appears that finally E is at last converging to an absolute minimum at around -2.760 Ry. The instabilities in the overlap matrix mentioned above are largely controlled, but there is a small uncertainty in the absolute energy on the order of 1 mRy. It is clear even from inspection of Curve "1", that a pure APW basis would require an extremely large number of PWs to reach this level of convergence.
LDA+U. To get started, LDA+U needs in addition to the LDA input two modifications:
IDU= 0 0 2 2 UH= 0 0 0.1 0.632 JH= 0 0 0 0.055The IDU token tells lmf that no U is to be added to the s or p channels, but that a U is to be added to the d and f channels. IDU=2 specifies LDA+U functional style 2; this is the "Fully Localized Limit" described in Liechtenstein, PRB 52, R5467 (1995)). IDU=1 specifies the "Around Mean Field" functional (Petukhov, PRB 67, 153106 (2003)). U=0.1 Ry is included on the d orbital, and U=0.632 is included on the f orbital. Additionally J=0.055 is put on the f orbital.
The density-matrix is read from and written to a file dmats.ext. Two density-matrices (1 for each spin) are written to this file in a (2l+1) by (2l+1) block for each l block for which a U is defined. dmats.ext is an ASCII file which you can read, and it's quite useful to interpret what's going on. The diagonal parts are the occupation numbers and are the most important. Note that the file may be stored in either spherical harmonics or real harmonics, depending on how SHARM= is set in the OPTIONS category.
Usually it is too tedious to supply dmats.ext itself, especially since you typically won't know what to choose. Instead, you can supply an "occupation numbers" file occnum.ext, which is a starting guess for the density-matrix (its diagonal part). occnum.ext has one line of (2l+1) numbers for the occupation numbers of the first spin, following by a line with the occupation numbers for the second spin. In this test:
fp/test/test.fp erasthe script assumes a particular starting spin configuration through the occupation number file occnum.eras it uses. Er has 11 f electrons, 7 of which go into the majority channel and 4 into the minority channel. There is some choice in which m states to fill and which to keep empty. A key point is that the self-consistent solution you end up with will depend on this choice. The ErAs test uses the following input file for occnum.eras :
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 1 1 1 1 0 0The first and second lines are occupation numbers for the majority and minority d channel; the third corresponds to the majority f channel where all states are taken to be filled. The last line corresponds to the minority f channel. In this case, m=-2,-1,0,1 are filled and m=-3,2,3 are empty. As the script notes, different choices of starting occupation numbers lead to different self-consistent solutions. The one with the lowest energy is that which satisfies Hund's rule (m=0,1,2,3 filled and m=-3,-2,-1 empty).
The occupation numbers are by default correspond to spherical harmonic representations of Y_{lm}. If you want to define the occupation numbers in real harmonics, put
% realon the first line of occnum.ext.
Note that in the LDA case, complete information is contained in the density, stored in file rst.ext. In the LDA+U case, complete information is contained in the combination of rst.ext and dmats.ext.
The density-matrix is at present mixed independently of the charge, with linear mixing. (This will likely change in future) To do the mixing there is a special-purpose parameter in the ITER category
UMIX=## is a parameter between 0 and 1. It plays the same role for the density-matrix as the mixing beta plays for the mixing of the regular density.
There is additionally a tolerance parameter
TOLU=#that tells lmf to stops mixing dmats when its rms change falls below TOLU. Usually it's not necessary, and setting TOLU=0 (or leaving it out) means it plays no role.
Partially occupied core holes. Calculations involving partial core hole occupancy are useful in the context of Slater transition-state theory, which undoes most of the error in the LDA description of the core hole eigenvalue (see example J. Phys. Cond. Mat. 12, 729 (2000)). You specify which orbital in which species is to be treated as a partially occupied by core by adding a token C-HOLE= to the SPEC category. You also have to specify what the partial occupation is, which you do with token C-HOLE= . An example is the N 1s core. A core hole of 1 electron can be put in by
C-HOLE=1s C-HQ=-1Note the sign of the charge. The number refers to the excess electron charge. to put in a hole, use a negative charge. As an illustration in CrN, you can run the test case
fp/test/test.fp crn
Core holes are also useful as an approximate workaround in the LDA context to deal with (almost) nonbonding f electrons. In most 4f systems, the f states get shifted away from the Fermi level, even though the LDA typically is unable to do this (except for Gd), because it lacks a nonlocal exchange as in LDA+U. An approximate workaround is to treat the 4f electrons as core. For Gd in particular, the 7 majority states should be filled, while the 7 minority states empty. Thus a core hole of -7 is required, but it is necessary to further specify the spin polarization of the core. This is accomplished with a second argument. For the Gd case (4f core, -7 excess electrons, with the core magnetic moment +7), use
C-HOLE=4f C-HQ=-7,7As an illustration, try
fp/test/test.fp gdn
Restart file editor. Not documented.
Q: Running lmf, I faced an error like
Exit -1 problem in locpot -- possibly low LMXA, or orbital mismatch, species Rh CPU time: 0.246s Wall clock 0.286s 10:17:34 23.10.2013 onHow can I fix it?
A: Usually this happens because the input conditions have changed between runing lmfa and lmf. If for any reason the core configuration is not consistent with the atm file (if overlapping atomic densities) or the rst file this problem may occur.
Suppose for example, you have HAM_AUTOBAS_LOC set, which tells lmf to read information about local orbitals from basp.ext, and also tells lmfa to identify local orbitals to include in the basis. What basis information it generates, lmfa puts into basp0.ext .
Suppose you invoke lmfa and it finds a new local orbital which wasn't included in the before (e.g. missing from the basp.ext). If you now include this local orbital, e.g. copy basp0.ext to basp.ext, you are changing the core configuration which causes problems for lmf. To fix the problem run lmfa again with the new conditions (and remove the rst file if you have one). See this tutorial for an example.
Q: I want to generate the total DOS (SAVDOS=T, DOS=range) over a specified range.
But lmf automatically resets the range of DOS to cover all the occupied states.
When I do NiO (with deep core), the range is very large. Is there a way to avoid this problem?
A: This is a small bug, introduced to avoid problems in the sampling integration algorithm.
Use:
--no-fixef0and the range will not change.
GVLIST: gmax = 4.702 a.u. created 973 vectors of 2744 (35%) mesh has 14 x 14 x 14 divisions; length 0.588, 0.588, 0.668 SGVSYM: 170 symmetry stars found for 973 reciprocal lattice vectors ... sugcut: make orbital-dependent reciprocal vector cutoffs for tol= 1.00E-05 spec l rsm eh gmax last term cutoff Te 0* 1.30 -0.10 5.220 9.30E-05 973* Te 1* 1.30 -0.10 5.597 4.36E-04 973* Te 2 1.30 -0.10 5.976 2.04E-03 973* Te 0* 1.40 -1.00 4.847 2.11E-05 973* Te 1 1.40 -1.00 5.182 9.90E-05 973* E 0* 1.50 -0.30 4.524 1.06E-05 847 E 1 1.50 -0.30 4.823 2.01E-05 973*The table reflects how well each of the basis orbitals is converged in a PW expansion, which is used to make matrix elements of the interstitial potential. Only the E s orbital is converged below the specified tolerance (10^{-5} ). Convergence in the Te d orbital is particularly poor; however, the Te d is only minimally part of the bonding, and its poor convergence doesn't affect the energy much. In this case (3 Te atoms), the HF and HK total energies differ by ~10^{-4} Ry. Below is a table showing the convergence in HF energy as a function of the number of divisions, and also the difference between HK and HF energies.
n ehf ehf-ehk 14 -14561.140624 0.000089 15 -14561.133966 0.000017 16 -14561.132309 0.000005 20 -14561.131170 0.000000
Exit -1 ICLBSJ: sought atom no. 1 in class 7 but only 0 atoms existA: Usually it means that there is an improper synchronization between site and species data. (e.g. species are defined that have no corresponding sites)
Exit -1 symprj: no site mapped into first under op #Why does it not detect this problem earlier?
SYMCRY: crystal invariant under 8 symmetry operations for tol=1e-5and can it be set somewhere to control how symmetry routines work?
--fixpos:#where # is the tolerance you specify for sites being `close enough' to coincident.