The ABINIT Users Mailing List ( CLOSED )

[Deyu Lu <dylu@ucdavis.edu>]

Eventually I got everything right. The fact is that the results of
epsilon^-1 depend on the choice of origin. Hybertsen and Louie set the
origin in the middle of the two silicon atoms in the unit cell allowing
inversion symmetry, while I followed the tutorial to set the origin at
one of the silicon atoms. Given the displacement tau, the difference is
a phase factor exp(i (G-G') . tau), which solves all the disagreement I
had before. Enclosed is a summary of the comparison, and the second last
column of epsilon^-1 corresponds to the raw Abinit output times the
phase factor. The overall agreement is beyond 1.e-3.

Deyu Lu

On Sun, 2006-06-11 at 19:20 -0700, Deyu Lu wrote:
> 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.
> 
        Static dielectric matrix calculation for bulk silicon
-------------------------------------------------------------------------------
1. Pseudopotential: Troullier-Martins pseudopotentail provided by Abinit
packageConverged lattice constant is 10.217 a.u, experiment value is 10.259 
a.u.

2. Lattice constant: use the same value of 10.260 a.u. as in Hybertsen's paper

3. SCF calculation has a convergence criteria of 1.0d-10 in total energy

4. Cutoff energy: from 6.5 to 8.5 Hartree (8.5 Ha corresponding to 331 plane
waves)

5. 60 k points are used for BZ sampling

6. The dielectric matrix has the same energy cutoff as wavefunc, i.e., 331*331
 
7. A total number of 200 bands are used to construct chi0 and epsilon 

8. Computor time: 15.1 hours on 3.0 G Intel processor for 2 omega and 3 q
points (a-c ecut=8.5 Ha).

9. The abinit 5.1.2 source code is modified to use double precision data type 
for complex variables.

-----------------------------------------------------------------------------------------
RPA dielectric constant:
a) With PSP: 14si.pspnc, TM, rcs=rcp=rcd=2.0872718, lmax=2, llocal=2
                             ecut=6.5  ecut=8.5  ecut=10.5   Ref       EXP
dielectric constant:         12.6068   12.6365   12.6950     12.2      11.7
dielectric constant NoLF:    14.0443   14.0755   14.1410     13.61

b) With PSP: Si.d.cpi generated from fhi98pp, 
TM, rcs=1.7039185, rcp=1.8786063, rcd=2.0212776, lmax=2, llocal=2
                             ecut=6.5  ecut=8.5  ecut=10.5   Ref       EXP
dielectric constant:         12.3043   12.4104               12.2      11.7
dielectric constant NoLF:    13.7126   13.8276               13.61

c) With PSP: Si.s.cpi generated from fhi98pp, 
TM, rcs=1.7039185, rcp=1.8786063, rcd=2.0212776, lmax=2, llocal=0
                             ecut=6.5  ecut=8.5  ecut=10.5   Ref       EXP
dielectric constant:         12.2280                         12.2      11.7
dielectric constant NoLF:    13.6275                         13.61
----------------------------------------------------------------------------------------
With PSP: Si.s.cpi generated from fhi98pp, llocal=0

a) TEST CHARGE epsilon^-1(q->0), q=5.e-6 b2 + 5.e-6 b3 = (1.e-5,0,0) 2\pi/a
               G              G'         ecut=6.5(complex)       real      Ref
     1    0    0    0    0    0    0    ( 0.077328  0.000000)  0.077328   
0.077
     2    1    1    1    0    0    0    (-0.015598 -0.015598)  0.022059   
0.022
     3    1    1    1    1    1    1    ( 0.522503  0.000000)  0.522503   
0.522
     4   -1    1    1    1    1    1    ( 0.000000  0.011301)  0.011301   
0.011
     5    1   -1    1    1    1    1    ( 0.000000 -0.001285) -0.001285  
-0.001
     6    1   -1    1   -1    1    1    (-0.006559  0.000000) -0.006559  
-0.006
     7    1    1   -1    1   -1    1    ( 0.006027  0.000000)  0.006027   
0.007
     8    1   -1   -1   -1    1    1    ( 0.000000  0.065146)  0.065146   
0.066

b) TEST CHARGE epsilon^-1(X)            
     1    0    0    0    0    0    0    ( 0.274256  0.000000)  0.274256   
0.2749 
     2    0    0   -2    0    0    0    ( 0.000000  0.000000)  0.000000   0.0
     3    1    1   -1    0    0    0    (-0.018502  0.018502) -0.026166  
-0.0259 
     4    1    1   -1    1    1   -1    ( 0.413871  0.000000)  0.413871   
0.4138
     5    1    1   -1    1   -1   -1    ( 0.000000  0.000056) -0.000056   0.0
     6    1   -1   -1   -1    1   -1    ( 0.044626  0.000000)  0.044626   
0.0456
c) TEST CHARGE epsilon^-1(L)
     1    0    0    0    0    0    0    ( 0.282390  0.000000)  0.282390   
0.2828
     2   -1   -1   -1    0    0    0    ( 0.050416 -0.050416) -0.071299  
-0.0719
     3   -1   -1    1    0    0    0    (-0.028626 -0.028626) -0.040483  
-0.0400 
     4   -1   -1    1   -1   -1   -1    ( 0.000000  0.001785) -0.001785  
-0.0019
     5   -1   -1    1   -1   -1    1    ( 0.482296  0.000000)  0.482296   
0.4822
     6   -1    1   -1   -1   -1    1    ( 0.008150  0.000000)  0.008150   
0.0089
     7   -2    0    0   -1   -1    1    (-0.037000 -0.037001) -0.052327  
-0.0520
     8   -2    0    0    1   -1   -1    (-0.001387 -0.001386) -0.001961  
-0.0013
-------------------------------------------------------------------------------
Ref: Phys. Rev. B, 35:5585, 1987.
EXP: Ref. 48 in Hybertsen and Louie's paper