The ABINIT Users Mailing List ( CLOSED )

[dylu@ucdavis.edu]

Hi, as a follow up of my previous post on epsilon^-1 calculation, I got some 
results close to Hybertsen & Louie's 1987 paper on bulk silicon, and
it is a quite long post.
------------------------------------------------------------------------------
With PSP: Si.s.cpi generated from fhi98pp, llocal=0
a) TEST CHARGE epsilon^-1(q->0)
               G              G'        ecut=6.5(A)ecut=6.5(B)           Ref
     1    0    0    0    0    0    0     0.077328  0.077328             0.077
     2    1    1    1    0    0    0    -0.020927 -0.015598             0.022
     3    1    1    1    1    1    1     0.527537  0.522503             0.522
     4   -1    1    1    1    1    1     0.000000  0.000000             0.011
     5    1   -1    1    1    1    1     0.000000  0.000000            -0.001
     6    1   -1    1   -1    1    1    -0.001524 -0.006559            -0.006
     7    1    1   -1    1   -1    1     0.003510  0.006027             0.007
     8    1   -1   -1   -1    1    1     0.000000  0.000000             0.066
A: q=1.e-5 b1 + 2.e-5 b2 + 3.e-5 b3
B: q=           5.e-6 b2 + 5.e-6 b3 = (1.e-5, 0 , 0) 2\pi/a
-----------------------------------------------------------------------------
The above results are obtained on parameters similar to the reference.

The G, G' values are defined in reciprocal coordinates, while they are in 
cartesian coordinates in the reference.

The pseudopotential are generated by fhi98pp with the following input:
14.00  3  2  8  0.00 : z  nc  nv iexc rnlc
    1  0   2.00      : n  l   f
    2  0   2.00
    2  1   6.00
    3  0   2.00
    3  1   2.00
2  t                 : lmax  s_pp_def

so that the cutoffs are  rcs=1.7039185, rcp=1.8786063, rcd=2.0212776 (a.u.),  
closer to the ones used in the reference, i.e., rcs=1.58, rcp=1.93, rcd=1.93 
than in 14si.pspnc. The 3s orbital is chosen as local channel, same as the 
reference.

To allow the output of test charge, not RPA, epsilon^-1, the following 
definition is made in screen.F90:
logical,parameter :: lda=.true.,rpa=.false.,testelectron=.false.

Finally, in order to approach Gama point from X axis as in the reference,
q(1) is reset to the following in findq.F90.
  do j=1,nq
   if(itst0(q(:,j))==1) then
!     q(1,j)=0.000010
     q(1,j)=0.d0
!     q(2,j)=0.000020
     q(2,j)=5.d-6
!     q(3,j)=0.000030
     q(3,j)=5.d-6
   end if

After all these modifications, I find that in table a elements 1, 3, 6, 7 are 
in very good agreement, better than 1.e-3. However, elements 2, 4, 5, 8 are 
not.
Element 2 even suffers an opposite sign.  Elements 4, 5, 8 are simply zero. 
If I output both real and imaginary parts of epsilon^-1 for element 4, 5, 8, 
I found that they are:
     4    1    0    0    1    1    1    0.000000    0.011301
     5    0    1    0    1    1    1    0.000000   -0.001285
     8   -1    0    0    1    0    0    0.000000    0.065146
Then the imaginary part actually corresponds to the number in the reference, 
again with an acurray no less than than 1.e-3. To make sure that the numbers 
I read out of the SCR file are correct, I double-check the log file. The 
numbers are identical. In the output of the log file, there is an interesting 
pattern:
Real (epsilon^-1(q->0, omega=0, G=2:5, G'=6:9))=0, which doesn't look quite 
right to me. To confirm this observation, I did a similar calculation in 
Abinit
4.6.5 with default parameters of q(1) for RPA epsilon^-1. Again, I found that
Real (epsilon^-1(q->0, omega=0, G=2:5, G'=6:9))=0.

I guess that it is not a coincidence, and I trust the result in Hybertsen and 
Louie's paper, since it was reproduced later by Engel and Farid in 1992 (PRB 
46: 15812) with a different numerical method.

Since I'm trying some method development based on GW part of Abinit, to get 
everything right in epsilon^-1 is critical. I'm digging into the code these 
days to see if there is anything wrong, but I would be really grateful if 
someone with more knowledge about the code could help me out of it.

To make the story complete, I list below the results for X and L. Still, some 
are close, and some are not. Like the case of Gamma point (element 2 in Table 
a), wing elements (element 2 in Table c) could have an opposite sign.

An input file I used is attached at the end.


Many thanks
Deyu Lu

b) TEST CHARGE epsilon^-1(X)
     1    0    0    0    0    0    0     0.274256                       0.2749
     2    0    0   -2    0    0    0     0.000000                       0.0
     3    1    1   -1    0    0    0    -0.018502                      -0.0259
     4    1    1   -1    1    1   -1     0.413871                       0.4138
     5    1    1   -1    1   -1   -1     0.000000                       0.0
     6    1   -1   -1   -1    1   -1     0.044626                       0.0456
c) TEST CHARGE epsilon^-1(L)
     1    0    0    0    0    0    0     0.282390                       0.2828
     2   -1   -1   -1    0    0    0     0.050416                      -0.0719
     3   -1   -1    1    0    0    0    -0.028626                      -0.0400
     4   -1   -1    1   -1   -1   -1     0.000000                      -0.0019
     5   -1   -1    1   -1   -1    1     0.482296                       0.4822
     6   -1    1   -1   -1   -1    1     0.008150                       0.0089
     7   -2    0    0   -1   -1    1    -0.037000                      -0.0520
     8   -2    0    0    1   -1   -1    -0.001387                      -0.0013
-------------------------------------------------------------------------------
ndtset     2

# Definition of parameters for the calculation of the KSS file
nbandkss1   -1         # Number of bands in KSS file (-1 means the maximum 
possible)
nband1      9         # Number of (occ and empty) bands to be computed
istwfk1     60*1

# Calculation of the screening (epsilon^-1 matrix)
optdriver2  3        # Screening calculation
getkss2     -1       # Obtain KSS file from previous dataset
nband2      200       # Bands to be used in the screening calculation
ecutwfn2    6.5      # Cut-off energy of the planewave set to represent the 
wavefunctions
ecuteps2    6.5      # Cut-off energy of the planewave set to represent the 
dielectric matrix
ppmfrq2    16.7 eV  # Imaginary frequency where to calculate the screening
nqptdm2    3
qptdm2     0.000000    0.000005    0.000005      #Gama
           0.500000    0.500000    0.500000      #L
          -0.500000   -0.500000    0.000000      #X

# Data common to the three different datasets

# Definition of the unit cell: fcc
acell  3*10.260        # This is equivalent to   10.260 10.260 10.260
rprim  0.0  0.5  0.5   # FCC primitive vectors (to be scaled by acell)
       0.5  0.0  0.5
       0.5  0.5  0.0

# Definition of the atom types
ntypat  1         # There is only one type of atom
znucl 14          # 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, the only type is Silicon.

# Definition of the atoms
natom 2           # There are two atoms
typat  1 1        # They both are of type 1, that is, Silicon.
xred              # Reduced coordinate of atoms
      0.0  0.0  0.0
      0.25 0.25 0.25


# Definition of the k-point grid
kptopt  1            # Option for the automatic generation of k points,
nkpt    60
ngkpt   8 8 8
nshiftk 4
shiftk  0.5 0.5 0.5  # These shifts will be the same for all grids
        0.5 0.0 0.0
        0.0 0.5 0.0
        0.0 0.0 0.5

# Use only symmorphic operations
symmorphi 0

# Definition of the planewave basis set (at convergence 16 Rydberg 8 Hartree)
ecut 6.5          # Maximal kinetic energy cut-off, in Hartree

# Definition of the SCF procedure
nstep   10        # Maximal number of SCF cycles
toldfe  1.0d-10   # Will stop when this tolerance is achieved on total energy
diemac  12.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.
                  # Here, we follow the prescription for bulk silicon.