Merge pull request #49 from antelmor/master

Merge algorithm for non-collinear systems

TB2J/rotate_atoms.py

@@ -5,19 +5,24 @@
 from TB2J.tensor_rotate import Rxx, Rxy, Rxz, Ryx, Ryy, Ryz, Rzx, Rzy, Rzz
 
 
-def rotate_atom_xyz(atoms):
+def rotate_atom_xyz(atoms, noncollinear=False):
     """
-    given a atoms, return:
-    - 'z'->'x'
-    - 'z'->'y'
-    - 'z'->'z'
+    given a atoms, return rotated atoms:
+    atoms_1, ..., atoms_n,
+    where we considered n diffeerent roation axes.
+
+    When noncollinear == True, more rotated structures
+    will be generated.
     """
-    atoms_x = copy.deepcopy(atoms)
-    atoms_x.rotate(90, "y", rotate_cell=True)
-    atoms_y = copy.deepcopy(atoms)
-    atoms_y.rotate(90, "x", rotate_cell=True)
-    atoms_z = atoms
-    return atoms_x, atoms_y, atoms_z
+
+    rotation_axes = [(1, 0, 0), (0, 1, 0)]
+    if noncollinear:
+        rotation_axes += [(1, 1, 0), (1, 0, 1), (0, 1, 1)]
+
+    for axis in rotation_axes:
+        rotated_atoms = copy.deepcopy(atoms)
+        rotated_atoms.rotate(90, axis, rotate_cell=True)
+        yield rotated_atoms
 
 
 def rotate_atom_spin_one_rotation(atoms, Rotation):
@@ -96,18 +101,15 @@ def check_ftype(ftype):
         print("=" * 40)
 
 
-def rotate_xyz(fname, ftype="xyz"):
+def rotate_xyz(fname, ftype="xyz", noncollinear=False):
     check_ftype(ftype)
     atoms = read(fname)
     atoms.set_pbc(True)
 
-    atoms_x, atoms_y, atoms_z = rotate_atom_xyz(atoms)
+    rotated = rotate_atom_xyz(atoms, noncollinear=noncollinear)
 
-    fname_x = f"atoms_x.{ftype}"
-    fname_y = f"atoms_y.{ftype}"
-    fname_z = f"atoms_z.{ftype}"
+    for i, rotated_atoms in enumerate(rotated):
+        write(f"atoms_{i+1}.{ftype}", rotated_atoms)
+    write(f"atoms_0.{ftype}", atoms)
 
-    write(fname_x, atoms_x)
-    write(fname_y, atoms_y)
-    write(fname_z, atoms_z)
-    print(f"The output has been written to {fname_x}, {fname_y}, {fname_z}")
+    print(f"The output has been written to the atoms_i.{ftype} files. atoms_0.{ftype} contains the reference structure.")

TB2J/rotate_siestaDM.py

@@ -0,0 +1,36 @@
+import sisl
+
+def rotate_siesta_DM(DM, noncollinear=False):
+
+    angles_list = [ [0.0, 90.0, 0.0], [0.0, 90.0, 90.0] ]
+    if noncollinear:
+        angles_list += [[0.0, 45.0, 0.0], [0.0, 90.0, 45.0], [0.0, 45.0, 90.0]]
+
+    for angles in angles_list:
+        yield DM.spin_rotate(angles)
+
+def read_label(fdf_fname):
+
+    label = 'siesta'
+    with open(fdf_fname, 'r') as File:
+        for line in File:
+            corrected_line = line.lower().replace('.', '').replace('-', '')
+            if 'systemlabel' in corrected_line:
+                label = line.split()[1]
+                break
+
+    return label
+
+def rotate_DM(fdf_fname, noncollinear=False):
+
+    fdf = sisl.get_sile(fdf_fname)
+    DM = fdf.read_density_matrix()
+    label = read_label(fdf_fname)
+
+    rotated = rotate_siesta_DM(DM, noncollinear=noncollinear)
+
+    for i, rotated_DM in enumerate(rotated):
+        rotated_DM.write(f"{label}_{i+1}.DM")
+    DM.write(f"{label}_0.DM")
+
+    print(f"The output has been written to the {label}_i.DM files. {label}_0.DM contains the reference density matrix.")

docs/src/rotate_and_merge.rst

@@ -3,32 +3,40 @@
 Averaging multiple parameters
 ===============================
 
-As discussed in the previous section, the :math:`z` component of the DMI, and the :math:`xz`, :math:`yz`, :math:`zz`, :math:`zx`, :math:`zy`, components of the anisotropic exchanges are non-physical, and an :math:`xyz` average is needed to get the full set of magnetic interaction parameters if the spins are all along :math:`z`. In this case, scripts to rotate the structure and merge the results are provided, they are named TB2J\_rotate.py and TB2J\_merge.py. The TB2J\_rotate.py reads the structure file and generates three files containing the :math:`z\rightarrow x`, :math:`z\rightarrow y`,and the non-rotated structures. The output files are named atoms\_x, atoms\_y, atoms\_z. A large number of output file formats is supported thanks to the ASE library and the format of the output structure files is provided using the --format parameter. An example for using the rotate file is:
+When the spins of sites :math:`i` and :math:`j` are along the directions :math:`\hat{\mathbf{m}}_i` and :math:`\hat{\mathbf{m}}_j`, respectively, the components of :math:`\mathbf{J}^{ani}_{ij}` and :math:`\mathbf{D}_{ij}` along those directions will be unphysical. In other words, if :math:`\hat{\mathbf{u}}` is a unit vector orthogonal to both :math:`\hat{\mathbf{m}}_i` and :math:`\hat{\mathbf{m}}_j`, we can only obtain the projections :math:`\hat{\mathbf{u}}^T \mathbf{J}^{ani}_{ij} \hat{\mathbf{u}}` and :math:`\hat{\mathbf{u}}^T \mathbf{D}_{ij} \hat{\mathbf{u}}`. Notice that for collinear systems, there will be two orthonormal vectors :math:`\hat{\mathbf{u}}` and :math:`\hat{\mathbf{v}}` that are also orthogonal to :math:`\hat{\mathbf{m}}_i` and :math:`\hat{\mathbf{m}}_j`. 
+
+The projection for :math:`\mathbf{J}^{ani}_{ij}` can be written as
+
+:math:`\hat{\mathbf{u}}^T \mathbf{J}^{ani}_{ij} \hat{\mathbf{u}} = \hat{J}_{ij}^{xx} u_x^2 + \hat{J}_{ij}^{yy} u_y^2 + \hat{J}_{ij}^{zz} u_z^2 + 2\hat{J}_{ij}^{xy} u_x u_y + 2\hat{J}_{ij}^{yz} u_y u_z + 2\hat{J}_{ij}^{zx} u_z u_x,`
+
+where we considered :math:`\mathbf{J}^{ani}_{ij}` to be symmetric. This equation gives us a way of reconstructing :math:`\mathbf{J}^{ani}_{ij}` by performing TB2J calculations on rotated spin configurations. If we perform six calculations such that :math:`\hat{\mathbf{u}}` lies along six different directions, we obtain six linear equations that can be solved for the six independent components of :math:`\mathbf{J}^{ani}_{ij}`. We can also reconstruct the :math:`\mathbf{D}_{ij}` tensor in a similar way. Moreover, if the system is collinear then only three different calculations are needed.
+
+To account for this, TB2J provides scripts to rotate the structure and merge the results; they are named TB2J\_rotate.py and TB2J\_merge.py. The TB2J\_rotate.py reads the structue file and generates three(six) files containing the rotated structures whenever the system is collinear (non-collinear). The --noncollinear parameters is used to specify wheter the system is noncollinear. The output files are named atoms\_i (i = 0, ..., 5), where atoms\_0 contains the unrotated structure. A large number of file formats is supported thanks to the ASE library and the output structure files format is provided through the --format parameter. An example for using the rotate file with a collinear system is:
 
 .. code-block:: shell
 
    TB2J_rotate.py BiFeO3.vasp --ftype vasp
 
+If te system is noncollinear, then we run the following instead:
 
-.. note::
-
-    Some file format does not include the cell orientation, e.g. the cif file. They should not be used as the output format.
-
+.. code-block:: shell
 
-The user has to perform DFT single point energy calculations for these three structures in different directories, keeping the spins along the $z$ direction, and run TB2J on each of them. After producing the TB2J results for the three rotated structures, we can merge the DMI results with the following command by providing the paths to the TB2J results of the three cases:
+   TB2J_rotate.py BiFeO3.vasp --ftype vasp --noncollinear   
 
-::
+The user has to perform DFT single point energy calculations for the generated structures in different directories, keeping the spins along the $z$ direction, and run TB2J on each of them. After producing the TB2J results for the rotated structures, we can merge the DMI results with the following command by providing the paths to the TB2J results of the three cases:
 
-   TB2J_merge.py BiFeO3_x BiFeO3_y BiFeO3_z --type structure
+::
 
+   TB2J_merge.py BiFeO3_1 BiFeO3_2 BiFeO3_0
 
-Note that the whole structure are rotated w.r.t. the laboratory axis but not to the cell axis. Therefore, the k-points should not be changed in both the DFT calculation and the TB2J calculation. 
+Here the last directory will be taken as the reference structure. Note that the whole structure are rotated w.r.t. the laboratory axis but not to the cell axis. Therefore, the k-points should not be changed in both the DFT calculation and the TB2J calculation. 
 
 A new TB2J\_results directory is then made which contains the merged final results. 
 
-Another method is to do the DFT calculation with spins along the :math:`x`, :math:`y` and :math:`z` axis, respectively, and then merge the result with:
+Another method is to do the DFT calculation with spins rotated globally. That is they are rotated with respect to an axis, but their relative orientations remain the same. This can be specified in the initial magnetic moments from a DFT calculation. For calculations done with SIESTA, there is a script that rotates the density matrix file along different directions. We can then use these density matrix files to run single point calculations to obtain the required rotated magnetic configurations. An example is:
 
 ::
 
-   TB2J_merge.py BiFeO3_x BiFeO3_y BiFeO3_z --type spin
+   TB2J_rotateDM.py --fdf_fname /location/of/the/siesta/*.fdf/file
 
+As in the previous case, we can use the --noncollinear parameter to generate more configurations. The merging process is performed in the same way.

scripts/TB2J_merge.py

@@ -2,7 +2,7 @@
 import argparse
 import os
 import sys
-from TB2J.io_merge import merge, merge2
+from TB2J.io_merge import merge
 
 
 def main():
@@ -28,11 +28,18 @@ def main():
         type=str,
         default="TB2J_results",
     )
+    parser.add_argument(
+        "--main_path",
+        help="The path containning the reference structure.",
+        type=str,
+        default=None
+    )
 
     args = parser.parse_args()
     # merge(*(args.directories), args.type.strip().lower(), path=args.output_path)
     # merge(*(args.directories), method=args.type.strip().lower(), path=args.output_path)
-    merge2(args.directories, args.type.strip().lower(), path=args.output_path)
+    #merge2(args.directories, args.type.strip().lower(), path=args.output_path)
+    merge(*args.directories, main_path=args.main_path, write_path=args.output_path)
 
 
 main()

scripts/TB2J_rotate.py

@@ -14,9 +14,14 @@ def main():
         default="xyz",
         type=str,
     )
+    parser.add_argument(
+        "--noncollinear",
+        action="store_true",
+        help="If present, six different configurations will be generated. These are required for non-collinear systems."
+    )
 
     args = parser.parse_args()
-    rotate_xyz(args.fname, ftype=args.ftype)
+    rotate_xyz(args.fname, ftype=args.ftype, noncollinear=args.noncollinear)
 
 
 if __name__ == "__main__":

scripts/TB2J_rotateDM.py

@@ -0,0 +1,21 @@
+#!/usr/bin/env python3
+import argparse
+from TB2J.rotate_siestaDM import rotate_DM
+
+def main():
+    parser = argparse.ArgumentParser(description="")
+    parser.add_argument(
+        "--fdf_fname", help="Name of the *.fdf siesta file."
+    )
+    parser.add_argument(
+        "--noncollinear",
+        action="store_true",
+        help="If present, six different configurations will be generated. These are required for non-collinear systems."
+    )
+
+    args = parser.parse_args()
+    rotate_DM(args.fdf_fname, noncollinear=args.noncollinear)
+
+
+if __name__ == "__main__":
+    main()

setup.py

@@ -19,6 +19,7 @@
         "scripts/siesta2J.py",
         "scripts/abacus2J.py",
         "scripts/TB2J_rotate.py",
+        "scripts/TB2J_rotateDM.py",
         "scripts/TB2J_merge.py",
         "scripts/TB2J_magnon.py",
         "scripts/TB2J_magnon_dos.py",