import pylab

# Co data from Table I, Pajda et at. PRB 64, 174402, (2001)
Co_data = numpy.array([
    1.085,
    0.110,
    0.116,
    -0.090,
    0.026,
    0.043,
    -0.024,
    0.012,
    0.026,
    0.006])


filename = 'Co-FCC.hdf5'

# Load the HeisenbergExchange object
heisenberg_exchange = nlread(filename, HeisenbergExchange)[0]

# Get the distances, unique coupling elements, and multiplicity.
distances, Jij, m = heisenberg_exchange.uniqueCouplingMatrixElementsAndDistances()

# Get the coupling elements in units of mRy
Jij = Jij.inUnitsOf(Rydberg)*1e3

# Plot the coupling elements J_ij vs. distance r_ij
pylab.figure()
pylab.plot(distances, Jij, 'ko-', label='QATK')
pylab.plot(distances[:len(Co_data)], Co_data, 'rs', label='Ref.')
pylab.xlabel('$r_{ij}$ (Ang.)')
pylab.ylabel('$J_{ij}$ (mRy)')
pylab.grid('on')
pylab.xlim([0, 12])
pylab.legend(loc=0)
pylab.title(filename.split('.')[0])
pylab.show()