Hello,
I am trying to calculate mean field temperature of a structure using HeisenbergMapper, but I have some basic questions about how it is implemented.
Here are the some lines from pymatgen/heisenberg.py at 489325fac44bf45e629907625587f551972adc32 · materialsproject/pymatgen · GitHub starting at line number 551.
else: # multiple magnetic sublattices
omega = np.zeros((num_sublattices, num_sublattices))
ex_params = self.ex_params
ex_params = {k: v for (k, v) in ex_params.items() if k != "E0"} # ignore E0
for k in ex_params:
# split into i, j unique site identifiers
sites = k.split("-")
sites = [int(num) for num in sites[:2]] # cut 'nn' identifier
i, j = sites[0], sites[1]
omega[i, j] += ex_params[k]
omega[j, i] += ex_params[k]
omega = omega * 2 / 3 / k_boltzmann
eigenvals, eigenvecs = np.linalg.eig(omega)
mft_t = max(eigenvals)
The part that doesnt make sense to me is that in the for loop, it loops over unique sites, but there is no information on the coordination number of those sites with each other. Ashcroft Mermin eq 33.58 (pg. 715) says that for a FM model, all identical sites, you should just sum all the J® in your repeat unit. I dont think that it is the case in the code above because it only loops over unique sites, but does not take into account how they are coordinated. I would appreciate if you could share your feedback on this. This is not a bug, thats why I didnt post it on github, but I would appreciate if we could discuss.
Additionally, I have some code here GitHub - kayahans/heisenberg_mapper: Simple extension of the HeisenbergMapper module in pymatgen that allows any nearest neighbors that I have built based on your code to include beyond nearest neighbor interactions. I am not sure why they are not done in your version yet. Is there a physical reason that I am missing which makes things more complicated? Besides, with the simple routine I have in that code, you can use all the possible magnetic orderings from MagneticOrderingsWF from atomate to fit a hamiltonian (you can see the fit in that page). In your version of the code, you decided to select only the lowest energy orderings, and solve the hamiltonian exactly. I found that a bit unreliable especially when you are looking at weak couplings. The final J values seem to depend critically how the routine selects specific magnetic orderings for weak interactions, because some low energy orderings seem to be excluded unfortunately because the hamiltonian becomes uninvertible and then they are substituted with a higher energy ordering. But for stronger interactions, I found that both the fitting I made and your exact method seems to work closely.
Thanks,
Kayahan