DFPT+U doesnt work for q != Gamma

tyst3273 · June 21, 2023, 3:58pm

I am trying to calculate phonon dispersions with DFPT and DFT+U. It seems like this isn’t allowed… however, this old thread seems to say it works (LDA+U for phonon and electron-phonon calculations - ABINIT Discussion Forums). I tried both LDA and GGA PAW files and, for both, I get the message:

chkdpr: ERROR -
Context : the value of the variable usepawu is 1.
The value of the input variable norm(qpt) is 6.250000000000E-02,
while it must be equal to 0.000000000000E+00.
Action: you should change the input variables norm(qpt) or usepawu.

Am I doing something wrong, or is DFT+U disabled for phonons at q != Gamma.

Thanks!

bguster · June 22, 2023, 2:12pm

Hi tyst3273,

Can you provide a little bit more context to your calculation, such a k/q-grid and the type of perturbations you are performing?
Thanks!
Bogdan

tyst3273 · June 23, 2023, 7:16pm

Thanks for the reply! My replies keep getting flagged as spam… I guess because of the attachments! So I will just paste the input here.

First I run an SCF calculation to get converged density then an NSCF calculation to converge the wavefunctions. Then I do E-field and then phonon perturbations at Gamma (q=0,0,0). All 4 of these calculations run with no errors or warnings! See the input below:

general

ndtset 4
autoparal 1

file I/O

prtwf 0
prtden 1
getwfk 2
getden 1

pp’s

pp_dirpath “/home/ty/pseudo/abi/paw_lda_JTH”
pseudos “Ti.LDA_PW-JTH.xml,O.LDA_PW-JTH.xml”

lattice

ntypat 2
natom 6
typat 1 1 2 2 2 2
znucl 22 8
acell 4.588027812745 4.588027812745 3.003952434907 Angstrom
rprim 1.000000000000 0.000000000000 0.000000000000
0.000000000000 1.000000000000 0.000000000000
0.000000000000 0.000000000000 1.000000000000
xred 0.500000000000 0.500000000000 0.500000000000
0.000000000000 0.000000000000 0.000000000000
0.195467607839 0.804532392161 0.500000000000
0.804532392161 0.195467607839 0.500000000000
0.304532392161 0.304532392161 -0.000000000000
0.695467607839 0.695467607839 0.000000000000

scf control

pawxcdev 0
nstep 200
ecut 1000 eV
pawecutdg 1500 eV
occopt 7
tsmear 0.05 eV
tolvrs 1e-16

mixing

diemac 10.0
diemix 0.7
dielng 1.0

kpoints

ngkpt 4 4 6
nshiftk 1
shiftk 0 0 0
kptopt 3

DFT+U

usepawu 1
lpawu 2 1
upawu 8.00 6.00 eV
jpawu 0.00 0.00 eV

electrons

nband 36
nbdbuf 6

ffts

ngfft 48 48 32
ngfftdg 60 60 40

qpoints – 4x4x6 = 24 qpts

ngqpt 4 4 6
nshiftq 1
shiftq 0 0 0
qptopt 1

respfn

nqpt 1
rfasr 1
rfphon 1
rfatpol 1 6
rfdir 1 1 1

======== ground state ========

gs

kptopt1 1
prtwf1 1
getden1 0
getwfk1 0
rfphon1 0
nqpt1 0

tighten convergence on gs wfk

tolwfr2 1e-20
kptopt2 1
iscf2 -2
prtwf2 1
getwfk2 0
rfphon2 0
nqpt2 0

======== phonons ===========

d/dk at q == 0,0,0

tolwfr3 1e-20
iscf3 -3
prtwf3 1
rfelfd3 2
rfphon3 0
iqpt3 1

pert at q == 0,0,0

getddk4 -1
rfelfd4 3
iqpt4 1

Then I try to run an NSCF calculation at k+q followed by a phonon perturbation at q != 0,0,0. abinit reads the input file and when it reads dataset 2 (the phonon pertubation) it crashes with the error message. Here is the input (it reads GS density and wavefunctions from the Gamma calculation before):

general

ndtset 2
autoparal 1

file I/O

prtwf 0
prtden 1
getwfk_filepath “…/gamma/runo_DS2_WFK”
getden_filepath “…/gamma/runo_DS1_DEN”

pp’s

pp_dirpath “/home/ty/pseudo/abi/paw_lda_JTH”
pseudos “Ti.LDA_PW-JTH.xml,O.LDA_PW-JTH.xml”

lattice

ntypat 2
natom 6
typat 1 1 2 2 2 2
znucl 22 8
acell 4.588027812745 4.588027812745 3.003952434907 Angstrom
rprim 1.000000000000 0.000000000000 0.000000000000
0.000000000000 1.000000000000 0.000000000000
0.000000000000 0.000000000000 1.000000000000
xred 0.500000000000 0.500000000000 0.500000000000
0.000000000000 0.000000000000 0.000000000000
0.195467607839 0.804532392161 0.500000000000
0.804532392161 0.195467607839 0.500000000000
0.304532392161 0.304532392161 -0.000000000000
0.695467607839 0.695467607839 0.000000000000

scf control

pawxcdev 0
nstep 200
ecut 1000 eV
pawecutdg 1500 eV
occopt 7
tsmear 0.05 eV
tolvrs 1e-16

mixing

diemac 10.0
diemix 0.7
dielng 1.0

kpoints

ngkpt 4 4 6
nshiftk 1
shiftk 0 0 0
kptopt 3

DFT+U

usepawu 1
lpawu 2 1
upawu 8.00 6.00 eV
jpawu 0.00 0.00 eV

electrons

nband 36
nbdbuf 6

ffts

ngfft 48 48 32
ngfftdg 60 60 40

qpoints – 4x4x6 = 24 qpts

ngqpt 4 4 6
nshiftq 1
shiftq 0 0 0
qptopt 1

respfn

nqpt 1
rfphon 1
rfatpol 1 6
rfdir 1 1 1

======== phonons ===========

nscf at q != 0,0,0

tolwfr1 1e-20
iscf1 -2
prtwf1 1
rfphon1 0
iqpt1 2

pert at q != 0,0,0

getwfq2 -1
iqpt2 2

The *.abo reports the message:

…
chkinp: Checking input parameters for consistency, jdtset= 1.

chkinp: Checking input parameters for consistency, jdtset= 2.

chkdpr: ERROR -
Context : the value of the variable usepawu is 1.
The value of the input variable norm(qpt) is 6.250000000000E-02,
while it must be equal to 0.000000000000E+00.
Action: you should change the input variables norm(qpt) or usepawu.

Checking consistency of input data against itself revealed some problem(s).
So, stopping. The details of the problem(s) are given in the error file or the standard output file (= “log” file).
In parallel, the details might not even be printed there. Then, try running in sequential to see the details.

ebousquet · June 28, 2023, 10:46am

Dear tyst3273,
Can you also tell us which version of Abinit are you using?
Best,
Eric

tyst3273 · June 28, 2023, 3:46pm

Hopefully this helps!

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

=== Build Information ===
Version : 9.8.4
Build target : x86_64_linux_gnu9.4
Build date : 20230620

=== Compiler Suite ===
C compiler : gnu9.4
C++ compiler : gnu9.4
Fortran compiler : gnu9.4
CFLAGS : -g -O2
CXXFLAGS : -g -O3 -mtune=native -march=native
FCFLAGS : -g -ffree-line-length-none -I/usr/include/mkl
FC_LDFLAGS :

=== Optimizations ===
Debug level : basic
Optimization level : aggressive
Architecture : intel_xeon

=== Multicore ===
Parallel build : yes
Parallel I/O :
openMP support :
GPU support :

=== Connectors / Fallbacks ===
LINALG flavor : mkl
FFT flavor : dfti
HDF5 : yes
NetCDF : yes
NetCDF Fortran : yes
LibXC : yes
Wannier90 : yes

=== Experimental features ===
Exports :
GW double-precision :

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Default optimizations:
-O3 -mtune=native -march=native -funroll-loops -faggressive-function-elimination

Optimizations for 43_ptgroups:
-O0

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
CPP options activated during the build:

                CC_GNU                   CXX_GNU                    FC_GNU

             HAVE_DFTI HAVE_FC_ALLOCATABLE_DT...             HAVE_FC_ASYNC

     HAVE_FC_BACKTRACE  HAVE_FC_COMMAND_ARGUMENT      HAVE_FC_COMMAND_LINE

    HAVE_FC_CONTIGUOUS           HAVE_FC_CPUTIME              HAVE_FC_EXIT

         HAVE_FC_FLUSH             HAVE_FC_GAMMA            HAVE_FC_GETENV

HAVE_FC_IEEE_ARITHMETIC HAVE_FC_IEEE_EXCEPTIONS HAVE_FC_INT_QUAD

         HAVE_FC_IOMSG     HAVE_FC_ISO_C_BINDING  HAVE_FC_ISO_FORTRAN_2008

    HAVE_FC_LONG_LINES        HAVE_FC_MOVE_ALLOC  HAVE_FC_ON_THE_FLY_SHAPE

       HAVE_FC_PRIVATE         HAVE_FC_PROTECTED           HAVE_FC_SHIFTLR

     HAVE_FC_STREAM_IO            HAVE_FC_SYSTEM          HAVE_FORTRAN2003

             HAVE_HDF5        HAVE_LIBPAW_ABINIT      HAVE_LIBTETRA_ABINIT

            HAVE_LIBXC         HAVE_LINALG_AXPBY        HAVE_LINALG_GEMM3M

     HAVE_LINALG_GEMMT  HAVE_LINALG_MKL_IMATCOPY   HAVE_LINALG_MKL_OMATADD

HAVE_LINALG_MKL_OMATCOPY HAVE_LINALG_MKL_THREADS HAVE_LINALG_SCALAPACK

              HAVE_MPI                 HAVE_MPI2       HAVE_MPI_IALLGATHER

   HAVE_MPI_IALLREDUCE        HAVE_MPI_IALLTOALL       HAVE_MPI_IALLTOALLV

       HAVE_MPI_IBCAST         HAVE_MPI_IGATHERV        HAVE_MPI_INTEGER16

           HAVE_MPI_IO HAVE_MPI_TYPE_CREATE_S...               HAVE_NETCDF

   HAVE_NETCDF_FORTRAN                HAVE_NUMPY             HAVE_OS_LINUX

     HAVE_TIMER_ABINIT            HAVE_WANNIER90

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

bguster · July 11, 2023, 11:37am

Indeed your post have been flagged as spam.
Nonetheless, your error might be coming your WFK on the k-point grid, as in your qpt is NOT the difference of 2 k-points on the grid. this is to say that you need to use a q-grid commensurate with the k-point grid from your GS or NSCF calculation or you need to interpolate your WFK file in order to include specifically that qpt.

From what I can tell you are using a 4x4x6 kpt and qpt grid while you are requesting a qpt = 1/16 which is not a point on your grid, so you need to enhance your qpt grid with one of the two solutions. If 16 kpoints are way too much for your calculations, look into interpolating your WFK.
Let me know if this works for you!
Bogdan

tyst3273 · August 29, 2023, 4:51am

Hi, sorry I have been very busy with other things for a while. I think I don’t understand then. Please excuse my ignorance if I am wrong in what follows!

(i) Doesn’t the nscf calculation at k+q (dataset 1) produce the wavefunctions at k+q? The wavefunctions at k (q=0) are read from a previous run. I came up with this work flow following the abinit help file (and trial and error!)

From the help file:

When one considers the response to an atomic displacement for a general q point, the following procedure is suggested:
first, a self-consistent ground-state computation with the restricted set of k-points in the Irreducible Brillouin Zone (with kptopt=1)*
second, a non-self-consistent ground-state run with the set of k+q points, that might be reduced thanks to symmetries (with kptopt=1)*
third, a self-consistent response-function computation of the atomic displacement perturbation, with the full set of k-points (with kptopt=3)*

I think this is what I am doing.

(ii) I choose a qpt that is on the grid… the qpts are generated automatically according to ngqpt and I select each for a respfn calculation using the iqpt variable. Please let me know if this is wrong!

(iii) The calculation runs perfectly if I turn of DFT+U, e.g. usepaw=0. Of course I expect the results to be very different. Similarly, using the exact same “template” for DFPT phonon calculations, I can do both PAW and NCPP calculations without DFT+U and I get results that are exactly what I expect. i.e. I don’t think I am doing things incorrectly …

Sorry again for my confusion. I would just really like electron-phonon matrix elements with DFT+U so am still trying to get this to work.