The ABINIT Users Mailing List ( CLOSED )

["Fernando D. Vila" <fer@tiziano.phys.washington.edu>]

Hello everyone..

  Some time ago I asked how to get the full dynamical matrix in real 
space. Thanks to the very nice tutorial that was recently published for 
the 4.6 version, I've been able to get all the values I needed. My 
interest has been in the calculation of EXAFS Debye-Waller factors. I've 
ben able to interface abinit's output with a small code we have to 
calculate them. All the resutls have been quite decent until I tried bcc 
Fe. I have calculated fcc metals (Cu, Ag, Au, Pb and Pt) and 
insulators/semiconductors like C, Si, Ge, GaAs. All this systems are 
essentially fcc, so I didn't have to modify the tutorial very much. Beside 
some problems with the results with the heavier metals (where the DW 
factors are a bit too small, meaning the force constants are a bit too 
large), all the values are within 10% of experiment.

Now, for Fe I had to do some changes in the input and my DW factors are a 
factor of 2.2 too large. I have traced back the problem to the second 
derivatives, ruling out a problem with my codes and with mrgddb and 
anaddb. I know that the highest freq for Fe is about 9.5 THz. I'm getting 
a value that is about 1.5 times smaller, which agrees well with the factor 
of 2.2 in the 2nd derivatives. Since I'm at a complete loss as to what 
might be wrong, I'm sending the input file to see if anyone can spot my 
(probably obvious) mistake. Finally, I'm using the fhi psp from the abinit 
website.

Best regards and thanks in advance..

PS: Any idea why the the heavier metals are stiffer?

Here is the input:
---------------------------------------------------------------------------
# Crystalline AlAs : computation of the total energy
#
# WARNING!!! Modified by Fer to do Fe
#            This files also include keywords from tutorial 4 for Al

#Definition of occupation numbers
occopt 4
tsmear 0.04

   ndtset  10

#Set 1 : ground state self-consistency

  getwfk1   0            # Cancel default
  kptopt1   1            # Automatic generation of k points, taking
                         # into account the symmetry
    nqpt1   0            # Cancel default
  tolvrs1   1.0d-18      # SCF stopping criterion (modify default)
  rfphon1   0            # Cancel default

#Q vectors for all datasets

#Complete set of symmetry-inequivalent qpt chosen to be commensurate
# with kpt mesh so that only one set of GS wave functions is needed.
#Generated automatically by running GS calculation with kptopt=1,
# nshift=0, shiftk=0 0 0 (to include gamma) and taking output kpt set
# file as qpt set. Set nstep=1 so only one iteration runs.

     nqpt   1            # One qpt for each dataset (only 0 or 1 allowed)
                         # This is the default for all datasets and must
                         #  be explicitly turned off for dataset 1.

     qpt2   0.00000000E+00  0.00000000E+00  0.00000000E+00
     qpt3   0.00000000E+00  0.00000000E+00  0.00000000E+00
     qpt4   2.50000000E-01  0.00000000E+00  0.00000000E+00
     qpt5   5.00000000E-01  0.00000000E+00  0.00000000E+00
     qpt6   2.50000000E-01  2.50000000E-01  0.00000000E+00
     qpt7   2.50000000E-01  2.50000000E-01  2.50000000E-01
     qpt8  -2.50000000E-01  2.50000000E-01  2.50000000E-01
     qpt9   5.00000000E-01  5.00000000E-01  2.50000000E-01
     qpt10  5.00000000E-01  5.00000000E-01  5.00000000E-01

#Set 2 : Response function calculation of d/dk wave function

    iscf2   -3         # Need this non-self-consistent option for d/dk
  kptopt2   2          # Modify default to use time-reversal symmetry
  rfphon2   0          # Cancel default
  rfelfd2   2          # Calculate d/dk wave function only
  tolvrs2   0.0        # Cancel default for d/dk
  tolwfr2   1.0d-22    # Use wave function residual criterion instead

#Set 3 : Response function calculation of Q=0 phonons and electric field pert.

  getddk3   2          # d/dk wave functions from last dataset
  kptopt3   2          # Modify default to use time-reversal symmetry
  rfelfd3   3          # Electric-field perturbation response only

#Sets 4-10 : Finite-wave-vector phonon calculations (defaults for all 
datasets)

   getwfk   1          # Use GS wave functions from dataset1
   kptopt   3          # Need full k-point set for finite-Q response
   rfphon   1          # Do phonon response
  rfatpol   1 1        # Treat displacements of all atoms
    rfdir   1 1 1      # Do all directions (symmetry will be used)
   tolvrs   1.0d-8     # This default is active for sets 3-10

#######################################################################
#Common input variables

#Definition of the unit cell
acell 3*5.42351397873  # This is equivalent to   10.61 10.61 10.61
rprim -0.5  0.5  0.5   # In lessons 1 and 2, these primitive vectors
       0.5 -0.5  0.5   # (to be scaled by acell) were 1 0 0  0 1 0  0 0 1
       0.5  0.5 -0.5   # that is, the default.

#chkprim 0

#Definition of the atom types
ntypat 1               # There are two types of atom
znucl 26               # The keyword "znucl" refers to the atomic number of 
the
                       # possible type(s) of atom. The pseudopotential(s)
                       # mentioned in the "files" file must correspond
                       # to the type(s) of atom. Here, type 1 is the Aluminum,
                       # type 2 is the Arsenic.

#Definition of the atoms
natom 1               # There are two atoms
typat                  # The first is of type 1 (Al), the second is of type 2 
(As).
 1
xred                   # This keyword indicate that the location of the atoms
                       # will follow, one triplet of number for each atom
   0.0  0.0  0.0       # Triplet giving the REDUCED coordinate of atom 1.

#Gives the number of band, explicitely (do not take the default)
nband 12             # For an insulator (if described correctly as an 
insulator
                     # by DFT), there is no need to include conduction bands
                     # in response-function calculations

#Exchange-correlation functional
ixc 11            # PBE GGA

#Definition of the planewave basis set
ecut   60.0           # Maximal kinetic energy cut-off, in Hartree

#Definition of the k-point grid
ngkpt   4  4  4
nshiftk 2
shiftk
        0.25  0.25  0.25
       -0.25 -0.25 -0.25

#Definition of the SCF procedure
iscf     5
nstep   80             # Maximal number of SCF cycles
#diemac 16.0           # Although this is not mandatory, it is worth to
                       # precondition the SCF cycle. The model dielectric
                       # function used as the standard preconditioner
                       # is described in the "dielng" input variable section.
                       # The dielectric constant of AlAs is smaller that the
                       # one of Si (=12).

Ubi dubium ibi libertas.
*******************************************************************************
Fernando D. Vila                Voice    (206)543-9697
Department of Physics           Fax      (206)685-0635
University of Washington        E-mail   fdv@u.washington.edu
Seattle, WA 98195, USA          WWW      http://faculty.washington.edu/fdv
*******************************************************************************