Allowing to pack Hermitian atomic matrices?
Firstly, pack does not pack Hermitian (but complex symmetric matrices). In other words, it does exactly what is said in the documentation
Pack a 2D array to 1D, adding offdiagonal terms.
The matrix::
/ a00 a01 a02 \
A = | a10 a11 a12 |
\ a20 a21 a22 /
is transformed to the vector::
(a00, a01 + a10, a02 + a20, a11, a12 + a21, a22)
but typically this is used to pack the atomic density matrix which is Hermitian, so I think this is not desired behavior
So when calculating the atomic density matrices, one neglects the imaginary part of density matrix first here
D_sii[kpt.s] += np.dot(P_ni.T.conj() * f_n, P_ni).real
and then for a second time here
D_sp[:] = [pack(D_ii) for D_ii in D_sii]
So, atomic density matrices may only be used to calculate expectation values of observables with purely real representation in the partial wave basis. For example, not current.
If one would allow the packed atomic density matrix to be complex, the traces should be changed from np.sum(D_p * H_p).real to np.sum(D_p.imag * H_p.imag) + np.sum(D_p.real * H_p.real),
where D_p is packed with pack complex conjugating the other offdiagonal term and H_p is packed with pack2.
I don't know how much work this would be, but it would help introducing the gauge including PAW structure.
from gpaw.utilities import pack, pack2, unpack, unpack2
import numpy as np
H = np.array([ [0, 1j], [-1j,0]])
print pack(H)
# [ 0.+0.j 0.+0.j 0.+0.j]
print pack2(H)
#[ 0.+0.j 0.+1.j 0.+0.j]
print unpack2(pack(H)) # Packing and unpacking a hermitian matrix returns just zeros!!!
#[[ 0.+0.j 0.+0.j]
# [ 0.+0.j 0.+0.j]]
print unpack(pack2(H))
#[[ 0.-0.j 0.+1.j]
# [ 0.-1.j 0.-0.j]]