Perovskite to Pseudospin: A Julia Challenge Part 2

Continued from the previous post. Mid way through the project, I realized that I had written the expression for PBC a tiny bit wrong and that led to disastrous calculation of U value. So in this post, let us rectify that and add dipole terms to the system.

Corrected Periodic Boundary Conditions

PBC for a given index turned out to be (after hours of breaking the nerves in my fried brain):

function pbc_1(system, index)
    eval_1(i,n) = ((abs(i)-1) ÷ n) *
        ((sign(i)+1) ÷ 2) -
        ((sign(i-1)-1) ÷ 2) * ((i-n) ÷ n)
    index_1(i,n) = ((abs(i+n-1) % n)+1)
    unit_size = size(system)[1]
    a = (-system[1,1,1].ge +
         system[1,1,2].ge)[3] * unit_size
    uc = eval_1.(index, unit_size) * a
    index_mod = index_1.(index, unit_size)
    return index_mod, uc
end

Dipole-Dipole Interaction

The dipole interaction energy is derived and given by

Where is the dipole vector and is the vector connecting the dipoles. For simplicity, we will ignore the pre-factors and code just the form of it. This doesn't remove any physics as the strength is arbitrary and controlled by the coupling between dipole energy and "bonding" energy.

The total energy is:

The beauty of Julia is that it is both a Lisp language as well as a JIT one, thus we can plot beautiful objects that are as fast as FORTRAN. Rotating the tetrahedron along the direction (in which it is oriented) should not change the dipole, but because of the symmetry of the tetrahedron, one should see a 3-fold symmetry in nearest-neighbor bonding energy, which is exactly what we observe.

Next stop: the full Monte Carlo simulation. The GitHub repo has the full code.