DFT calculations for atoms

Paolo Giannozzi has developped a code that allows you to do density functional calculations for atoms (actually it is meant to produce pseudopotentials). Since this is a well defined, small problem it is ideal for getting an idea of how DFT calculations are done.
  • Download code (about 200 kB) and unpack in a new directory (needs about 1 MB)
  • edit Src/Makefile: change Fortran compiler (FC) to g77 and, if you don't have the BLAS and LAPACK libraries installed, comment LIBS= and uncomment LAPACK=. cd .. and make ld1 wfcgraf
  • to check if everything went ok, run Src/ld1 < Doc/ld1.in | tee ld1.out and diff ld1.out Doc/ld1.out. The input to ld1 is explained at the beginning of Src/ld1.f (see also below).
  • The file vs2p0d3.wfc contains the wavefunctions, which can be plotted with Src/wfcgraf with the input
        filwf > vs2p0d3.wfc
        filwfp >
          wavefunction #  1 : label > 1s
          wavefunction #  2 : label > 2s
          wavefunction #  3 : label > 2p
          wavefunction #  4 : label > 3s
          wavefunction #  5 : label > 3p
          wavefunction #  6 : label > 3d
          wavefunction #  7 : label > 4s
          wavefunction #  8 : label > 4p
          wavefunction #  9 : label >
        show wavefunction () or rho (else) ? >
        r max, r tic  > 5 0.5
          output file (postscript) > wf.ps
          output file (xmgr) >
    or, to just see, e.g., the 4s wavefunction, enter
        filwf > vs2p0d3.wfc
        filwfp >
          wavefunction #  1 : label > 4s
          wavefunction #  2 : label >
        show wavefunction () or rho (else) ? >
        r max, r tic  > 10 0.5
          output file (postscript) > wf4s.ps
          output file (xmgr) >
  • Now you are ready to run your own calculations. The input to ld1 is explained in Src/ld1.f:
          program ld1
    c-    Self-consistent calculation of electronic states for atoms
    c     (all-electron, AE) and pseudo-atoms (pseudopotential, PP)
    c     Atomic Rydberg units are used : e^2=2, m=1/2, hbar=1
    c     Atomic wavefunctions: Psi(r)=(chi(r)/r)Y_lm(r)
    c     The charge density rho(r) has spherical symmetry:
    c     rho(r) = \sum oc_nl*chi_nl(r)^2/(4\pi r^2)
    c     Gradient Correction added by Andrea Dal Corso
    c namelist "input":
    c          atom  = atomic symbol
    c          xmin  = mesh parameter                        default: -8.0
    c          dx    = mesh parameter                        default: 0.0500
    c          rmax  = outermost mesh point                  default: 60.0
    c                  The mesh is: r(i+1) = exp(xmin+i*dx)/zed (i=0,...,mesh)
    c                  Larger negative xmin and smaller dx yield denser and
    c                  more accurate grids.
    c          rel   = 1 (0): (non) scalar-relativistic calculation
    c                  default: non-relativistic if Z < 19 or for PP
    c          beta  = parameter for potential mixing (default: 0.5)
    c          rcut  = cut the exch-corr potential for r < rcut  (default:0.0)
    c                  may be useful for gradient-corrected functionals 
    c          config= electronic configuration. Acceptable format:
    c                  '1s2 2s2 2p6 3s0.5' or '[Ne] 3s0.5'
    c                  For PP calculations, the labels should be the
    c                  labels of corresponding AE wavefunctions:
    c                  the main quantum number is "renormalized" to the
    c                  correct value (1s for the first s state and so on)
    c          oldps = 0: new (PWSCF) PP format   (default)
    c                = 1: old PP format
    c          dft =   a string containing any of the following keywords
    c                (case-insensitive, in any nonconflicting arrangement):
    c                "lsd"    LSD calculation of the energy    (default:no)
    c                "sic"    Self-Interaction Correction      (default:no)
    c                "sla"    Slater  Exchange                 (default)
    c                "noexc"  No Exchange
    c                "nocor"  No Correlation
    c                "pz"     Perdew-Zunger correlation        (default)
    c                "gl"     Gunnarson-Lunqvist correlation
    c                "hl"     Hedin-Lunqvist correlation
    c                "pw"     Perdew-Wang correlation
    c                "vwn"    Vosko-Wilk-Nusair correlation
    c                "wig"    Wigner correlation
    c                "lyp"    Lee-Yang-Parr correlation
    c                "obz"    Ortiz-Ballone form for Perdew-Zunger corr.
    c                "obw"    Ortiz-Ballone form for Perdew-Wang corr.
    c                "nogcx"  No Gradient Correction on exchange    (default)
    c                "nogcc"  No Gradient Correction on correlation (default)
    c                "b88"    Becke88  grad-corr exchange (beta=0.0042)
    c                "p86"    Perdew86 grad-corr correlation
    c                "bp"     Becke88 + Perdew86
    c                "pw91"   Perdew-Wang 91 GGA
    c                "blyp"   Becke88 + Lee-Yang-Parr
    c                "pbe"    Improved GGA by Perdew-Burke-Erzenhof
    c                "revPBE" Revised PBE (Zhang and Yang)
    c                "hcth"   Cambridge functional after Handy & Co. (HCTH/120)
    c                "optx"   Handy's exchange functional (to be used with LYP)
    c          file_pseudo         input PP file             (default: ' ')
    c          file_wavefunctions  output file for wavefcts  (default: ' ')
    c          file_logder         output file for logarithmic derivatives
    c                                                        (default: 'logder')
    c                              Format readable by plotting program "xmgr"
    c                              A suffix is added to file_logder:
    c                              '.ae' for an AE calculation
    c                              '.pp' and '.kb' for a PP calculation
    c                              (the first is for semilocal PP, the second
    c                               for PP in the Kleinman-Bylander scheme.
    c                               The local part is read from the PP file)
    c          emin, emax = energy range for logarithmic derivatives calculation
    c                      (default: emin = emax = 0)
    c          ne         = number of energy steps for log. der. calculations
    c                       (default: 100)
    c          rd         = radius used to calculated log.der. (default: 0.0)
    c          ibound= 0 : free boundary conditions         (default)
    c                = 1 : atom in a box                    (obsolete)
    c          latt= 1(0): (no) Latter correction           (default: no)
    You should compare your calculations to the Atomic Reference Data
  • Here you can find the documentation of how to use the codes and on the background (all relevant formulas are in appedix A)

Erik Koch