Caluculation of exchange coupling coefficients using the four-state method and \(T_c\) using the Monte Carlo simulations

Author: Yu Qisheng

---

In this tutorial, I will show how to calculate exchange coupling coefficients using the four-state method with PASP (Property analysis and simulation package for materials). The example system is VHfO\(_4\), where the magnetic atoms is vanadium, DFT+U method is applied in some self-consistent energy calculations with VASP. Then calculate \(T_c\) using the Monte Carlo simulations.

The principle of the four-state method is very easy: \(\(J_{ij}=\frac{E_{ij,\uparrow\uparrow}+E_{ij,\downarrow\downarrow}-E_{ij,\uparrow\downarrow}-E_{ij,\downarrow\uparrow}}{4}\)\)

So what we need to do is find the environmentally inequivalent magnetic atom pairs and calculation the energies in different magnetic moments.

The first step is finding the environmentally inequivalent magnetic atom pairs \(A_{ij}\).

To begin, prepare an optimized unitcell structure file(POSCAR) with corresponding pseudopotential file(POTCAR). Create a new directory fourstates, using the post_VASP script to construct a 222 supercell(usage: ./post_VASP \(\rightarrow\) 1 \(\rightarrow\) 2 2 2), copy the generated POSCAR.multi file as the new POSCAR in a new created subdirectory input1, and the original POSCAR as POSCAR_UC. Then modify the PASP.input file:

Model find_J_path_sym

#Read the prepared POSCAR and POSCAR_UC file
READ_POSCAR T
READ_POSCAR_UC T

#The scale of supercell
%block SuperCell
2 0 0
0 2 0
0 0 2
%endblock SuperCell

#The atomic number of the different ions, the order correspond to POSCAR file
%block Z_values
23 72 8
%endblock Z_values

#The magnetic moments value of the different ions in ${\mu}B$
%block Spin_values
1 0 0
%endblock Spin_values

#1 for magnetic atoms; 2 for ligand atoms; 0 for others
%block magnetic_or_ligand
1 0 2
%endblock magnetic_or_ligand

#Cutoff distance in Bohr, the distance of $A_{ij}$ should be in this range
%block fig.mdist
1 10.5
%endblock fig.mdist

#The index of elements whose magnetic interactions need to be calculated(count start at 1)
%block isaAB
1 1
%endblock isaAB

And prepare a spin3.in file to create an inital spin configuration, each row corresponds to the magnetic moment of the corresponding atom in POSCAR.

Before you sbatch pasp.slurm, you should check if the files correct by ls:

PASP.input pasp.slurm POSCAR POSCAR_UC spin3.in

After step1 calculation(often takes few seconds), if PASP is finished gracefully, copy the output files cell.str, unit_cell.str,pair_POS.*,J.POS and pasp.out(change the name to log.out) to the new created subdirectory input2, and prepare the VASP SCF calculation files except POSCAR. Modify the run_DFT.sh script so that the parameters match your needs, then run this script. Sometimes the SCF calculations may meet some problems such as convergence or compute node issues, so check each SCF calculations before you start the next step. After checking, modify the PASP.input file: add a line cal_spin_exchange T after Model find_J_path_symi, change the READ_POSCAR T and READ_POSCAR_UC T to %include cell.str and %include unit_cell.str.

Again, before you start the calculation, check the files by ls:

cell.str INCAR J1 J2 J3 J4 J5 J.POS KPOINTS log.out vasp.sbatch pair_POS.001 pair_POS.002 pair_POS.003 pair_POS.004 pair_POS.005 PASP.input pasp.slurm POTCAR run_DFT.sh unit_cell.str

The number of J*directories and pair_POS.* files shoule be equal to the number of environmentally inequivalent magnetic atom pairs.

Before do the next PASP calculation, do check all the scf calculations are correct! e.g. all the magnetic moments are maintained at the set value. If some moments changed, you can using the specially complied VASP to fix the moments. The example sbatch script follows:

#!/bin/bash
#SBATCH --job-name=vasp
#SBATCH --nodes=1
#SBATCH --cpus-per-task=1
#SBATCH --ntasks-per-node=64
#SBATCH --mem=489789M
#SBATCH --qos=huge

    module load vasp/5.4.4-MCONSTR2-fast
    readonly map_opt='ppr:8:l3cache:pe=1'

export OMP_NUM_THREADS=1

mpirun -np 64 -map-by $map_opt vasp_ncl

The ncl VASP is necessary, do not change it! Then add or change these lines in your INCAR file:

COL_CONSTRAINED_M = 1/2
M_CONSTR = * * * *
LAMBDA = $\lambda$
MAGMOM = * * * *
RWIGS = a b c

This part of the settings applies to collinear magnetic moments. The value of COL_CONSTRAINED_M corresponds to two methods, 1 for methond [1] and 2 for method [2]; the values of M_CONSTR are the moments you want to maintain; the value of LAMBDA should be larger much smaller than 1 for method1 or larger than 5 for method2; the values of MAGMOM are the initial magnetic moments; the values of RWIGS are Wigner-Seitz radius, you can find them in POTCAR.

For non-collinear magnetic moments, new options 5 and 7 have been added to the existing I_CONSTRAINED_M schemes (1, 2, 3). Option 5 implements method [1], while option 7 implements method [2]. When using I_CONSTRAINED_M = 5, in addition to the parameters required for I_CONSTRAINED_M = 1 (such as LAMBDA and M_CONSTR), you must also set H_LAMBDA to specify the magnitude of the applied magnetic field (typically recommended to be much less than 1) and H_CONSTR to define its direction and magnitude (without normalization, following the same format as M_CONSTR). In contrast, I_CONSTRAINED_M = 7 requires the same parameter setup as I_CONSTRAINED_M = 1.

If you check all the scf calculations are correct, you can do the next PASP calculation by sbatch pasp.slurm.

After few seconds, if PASP is finished gracefully, check the spin_exchange.dat file to get the information about exchange coupling coefficients:

5 # units in eV and Ang
      5.370859943590      0.000000000000      0.000000000000
      0.000000000000      4.958931042774      0.000000000000
     -0.508019826630      0.000000000000      4.938146270730
1      -0.013740000000  1  1      0.000000000000      0.000000000000      0.000000000000      2.685429971795     -2.276101291285      0.000000000000
2       0.002000000000  1  1      0.000000000000      0.000000000000      0.000000000000      2.685429971795      2.682829751489      0.000000000000
3      -0.025000000000  1  1      0.000000000000      0.000000000000      0.000000000000      0.000000000000      4.958931042774      0.000000000000
4       0.000180000000  1  1      0.000000000000      0.000000000000      0.000000000000     -0.508019826630      0.000000000000      4.938146270730
5      -0.029620000000  1  1      0.000000000000      0.000000000000      0.000000000000      5.370859943590      0.000000000000      0.000000000000
0 # no dJ/du, i.e., no spin-lattice interaction
0 # single site anisotropy
0 # DM
0 # higher order interaction

The values of the second column are the exchange coupling coefficients between the atoms with the position following columns.

Copy the spin_exchange.dat to the new created subdirectory MC, and modify the PASP.input file(compared to input1): add a line %include MC.input at the head, change the Model find_J_path_sym to Model Spin_lattice_interaction. And prepare a MC.input file:

#number of intervals
PT.ninterval 1

#Whether to enable the MC function
MC T

#Under the premise that the MC function is enabled, whether to enable the PTMC function
Parallel_Tempering T

#number of replicas
PT.M 64 

#lowest T (in Hartree)
PT.low_temp  0.00000317 #1K

#highest T (in Hartree)
PT.high_temp  0.001585  #500K

#PT.Tmesh
# 1: T are evenly distributed according to 1/T
# 2: T are evenly distributed according to T
PT.Tmesh  2

# One PT iteration contains nblk_T exchange steps
# One exchange step contains nblk MCS steps
# One MCS step contains ns_blk (set to number of magnetic sites in the supercell) local moves

PT.niter   1  #number of PT iteration
PT.nblk_T  4000

nblk  500
#nbeg_blk: number of exchange steps to reach equilibrium, after that we will do measurement
nbeg_blk 2000

For Monte Carlo simulations, prepare a larger supercell(such as 666,888) by using the post_VASP script mentioned before as the POSCAR file, also prepare the POSCAR_UC file.

Checking the files by ls:

MC.input PASP.input POSCAR POSCAR_UC pasp.slurm spin_exchange.dat

After few hours(maybe longer, depends on your system), if PASP is finished gracefully, check the C_PTMC_K files, the first column is temprature and the second column is heat capacity, and the temprature correspond to the maximum heat capacity is the \(T_c\) we want to calculate. Further, you can open the lowest_onlyspin.xsf file in VESTA to see the magnetic ordering.

References [1] Y. Chen, Y. Yang, C. Xu, and H. Xiang, Constraining spin directions in density functional theory calculations by imposing a local magnetic field, Phys. Rev. B 107, 214439 (2023). [2] P.W. Ma and S. L. Dudarev, Constrained density functional for noncollinear magnetism, Phys. Rev. B 91,054420 (2015).