The ABINIT Users Mailing List ( CLOSED )

[fer@tiziano.phys.washington.edu]

Some time ago I sent a message asking for some help with a calculation of 
phonons in Fe. The problem was that the results were pretty bad. As I said at 
the time, I am a "molecules kind of guy" and completely forgot that Fe is 
ferromagnetic.. :-P After I realized that I started trying to do 
spin-polarized calculations and was able to get the ground state of Fe to 
converge correctly. Now I'm trying to do calculations for each of the 4x4x4 
q-points required to get the whole dynamical matrix and I'm having some 
convergence problems.

For some of the q-points the second derivatives just "explode" and I can't 
figure out how to get the whole set to converge. I'm using T-M FHI 
pseudopotential as provided in the ABINIT website. Here is a typical 
"explosion":

     iter   2DEtotal(Ha)       deltaE(Ha) residm    vres2
-ETOT  1   81763.971689883     7.647E+04 1.590E+01 6.956E+09
 ETOT  2   2062994.7382198     1.981E+06 3.385E+01 1.630E+11
 ETOT  3   18402.974058305    -2.045E+06 2.871E+01 1.375E+09
 ETOT  4   171.52011047491    -1.823E+04 2.986E-01 1.301E+07
 ETOT  5   2065.5006653301     1.894E+03 3.534E-02 1.655E+08
 ETOT  6   8421.5060419704     6.356E+03 2.401E-01 6.563E+08
 ETOT  7   82.544911882740    -8.339E+03 1.523E-01 8.027E+06
 ETOT  8  -9.3896549601841    -9.193E+01 1.350E-02 2.633E+06
 ETOT  9  -53.435097572993    -4.405E+01 8.114E-03 6.720E+05
 ETOT 10  -111.41793571227    -5.798E+01 6.568E-02 2.150E+06
 ETOT 11  -56.621827824279     5.480E+01 1.006E-03 5.641E+06
 ETOT 12  -254.30296451650    -1.977E+02 1.324E-01 2.588E+06
 ETOT 13   75.421964034640     3.297E+02 3.166E-02 3.569E+07
 ETOT 14   216.68416962394     1.413E+02 1.148E-01 5.954E+07
 ETOT 15   7081.0318423753     6.864E+03 3.502E-01 6.195E+08
 ETOT 16   11951.897429366     4.871E+03 7.451E-01 1.027E+09
 ETOT 17   123917.98415428     1.120E+05 4.464E+00 9.869E+09
 ETOT 18   198390.86476296     7.447E+04 6.242E+00 1.576E+10
 ETOT 19   1963683.9673588     1.765E+06 3.848E+01 1.546E+11
 ETOT 20   3112599.3949469     1.149E+06 7.504E+01 2.450E+11
 ETOT 21   30417037.058737     2.730E+07 5.228E+02 2.391E+12
 ETOT 22   48329453.803039     1.791E+07 1.099E+03 3.799E+12
 ETOT 23   471752329.56943     4.234E+08 7.773E+03 3.708E+13
 ETOT 24   748810687.23925     2.771E+08 1.714E+04 5.885E+13
 ETOT 25   7311483084.4625     6.563E+09 1.195E+05 5.746E+14
 ETOT 26   11606974126.738     4.295E+09 2.650E+05 9.122E+14
 ETOT 27   113305917642.26     1.017E+11 1.858E+06 8.904E+15
 ETOT 28   179877811467.53     6.657E+10 4.105E+06 1.414E+16
 ETOT 29   1755990209022.8     1.576E+12 2.875E+07 1.380E+17
 ETOT 30   2787721847839.4     1.032E+12 6.365E+07 2.191E+17
 ETOT 31   27213805885607.     2.443E+13 4.457E+08 2.139E+18
 ETOT 32   43204197949368.     1.599E+13 9.862E+08 3.395E+18
 ETOT 33  4.21752761389320E+14 3.785E+14 6.908E+09 3.314E+19
 ETOT 34  6.69578444121069E+14 2.478E+14 1.529E+10 5.262E+19
 ETOT 35  6.53625338303775E+15 5.867E+15 1.070E+11 5.137E+20
 ETOT 36  1.03772000684718E+16 3.841E+15 2.369E+11 8.155E+20
 ETOT 37  1.01298240614436E+17 9.092E+16 1.659E+12 7.961E+21
 ETOT 38  1.60827627813465E+17 5.953E+16 3.671E+12 1.264E+22
 ETOT 39  1.56991748691401E+18 1.409E+18 2.571E+13 1.234E+23


and so on until ~1E+42.

Here is my input file (please don't pay to much attention to the comments 
since they correspond to the tutorial file, I kept them just for reference):

###########################################################################
occopt 4
tsmear 0.01
spinat   0.0 0.0 4.0
nsppol   2
                                                                              
  
   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     3
nstep   80             # Maximal number of SCF cycles
#diemac 1000000.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).
###############################################################################

Thanks in advance for any help you can give me.

Fer.

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