Here I'll explain the creation of this result: https://twitter.com/Jules_Oppenheim/status/1576596723933200384
In a recent paper (https://doi.org/10.1002/anie.202213960), Kampouri et al. synthesize a metal-organic framework from 3,3’,4,4’,5,5’-hexahydroxybiphenyl and FeCl2. The MOF crystallizes into a net with a two-dimensional secondary building unit conisting of alternating octahedral and square planar iron atoms, bridged by the oxygen of the hexahydroxybiphenyl linkers. Topologically, this SBU is equivalent to a cpc net (http://rcsr.net/layers/cpc). In this net, there are two unique vertices (corresponding to a Fe connected to 3 other Fe atoms and a Fe connected to 4 other Fe atoms) and two unique edges (corresponding to a 3-connected Fe : 4-connected Fe linkage and a 4-connected Fe : 4-connected Fe linkage).
Slice of the crystal struccture in the bc-plane, depicting the 2D SBU
Topological simplification of the 2D SBU, depicting the connectivity of Fe bridges by oxygens
During the synthesis of this MOF, the iron atoms undergo a partial oxidation to achieve a final composition of about (Fe2+,Fe3+2). For the sake of a simple calculation of the magnetism, let's assume that the MOF is fully oxidized to Fe3+ and that the lattice of high spin S = 5/2 atoms behave similarly to a lattice of S = 1/2 atoms, with only nearest neighbor interactions.
There are multiple ways to calculate the magnetic susceptibility for an lattice of spins (in the absence of an external magnetic field). If the system is small or simple, the result can be calculated analytically. However, for an arbitrary system, there is no analytic solution. In these cases, it makes more sense to perform a Monte Carlo simulation to numerically calculate the result.
Here, the cpc net can be written as an n x n matrix where n is the number of included vertices. As n gets large, the simulated result should converge to the correct value
Each of these Hamiltonian graphs, H, is encoded into matrix form. For example, here is the n = 4 matrix (the 32 x 32 matrices take some time to prepare...)
| 0 | -J2 | -J1 | -2J2 |
| -J2 | 0 | -2J2 | 0 |
| -J1 | -2J2 | 0 | -J2 |
| -2J2 | 0 | -J2 | 0 |
To calculate the energy of a certain state A (i.e. a n x 1 array where each index is 5/2 or -5/2 depicting spin up or down), the dot product of A.H.A is computed.
Now to perform a Monte Carlo simulation, we can run the following code
# count is set to 1,000,000# grid is the array A# ham is the matrix Hfor i in xrange(count): for q in xrange(4): # Choose random point on grid and flip x = int(random.random()*length) z = random.random() grid2 = np.copy(grid) grid2[x] *= -1 # Calculate energy the smart way diff = np.dot(grid2,np.dot(ham,grid2)) - e # If energy is less, accept, else, accept with some probability if diff < -1: p = 1 else: p = min(1.0, math.exp(-(diff)/(0.695034800*t))) num = int(z + (1-p)) # 0 if accept and 1 if reject if num == 0: grid = grid2 e += (num^1) * diff m = np.sum(grid) * 5/2
And the magnetic susceptibility is calculated by the difference of the average of the square of magnetization and the average of the magnetization squared, divided by the temperature.
Using this code for one million cycles per data point and calculating the magnetic susceptibility every 3 degrees from 0 K to 300 K takes about 2 hours of computational time on a laptop
The result tends to converge to at n = 32
The result for n = 32 as a function of the J1/J2 ratio is as follows
When both J1 and J2 are negative, the system is strongly anti-ferromagnetic. As J2 becomes larger, we start to see the increase in chiT.
When J2 becomes positive, ChiT becomes even larger, reaching a critical point at J2 = 1, where the system begins to behave as a ferromagnet at low T.
There are a few special cases here.
- When J2 = 0, the system reduces down to a collection of independent dimers with coupling constant J1.
- When J2 = 0.5 J1 and J2 < 0, there are two degenerate ground states (depicting spin up with purple and spin down with yellow circles)
We can observe these special points by looking at the eigenvalues of the Hamiltonian for n = 16. Notice the degeneracy in the ground state when J2 = 0.5 J1 and when J2 = -0.5 J1
