Hello Ponce,
Thank you very much for the response. I understand that this topic isn’t recent but I would still appreciate if you could consider the following:
 Considering
temperature_final.py
, I don’t think the SOC would not work by default:
This segment in rf_final.py
code:
occtmp = EIGR2D.nsppol*EIGR2D.occ[:,:,:]/2 # jband # should be 1 !
delta_E_ddw = N.einsum('lij,k>lijk',eig0[:,:,:].real,N.ones(EIGR2D.nband))  \
N.einsum('lij,k>likj',eig0[:,:,:].real,N.ones(EIGR2D.nband))  \
N.einsum('i,ljk>ljik',N.ones((EIGR2D.nband)),(2*occtmp1))*smearing*1j # spin,ikpt,iband,jband
ddw_tmp = N.einsum('ijkln,lm>mijkn',ddw_addQ,2*bose+1.0) # itemp,ikpt,iband,jband,ispin
ddw_add = N.einsum('ijklm,mjkl>imjk',ddw_tmp,1.0/delta_E_ddw) # temp,spin,ikpt,iband
delta_E = N.einsum('lij,k>lijk',eig0[:,:,:].real,N.ones(EIGR2D.nband))  \
N.einsum('lij,k>likj',eigq.EIG[:,:,:].real,N.ones(EIGR2D.nband)) # spin,ikpt,iband,jband
delta_E_sm = N.einsum('i,ljk>ljik',N.ones((EIGR2D.nband)),(2*occtmp1))*smearing*1j # spin,ikpt,iband,jband
num1 = N.einsum('ij,mkl>mkijl',bose,N.ones((EIGR2D.nsppol,EIGR2D.nkpt,EIGR2D.nband))) +1.0 \
 N.einsum('ij,mkl>mkijl',N.ones((3*EIGR2D.natom,ntemp)),occtmp) # spin,k,mod,temp,band # bef was (imode,tmp,band)
deno1 = N.einsum('mijk,l>mijkl',delta_E,N.ones(3*EIGR2D.natom),dtype=N.complex)
Has a problem with SOC, since occupations are 1 per band (instead of 2) so this array occtmp
gets filled with 1/2 which results division with 0 in ddw_add
i.e. (2*occtmp1)
in delta_E_ddw
becomes zero.
I suspect that line is related to eq (15) and (16) in Ponce2015. There the 1j*smearing
is the adhoc parameter of numerical nature i*delta
inferred as finite lifetime. So, the occtmp
variable is an array over states that takes 1 for occupied and 0 for unoccupied in nonSOC case (0.5 and 0 for SOC) so the delta_E_ddw
would suggest that i*delta
changes sign in equations (15) and (16) from + to  for unoccupied states. I am confused why is that the case?
Anyhow, changing (2*occtmp1)
to N.ones(occtmp.shape)
in delta_E_ddw
and other places where it occurs can remove this division with 0, when dealing with SOC case. Since I am only interested in temperature dependence, that takes only the real part of temperature dependent renormalization delta Epsilon
for adiabatic/dynamic case of R.I.A. , I guess the sign of i*delta
shouldn’t be an issue?
There is a part of rf_final.py
the code that calculates average renormalization for degenerate states def make_average
that for SOC should consider 4fold degenerate states, that part I edited too but the results of my renormalization do not look promising.

Another question that I have is about III B section of Ponce2015, Restoration of charge neutrality:
Since polar materials already exibit nonzero Born effective charges, will the restoration of charge neutrality help the convergence over kpoints for them like it does for nonpolar materials? Is it even correct to impose charge neutrality in that case? If so, why? Does the code add the longrange contribution to GKK when getddb
is activated for elph interaction ( phonon=1
and ieig2rf=5
).

I use the Legacy code since I was under the impression that SOC is not supported in the new EPH code? Is that still the case or has the support been added, because I couldn’t find it in release notes? In the new EPH tutorial, ephintro, long range dipolar fields are adressed as in Verdi2015, so does the code know when it is dealing with polar or nonpolar materials?
Kind regards,
Milan Jocic