Density functional perturbation theory for lattice dynamics with fully relativistic ultrasoft pseudopotentials: the magnetic case
Andrea Urru Andrea Dal Corso





International School for Advanced Studies (SISSA), 

Via Bonomea 265, 34136 Trieste (Italy).

International School for Advanced Studies (SISSA), 

Via Bonomea 265, 34136 Trieste (Italy).
DEMOCRITOS IOM-CNR Trieste (Italy).




We extend density functional perturbation theory for lattice dynamics with fully 
relativistic ultrasoft pseudopotentials to magnetic materials.
Our approach is based on the application of the time-reversal operator to the Sternheimer
linear system and to its self-consistent solutions.
Moreover, we discuss how to include in the formalism the symmetry operations of the magnetic 
point group which require the time-reversal operator. We validate our implementation by comparison 
with the frozen phonon method in fcc Ni and in a monatomic ferromagnetic Pt wire.



 Density functional perturbation theory for lattice dynamics with fully relativistic ultrasoft pseudopotentials: the                                                                                    
magnetic case
    Andrea Dal Corso
    December 30, 2023
======================================================================================================================================================================================================================




§ INTRODUCTION

Density Functional Perturbation Theory (DFPT) is widely used for the computation of the 
linear response properties of solids, and in particular for the study of their lattice dynamics. <cit.>
Some years ago, one of us applied DFPT <cit.> to a scheme based on plane waves and norm conserving (NC) or 
ultrasoft (US) <cit.> pseudopotentials (PPs), that allow the introduction of spin-orbit effects within a 
fully relativistic (FR) density functional formalism and can be written in a form very similar to the scalar 
relativistic (SR) one. 
However, the theory presented in Ref. us_fr_dfpt was implemented only for time-reversal invariant systems, 
and therefore applications that include spin-orbit so far have been limited to non-magnetic solids. <cit.> 

In this work we extend this theory to the case of magnetic systems, by explicitly considering the presence 
of an exchange-correlation magnetic field in the Hamiltonian. DFPT equations in presence of a magnetic 
field have been recently written to calculate magnons with NC PPs in Refs. Giustino 
and magnon. In Ref. Giustino, the charge density induced by 
a periodic perturbation was computed by using the response to a perturbation at wave vector q and the response to 
a perturbation at -q, 
while in Ref. magnon the problem at -q was not solved, but the time-reversal operator was used 
to obtain a second Sternheimer equation with a reversed magnetic field. The two formulations are equivalent. 
We generalize the theory of Ref. magnon to a phonon perturbation, avoiding the study of the 
response at -q, and write it in a form applicable to both NC and US PPs. 

In presence of a magnetic field, the solid is invariant upon the symmetry operations of the magnetic space 
group. Some of these operations require the time-reversal operator. We discuss how to exploit these symmetries 
for the symmetrization of the induced charge and magnetization densities and for the dynamical matrix. 

Finally, we validate our method in ferromagnetic fcc Ni first computing the phonon frequencies at the X point in 
the Brillouin zone (BZ) and comparing with DFPT results, and then by computing the phonon dispersions.  
Moreover, we apply our method to a monatomic ferromagnetic Pt nanowire and compare its 
vibrational properties when the magnetization is parallel or perpendicular to the wire. Also for this case 
we compare the DFPT results to the frozen phonon method for a phonon wavevector q = π/a and q = π/2a, and then 
we compute by DFPT the phonon dispersion in the 1D BZ.



§ DFPT WITH FULLY RELATIVISTIC US-PPS

In the Density Functional Theory (DFT) with FR US PPs, which accounts for spin-orbit effects, 
the minimization of the total energy functional leads to the Kohn-Sham <cit.> equations for the two-component spinor wave functions
<cit.>: 


    ∑_σ ' H^σσ'|Ψ_i^σ'⟩ = ϵ_i ∑_σ' S^σσ'|Ψ_i^σ'⟩,


where S^σσ' is the overlap matrix needed in the US scheme, and the Hamiltonian H^σσ' is:


    H^σσ' = - 1/2∇^2 δ_σσ' + V_KS^σσ'.


V_KS^σσ' is the total Kohn-Sham potential:


    V_KS^σσ' = V_NL^σσ' + ∑_σ_1 σ_2∫ d^3r V_LOC^σ_1 σ_2 K_σσ'^σ_1 σ_2,


where: 


    K_σσ'^σ_1 σ_2 (r, r_1, r_2)     = δ(r - r_1)  δ(r - r_2)  δ_σ_1 σ δ_σ_2 σ'
        
    + ∑_I m n∑_m_1 n_1 Q_m n^I f_m_1 m^σσ_1β_m^I (r_1) 
       × f_n n_1^σ_2 σ'β_n_1^* I (r_2),
 

where I = {ρ, s”}, while f_m_1 m^σ_1 σ, f_n n_1^σ' σ_2, as well as the indeces m, 
n, m_1, and n_1 are defined in Eq. (5) of Ref. us_fr_dfpt.
In particular, in Eq. (<ref>) V_LOC^σσ' = V_effδ_σσ' - μ_B B_xc·σ^σσ', and V_eff = V_loc + V_H + V_xc, is the sum of local, Hartree, and exchange and 
correlation potential, and V_NL is the bare non local potential: in magnetic systems the spin, represented 
by the Pauli matrices, is coupled to the exchange-correlation magnetic field, 
B_xc, defined as B_xc = -δ E_xc / δm. This term breaks the time-reversal 
symmetry: indeed, introducing the time-reversal operator, T = σ_y K, where 
K is the complex-conjugation operator and σ_y is the Pauli matrix, the following relationship holds:


    T H^[B_xc]T^† = H^[-B_xc].



We first exploit the time-reversal operator to rewrite the induced spin density. Following 
the notation of Ref. us_fr_dfpt, we consider both the metallic and the insulating case. The change of the 
spin density induced by the variation of an external parameter μ (Eq. (10) of Ref. us_fr_dfpt) may 
be written as:


    dn^σσ'/d μ =    ∑_i∑_σ_1 σ_2[ ⟨Ψ_i^σ_1| 
    K_σ_1 σ_2^σσ'|Δ^μΨ_i^σ_2⟩
        + ∑_σ_3 σ_4⟨(TΨ_i)^σ_1|T_σ_1 σ_3  K_σ_3 σ_4^σ' σ T_σ_4 σ_2^†
       ×|( TΔ^μΨ_i )^σ_2⟩] + Δ^μ n^σσ'.


Δ^μ n^σσ' is defined as in the non-magnetic case and corresponds to 
the last two terms of Eq. (10) of Ref. us_fr_dfpt.

The same idea can be applied to the second-order derivatives of the total energy. Only the term 
d^2 E_tot^(2)/d μ dλ (Eq. (19) of Ref. us_fr_dfpt) needs to be rewritten 
by using T:



    d^2 E_tot^(2)/d μ dλ = ∑_i∑_σ_1 σ_2⟨Ψ_i^σ_1|∂ V_KS^[B] 
    σ_1 σ_2/∂λ - ϵ_i ∂ S^σ_1 σ_2/∂λ|Δ^μΨ_i^σ_2⟩
        + ∑_i∑_σ_1 σ_2⟨(TΨ_i)^σ_1|∂ V_KS^[-B] σ_1 σ_2/∂λ - ϵ_i ∂ S^σ_1 σ_2/∂λ|( TΔ^μΨ_i )^σ_2⟩,



while the other contributions can be kept in their original form. Both Eq. (<ref>) and Eq. (<ref>) contain 
two unknown terms, namely |Δ^μΨ_i^σ⟩ and |(TΔ^μΨ_i)^σ⟩. 
The first can be computed by means of the Sternheimer linear system (Eqs. 13 and 14 of Ref. us_fr_dfpt), 
while the second is obtained by solving the following linear system:


    ∑_σ_2[ H^[-𝐁] σ_1 σ_2 - ϵ_i S^σ_1 σ_2] |( TΔ^μΨ_i )^σ_2⟩ = 
        - ∑_σ_2 σ_3Π_c,i^†σ_1 σ_2[ 
    dV_KS^[-𝐁] σ_2 σ_3/d μ - ϵ_i ∂ S^σ_2 σ_3/∂μ]
    |(TΨ_i)^σ_3⟩,


obtained by applying T to both sides of the Sternheimer linear system (Eq. (13) of Ref. us_fr_dfpt) and 
using the fact that T  dV_KS^[B]/d μ T^† = dV_KS^[-B]/d μ. In particular, here we 
introduced the time-reversed projector on the conduction manifold, namely 
Π_c,i^†σ_1 σ_2 = ∑_σσ'T_σ_1 σ P_c,i^†σσ'T^†_σ' σ_2, 
similarly to what proposed in Ref. magnon for the calculation of magnons.
Eqs. (<ref>), (<ref>), and (<ref>) are valid for the US PPs scheme, the 
NC formulation can be obtained as a particular case by writing 
K_σ_1 σ_2^σσ'(r,r_1,r_2) = δ(r - r_1)  δ(r - r_2)  δ_σσ_1 δ_σ' σ_2 
and S^σ_1 σ_2 = δ_σ_1 σ_2. Moreover, the insulating case can be dealt with by putting 
θ_F,i = 1 if the state is occupied or 0 if the state is empty (see Ref. us_dfpt for the 
definition of θ_F,i).



§ PHONONS IN PERIODIC SOLIDS

In this section, we consider a phonon perturbation with a wavevector q perturbing a periodic solid, for which 
the wave functions Ψ_i^σ may be written in the Bloch form, Ψ_k v^σ= e^k·r u_k v^σ, where u_k v^σ is lattice periodic. 
Following the discussion reported in Appendix A of Ref. us_dfpt, we introduce the variation of the density 
and of the wave functions, induced by a phonon perturbation u_ν s' β = 1/√(M_s')( u_s' β 
(q) e^q·R_ν), where  indicates the real part, and define:


    dn^σσ'/d u_s' β (q)    = ∑_ν e^q·R_νdn^σσ'/d u_ν s' β, 
    |Δ^u_s' β (q)Ψ_kv^σ⟩    = ∑_ν e^q·R_ν|Δ^u_ν s' βΨ_kv^σ⟩.


Eq. (<ref>) then becomes:  


    dn^σσ'/d u_s' β (q)    = ∑_k v∑_σ_1 σ_2[ ⟨Ψ_k v^σ_1| K_σ_1 σ_2^σσ'|Δ^u_s' β (q)Ψ_k v^σ_2⟩
        + ∑_σ_3 σ_4⟨(TΨ_- k v)^σ_1| T_σ_1 σ_3 K_σ_3 σ_4^σ' σ T_σ_4 σ_2^†
       ×|( TΔ^u_s' β (- q)Ψ_-k v)^σ_2⟩] + Δ^u_s' β (q) n^σσ'.


In Eq. (<ref>) the second term is identical to the first and is not explicitly computed in time-reversal invariant systems. The same holds for 
(<ref>) below (for the dynamical matrix). Instead, for magnetic systems the two terms must be computed separately. In particular, the 
time-reversed response of the wave functions can be computed by solving the linear system Eq. (<ref>), which becomes: 


    ∑_σ_2[ H^[-B] σ_1 σ_2 - ϵ_-k v S^σ_1 σ_2] |(TΔ^u_s' β (- q)Ψ_- k v)^σ_2⟩ = 
        - ∑_σ_2Π_c, -kv^†σ_1 σ_2[ |ϕ_T -kv^u_s' β (q) [-B]σ_2⟩
        + ∑_σ_3∑_σ_4 σ_5∫ d^3 r d V_H,xc^[-B] σ_4 σ_5/d u_s' β (q) K^σ_4 σ_5_σ_2 σ_3 (r) |(TΨ_- k v)^σ_3⟩],


where, similarly to Ref. us_dfpt we defined:


    |ϕ_T -k v^u_s' β (q) [-B] σ_2⟩ = ∑_σ_3( ∂ V_KS^[-B] σ_2 σ_3/∂ u_s' β (q) - ϵ_-k v∂ S^σ_2 σ_3/∂ u_s' β (q)) |(TΨ_-k v)^σ_3⟩,


in which ∂ V_KS^[-B] σ_2 σ_3 / ∂ u_s' β (q) and ∂ S^σ_2 σ_3 / ∂ u_s' β (q)
are defined similarly to Eq. (<ref>). The action of the time-reversal operator on the linear system changes the sign of the exchange 
and correlation magnetic field, which enters in the Hamiltonian, in d V_H, xc^[B] σ_4 σ_5 / d u_s' β (q), 
and in |ϕ^u_s' β (q) [B] σ_2_k v⟩ through the third term in Eq. (9) of Ref. us_fr_dfpt.
We can then write the contribution to the dynamical matrix coming from d^2 E_tot^(2)/d u_μ s α d u_ν s' β 
in the following way:


    Φ^(2)_s α
     s' β (q)     = 1/N∑_k v∑_σ[ ⟨ϕ^u_s α (q) [B] σ_k v|Δ^u_s' β (q)Ψ_k v^σ⟩
        + ⟨ϕ^u_s α (q) [-B] σ_T -k v|( TΔ^u_s' β (-q)Ψ_-k v) ^σ⟩],


while the other contributions may be kept in their original form, discussed in Ref. us_dfpt. Here, N is the 
number of cells in the solid.

Eq. (<ref>) may be further manipulated by writing explicitly K_σ_1 σ_2^σσ' 
(Eq. (<ref>)). Introducing the periodic parts of the Bloch functions and of the responses of the wave functions, 
we obtain the periodic part of the induced spin density (indicated with a tilde (∼)):


    dn^σσ'/d u_s' β (q)    = ∑_k v[ u_k v^*  σΔ^u_s' β (q) u_k v^σ' + ∑_σ_1 σ_2 U_σ' σ_1
       × (T u_- k v)^*  σ_1 (TΔ^u_s' β (- q) u_- k v)^σ_2 U_σ_2 σ^†] 
        + ∑_s” m   n( Q_m n^s”qΔ^u_s' β (q)ρ_m n^s”σσ') 
        + Δ^u_s' β (q) n^σσ',


where we defined the quantities Q_m n^s”q and Δ^u_s' β (q)ρ_m n^s”σσ' as: 


    Q_m n^s”q = e^-q·r∑_ρ e^q·R_ρ Q_m n^I ,




    Δ^u_s' β (q)ρ_m n^s”σσ'    = ∑_m_1 n_1∑_σ_1 σ_2∑_k v( β_k v^*   s” m_1 σ_1  f_m_1 m^σ_1 σ f_n n_1^σ' σ_2Δ^u_s' β (q)β_k v^s” n_1 σ_2
        + ∑_σ”σ”'β_T   -k v^*   s” m_1 σ_1 f_m_1 m^σ_1 σ” U_σ”σ' U_σσ”'^† f_n n_1^σ”' σ_2Δ^u_s' β (- q)β_T   -k v^s” n_1 σ_2),



where β_T   -k v^s” m_1 σ_1 = e^- k·R_ρ ⟨β^I_m_1| 
(TΨ_- k v)^σ_1⟩, Δ^u_s' β (- q)β_T   -k v^s” n_1 σ_2 = 
e^-  (k+q) ·R_ρ⟨β^I_n_1| (TΔ^u_s' β (- q)Ψ_- k v)^σ_2⟩, 
and we used the fact that: 


    ∑_σ_3 σ_4T_σ_1 σ_3 K_σ_3 σ_4^σ' σT_σ_4 σ_2^† = ∑_σ”σ”' U_σ' σ” K_σ_1 σ_2^σ”σ”' U_σ”' σ^†,


where U = σ_y is the unitary part of the time-reversal operator.
The induced charge and magnetization densities can be computed from the induced spin density (Eq. (<ref>)) as: 


    d n/d u_s α(q)    = ∑_σdn^σσ/d u_s α (q), 
    d m_β/d u_s α(q)    = μ_B ∑_σ σ'dn^σσ'/d u_s α (q)σ_β^σσ'.


In particular, for the induced charge density we use the fact that ∑_σ U_σ_2 σ^† U_σσ_1 = δ_σ_1 σ_2 
and for the induced magnetization density the fact that 
∑_σσ' U_σ_2 σ^†σ_α^σσ' U_σ' σ_1 = - σ_α^σ_1 σ_2, 
so that the terms of the induced spin density that contain the time-reversed wave functions are subtracted in Eq. (<ref>).

The linear system (<ref>) may be written in terms of lattice periodic functions: 

    ∑_σ_2( H_k + q^[-B] σ_1 σ_2 - ϵ_-k v S_k + q^σ_1 σ_2) |(TΔ^u_s' β (- q) u_- k v)^σ_2⟩
        = - ∑_σ_2Π_c, -k v^†σ_1 σ_2, -k - q[ |ϕ_T - k v^u_s' β (q) [-B] σ_2⟩ + ∑_σ_3d V_H, xc^[- B] σ_2 σ_3/d u_s' β (q)
       ×|(T u_k v)^σ_3⟩ + ∑_σ_3∑_s” m_1 n_1^3 I_s” m_1 n_1^u_s' β (q) [-B] σ_2 σ_3|β_m_1 k + q^s”⟩
       ×β_T   -k v^s” n_1 σ_3],


where we defined: 


    ⟨r|β_m_1 k + q^s”⟩ = e^- (k + q) ·r∑_ρ e^ (k + q) ·R_ρβ_m_1^s” (r - R_I)




    ^3 I_s” m_1 n_1^u_s' β (q) [-B] σ_2 σ_3 = ∑_m n∑_σ_4 σ_5 f_m_1 m^σ_2 σ_4 f_n n_1^σ_5 σ_3∑_α A_α^σ_4 σ_5^3 I_s” m n α^u_s' β (q) [-B],


and:

    ^3 I_s” m n α^u_s' β (q) [-B] = ∫ Q_m n^γ(s”) (r - d_s”)  d C_α^[-B]/d u_s' β (q)  d^3 r,


where α = 1, …, 4, A = (1, σ_x, σ_y, σ_z), similarly to Ref. us_fr_pseudo 
(1 is the 2 × 2 identity matrix), and C^[B] = (V_H,xc, -μ_B B_xc, x, -μ_B B_xc, y, -μ_B B_xc, z). 



§ SYMMETRIZATION


We indicate with {S|f} the symmetry operations of the space group of the crystal, where S 
is a rotation (proper or improper) and f is a translation. In a magnetic crystal, we have to 
consider also the operations S such that {TS|f} is a symmetry of the crystal.

Since, for a phonon perturbation, the charge (and magnetization) density response and the dynamical matrix 
are computed at a given finite wave vector q, we use as symmetry operations only those N_S operations 
of the antiunitary small space group of q, the subgroup of the antiunitary space group of the crystal, which 
contains the symmetry operations {S|f} such that:


    Sq = q + G_S,


if {S|f} is a symmetry of the crystal, or:


    Sq = - q + G_S,


if {TS|f} is a symmetry of the crystal. Here, G_S is a reciprocal lattice vector 
that might appear when q is at zone border. In order to distinguish the two cases we introduce a variable 
τ (S) which may take the values τ = 0 or τ = 1 if Eq. (<ref>) or Eq. (<ref>) holds, respectively.
We compute the unsymmetrized induced spin density by summing over the Irreducible Brillouin Zone (IBZ) in Eqs. 
(<ref>) and (<ref>), introducing a weight proportional to the number of elements in the star of q. 
Then, we calculate the unsymmetrized induced charge and magnetization densities d n^NS/d u_s' β(q) 
and d m_δ^NS/d u_s' β(q) using Eqs. (19) and (<ref>). Finally, the complete 
responses are obtained through the following relationships:


    d n /d u_s' β(q) =  1/N_S∑_{S|f}O_τ(S)   [ ∑_γ S_γβd n^NS/d u_s̅'̅γ(q)
       × e^G_S^-1·r e^- q·R^S_s'],



    d m_δ/d u_s' β(q)    =  1/N_S∑_{S|f} (-1)^τ(S)O_τ(S)[ ∑_γηS̃^-1_δη S_γβ
       ×d m_η^NS/d u_s̅'̅γ(q) e^G_S^-1·r e^- q·R^S_s'],


where S̃ is the proper part of S, O_τ(S) is the identity if τ(S)=0, 
or O_τ(S) = K if τ(S)=1. Moreover, R^S_s' = 
Sd_s' - d_s̅'̅, where d_s' identifies the position of the atom s' with 
respect to the origin of its primitive cell, while d_s̅'̅ is obtained by applying the 
rotation S to the atom s' ({S|f} (R_ν + d_s') = R_ν̅ + d_s̅'̅). 
Similarly, the dynamical matrix becomes:


    Φ_s α
     s' β (q) = 1/N_S∑_{S|f}O_τ(S)   [ ∑_γδ S_γα  S_δβ Φ_s̅γ
    s̅'̅δ^NS (q) 
       × e^q· (R^S_s - R^S_s')],


where Φ_s̅γ
s̅'̅δ^NS (q) is obtained summing over the IBZ in 
Eq. (<ref>) and including the terms coming from d^2 E_tot^(1)/d u_μ s α d u_ν s' β, d^2 E_tot^(3)/d u_μ s α 
d u_ν s' β, and d^2 E_tot^(4)/d u_μ s α d u_ν s' β, defined in Ref. us_fr_dfpt.



§ APPLICATIONS

In this section we use the theory described above to compute the phonon dispersions of ferromagnetic fcc Ni 
and of a monatomic ferromagnetic Pt nanowire. We validate the theory by comparing the phonon frequencies 
obtained by diagonalizing the dynamical matrix (Eq. (<ref>)) with those obtained by the frozen phonon method.



 §.§ Computational details

First-principle calculations were performed within the Local Density Approximation (LDA) <cit.> and the Perdew-Burke-Ernzerhof 
(PBE) <cit.> schemes, as implemented in the Quantum ESPRESSO <cit.> and 
 <cit.> packages. The atoms are described by FR 
US PPs <cit.>, with 4s and 3d valence electrons for Ni 
(PPs  and  from pslibrary 0.1) and with 
6s and 5d valence electrons for Pt ( from pslibrary 1.0.0 
<cit.>).

DFPT calculations on ferromagnetic fcc Ni are at the theoretical LDA and PBE lattice constants, a = 6.483 a.u. 
and a=6.658 a.u., which are 2.6% and 0.02% smaller than experiment <cit.> (a=6.659 a.u.), 
respectively. The pseudo wavefunctions 
(charge density) are expanded in a plane waves basis set with a kinetic energy cut-off of 120 (600) Ry. The 
BZ integrations were done using a shifted uniform Monkhorst-Pack <cit.> k-point 
mesh of 28 × 28 × 28 points for the phonon calculations at a single wave vector q. 
The same computational parameters, except the k-point mesh which has been reduced to 18 × 18 × 18 
points, have been used for the phonon dispersions. The dynamical matrices have 
been computed by DFPT on a 6 × 6 × 6 q-point mesh, and Fourier interpolated to obtain the 
complete dispersions. 
Phonon frequencies of ferromagnetic Ni with the frozen phonon method, were calculated with 
a simple cubic supercell with 4 Ni atoms. The kinetic energy cut-offs used are the same as for the DFPT 
calculations, while the BZ integrations were performed on a k-point mesh of 24 × 24 × 24 points. 
The presence of a Fermi surface has been dealt with by the Methfessel-Paxton smearing method <cit.> with a 
smearing parameter σ=0.02 Ry. 

DFPT calculations on monatomic ferromagnetic Pt nanowire were done at a stretched geometry with interatomic distance 
d = 4.927 a.u.. The wire replicas have been separated by a vacuum space of 20 a.u.. We have checked that 
by increasing the vacuum space the computed frequencies do not change more than 0.2 cm^-1. The system has been 
studied in a ferromagnetic configuration, with magnetization either parallel or perpendicular to the wire. 
The kinetic energy cut-off was 60 (400) Ry for the wave functions (charge density). The k-point 
mesh is a shifted uniform Monkhorst-Pack mesh of 300 points. Frozen phonon calculations 
were performed with supercells with 2 and 4 Pt atoms, and Monkhorst-Pack meshes of 150 
and 75 k-points, respectively. The smearing parameter was σ = 0.002 Ry.
 


 §.§ Fcc Ni

We start our discussion from the computation of the phonon frequencies of ferromagnetic fcc Ni 
with the magnetization along [001] (and with a magnitude that turns out to be 0.62  μ_B per atom), 
and compare the DFPT and the frozen phonon method at the Y and Z points. The results obtained are reported 
in Table <ref>. The frequencies of the transverse modes at q = (0, 0, 2 π / a) (Z) are degenerate with 
both methods, as a consequence of the tetragonal magnetic symmetry (D_4h (C_4h)): indeed both transverse modes 
have atomic displacements perpendicular to the magnetization. Instead, the transverse modes at 
q = (0, 2 π / a, 0) (Y) show a small splitting of 0.04  cm^-1. The two modes are actually 
different because the atomic displacements are either parallel or perpendicular to the magnetization. A frequency 
splitting arises as a consequence of spin-orbit coupling. The DFPT and frozen phonon methods agree within 
0.3  cm^-1. The DFPT and the frozen phonon method predict the same splitting, which however is 
small compared to the agreement of the absolute values of the frequencies obtained with the two methods, hence 
it is not possible to give an accurate quantitative prediction, but only an order of magnitude. With the kinetic 
energy cut-offs and k-point mesh used, the frequencies obtained are converged 
within 5 × 10^-3 cm^-1, the same order of magnitude as the errorbar reported in Table <ref>
and due to the fit.




In Fig. <ref> we show the complete phonon dispersion of fcc Ni obtained by DFPT. Both LDA and PBE theoretical 
dispersions are shown, together with inelastic neutron scattering data <cit.>. The agreement between 
the LDA result and the experiment is poor, mainly because LDA underestimates the lattice constant: the 
highest frequencies of the dispersion (e.g. at the X and L points) are about 30  cm^-1 higher 
than the experiment. On the other hand, the PBE phonon dispersions are in excellent agreement with the experiment. 
Note however that this agreement is slightly worsened by temperature effects <cit.> not included 
in the present study.





 §.§ Pt monatomic wire





In this section we consider a monatomic Pt nanowire, a metal with ferromagnetic ordering. It has been 
shown <cit.> that at its equilibrium geometry (atomic distance d = 4.441 a.u.) 
the system shows a colossal magnetic anisotropy, since the preferred orientation of the magnetization is 
parallel to the wire and the magnetization vanishes when forced to be perpendicular to the wire. Instead, 
for stretched geometries with atomic distance higher than 4.913 a.u. a non-zero magnetization perpendicular 
to the wire is allowed. Here we consider a stretched geometry with d = 4.927 a.u. and compute the phonon 
dispersions with both a magnetization parallel and perpendicular to the wire.
In the following the nanowire is along the z direction. 
In Table <ref> we compare the phonon frequencies, at q = π/a and q = π/2a with 
m∥ x and m∥ z, computed by the DFPT and with the frozen phonon method. 
With a magnetization m∥ z (m = 0.65  μ_B per atom), the frequencies of the transverse modes 
are degenerate, while with m∥ x (m = 0.13  μ_B per atom) at q = π/a the two transverse 
modes show a splitting of about 0.5  cm^-1, which is of the same order of magnitude 
as the overall agreement of the two methods. At q = π/2a this splitting is about 
7  cm^-1, one order of magnitude larger than at q = π/a. In both cases the polarization of 
the transverse mode with higher frequency is perpendicular to the magnetization. As discussed above for fcc Ni, 
the two transverse modes are not equivalent due to the 
presence of the magnetization and of spin-orbit coupling. Pt atoms are heavier than Ni and 
show a stronger spin-orbit interaction: indeed, the splitting reported for Pt is 1-2 orders of magnitudes 
higher than in Ni. The DFPT and frozen phonon results agree within 0.4  cm^-1 on average. 
As before, the errorbars reported in Table <ref> come from the linear fit. With the kinetic energy cut-offs 
and the k-point mesh used all the frequencies reported are converged within 0.03  cm^-1 

In Fig. <ref> we show the phonon branches along Γ - Z for a ferromagnetic wire with magnetization 
parallel (left panel) or perpendicular to the wire (right panel). The two dispersions show evident differences: 
at q = π/a, the longitudinal mode for the wire with m∥ z is lower in frequency than for the wire 
with m∥ x, while the transverse modes are higher in frequency. 
In the central part of the BZ, around q = π/2a, the longitudinal mode of the 
wire with m∥ z has a higher frequency at the Kohn anomaly than the wire with m∥ x, 
while the transverse modes show a Kohn anomaly only for m∥ x. We remark that at the stretched 
geometry studied (d = 4.927 a.u.) the phonon modes are still stable, but the range of atomic distances at which 
both modes are stable is quite narrow. 






§ CONCLUSIONS

We extended the DFPT for lattice dynamics with FR NC 
and US PPs to the magnetic case. 
We validated the theory by comparing the DFPT to the frozen phonon method for ferromagnetic fcc Ni and for a 
monatomic ferromagnetic Pt nanowire. The agreement between the two methods is within 0.5  cm^-1. For 
both systems, we computed by DFPT also the complete phonon dispersions and discussed their features, 
showing that magnetism together with spin-orbit coupling may lift the degeneracy of some phonon modes. 
For our systems these splittings range from 10^-2 cm^-1 (in Ni) to a few cm^-1 
(in Pt nanowire). 


§ ACKNOWLEDGMENTS


Computational facilities have been provided by SISSA through its Linux Cluster and ITCS and by CINECA through 
the SISSA-CINECA 2018-2019 Agreement.

00







dfpt_review S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001)
us_fr_dfpt A. Dal Corso, Phys. Rev. B 76, 054308 (2007).
us_vanderbilt D. Vanderbilt, Phys Rev. B 41, 7892 (1990).
Pb_phonon A. Dal Corso, J. Phys. Condens. Matter 20, 445202 (2008).  
Giustino K. Cao, H. Lambert, P. G. Radaelli, and F. Giustino, Phys. Rev. B 97, 024420 (2018)
magnon T. Gorni, I. Timrov, and S. Baroni, Eur. Phys. J. B 91, 249 (2018)
KS W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
us_dfpt A. Dal Corso, Phys. Rev. B 64, 235118 (2001).
PZ J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
PBE J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
QE P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al., J. Phys. Condens. Matter 21, 395502 (2009) (See <http://www.quantum-espresso.org>).
QE_2 P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, et al., J. Phys. Condens. Matter 29, 465901 (2017).
thermo_pw  is an extension of the Quantum ESPRESSO (QE) package which provides an alternative
organization of the QE work-flow for the most common tasks. For more information see <https://dalcorso.github.io/thermo_pw>.
us_fr_pseudo A. Dal Corso and A. Mosca Conte, Phys. Rev. B 71, 115106 (2005). 
pslibrary A. Dal Corso, Comp. Mat. Sci. 95, 337 (2014).
pslibrary_2 See <https://dalcorso.github.io/pslibrary>.
COD R. W. G. Wyckoff, Crystal Structures, John Wiley & Sons, Inc., (Wiley, New York, 1963), Volume 1, 7-83.
k_grid J. Monkhorst and J.D. Pack, Phys. Rev. B 13, 5188 (1976).
MP M. Methfessel and A.T. Paxton, Phys. Rev. B 40, 3616 (1989).
exp_phonon_Ni R. J. Birgeneau, J. Cordes, G. Dolling, and A. D. B. Woods, Phys. Rev. 136, A1359 (1964).
Ni_temp_phonon A. Dal Corso, J. Phys. Condens. Matter 25, 145401 (2013).
Pt_nature A. Smogunov, A. Dal Corso, A. Delin, R. Weht, and E. Tosatti, Nature Nanotechnology 3, 22 (2008). 
Pt_PRB A. Smogunov, A. Dal Corso, and E. Tosatti, Phys. Rev. B 78 014423 (2008).


