Printing dynamical matrix with LO-TO splitting

Running Abinit Response calculations
1

Hi everyone,

I am trying to look at the dynamical matrix in the DFPT Gamma point calculation with LO-TO splitting. The output printed in the .abo file is the dynamical matrix with no LO-TO splitting. Is there a way to get the dynamical matrix along the different directions with the non-analytic core corrections?

asr 1
eivec 1
nph1l 1
qph1l 0.0 0.0 0.0 1.0

nph2l 1
qph2l 0.0 0.0 1.0 0.0

This post processing gives the phonon frequency and the eigen vectors corresponding to the LO-TO splitting, but is there a way to get the dynamical matrix from which these were obtained in ABINIT?

Thanks
Dominic

2

Hi Dominic,

  • either you recompose it from the eigenvectors (mind the mass normalizations) and values

  • or you hack the code to print out the corrected dynamical matrix:

77_ddb/m_phonons.F90 around line 1933, print d2cart after gtdyn9. You will get output for every q point calculation, so check for qphnrm(1) being = 0, and you have to distinguish which printing comes from which q→ limit

3

Hi,

Thanks for the reply. I tried changing the code in mkphbs subroutine.

    call gtdyn9(ddb%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart,Crystal%gmet,ddb%gprim,ddb%mpert,natom, &
          Ifc%nrpt,qphnrm(1),qphon,Crystal%rmet,ddb%rprim,Ifc%rpt,Ifc%trans,Crystal%ucvol,Ifc%wghatm,Crystal%xred,ifc%zeff,&
          ifc%qdrp_cart,ifc%ewald_option,xmpi_comm_self,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad)
      
     ! ========== ADD YOUR CODE HERE (START) ==========
     if (abs(qphnrm(1)) < tol8) then
       write(std_out,*) '==== LO-TO CORRECTED DYNAMICAL MATRIX ===='
       write(std_out,*) 'Direction:', qphon(:)
       write(std_out,*) 'qphnrm:', qphnrm(1)
       do ii=1,3*natom
         do jj=1,3*natom
           write(std_out,'(2I5,2ES22.12)') ii, jj, d2cart(1,ii+(jj-1)*3*natom), d2cart(2,ii+(jj-1)*3*natom)
         end do
       end do
       write(std_out,*) '=========================================='
     end if
     ! ========== ADD YOUR CODE HERE (END) ==========

I compiled this and was not printing the dynamical matrix with the anaddb code. Is this the right way to change the code?

Thanks
Dominic

4

Hi,

I tried the first suggestion to reconstruct from eigen values and eigen vectors. Initially, I used anaddb to obtain the phonon frequencies and eigen-displacements along the ⟨110⟩ direction, with the following setup:

asr 1
eivec 1
nph1l 1
qph1l 0.0 0.0 0.0    1.0

nph2l  1
qph2l  1.0 1.0 0.0 0.0

Using the resulting phonon energies and eigenvectors, I applied mass normalization and then reconstructed the dynamical matrix using and then did D = V E V^T to get the dynamical matrix (attaching the code here).

I then enforced symmetry by explicitly symmetrizing the Hessian as H = 0.5*(H + H^T).

However, I am encountering an issue in the reconstructed matrix. Even after symmetrization, certain symmetry-equivalent force-constant elements are not exactly degenerate. For example, in a 10-atom SrTiO_3 cell (atoms 1–2 = Sr, 3–4 = Ti, 5–10 = O), terms such as Sr1x–Ti3x and Sr1x–Ti4x differ slightly:

1 1 1 3 -0.0158635170 0.0000000000 vs 1 1 1 4 -0.0158637861 0.0000000000

Then for terms: Sr1x-Ti3y vs Sr1x-Ti4y,
1 1 2 3 0.0019849360 0.0000000000 vs 1 1 2 4 -0.0019850831 0.0000000000

and for the Sr1x-O7x vs O8x vs O9x vs O10x, we have,
1 1 1 7 -0.0031236763 0.0000000000 vs 1 1 1 8 -0.0031239747 0.0000000000 vs 1 1 1 9 -0.0031236763 0.0000000000 vs 1 1 1 10 -0.0031239747 0.0000000000

This small error in the dynamical matrix gets transferred with the phonon frequency computation as well (loss in the degeneracy of the modes).

Is there any mistake/missing in my reconstruction of the dynamical matrix process?

Thanks
Dominic
dynamical.py.txt (37.9 KB)

5

sounds correct, there may be different sources for symmetry breaking noise:

  • in the DFPT the derivatives are all calculated explicitly, often with non stationary expressions for the off diagonal terms. These will have some noise, and larger than the diagonal elements.
  • in anaddb there are routines to symmetrize the dynamical matrix using the symops of the space group. These should be activated, perhaps your space group was not correctly recognized? Anyhow, degeneracies are usually ensured at this stage, and I have not had problems on this level.
  • In some cases imposing the symmetries, charge neutrality, acoustic sum rule, and hermiticity can be contradictory to some extent. The order you do them in can break the previous symmetrizations using the last one. Abinit does it sequentially (I forgot the exact order, should be in m_phonon.F90). The clean option (to be implemented) is that of the TDEP code, which is to optimize under all constraints at once (eg least means squared), with a strong weight on these exact conditions
  • It could be that the diagonalization algorithm (lapack I think) is not precise enough in the eigenvector outputs. I have had similar cases with only 5-6 digits of precision for degenerate modes. This would propagate as you mention. If you dump the dynmat from the code you should be able to see how many digits you have on the input matrix
6

I fixed the issue by setting a tolerance factor. In the anaddb output, it mentions that values less than 1e-7 would be set to zero, so I set that and now diagonalising the dynamical matrix gives me proper degeneracy of the mode. Thanks for the suggestions.

Dominic