Is it possible for 3D structures?

[AdityaRoy-1996]

I am trying to solve the spin wave neutron simulations for YCrO3 using the exchange parameters as given in this paper. But it seems like for 3d structures as well, it is showing flat bands. Can it only solve for 1d structures as of now?

A quick reply will be appreciated as It is an urgent task.

Thanks, Aditya.

Bajaj_2021_JPCM_Magnetoelastic coupling and spin contributions to entropy and thermal transport in biferroic yttrium orthochromite.pdf

[bingli621]

This code can give you dispersion along all three reciprocal space dimensions. If you’re seeing flat bands along given directions, it might mean that the interactions are not set properly. It you post the code you used to set up the cell I can take a look here.

[AdityaRoy-1996]

This code can give you dispersion along all three reciprocal space dimensions. If you’re seeing flat bands along given directions, it might mean that the interactions are not set properly. It you post the code you used to set up the cell I can take a look here.

@bingli621 This is the following code for YCrO3. Please help me solve the dispersion. The powder spectra and exchange interaction values are present in the paper.

Thanks a lot!

-- coding: utf-8 --

""" Created on Sun Nov 19 14:22:57 2023

@author: Adi """

-- coding: utf-8 --

import sys from pathlib import Path

sys.path.append(str(Path(file).parent.parent) + "/pyLiSW") from Atoms import Atoms from Bonds import Bonds from Sample import Sample from LSWT import LSWT from utils import gamma_fnc from matplotlib import pyplot as plt import numpy as np

-------------------------------------------------------------

1d AFM chain, Neel type, spin along z

-------------------------------------------------------------

if name == "main":

lattice parameters in Angstrom

# a = 3
# # determin the effective lattice parameters
lat_params  =  [5.5225, 7.5474, 5.2521]
# propagation vector
tau = (0, 0, 0)
# vector perpendicular to the plane of rotation
n = (0, 0, 1)
# temperature
te = 0
afm_chain = Sample(
    lat_params, tau, n, te,  gamma_fnc=gamma_fnc
)

# -------------------------------------------------------------
# Add atoms
# -------------------------------------------------------------
s1 = 3 / 2
aniso = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
# atom postions with effective lattice parameters
positions  =  np.array([[0.5, 0.5, 0.0],
                        [0.0, 0.5, 0.5],
                        [0.0, 0.0, 0.5],
                        [0.5, 0.0, 0.0],
                        ])
atoms  =  []
for i in range(positions.shape[0]):
    atoms.append(Atoms(t=positions[i,:], ion="Cr3", spin=s1, aniso=aniso))
afm_chain.add_atoms(atoms)

# -------------------------------------------------------------
# Add bonds
# -------------------------------------------------------------
j1 = -2.725   # meV  # AFM is positive, FM is negative
bonds = []
my_bonds  =   np.array([[1,2], [0,3], [1,0], [2,3]])
for i in range(my_bonds.shape[0]):
    bonds.append(Bonds(afm_chain, idx0 = my_bonds[i,0],
                      idx1 = my_bonds[i,1],
                      j=j1))

j2 = -0.1   # meV
my_bonds  =   np.array([[0,0], [1,1], [2,2], [3,3], [3,1], [0,2]])
for i in range(my_bonds.shape[0]):
    bonds.append(Bonds(afm_chain, idx0 = my_bonds[i,0],
                      idx1 = my_bonds[i,1],
                      j=j2))

afm_chain.add_bonds(bonds)

# -------------------------------------------------------------
# Simulate dispersion
# -------------------------------------------------------------
qe_range = (
    [0.00, 1.00, 0.01],
    [0.00, 0.01, 0.01],
    [0.00, 0.01, 0.01],
    [-20, 20, 0.01],
)
# Qlab  = {'[001]'  '\Gamma'   '[010]' '\Gamma' '[100]' };

sim_qespace = LSWT(qe_range, afm_chain)
sim_qespace.dispersion_calc()
fig, ax  =  sim_qespace.plot_disp("x")
ax.get_legend().remove()
fig.tight_layout()

# # -------------------------------------------------------------
# # Simulate intensities
# # -------------------------------------------------------------
# slice_range = (
#     [0, 3.01, 0.005],
#     [0.00, 0.01],
#     [0.00, 0.01],
#     [-6, 6.01, 0.01],
# )
# sim_qespace.inten_calc()
# sim_qespace.slice(Qlist, plot_axes=(0, 3), SIM=True, vmin=0, vmax=20)
# # -------------------------------------------------------------
# # Making cuts
# # -------------------------------------------------------------
# # cut_range = (
# #     [0.5, 0.52],
# #     [0.00, 0.01],
# #     [0.00, 0.01],
# #     [-0.5, 0.51, 0.01],
# # )
# # sim_qespace.cut(cut_range, plot_axis=3, SIM=True)

# plt.show()

#% Clear all memory
# %reset -sf

[bingli621]

Hi @AdityaRoy-1996, there are a couple of things you should pay attention to while setting up the model. 1) Magnetic ground state. You need to rotate some spins 180 deg because the group state is AF. This is reflected by the spin S' values in the terminal output. 2) Anisotropy. Introducing a very small uniaxial single-ion anisotropy is helpful to stabilize the calculation. In this case, I recommend having the easy-axis anisotropy along the c-axis, Dz=[[0,0,0],[0,0,0],[0,0,-0.01]], since the spins point along the c-axis. 3) Interactions. This code doesn't apply symmetry, so all interactions have to be introduced manually, and sometimes it can be difficult to track down all of them. The output of bond length dij in the terminal can be helpful.

I've pasted the proper setup file below with some results for the spin wave dispersion along (H00) and (HH0) directions.

[bingli621]

import sys
from pathlib import Path

sys.path.append(str(Path().parent) + "/pyLiSW")
from Atoms import Atoms
from Bonds import Bonds
from Sample import Sample
from LSWT import LSWT
from utils import gamma_fnc
from matplotlib import pyplot as plt
import numpy as np

if __name__ == "__main__":
    # lattice parameters in Angstrom
    # a = 3
    # # determin the effective lattice parameters
    lat_params = [5.5225, 7.5474, 5.2521]
    # propagation vector
    tau = (0, 0, 0)
    # vector perpendicular to the plane of rotation
    n = (1, 0, 0)
    # temperature
    te = 0
    ycro3 = Sample(lat_params, tau, n, te, gamma_fnc=gamma_fnc)

    # -------------------------------------------------------------
    # Add atoms
    # -------------------------------------------------------------
    s1 = 3 / 2
    aniso = [[0, 0, 0], [0, 0, 0], [0, 0, -0.01]]  # force moments along c
    # atom postions with effective lattice parameters
    atoms = [
        Atoms(
            t=(0.5, 0.5, 0.0),
            ion="Cr3",
            spin=s1,
            aniso=aniso,
            theta=0,
            n_Rp=(1, 0, 0),
        ),  # moments along positive c-axis
        Atoms(
            t=(0.0, 0.5, 0.5),
            ion="Cr3",
            spin=s1,
            aniso=aniso,
            theta=np.pi,
            n_Rp=(1, 0, 0),
        ),  # moments along negative c-axis, rotation needed
        Atoms(
            t=(0.0, 0.0, 0.5),
            ion="Cr3",
            spin=s1,
            aniso=aniso,
            theta=0,
            n_Rp=(1, 0, 0),
        ),  # moments along positive c-axis
        Atoms(
            t=(0.5, 0.0, 0.0),
            ion="Cr3",
            spin=s1,
            aniso=aniso,
            theta=np.pi,
            n_Rp=(1, 0, 0),
        ),  # moments along negative c-axis, rotation needed
    ]

    ycro3.add_atoms(atoms)

    # -------------------------------------------------------------
    # Add bonds
    # -------------------------------------------------------------
    j1 = 2.725  # meV  # AFM is positive, FM is negative
    # my_bonds = np.array([[1, 2], [0, 3], [1, 0], [2, 3]])
    j2 = 0.1  # meV
    # my_bonds = np.array([[0, 0], [1, 1], [2, 2], [3, 3], [3, 1], [0, 2]])
    bonds = [
        # j1_1 bonds
        Bonds(ycro3, idx0=1, idx1=2, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=1, idx1=2, r0=(0, 0, 0), r1=(0, 1, 0), j=j1),
        Bonds(ycro3, idx0=0, idx1=3, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=0, idx1=3, r0=(0, 0, 0), r1=(0, 1, 0), j=j1),
        # j1_2 bonds
        Bonds(ycro3, idx0=1, idx1=0, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=1, idx1=0, r0=(0, 0, 0), r1=(0, 0, 1), j=j1),
        Bonds(ycro3, idx0=0, idx1=1, r0=(0, 0, 0), r1=(1, 0, 0), j=j1),
        Bonds(ycro3, idx0=0, idx1=1, r0=(0, 0, 1), r1=(1, 0, 0), j=j1),
        Bonds(ycro3, idx0=2, idx1=3, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=2, idx1=3, r0=(0, 0, 0), r1=(0, 0, 1), j=j1),
        Bonds(ycro3, idx0=2, idx1=3, r0=(1, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=2, idx1=3, r0=(1, 0, 0), r1=(0, 0, 1), j=j1),
        # j2_1 bonds
        Bonds(ycro3, idx0=0, idx1=0, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
        Bonds(ycro3, idx0=0, idx1=0, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=1, idx1=1, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
        Bonds(ycro3, idx0=1, idx1=1, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=2, idx1=2, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
        Bonds(ycro3, idx0=2, idx1=2, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=3, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
        Bonds(ycro3, idx0=3, idx1=3, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        # j2_2 bonds
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 0), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 0), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 1), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 1), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 1), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 1), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(0, 1, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(0, 1, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(1, 1, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(1, 1, 0), j=j2),
    ]
    ycro3.add_bonds(bonds)

    # -------------------------------------------------------------
    # Simulate dispersion along (H00)
    # -------------------------------------------------------------
    proj = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    axes = ["(H,0,0)", "(0,K,0)", "(0,0,L)"]
    qe_range = (
        [0.00, 2.01, 0.01],
        [0.00, 0.01, 0.01],
        [0.00, 0.01, 0.01],
        [-30, 30, 0.01],
    )
    sim_qespace = LSWT(qe_range, ycro3, proj_axes=proj, axes=axes)
    sim_qespace.dispersion_calc()
    sim_qespace.plot_disp("x")

    slice_range = (
        [0, 2.01, 0.01],
        [0.00, 0.01],
        [0.00, 0.01],
        [-30, 30.01, 0.01],
    )
    sim_qespace.inten_calc()
    sim_qespace.slice(slice_range, plot_axes=(0, 3), vmin=0, vmax=5)

    # -------------------------------------------------------------
    # Simulate dispersion along (HH0)
    # -------------------------------------------------------------
    proj = [[1, 1, 0], [-1, 1, 0], [0, 0, 1]]
    axes = ["(H,H,0)", "(-K,K,0)", "(0,0,L)"]
    qe_range_2 = (
        [0.00, 2.01, 0.01],
        [0.00, 0.01, 0.01],
        [0.00, 0.01, 0.01],
        [-30, 30, 0.01],
    )
    # Qlab  = {'[001]'  '\Gamma'   '[010]' '\Gamma' '[100]' };

    sim_qespace_2 = LSWT(qe_range_2, ycro3, proj_axes=proj, axes=axes)
    sim_qespace_2.dispersion_calc()
    sim_qespace_2.plot_disp("x")

    slice_range_2 = (
        [0, 2.01, 0.01],
        [0.00, 0.01],
        [0.00, 0.01],
        [-30, 30.01, 0.01],
    )
    sim_qespace_2.inten_calc()
    sim_qespace_2.slice(slice_range_2, plot_axes=(0, 3), vmin=0, vmax=5)

    plt.show()

    # % Clear all memory
    # %reset -sf

[bingli621]

[AdityaRoy-1996]

import sys
from pathlib import Path

sys.path.append(str(Path().parent) + "/pyLiSW")
from Atoms import Atoms
from Bonds import Bonds
from Sample import Sample
from LSWT import LSWT
from utils import gamma_fnc
from matplotlib import pyplot as plt
import numpy as np

if __name__ == "__main__":
    # lattice parameters in Angstrom
    # a = 3
    # # determin the effective lattice parameters
    lat_params = [5.5225, 7.5474, 5.2521]
    # propagation vector
    tau = (0, 0, 0)
    # vector perpendicular to the plane of rotation
    n = (1, 0, 0)
    # temperature
    te = 0
    ycro3 = Sample(lat_params, tau, n, te, gamma_fnc=gamma_fnc)

    # -------------------------------------------------------------
    # Add atoms
    # -------------------------------------------------------------
    s1 = 3 / 2
    aniso = [[0, 0, 0], [0, 0, 0], [0, 0, -0.01]]  # force moments along c
    # atom postions with effective lattice parameters
    atoms = [
        Atoms(
            t=(0.5, 0.5, 0.0),
            ion="Cr3",
            spin=s1,
            aniso=aniso,
            theta=0,
            n_Rp=(1, 0, 0),
        ),  # moments along positive c-axis
        Atoms(
            t=(0.0, 0.5, 0.5),
            ion="Cr3",
            spin=s1,
            aniso=aniso,
            theta=np.pi,
            n_Rp=(1, 0, 0),
        ),  # moments along negative c-axis, rotation needed
        Atoms(
            t=(0.0, 0.0, 0.5),
            ion="Cr3",
            spin=s1,
            aniso=aniso,
            theta=0,
            n_Rp=(1, 0, 0),
        ),  # moments along positive c-axis
        Atoms(
            t=(0.5, 0.0, 0.0),
            ion="Cr3",
            spin=s1,
            aniso=aniso,
            theta=np.pi,
            n_Rp=(1, 0, 0),
        ),  # moments along negative c-axis, rotation needed
    ]

    ycro3.add_atoms(atoms)

    # -------------------------------------------------------------
    # Add bonds
    # -------------------------------------------------------------
    j1 = 2.725  # meV  # AFM is positive, FM is negative
    # my_bonds = np.array([[1, 2], [0, 3], [1, 0], [2, 3]])
    j2 = 0.1  # meV
    # my_bonds = np.array([[0, 0], [1, 1], [2, 2], [3, 3], [3, 1], [0, 2]])
    bonds = [
        # j1_1 bonds
        Bonds(ycro3, idx0=1, idx1=2, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=1, idx1=2, r0=(0, 0, 0), r1=(0, 1, 0), j=j1),
        Bonds(ycro3, idx0=0, idx1=3, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=0, idx1=3, r0=(0, 0, 0), r1=(0, 1, 0), j=j1),
        # j1_2 bonds
        Bonds(ycro3, idx0=1, idx1=0, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=1, idx1=0, r0=(0, 0, 0), r1=(0, 0, 1), j=j1),
        Bonds(ycro3, idx0=0, idx1=1, r0=(0, 0, 0), r1=(1, 0, 0), j=j1),
        Bonds(ycro3, idx0=0, idx1=1, r0=(0, 0, 1), r1=(1, 0, 0), j=j1),
        Bonds(ycro3, idx0=2, idx1=3, r0=(0, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=2, idx1=3, r0=(0, 0, 0), r1=(0, 0, 1), j=j1),
        Bonds(ycro3, idx0=2, idx1=3, r0=(1, 0, 0), r1=(0, 0, 0), j=j1),
        Bonds(ycro3, idx0=2, idx1=3, r0=(1, 0, 0), r1=(0, 0, 1), j=j1),
        # j2_1 bonds
        Bonds(ycro3, idx0=0, idx1=0, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
        Bonds(ycro3, idx0=0, idx1=0, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=1, idx1=1, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
        Bonds(ycro3, idx0=1, idx1=1, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=2, idx1=2, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
        Bonds(ycro3, idx0=2, idx1=2, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=3, r0=(0, 0, 0), r1=(0, 0, 1), j=j2),
        Bonds(ycro3, idx0=3, idx1=3, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        # j2_2 bonds
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 0), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 0), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 1), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 1), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 0, 1), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=3, idx1=1, r0=(0, 1, 1), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(0, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(1, 0, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(0, 1, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(0, 1, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 0), r1=(1, 1, 0), j=j2),
        Bonds(ycro3, idx0=0, idx1=2, r0=(0, 0, 1), r1=(1, 1, 0), j=j2),
    ]
    ycro3.add_bonds(bonds)

    # -------------------------------------------------------------
    # Simulate dispersion along (H00)
    # -------------------------------------------------------------
    proj = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    axes = ["(H,0,0)", "(0,K,0)", "(0,0,L)"]
    qe_range = (
        [0.00, 2.01, 0.01],
        [0.00, 0.01, 0.01],
        [0.00, 0.01, 0.01],
        [-30, 30, 0.01],
    )
    sim_qespace = LSWT(qe_range, ycro3, proj_axes=proj, axes=axes)
    sim_qespace.dispersion_calc()
    sim_qespace.plot_disp("x")

    slice_range = (
        [0, 2.01, 0.01],
        [0.00, 0.01],
        [0.00, 0.01],
        [-30, 30.01, 0.01],
    )
    sim_qespace.inten_calc()
    sim_qespace.slice(slice_range, plot_axes=(0, 3), vmin=0, vmax=5)

    # -------------------------------------------------------------
    # Simulate dispersion along (HH0)
    # -------------------------------------------------------------
    proj = [[1, 1, 0], [-1, 1, 0], [0, 0, 1]]
    axes = ["(H,H,0)", "(-K,K,0)", "(0,0,L)"]
    qe_range_2 = (
        [0.00, 2.01, 0.01],
        [0.00, 0.01, 0.01],
        [0.00, 0.01, 0.01],
        [-30, 30, 0.01],
    )
    # Qlab  = {'[001]'  '\Gamma'   '[010]' '\Gamma' '[100]' };

    sim_qespace_2 = LSWT(qe_range_2, ycro3, proj_axes=proj, axes=axes)
    sim_qespace_2.dispersion_calc()
    sim_qespace_2.plot_disp("x")

    slice_range_2 = (
        [0, 2.01, 0.01],
        [0.00, 0.01],
        [0.00, 0.01],
        [-30, 30.01, 0.01],
    )
    sim_qespace_2.inten_calc()
    sim_qespace_2.slice(slice_range_2, plot_axes=(0, 3), vmin=0, vmax=5)

    plt.show()

    # % Clear all memory
    # %reset -sf

@bingli621 Thanks a lot to help me out with the code.

Just an extension to my previous question, is there any chance you extend your code to implement the symmetry (searching equivalent bonds will be a lot more easier, you can see the spinW code 'gencoupling' library for this) and finding powder spectra by direction averaging the single crystal spectra?

Again thanks a lot for your time.

[bingli621]

To have the symmetry analysis implemented would be great. The full-on powder-averaging should be straightforward but time-consuming. I was thinking about implementing Monte-Carlo sampling to reduce the time. But unfortunately, I don't have any specific plans for major enhancement at this time. If you're interested in contributing to this code and your analysis requires more features, I'd be more than happy to work with you.

[AdityaRoy-1996]

To have the symmetry analysis implemented would be great. The full-on powder-averaging should be straightforward but time-consuming. I was thinking about implementing Monte-Carlo sampling to reduce the time. But unfortunately, I don't have any specific plans for major enhancement at this time. If you're interested in contributing to this code and your analysis requires more features, I'd be more than happy to work with you.

@bingli621 My work is extensively in solving the powder spectra, since we work on powders. I was previously using SpinW, but its too memory heavy and crashes my laptop. Your code as of now is working like charm. I am going through code to find the powder dispersion. In case i succeed, I will definitely send you the modified code. Also i will try the symmetry analysis later.

Thanks a lot!

[bingli621]

Let me see what I can do about the powder averaging part. It might be easier than I thought.