Discrepancy in force constants spatial decay: Flexoelectric vs. Phonon DDBs

Running Abinit Response calculations
1

Hi everyone,

I am analyzing the spatial decay of longitudinal interatomic force constants (IFCs) using the abipy package, via the ifc.plot_longitudinal_ifc() function. I have observed a severe truncation in the real-space IFCs when extracting them from a flexoelectric DFPT calculation compared to a GM point phonon calculation.

Flexoelectric Run: I merged the DDB files from my flexoelectric calculations using mrgddb. When plotting the longitudinal IFCs via abipy, the output yields only a single data point at roughly 5 Bohr, abruptly truncating the long-range interactions.

640×480 11.4 KB

GM Phonon Run: I generated a DDB from a phonon run on the same system. Running the exact same abipy plotting function yields the expected spatial decay curve, extending out to ~20 Bohr with 5-6 distinct neighbor shells.

640×480 12.8 KB

  • Does the flexoelectric DFPT implementation in ABINIT inherently ignore or truncate the long-range force constants during the calculation?

The phonons were obtained with script:

ndtset  3

# Ground state calculation
  kptopt1   1             # Automatic generation of k points, taking symmetry into account
  tolvrs1   1.0d-18       # SCF stopping criterion

# Response Function calculation: d/dk
  rfelfd2   2             # Activate the calculation of the d/dk perturbation
  nqpt2     1
  qpt2      0.0 0.0 0.0   # This is a calculation at the Gamma point
  getwfk2   -1            # Use the output wf of the previous dataset as input
  kptopt2   2             # Automatic generation of k points, using time-reversal symmetry
  iscf2     -3            # Non-self-consistent d/dk perturbation
  tolwfr2   1.0d-22       # Very low tolwfr for non-self-consistent calculations

# Response Function calculation: electric field perturbation and phonons
  rfphon3   1             # Activate the calculation of the atomic displacement perturbations
  rfelfd3   3             # Activate the calculation of the electric field perturbation
  nqpt3     1
  qpt3      0.0 0.0 0.0   # This is a calculation at the Gamma point
  getwfk3   -2            # Use output wfs from dataset 1
  getddk3   -1            # Use output ddk wfs from dataset 2
  kptopt3   2             # Automatic generation of k points, using time-reversal symmetry
  tolvrs3   1.0d-18      # Very low tolwfr for non-self-consistent calculations

acell   7.3609984220E+00  7.3609984220E+00  7.3609984220E+00

rprim   0.0000000000E+00  7.0710678119E-01  7.0710678119E-01
                       7.0710678119E-01  0.0000000000E+00  7.0710678119E-01
                       7.0710678119E-01  7.0710678119E-01  0.0000000000E+00
       
natom    2  
ntypat   2  
typat    1 2  
znucl    56 8  
nband    11
nbdbuf   3
xred     0.000000000000000   0.000000000000000   0.000000000000000 
         0.500000000000000   0.500000000000000   0.500000000000000 

pp_dirpath "~/pseudo/PBEsol_NCPP/"      # Path to the pseudopotential directory
pseudos "Ba.oncvpsp.psp8,O.oncvpsp.psp8"

ecut 150         # Maximal kinetic energy cut-off, in Hartree
ecutsm 0.5

# Definition of the k-point grid
ngkpt  10 10 10
nshiftk   1
shiftk  0.0 0.0 0.0 

nstep 500       
diemac 4.0

while the flexoelectric calculations were done using the script:

ndtset 5

#Set 1 : ground state self-consistency
#*************************************
getwfk1   0          
kptopt1   1         
tolvrs1   1.0d-18 

#Set 2: Response function calculation of d/dk wave function
#**********************************************************
iscf2    -3       
rfelfd2   2      
tolwfr2   1.0d-20 

#Set 3: Response function calculation of d2/dkdk wave function
#*************************************************************
getddk3   2 
iscf3    -3 
rf2_dkdk3 3 
tolwfr3   1.0d-20  

#Set 4 : Response function calculation to q=0 phonons, electric field and strain
#*******************************************************************************
getddk4   2          
rfelfd4   3         
rfphon4   1        
rfstrs4   3      
rfstrs_ref4   1      
tolvrs4   1.0d-16
prepalw4  1 

#Set 5 : Long-wave magnitudes calculation
#****************************************
optdriver5 10        
get1wf5   4         
get1den5  4        
getddk5   2       
getdkdk5  3      
lw_flexo5 1 

getwfk  1  
useylm  1
kptopt 2

acell   7.3609984220E+00  7.3609984220E+00  7.3609984220E+00

rprim   0.0000000000E+00  7.0710678119E-01  7.0710678119E-01
                       7.0710678119E-01  0.0000000000E+00  7.0710678119E-01
                       7.0710678119E-01  7.0710678119E-01  0.0000000000E+00
       
natom    2  
ntypat   2  
typat    1 2  
znucl    56 8  
nband    11
nbdbuf   3
xred     0.000000000000000   0.000000000000000   0.000000000000000 
         0.500000000000000   0.500000000000000   0.500000000000000 

pp_dirpath "~/pseudo/PBEsol_NCPP/"      # Path to the pseudopotential directory
pseudos "Ba.oncvpsp.psp8,O.oncvpsp.psp8"

ecut 150         # Maximal kinetic energy cut-off, in Hartree
ecutsm 0.5

ngkpt  10 10 10
nshiftk   1
shiftk  0.0 0.0 0.0 

nstep 500      
diemac 4.0

The k-grids are the same in both cases, then why are the force constants so different with the flexoelectric calculations?
Are the long-range interactions not considered in the flexoelectric tensor calculation?

Thanks
Dominic

2

Hi Dominic,

I’m not an expert of the calculations themselves, but what the LW formalism does is a perturbation wrt q, so it gives you a kind of Taylor expansion in wavelengths. The question is what is abipy plotting? I suspect it just gives you two points to show the affine function which is effectively calculated by the LW. The point at finite r could be anywhere, as abinit is giving you the slope.

In the phonon case, you calculate the dynamical matrix, and even for q=0 in a multi atom system you have some distance dependence (atom 1-2 1-3…) within the primitive cell. If you do finite q then you get full distance dependent ifc(r) with a Fourier transform.

We might get confirmation from @gmatteo who coded most of abipy, or @mroyo who worked on the LW. or you can have a look at the abipy sources.

3

Dear Dominic,
I don’t know about abipy, but in anaddb the default value of the variables dipquad and quadquad is 1. This means that, if anaddb finds the quadrupoles in the DDB file -which are calculated in the flexoelectricity run, it not only removes the dipole-dipole LR contributions from the IFCs but also the dipole-quadrupole and quadrupole-quadrupole ones. This usually makes the IFCs to decay faster.

However, I’m not sure if this is the cause of your problem because you are using IFCs calculated at the Gamma point only (right?) whereas the DQ and QQ interactions appear at first and second order in q, so they should vanish at Gamma.

4

Thanks @mroyo and @mverstra for the reply.

I was trying to understand how the individual clamped-ion flexoelectric force response tensor components is estimated. Since these terms would be the second order moment of the real space force constants, I was looking at the IFC decay.

So does this mean the clamped-ion flexoelectric force response tensor, does not include the LR interactions or is it not included in the IFC plot from abipy (because of the default settings of dipquad=1 and quadquad=1). I was interested in seeing the long range dependence of the CI-flexoelectric force response tensor.

5

Hi Dominic,
the CI flexoelectric force-response tensor is defined, and calculated in abinit, under short-circuit electrical boundary conditions, hence it does not include LR interactions at all.

A direction dependent nonanalytical force-response tensor can be defined by incorporating the LR contributions, see Eq. (S19) of this reference: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.125.217602 . However, such a calculation is not implemented currently in abinit or anaddb, although the necessary ingredients are there.

6

Thanks for the explanation @mroyo. So the CI flexo force response tensor does not depend on the Non-analytical corrections.

By long range earlier, I was referring to the distance cutoff to which the force constants are estimated accurately. As the force response tensor can be described a second order moment in the real space, the distance to which we have the force constants estimated has a crucial role with flexoelectric coefficient calculation (I believe).

Is there a way to know the distance cutoff/extend of neighbors considered in the Flexoelectric calculation for the force response tensor?

7

Hi Dominic,

Is there a way to know the distance cutoff/extend of neighbors considered in the Flexoelectric >calculation for the force response tensor?

Since it is calculated under short-circuit conditions, this implies a fast decay, of the order of the inter-atomic distance, for the IFCs.

To exactly estimate the decay one could try to calculate the force constants over a regular grid of q points, while imposing short-circuit conditions. The Fourier transform of the IFCs into real space should serve to estimate the decay, if any, beyond the original primitive cell.

The possibility to run finte-q response functions calculations under short-circuit boundary conditions is not yet available in the public version of ABINIT. Most surely it will appear in the next release.

8

Hi Miquel,

Just to confirm - by the current method in which the flexoelectric force response tensor is computed, it does not account for the non-analytic long-range electrostatic corrections (short-circuit condition).

The Tensor Calculation (Reciprocal Space): Because the tensor is calculated via analytical q-derivatives at the Gamma point using the LW method, it implicitly captures the full, mathematically converged sum of these short-range interactions across all neighbor shells (i.e., the full \sum \Phi R^2). It does not stop at the nearest neighbor.

The Real-Space Output Artifact: However, if we try to explicitly extract the real-space force constants (\Phi) from this Gamma-only DDB file, the mathematical limitation of inverse-Fourier-transforming a single q-point artificially truncates the data at the primitive cell boundary. This perfectly explains why my real-space output files only show the inside-cell Ba-O interactions and immediately zero out everything outside the unit cell.

And the final CI-flexo force response tensor value we see from the anaddb post-processing of the flexoelectric calculations have the correct values as all IFCs in the reciprocal space are considered in estimating the response.

Thanks
Dominic

9

That’s how I see it, yes.

Miquel