RF2 tutorial calculation with monolayer Graphene(PAW+GGA)

Running Abinit Ground state
1

Hi everyone,
I recently encountered a SCF issue when I tried to conduct a RF2 tutorial calculation with monolayer Graphene using PAW+GGA method. Here’s the output file result:

— !BeginCycle
iteration_state: {dtset: 1, }
solver: {iscf: 7, nstep: 100, nline: 4, wfoptalg: 10, }
tolerances: {toldfe: 1.00E-10, }
…

 iter   2DEtotal(Ha)        deltaE(Ha) residm    vres2

-ETOT 1 -15963898041793. -1.596E+13 3.894E+09 4.319E+30
ETOT 2 -1.65815441111300E+20 -1.658E+20 6.371E+15 1.723E+37
ETOT 3 -2.88800694657322E+20 -1.230E+20 7.119E+16 2.793E+36
ETOT 4 -8.62042204578756E+21 -8.332E+21 5.180E+16 3.489E+36
ETOT 5 -9.39735751796615E+21 -7.769E+20 7.861E+16 2.585E+36
ETOT 6 -1.30662913848526E+24 -1.297E+24 1.135E+18 2.683E+36
ETOT 7 -3.26168574721839E+22 1.274E+24 4.427E+18 4.032E+35
ETOT 8 -1.54237497000655E+22 1.719E+22 3.792E+17 1.107E+35
ETOT 9 -1.08364952902480E+22 4.587E+21 8.500E+16 2.897E+34
ETOT 10 -1.00337626489670E+29 -1.003E+29 1.888E+23 4.363E+40
ETOT 11 -2.60555374061286E+28 7.428E+28 8.096E+23 5.246E+40
ETOT 12 -6.07151260979212E+27 1.998E+28 1.336E+23 9.959E+40
ETOT 13 -9.65089236989867E+27 -3.579E+27 1.324E+23 2.135E+41
ETOT 14 -1.95027592176441E+28 -9.852E+27 3.384E+22 8.789E+40
ETOT 15 -7.53479189198434E+27 1.197E+28 9.007E+22 6.459E+40
ETOT 16 -1.12229421714361E+29 -1.047E+29 9.649E+22 7.904E+40
ETOT 17 -5.31687921968093E+28 5.906E+28 7.097E+23 8.962E+40
ETOT 18 -3.57434390752007E+34 -3.574E+34 8.961E+28 1.753E+46
ETOT 19 -3.13606942971243E+34 4.383E+33 1.396E+29 1.259E+46
ETOT 20 -5.04758763983600E+34 -1.912E+34 1.125E+29 1.950E+47

Here’s my input file:

Chksymbreak 0

occopt 4
tsmear 0.01

ndtset 4
prtden 0
prteig 2
#Q vectors for all datasets

#Complete set of symmetry-inequivalent qpt chosen to be commensurate

nqpt 1 # One qpt for each dataset (only 0 or 1 allowed)
qptopt 1 # activate determining qpts in IBZ with symmetry
ngqpt 4 4 1 # these variables mirror those used for the kpt
nshiftq 1 # mesh below, but with no shift. In this way,
shiftq 0 0 0 # qpts coherent with the kpts are constructed
# the number of resulting q pts is pretdetermined
# through use of the abitk tool

iqpt: 1
iqpt+ 1

#Set 1 : iqpt 1 is the gamma point, so Q=0 phonons and electric field pert.
getddk1 98 # d/dk wave functions
kptopt1 2 # Modify default to use time-reversal symmetry
rfelfd1 3 # Electric-field perturbation response
# (in addition to default phonon)

#Sets 2-8 : Finite-wave-vector phonon calculations (defaults for all datasets)

getwfk 99 # Use GS wave functions
kptopt 3 # Need full k-point set for finite-Q response
rfphon 1 # Do phonon response
toldfe 1.0d-10 # Converge on potential residual

prtwf 0
prtpot 0

#######################################################################
acell 4.6597825725E+00 4.6597825725E+00 1.2999393442E+01
rprim sqrt(0.75) -0.5 0.0 # FCC primitive vectors (to be scaled by acell)
sqrt(0.75) 0.5 0.0
0.0 0.0 1.0

#Definition of the atom types
ntypat 1 # There is only one type of atom
znucl 6 # 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.
pp_dirpath “$ABI_PSPDIR” # This is the path to the directory were
# pseudopotentials for tests are stored
pseudos “PAW/C.GGA_PBE-JTH.xml”
# Name and location of the pseudopotential

#Definition of the atoms
natom 2 # There are two atoms
typat 1 1 # They both are of type 1, that is, Silicon.
xred # This keyword indicate that the location of the atoms
# will follow, one triplet of number for each atom
1/3 1/3 0.0 # Triplet giving the REDUCED coordinate of atom 1.
2/3 2/3 0.0 # Triplet giving the REDUCED coordinate of atom 2.

#Definition of the planewave basis set
ecut 20.0 # Maximal kinetic energy cut-off, in Hartree
pawecutdg 40.0
pawovlp -1

              # from previous wavefunctions, transferred from the old
              # to the new k-points.

#Definition of the SCF procedure
nstep 100 # Maximal number of SCF cycles
# 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.

ngkpt 4 4 1
nshiftk 1
shiftk 0 0 0

pawxcdev 0

Here’s my input file generating ddk and wfk files:

occopt 4
tsmear 0.01

ndtset 2

kptopt1 1 # Automatic generation of k points, taking
# into account all symmetries
toldfe1 1.0d-10 # SCF stopping criterion

kptopt2 2 # DDK can use only time reveral symmetry
getwfk2 -1 # require ground state wavefunctions from previous run
rfelfd2 2 # activate DDK perturbation
iscf2 -3 # this is a non-self-consistent calculation
tolwfr2 1.0D-20 # tight convergence on wavefunction residuals

#Definition of the unit cell
acell 4.6597825725E+00 4.6597825725E+00 1.2999393442E+01
rprim sqrt(0.75) -0.5 0.0 # FCC primitive vectors (to be scaled by acell)
sqrt(0.75) 0.5 0.0
0.0 0.0 1.0

#Definition of the atom types
ntypat 1 # There is only one type of atom
znucl 6 # 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.
pp_dirpath “$ABI_PSPDIR” # This is the path to the directory were
# pseudopotentials for tests are stored
pseudos “PAW/C.GGA_PBE-JTH.xml”
# Name and location of the pseudopotential

#Definition of the atoms
natom 2 # There are two atoms
typat 1 1 # They both are of type 1, that is, Silicon.
xred # This keyword indicate that the location of the atoms
# will follow, one triplet of number for each atom
1/3 1/3 0.0 # Triplet giving the REDUCED coordinate of atom 1.
2/3 2/3 0.0 # Triplet giving the REDUCED coordinate of atom 2.

ecut 20.0 # Maximal kinetic energy cut-off, in Hartree
pawecutdg 40.0
pawovlp -1

              # from previous wavefunctions, transferred from the old
              # to the new k-points.

#Definition of the SCF procedure
nstep 100 # Maximal number of SCF cycles
# 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.

nshiftk 1
shiftk 0.0 0.0 0.0 # These shifts will be the same for all grids
ngkpt 4 4 1

pawxcdev 0

I have stucked in the same page for almost a month, so I sincerely appreciate the help!

Best and thank you guys,

Lucas Lu

2

Hi,

Is your system a metal? (zero gap)? from your choice of occopt etc it would appear so. DDK will converge slowly if at all for a metal, and anyway, the reason you’d run DDK is so that you can then run the electric field perturbation, which doesn’t really make sense to run on a metal. You might try the same calculation with ONLY the phonon perturbation.

3

I see tsmear in some of the inputs. If you include that and occopt >= 3 systematically, then ddk for metals actually can converge very quickly. It’s main use is to get the Fermi velocity in EPH, though for graphene the E field response is complex and non adiabatic in reality, due to its strange band structure.

Your k grid is super unconverged, however. Even with a massive smearing you will need 16x16 or 32x32 k. Please check the literature

Second check you can do is to run LDA: in some cases the GGA perturbed potentials are not well bounded in the vacuum and can get unphysically large. If this works tell us and send on the 2 input files for lda and gga - we’ll make a test out of it.

M.

4

Hi,

I tried a run with LDA pseudo, and etot becomes unphysically large when I try to run phonon response in dtset 3, 4, 5. Also the electrical field perturbation run, dtset 6 wasn’t converged at all:

Here’s the input file:

ndtset 5
#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
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.0000E+00 0.0000E+00 0.0000E+00 # wtk 6.2500E-02
qpt3 2.5000E-01 0.0000E+00 0.0000E+00 # wtk 3.7500E-01
qpt4 5.0000E-01 0.0000E+00 0.0000E+00 # wtk 1.8750E-01
qpt5 5.0000E-01 2.5000E-01 0.0000E+00 # wtk 3.7500E-01
qpt6 0 0 0

#Set 2: Response function for q=0, print out the ddk
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
getwfk2 1 # Calculate d/dk wave function only
#Set 6: Response function calculation of Q=0 phonons and electric field pert.

getwfk6 0
getddk6 2 # d/dk wave functions from last dataset
kptopt6 2 # Modify default to use time-reversal symmetry
rfelfd6 3 # Electric-field perturbation response only
rfphon6 0

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

getwfk3 1
kptopt3 3
rfphon3 1
#rfatpol3 1 2
#rfdir3 1 1 1

getwfk4 1
kptopt4 3
rfphon4 1
#rfatpol4 1 2
#rfdir4 1 1 1

getwfk5 1
kptopt5 3
rfphon5 1
#rfatpol5 1 2
#rfdir5 1 1 1

toldfe 1.0d-10 # This default is active for sets 3-10
#######################################################################
#Common input variables
#Definition of the planewave basis set
acell 4.6597825725E+00 4.6597825725E+00 1.2999393442E+01
rprim sqrt(0.75) -0.5 0.0 # FCC primitive vectors (to be scaled by acell)
sqrt(0.75) 0.5 0.0
0.0 0.0 1.0

#Definition of the atom types
ntypat 1 # There is only one type of atom
znucl 6 # 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.
pp_dirpath “$ABI_PSPDIR” # This is the path to the directory were
# pseudopotentials for tests are stored
pseudos “LDA/C.LDA_PW-JTH.xml”
# Name and location of the pseudopotential

#Definition of the atoms
natom 2 # There are two atoms
typat 1 1 # They both are of type 1, that is, Silicon.
xred # This keyword indicate that the location of the atoms
# will follow, one triplet of number for each atom
1/3 1/3 0.0 # Triplet giving the REDUCED coordinate of atom 1.
2/3 2/3 0.0 # Triplet giving the REDUCED coordinate of atom 2.

#Definition of the planewave basis set
ecut 20.0 # Maximal kinetic energy cut-off, in Hartree
pawecutdg 40.0
pawovlp -1

              # from previous wavefunctions, transferred from the old
              # to the new k-points.

#Definition of the SCF procedure
nstep 50 # Maximal number of SCF cycles
diemac 12 # Although this is not mandatory, i
#GS calculation
ngkpt 4 4 1 # should be 8 8 1
nshiftk 1
shiftk 0.0 0.0 0.0 # This choice preserves the hexagonal symmetry of the grid,
pawxcdev 0

The output file:

== DATASET 3 ==================================================================

  • mpi_nproc: 1, omp_nthreads: -1 (-1 if OMP is not activated)

— !DatasetInfo
iteration_state: {dtset: 3, }
dimensions: {natom: 2, nkpt: 16, mband: 5, nsppol: 1, nspinor: 1, nspden: 1, mpw: 1047, }
cutoff_energies: {ecut: 20.0, pawecutdg: 40.0, }
electrons: {nelect: 8.00000000E+00, charge: 0.00000000E+00, occopt: 1.00000000E+00, tsmear: 1.00000000E-02, }
meta: {optdriver: 1, rfphon: 1, }
…

mkfilename : getwfk/=0, take file _WFK from output of DATASET 1.

Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):
R(1)= 4.0354901 -2.3298913 0.0000000 G(1)= 0.1239007 -0.2146023 0.0000000
R(2)= 4.0354901 2.3298913 0.0000000 G(2)= 0.1239007 0.2146023 0.0000000
R(3)= 0.0000000 0.0000000 12.9993934 G(3)= 0.0000000 0.0000000 0.0769267
Unit cell volume ucvol= 2.4444718E+02 bohr^3
Angles (23,13,12)= 9.00000000E+01 9.00000000E+01 6.00000000E+01 degrees
setup1 : take into account q-point for computing boxcut.

Coarse grid specifications (used for wave-functions):

getcut: wavevector= 0.2500 0.0000 0.0000 ngfft= 20 20 54
ecut(hartree)= 20.000 => boxcut(ratio)= 2.06435

Fine grid specifications (used for densities):

getcut: wavevector= 0.2500 0.0000 0.0000 ngfft= 30 30 80
ecut(hartree)= 40.000 => boxcut(ratio)= 2.16202

==> initialize data related to q vector <==

The list of irreducible perturbations for this q vector is:
1) idir= 1 ipert= 1
2) idir= 2 ipert= 1
3) idir= 3 ipert= 1

================================================================================

The perturbation idir= 1 ipert= 2 is
symmetric of a previously calculated perturbation.
So, its SCF calculation is not needed.

The perturbation idir= 2 ipert= 2 is
symmetric of a previously calculated perturbation.
So, its SCF calculation is not needed.

The perturbation idir= 3 ipert= 2 is
symmetric of a previously calculated perturbation.
So, its SCF calculation is not needed.

Perturbation wavevector (in red.coord.) 0.250000 0.000000 0.000000
Perturbation : displacement of atom 1 along direction 1
Found 2 symmetries that leave the perturbation invariant.
symkpt : not enough symmetry to change the number of k points.

Initialisation of the first-order wave-functions :
ireadwf= 0

— !BeginCycle
iteration_state: {dtset: 3, }
solver: {iscf: 7, nstep: 50, nline: 4, wfoptalg: 10, }
tolerances: {toldfe: 1.00E-10, }
…

 iter   2DEtotal(Ha)        deltaE(Ha) residm    vres2

-ETOT 1 -375249475043.77 -3.752E+11 1.193E+08 1.497E+23
ETOT 2 -1.23825221251966E+19 -1.238E+19 1.621E+15 4.942E+30
ETOT 3 -3.46047503906718E+18 8.922E+18 1.466E+15 1.925E+30
ETOT 4 -4.50256769154467E+17 3.010E+18 2.545E+14 7.200E+29
ETOT 5 -2.19757504197024E+17 2.305E+17 1.151E+14 3.213E+29
ETOT 6 -3.09189085164588E+16 1.888E+17 2.317E+13 2.906E+28
ETOT 7 -2.10394518903089E+16 9.879E+15 1.071E+13 1.468E+28
ETOT 8 -4.57356932783485E+15 1.647E+16 1.976E+12 3.970E+27
ETOT 9 -1.10918578148256E+16 -6.518E+15 1.598E+12 6.031E+27
ETOT 10 -6.31610758001888E+21 -6.316E+21 2.794E+18 2.520E+33
ETOT 11 -6.38313934304950E+21 -6.703E+19 1.871E+18 3.873E+33
ETOT 12 -1.59574033415398E+21 4.787E+21 6.360E+17 1.623E+33
ETOT 13 -1.83414714589283E+21 -2.384E+20 5.337E+17 1.013E+33
ETOT 14 -6.94228795675008E+20 1.140E+21 4.301E+17 5.481E+32
ETOT 15 -5.22011349845945E+21 -4.526E+21 8.938E+17 2.088E+33
ETOT 16 -1.39191901316062E+22 -8.699E+21 2.364E+18 6.169E+33
ETOT 17 -1.73959468100041E+23 -1.600E+23 9.458E+18 7.024E+34
ETOT 18 -5.48088496236819E+27 -5.481E+27 2.921E+23 2.187E+39
ETOT 19 -2.06584509246216E+27 3.415E+27 5.992E+23 7.246E+38
ETOT 20 -2.42310870572744E+27 -3.573E+26 3.070E+23 1.429E+39
ETOT 21 -3.03682668117395E+27 -6.137E+26 4.644E+23 1.387E+39
ETOT 22 -1.17836772121987E+28 -8.747E+27 2.079E+24 5.169E+39
ETOT 23 -1.20666930630534E+29 -1.089E+29 1.145E+25 5.058E+40
ETOT 24 -3.94438955558455E+29 -2.738E+29 3.200E+25 1.552E+41
ETOT 25 -6.60697459207140E+30 -6.213E+30 4.935E+26 2.625E+42
ETOT 26 -3.32093663898923E+34 -3.320E+34 1.900E+30 1.326E+46
ETOT 27 -2.70732903422640E+33 3.050E+34 1.467E+30 3.913E+45
ETOT 28 -3.28094153959660E+34 -3.010E+34 2.152E+30 1.390E+46
ETOT 29 -3.58850834877487E+34 -3.076E+33 6.262E+30 1.387E+46
ETOT 30 -3.63379950342011E+35 -3.275E+35 3.283E+31 1.511E+47
ETOT 31 -4.46770054852544E+36 -4.104E+36 3.215E+32 1.828E+48
ETOT 32 -1.41943688135983E+37 -9.727E+36 1.213E+33 5.490E+48
ETOT 33 -3.02813162778023E+38 -2.886E+38 2.011E+34 1.195E+50
922,19 33%

5

How about norm conserving potentials? This will narrow down the possible culprits…

6

Hi mverstra,

I just tried a run with NC+PBE, the calculation took less than a minute to give converged values! However, since I want to keep a pseudopotential consistency in my research(I used PAW+GGA+PBE before for band structure and relax studies), is there a possible way I can still use PAW+GGA+PBE pseudos for a 2D material’s phonon response study? (Ultimately I want to obtain a graph of phonon band structure for my 2D materials)

7

Hi again.

The structure and electronic bands will change only minimally: Use NC for everything if you just want the phonons. PAW is crucial only if you need +U or nuclear quantities.

Thanks for doing this very useful check - we have something to work from to debug.

@torrent do you have any ideas ? I bet this is due to xc potential terms in the vacuum.

@Lucas would you mind trying the paw dfpt dataset with ixc 0? This removes xc effects. The answer will be wrong but perhaps convergence will be quick.

Matthieu

8

Hi,

I was able to reproduce the issue.
But I wasn’t able to find a solution…
… except by changing the PAW dataset.
I tried a PAW dataset from GBRV pseudopotentials and then was able to converge (you have to redo the calculation from scratch).
I don’t know why the JTH dataset is problematic. And even if it is related to the dataset itself; it could be an indirect cause…
Use the xml dataset named c_pbe_v1.2_abinit.paw.xml

Marc

9

Hi everyone,

I just tried to use the pseudo @torrent has provided, and used the “ixc 0” for my dfpt dtsets: 3, 4, 5, and 6. And a Bug message shows up:

Input file:
#Crystalline AlAs computation of the phonon spectrum

ndtset 5
occopt 4
tsmear 0.01
#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
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.0000E+00 0.0000E+00 0.0000E+00 # wtk 6.2500E-02
qpt3 2.5000E-01 0.0000E+00 0.0000E+00 # wtk 3.7500E-01
qpt4 5.0000E-01 0.0000E+00 0.0000E+00 # wtk 1.8750E-01
qpt5 5.0000E-01 2.5000E-01 0.0000E+00 # wtk 3.7500E-01
qpt6 0 0 0

#Set 2: Response function for q=0, print out the ddk
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
getwfk2 1 # Calculate d/dk wave function only
#Set 6: Response function calculation of Q=0 phonons and electric field pert.

getwfk6 0
getddk6 2 # d/dk wave functions from last dataset
kptopt6 2 # Modify default to use time-reversal symmetry
rfelfd6 3 # Electric-field perturbation response only
rfphon6 0

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

getwfk3 1
kptopt3 3
rfphon3 1
#rfatpol3 1 2
#rfdir3 1 1 1

getwfk4 1
kptopt4 3
rfphon4 1
#rfatpol4 1 2
#rfdir4 1 1 1

getwfk5 1
kptopt5 3
rfphon5 1
#rfatpol5 1 2
#rfdir5 1 1 1

toldfe 1.0d-10 # This default is active for sets 3-10
#######################################################################
#Common input variables
#Definition of the planewave basis set
acell 4.6597825725E+00 4.6597825725E+00 1.2999393442E+01
rprim sqrt(0.75) -0.5 0.0 # FCC primitive vectors (to be scaled by acell)
sqrt(0.75) 0.5 0.0
0.0 0.0 1.0

#Definition of the atom types
ntypat 1 # There is only one type of atom
znucl 6 # 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.
pp_dirpath “$ABI_PSPDIR” # This is the path to the directory were
# pseudopotentials for tests are stored
pseudos “c_pbe_v1.2_abinit.paw.xml”

                            # Name and location of the pseudopotential

#Definition of the atoms
natom 2 # There are two atoms
typat 1 1 # They both are of type 1, that is, Silicon.
xred # This keyword indicate that the location of the atoms
# will follow, one triplet of number for each atom
1/3 1/3 0.0 # Triplet giving the REDUCED coordinate of atom 1.
2/3 2/3 0.0 # Triplet giving the REDUCED coordinate of atom 2.

#Definition of the planewave basis set
ecut 20.0 # Maximal kinetic energy cut-off, in Hartree
pawecutdg 40.0
pawovlp -1
# from previous wavefunctions, transferred from the old
# to the new k-points.
#Definition of the SCF procedure
nstep 50 # Maximal number of SCF cycles
diemac 12 # Although this is not mandatory, i
#GS calculation
ngkpt 4 4 1 # should be 8 8 1
nshiftk 1
shiftk 0.0 0.0 0.0 # This choice preserves the hexagonal symmetry of the grid,

ixc3 0
ixc4 0
ixc5 0
ixc6 0
pawxcdev 0

Output file:
== DATASET 3 ==================================================================

  • mpi_nproc: 1, omp_nthreads: -1 (-1 if OMP is not activated)

— !DatasetInfo
iteration_state: {dtset: 3, }
dimensions: {natom: 2, nkpt: 16, mband: 5, nsppol: 1, nspinor: 1, nspden: 1, mpw: 1047, }
cutoff_energies: {ecut: 20.0, pawecutdg: 40.0, }
electrons: {nelect: 8.00000000E+00, charge: 0.00000000E+00, occopt: 4.00000000E+00, tsmear: 1.00000000E-02, }
meta: {optdriver: 1, rfphon: 1, }
…

mkfilename : getwfk/=0, take file _WFK from output of DATASET 1.

Exchange-correlation functional for the present dataset will be:
No xc applied (usually for testing) - ixc=0

Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):
R(1)= 4.0354901 -2.3298913 0.0000000 G(1)= 0.1239007 -0.2146023 0.0000000
R(2)= 4.0354901 2.3298913 0.0000000 G(2)= 0.1239007 0.2146023 0.0000000
R(3)= 0.0000000 0.0000000 12.9993934 G(3)= 0.0000000 0.0000000 0.0769267
Unit cell volume ucvol= 2.4444718E+02 bohr^3
Angles (23,13,12)= 9.00000000E+01 9.00000000E+01 6.00000000E+01 degrees
setup1 : take into account q-point for computing boxcut.

Coarse grid specifications (used for wave-functions):

getcut: wavevector= 0.2500 0.0000 0.0000 ngfft= 20 20 54
ecut(hartree)= 20.000 => boxcut(ratio)= 2.06435

Fine grid specifications (used for densities):

getcut: wavevector= 0.2500 0.0000 0.0000 ngfft= 30 30 80
ecut(hartree)= 40.000 => boxcut(ratio)= 2.16202

— Pseudopotential description ------------------------------------------------

  • pspini: atom type 1 psp file is /home/lucas/abinit-10.0.5/tests//Psps_for_tests/c_pbe_v1.2_abinit.paw.xml
  • pspatm: opening atomic psp file /home/lucas/abinit-10.0.5/tests//Psps_for_tests/c_pbe_v1.2_abinit.paw.xml
  • pspatm : Reading pseudopotential header in XML form from /home/lucas/abinit-10.0.5/tests//Psps_for_tests/c_pbe_v1.2_abinit.paw.xml
    Pseudopotential format is: paw10
    basis_size (lnmax)= 4 (lmn_size= 8), orbitals= 0 0 1 1
    Spheres core radius: rc_sph= 1.25214716
    5 radial meshes are used:
    • mesh 1: r(i)=AA*[exp(BB*(i-1))-1], size= 479 , AA= 0.41313E-03 BB= 0.16949E-01
    • mesh 2: r(i)=AA*[exp(BB*(i-1))-1], size= 475 , AA= 0.41313E-03 BB= 0.16949E-01
    • mesh 3: r(i)=AA*[exp(BB*(i-1))-1], size= 508 , AA= 0.41313E-03 BB= 0.16949E-01
    • mesh 4: r(i)=AA*[exp(BB*(i-1))-1], size= 596 , AA= 0.41313E-03 BB= 0.16949E-01
    • mesh 5: r(i)=AA*[exp(BB*(i-1))-1], size= 618 , AA= 0.41313E-03 BB= 0.16949E-01
      Shapefunction is BESSEL type: shapef(r,l)=aa(1,l)*jl(q(1,l)*r)+aa(2,l)*jl(q(2,l)*r)
      Radius for shape functions = 1.25214716
      mmax= 479
      Radial grid used for partial waves is grid 1
      Radial grid used for projectors is grid 2
      Radial grid used for (t)core density is grid 3
      Radial grid used for Vloc is grid 4
      Radial grid used for pseudo valence density is grid 5
      Compensation charge density is taken into account in XC energy/potential
      pspatm: atomic psp has been read and splines computed

==> initialize data related to q vector <==

The list of irreducible perturbations for this q vector is:
1) idir= 1 ipert= 1
2) idir= 2 ipert= 1
3) idir= 3 ipert= 1

================================================================================

The perturbation idir= 1 ipert= 2 is
symmetric of a previously calculated perturbation.
So, its SCF calculation is not needed.

The perturbation idir= 2 ipert= 2 is
symmetric of a previously calculated perturbation.
So, its SCF calculation is not needed.

The perturbation idir= 3 ipert= 2 is
symmetric of a previously calculated perturbation.
So, its SCF calculation is not needed.

Perturbation wavevector (in red.coord.) 0.250000 0.000000 0.000000
Perturbation : displacement of atom 1 along direction 1
Found 2 symmetries that leave the perturbation invariant.
symkpt : not enough symmetry to change the number of k points.

dfpt_looppert : total number of electrons, from k and k+q
fully or partially occupied states are 8.000000E+00 and 8.000000E+00.
Initialisation of the first-order wave-functions :
ireadwf= 0

— !BeginCycle
iteration_state: {dtset: 3, }
solver: {iscf: 7, nstep: 50, nline: 4, wfoptalg: 10, }
tolerances: {toldfe: 1.00E-10, }
…

 iter   2DEtotal(Ha)        deltaE(Ha) residm    vres2

— !BUG
src_file: m_paw_denpot.F90
src_line: 246
mpi_rank: 0
message: |
XC kernels for ground state must be in memory!
…

Action: contact ABINIT group (please attach the output of abinit -b)

Thank you guys in advance for trying to help me, this really means a lot and sheds a light on my research which has put me in a bottleneck for two months.

10

Hi again Lucas,

so I think your last error arises because you have different ixc for GS and DFPT, and abinit can tell. Given Marc has reproduced the error, the last useful test would be a full run (GS and DFPT) with the JTH (problematic) paw dataset and ixc 0 at all stages. If that does not diverge, then we have a culprit and can debug further.

For your production runs (if you still want PAW), you should follow Marc’s recommendation with the GBRV potential and do everything consistently GS + DFPT etc. You should be fine

best

Matthieu