From ebe2e572f490e2e7f3b83ac0684adf5e2c5cabcb Mon Sep 17 00:00:00 2001 From: Stefan Stoll Date: Mon, 1 Mar 2021 20:55:32 -0800 Subject: [PATCH 1/8] add FFT over SO(3) --- docsrc/fftso3.html | 126 ++++++++++++++++++++++++++++++++++++ docsrc/funcsalphabet.html | 1 + docsrc/funcscategory.html | 1 + docsrc/releases.html | 1 + easyspin/fftso3.m | 109 +++++++++++++++++++++++++++++++ tests/fftso3_singlewigner.m | 26 ++++++++ tests/fftso3_value.m | 27 ++++++++ 7 files changed, 291 insertions(+) create mode 100644 docsrc/fftso3.html create mode 100644 easyspin/fftso3.m create mode 100644 tests/fftso3_singlewigner.m create mode 100644 tests/fftso3_value.m diff --git a/docsrc/fftso3.html b/docsrc/fftso3.html new file mode 100644 index 00000000..dfa1451a --- /dev/null +++ b/docsrc/fftso3.html @@ -0,0 +1,126 @@ + + + + + + + + + + fftso3 + + + + +
+ +
+ +
+ +
fftso3
+ +

+Fast Fourier transform over the rotation group SO(3). +

+ + +
Syntax
+ +
+c = fftso3(fcn,Lmax)
+c = fftso3(fcn,Lmax,nrm)
+
+ + +
Description
+ +

+fftso3 takes the function specified by the function handle fcn and calculates its Fourier coefficients over the space of rotations. +

+ +

+fcn is a function handle for a function fcn(alpha,beta,gamma defined over the three Euler angles alpha, beta, and gamma, in radians. +

+ +

+The Fourier expansion is in terms of Wigner functions +

+ +
+ [eqn] +
+ +

+Lmax is the maxmium value of L. +

+ +

+nrm is optional and can be false (default) or true. If true, the normalized Wigner functions are used in the above expansion instead of the standard un-normalized Wigner functions. +

+ + +

+c is a cell array with Lmax elements. c{L+1} contains the (2L+1)x(2L+1) matrix of coefficients for a given L, with both M and K in ascending order. In other words, [eqn] is given by c{L+1}(L+1+M,L+1+K). +

+ + +
Examples
+ +

+Here is a simple example that takes a Wigner function and returns the Fourier transform. Since the function is a single Wigner function, only a single coefficient is non-zero. +

+ +
+f = @(a,b,c)132*wignerd([1 1 -1],a,b,c);
+c = fftso3(f,2)
+c{2}
+
+
+ans =
+    0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
+    0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
+   10.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
+
+ + + +
Algorithm
+ +

+See Kostelec & Rockmore, FFTs on the Rotation Group, J. Fourier Anal. Appl., 2008, 14, 145-179. +

+ + +
See also
+ +

+erot, +eulang, +wignerd +

+ +
+
+ + + + + diff --git a/docsrc/funcsalphabet.html b/docsrc/funcsalphabet.html index b46e830a..2b5d3117 100644 --- a/docsrc/funcsalphabet.html +++ b/docsrc/funcsalphabet.html @@ -82,6 +82,7 @@ fastmotionFast-motion regime line width parameters from rotational correlation time. fdaxisFrequency domain axis for FFT fieldmodField modulation of EPR absorption spectra +fftso3Fast Fourier transform on the rotation group SO(3) garlicLiquid-state cw EPR spectra gaussianGaussian line shape gfreeg value of the free electron diff --git a/docsrc/funcscategory.html b/docsrc/funcscategory.html index 1976120b..3b401eb6 100644 --- a/docsrc/funcscategory.html +++ b/docsrc/funcscategory.html @@ -197,6 +197,7 @@ Angular momentum clebschgordanClebsch-Gordan coefficients +fftso3Fast Fourier transform on the rotation group SO(3) plegendreLegendre polynomials and Associated Legendre polynomials spherharmSpherical harmonics wigner3jWigner 3-j symbols diff --git a/docsrc/releases.html b/docsrc/releases.html index 35c78692..742ae657 100644 --- a/docsrc/releases.html +++ b/docsrc/releases.html @@ -50,6 +50,7 @@

Changes from release to release

  • The field Sys.aF in spin systems is no longer supported. Use Sys.B4 and Sys.B4Frame instead.
  • Interconvert between spherical and cartesian tensors with the new functions tensor_cart2sph and tensor_sph2cart.
  • rotview provides a graphical interface to visualize and explore Euler angles and their associated rotations.
    • +
  • fftso3 performs a Fast Fourier Transform over the space of rotations, SO(3).

    • Major bug fixes

    diff --git a/easyspin/fftso3.m b/easyspin/fftso3.m new file mode 100644 index 00000000..51ea851d --- /dev/null +++ b/easyspin/fftso3.m @@ -0,0 +1,109 @@ +% fftso3 Fast Fourier transform over the rotation group SO(3) +% +% c = fftso3(fcn,Lmax) +% c = fftso3(fcn,Lmax,nrm) +% +% Takes the function fcn(aplha,beta,gamma) defined over the Euler angles +% 0<=alpha<2*pi, 0<=beta<=pi, 0<=gamma<2*pi and calculates the coefficients c^L_MK +% of its truncated Fourier expansion in terms of Wigner D-functions D^L_MK: +% +% fcn(alpha,beta,gamma) = sum_(L,M,K) f^L_MK * D^L_MK(alpha,beta,gamma) +% +% where L = 0..Lmax, M = -L..L, K = -L..L. +% +% Input: +% fcn function handle for function fcn(alpha,beta,gamma), with angles +% in units of radians +% Lmax maximum L for Wigner expansion +% nrm true: use normalized Wigner functions sqrt((2L+1)/8pi^2)*D^L_MK +% false: use standard un-normalized Wigner functions D^L_MK +% (default) +% +% Output: +% c (Lmax+1)-element cell array of Fourier coefficients, L = 0 .. Lmax +% c{L+1} contains the (2L+1)x(2L+1) matrix of Fourier coefficients +% f^L_MK, with M and K ordered in ascending order: M,K = +% -L,-L+1,..,L + +% References: +% 1) P. J. Kostelec, D. N. Rockmore +% FFTs on the Rotation Group +% J.Fourier.Anal.Appl. 2008, 14, 145/179 +% https://doi.org/10.1007/s00041-008-9013-5 +% 2) D.-M. Lux, C. Wülker, G. S. Chirikjian, +% Parallelization of the FFT on SO(3) +% arXiv:1808.00896 [cs.DC] +% https://arxiv.org/abs/1808.00896 + +function c = fftso3(fcn,Lmax,nrm) + +if nargin==0 + help(mfilename); +end + +if nargin<3 + nrm = false; +end + +B = Lmax+1; % bandwidth (0 <=L < B) + +% Set up sampling grid over (alpha,beta,gamma) +k = 0:2*B-1; +alpha = 2*pi*k/(2*B); % radians +beta = pi*(2*k+1)/(4*B); % radians +gamma = 2*pi*k/(2*B); % radians + +% Calculate 2D inverse FFT along alpha and gamma, one for each beta +S = cell(1,2*B); +for j = 0:2*B-1 + f = fcn(alpha.',beta(j+1),gamma); % evaluate function over (alpha,gamma) grid + S_ = (2*B)^2*ifft2(f); % remove the prefactor included by ifft2() + S_ = fftshift(S_); % reorder to M1,M2 = -B, -B+1, ..., 0, ..., B-1 + S{j+1} = S_(2:end,2:end); % drop first row and column (M1,M2 = -B) +end + +% Calculate beta weights +i = 0:B-1; % summation index +wB = zeros(1,2*B); +for j = 0:2*B-1 + wB(j+1) = 2*pi/B^2 * sin(beta(j+1)) * sum(sin((2*i+1)*beta(j+1))./(2*i+1)); +end + +% Calculate all Wigner d-function values for all L and beta, using matrix exponential +d = cell(Lmax+1,2*B); +for L = 0:Lmax + v = sqrt((1:2*L).*(2*L:-1:1))/2; % off-diagonals of Jy matrix (without 1/i) + Jy = diag(v,-1) - diag(v,1); % assemble Jy matrix (M1,M2 = -Lmax,-Lmax+1,...,Lmax) + for j = 0:2*B-1 + d{L+1,j+1} = expm(-beta(j+1)*Jy); % (i dropped since 1/i dropped in Jy) + end +end + +% Calculate Fourier coefficients +c = cell(1,Lmax+1); +for L = 0:Lmax + c_ = zeros(2*L+1); + idx = (Lmax+1)+(-L:L); + for j = 0:2*B-1 % run over all beta values + c_ = c_ + wB(j+1)*d{L+1,j+1}.*S{j+1}(idx,idx); + end + c{L+1} = (2*L+1)/(8*pi*B)*c_; +end + +% Normalize +if nrm + for L = 0:Lmax + N = sqrt((2*L+1)/(8*pi^2)); + c{L+1} = c{L+1}/N; + end +end + +% Zero coefficients that are below threshold +maxcoeff = cellfun(@(x)max(abs(x),[],'all'),c); +threshold = 1e-13*max(maxcoeff); +for L = 0:Lmax + idx = abs(c{L+1}) Date: Thu, 4 Mar 2021 16:08:56 -0800 Subject: [PATCH 2/8] fftso3 -> ffteuler; improve output interface --- docsrc/eulerangles.html | 20 ++-- docsrc/{fftso3.html => ffteuler.html} | 29 ++--- docsrc/funcsalphabet.html | 2 +- docsrc/funcscategory.html | 2 +- docsrc/hamiltonian.html | 4 +- docsrc/releases.html | 2 +- easyspin/ffteuler.m | 128 +++++++++++++++++++++ easyspin/fftso3.m | 109 ------------------ tests/ffteuler_singlewigner.m | 22 ++++ tests/{fftso3_value.m => ffteuler_value.m} | 2 +- tests/fftso3_singlewigner.m | 26 ----- 11 files changed, 181 insertions(+), 165 deletions(-) rename docsrc/{fftso3.html => ffteuler.html} (70%) create mode 100644 easyspin/ffteuler.m delete mode 100644 easyspin/fftso3.m create mode 100644 tests/ffteuler_singlewigner.m rename tests/{fftso3_value.m => ffteuler_value.m} (95%) delete mode 100644 tests/fftso3_singlewigner.m diff --git a/docsrc/eulerangles.html b/docsrc/eulerangles.html index 03ab4f02..50e2168e 100644 --- a/docsrc/eulerangles.html +++ b/docsrc/eulerangles.html @@ -135,7 +135,7 @@

    Rotations and Euler angles

    - [eqn] and [eqn] and +i.e., interchange [eqn] and [eqn] and invert all the signs.

    Equivalent sets of Euler angles

    -Of course, the addition to any angle of an arbitrary multiple of 2[eqn] +Of course, the addition to any angle of an arbitrary multiple of 2[eqn] has no effect on the rotation matrix.

    -[eqn] to both α and γ. +(or subtract) [eqn] to both α and γ.

    @@ -260,7 +260,7 @@

    Rotations and Euler angles

    -[eqn]
    Description

    -fftso3 takes the function specified by the function handle fcn and calculates its Fourier coefficients over the space of rotations. +ffteuler takes the function specified by the function handle fcn and calculates its Fourier coefficients over the space of rotations.

    -fcn is a function handle for a function fcn(alpha,beta,gamma defined over the three Euler angles alpha, beta, and gamma, in radians. +fcn is a function handle for a function fcn(alpha,beta,gamma) defined over the three Euler angles alpha, beta, and gamma, in radians. fcn should be able to handle array inputs for alpha and gamma.

    - +

    The Fourier expansion is in terms of Wigner functions

    - [eqn] is given by c{L+1}(L+1+M,L+1+K).

    -

    -c is a cell array with Lmax elements. c{L+1} contains the (2L+1)x(2L+1) matrix of coefficients for a given L, with both M and K in ascending order. In other words, [eqn] is given by c{L+1}(L+1+M,L+1+K). +LMK is a list of (L,M,K) triplets of the coefficients with non-zero values, with one triplet per row. vals is the corresponding list of values.

    +
    Examples
    @@ -89,8 +90,8 @@ 

    -f = @(a,b,c)132*wignerd([1 1 -1],a,b,c);
-c = fftso3(f,2)
+f = @(a,b,c)10*wignerd([1 1 -1],a,b,c);
+c = ffteuler(f,2);
 c{2}
 
    @@ -105,7 +106,7 @@
 
    Algorithm
    
 
 
    
-See Kostelec & Rockmore, FFTs on the Rotation Group, J. Fourier Anal. Appl., 2008, 14, 145-179.
+See Kostelec & Rockmore, FFTs on the Rotation Group, J. Fourier Anal. Appl., 2008, 14, 145-179, https://doi.org/10.1007/s00041-008-9013-5.
 
    
 
 
diff --git a/docsrc/funcsalphabet.html b/docsrc/funcsalphabet.html
index 2b5d3117..62435e8c 100644
--- a/docsrc/funcsalphabet.html
+++ b/docsrc/funcsalphabet.html
@@ -82,7 +82,7 @@
 fastmotionFast-motion regime line width parameters from rotational correlation time.
 fdaxisFrequency domain axis for FFT
 fieldmodField modulation of EPR absorption spectra
-fftso3Fast Fourier transform on the rotation group SO(3)
+ffteulerFast Fourier transform on the rotation group SO(3)
 garlicLiquid-state cw EPR spectra
 gaussianGaussian line shape
 gfreeg value of the free electron
diff --git a/docsrc/funcscategory.html b/docsrc/funcscategory.html
index 3b401eb6..a5cbcba0 100644
--- a/docsrc/funcscategory.html
+++ b/docsrc/funcscategory.html
@@ -197,7 +197,7 @@
 
 Angular momentum
 clebschgordanClebsch-Gordan coefficients
-fftso3Fast Fourier transform on the rotation group SO(3)
+ffteulerFast Fourier transform on the rotation group SO(3)
 plegendreLegendre polynomials and Associated Legendre polynomials
 spherharmSpherical harmonics
 wigner3jWigner 3-j symbols
diff --git a/docsrc/hamiltonian.html b/docsrc/hamiltonian.html
index b48fda06..68c77e4c 100644
--- a/docsrc/hamiltonian.html
+++ b/docsrc/hamiltonian.html
@@ -667,10 +667,10 @@ 
    The Spin Hamiltonian
    
 It is usually given in barn (1 barn = 10-28 m2 = 10-24 cm2).
 
 
    
-The nuclear quadrupole tensor is related to the EFG tensor [eqn] via
+The nuclear quadrupole tensor is related to the EFG tensor [eqn] via
 
 
    
-[eqn]
 
    Description
    
 
 
    
-ffteuler takes the function specified by the function handle fcn and calculates its Fourier coefficients over the space of rotations.
+fftso3 takes the function specified by the function handle fcn and calculates its Fourier coefficients over the space of rotations.
 
    
 
 
    
@@ -57,7 +57,7 @@
 
    
 
 
    
-    [eqn] is given by c{L+1}(L+1+M,L+1+K).
+c is a cell array with Lmax elements. c{L+1} contains the (2L+1)x(2L+1) matrix of coefficients for a given L, with both M and K in ascending order. In other words, [eqn] is given by c{L+1}(L+1+M,L+1+K).
 
    
 
 
    
@@ -91,7 +91,7 @@
 
 
     f = @(a,b,c)10*wignerd([1 1 -1],a,b,c);
-c = ffteuler(f,2);
+c = fftso3(f,2);
 c{2}
 
    
 
    diff --git a/docsrc/funcsalphabet.html b/docsrc/funcsalphabet.html
index 62435e8c..2b5d3117 100644
--- a/docsrc/funcsalphabet.html
+++ b/docsrc/funcsalphabet.html
@@ -82,7 +82,7 @@
 fastmotionFast-motion regime line width parameters from rotational correlation time.
 fdaxisFrequency domain axis for FFT
 fieldmodField modulation of EPR absorption spectra
-ffteulerFast Fourier transform on the rotation group SO(3)
+fftso3Fast Fourier transform on the rotation group SO(3)
 garlicLiquid-state cw EPR spectra
 gaussianGaussian line shape
 gfreeg value of the free electron
diff --git a/docsrc/funcscategory.html b/docsrc/funcscategory.html
index a5cbcba0..3b401eb6 100644
--- a/docsrc/funcscategory.html
+++ b/docsrc/funcscategory.html
@@ -197,7 +197,7 @@
 
 Angular momentum
 clebschgordanClebsch-Gordan coefficients
-ffteulerFast Fourier transform on the rotation group SO(3)
+fftso3Fast Fourier transform on the rotation group SO(3)
 plegendreLegendre polynomials and Associated Legendre polynomials
 spherharmSpherical harmonics
 wigner3jWigner 3-j symbols
diff --git a/docsrc/releases.html b/docsrc/releases.html
index e75aa829..742ae657 100644
--- a/docsrc/releases.html
+++ b/docsrc/releases.html
@@ -50,7 +50,7 @@ 
    Changes from release to release
    
 
  • The field Sys.aF in spin systems is no longer supported. Use Sys.B4 and Sys.B4Frame instead.
    • 
 
  • Interconvert between spherical and cartesian tensors with the new functions tensor_cart2sph and tensor_sph2cart.
    • 
 
  • rotview provides a graphical interface to visualize and explore Euler angles and their associated rotations.
    • 
-
  • ffteuler performs a Fast Fourier Transform over the space of rotations, SO(3).
    • 
+
  • fftso3 performs a Fast Fourier Transform over the space of rotations, SO(3).
    • 
 
 
 
    Major bug fixes
    
diff --git a/easyspin/ffteuler.m b/easyspin/fftso3.m
similarity index 91%
rename from easyspin/ffteuler.m
rename to easyspin/fftso3.m
index 462cfbf2..6fd9006b 100644
--- a/easyspin/ffteuler.m
+++ b/easyspin/fftso3.m
@@ -1,8 +1,8 @@
-% ffteuler   Fast Fourier transform over the rotation group SO(3)
+% fftso3   Fast Fourier transform over the rotation group SO(3)
 %
-%   c = ffteuler(fcn,Lmax)
-%   [LMK,vals] = ffteuler(fcn,Lmax)
-%   [c,LMK,vals] = ffteuler(fcn,Lmax)
+%   c = fftso3(fcn,Lmax)
+%   [LMK,vals] = fftso3(fcn,Lmax)
+%   [c,LMK,vals] = fftso3(fcn,Lmax)
 %
 % Takes the function fcn(alpha,beta,gamma) defined over the Euler angles
 % 0<=alpha<2*pi, 0<=beta<=pi, 0<=gamma<2*pi and calculates the coefficients
@@ -36,7 +36,7 @@
 %    arXiv:1808.00896 [cs.DC]
 %    https://arxiv.org/abs/1808.00896
 
-function varargout = ffteuler(fcn,Lmax)
+function varargout = fftso3(fcn,Lmax)
 
 if nargin==0
   help(mfilename);
@@ -49,11 +49,13 @@
 B = Lmax+1; % bandwidth (0 <=L < B)
 
 % Set up sampling grid over (alpha,beta,gamma)
+k = 0:2*B-2;
+alpha = 2*pi*k/(2*B-1); % radians
+gamma = 2*pi*k/(2*B-1); % radians
+[alpha,gamma] = ndgrid(alpha,gamma); % 2D array for function evaluations
+
 k = 0:2*B-1;
-alpha = 2*pi*k/(2*B); % radians
 beta = pi*(2*k+1)/(4*B); % radians
-gamma = 2*pi*k/(2*B); % radians
-[alpha,gamma] = ndgrid(alpha,gamma); % 2D array for function evaluations
 nbeta = numel(beta);
 
 % Calculate 2D inverse FFT along alpha and gamma, one for each beta
@@ -89,7 +91,7 @@
 c = cell(Lmax+1,1);
 for L = 0:Lmax
   c_ = zeros(2*L+1);
-  idx = (Lmax+2)+(-L:L); % to access the range M,K = -L:L
+  idx = (Lmax+1)+(-L:L); % to access the range M,K = -L:L
   for j = 1:nbeta % run over all beta values
     c_ = c_ + wB(j)*d{L+1,j}.*S{j}(idx,idx);
   end
diff --git a/tests/fftso3_general.m b/tests/fftso3_general.m
new file mode 100644
index 00000000..dfc486fd
--- /dev/null
+++ b/tests/fftso3_general.m
@@ -0,0 +1,34 @@
+function ok = test()
+
+Lmax = 5;
+rng(623);
+
+% Generate a general bandlimited function from Wigner expansion
+idx = 0;
+N = (Lmax+1)*(2*Lmax+1)*(2*Lmax+3)/3;
+LMK = zeros(N,3);
+vals = zeros(N,1);
+for L = 0:Lmax
+  for K = -L:L
+    for M = -L:L
+      idx = idx + 1;
+      LMK(idx,:) = [L M K];
+      vals(idx) = rand + 1i*rand;
+    end
+  end
+end
+f = @(a,b,c)fcn(LMK,vals,a,b,c);
+
+% Calculate Fourier transform
+[LMK_,vals_] = fftso3(f,Lmax);
+
+ok = areequal(vals,vals_,1e-10,'abs') && areequal(LMK,LMK_);
+
+end
+
+function y = fcn(LMK,v,a,b,c)
+y = 0;
+for q = 1:numel(v)
+  y = y + v(q)*wignerd(LMK(q,:),a,b,c);
+end
+end
diff --git a/tests/ffteuler_largeorder.m b/tests/fftso3_largeorder.m
similarity index 85%
rename from tests/ffteuler_largeorder.m
rename to tests/fftso3_largeorder.m
index dd6ded9d..97ff93d9 100644
--- a/tests/ffteuler_largeorder.m
+++ b/tests/fftso3_largeorder.m
@@ -6,7 +6,7 @@
 
 f = @(a,b,c) A*wignerd(LMK,a,b,c);
 
-[LMK_,A_] = ffteuler(f,Lmax);
+[LMK_,A_] = fftso3(f,Lmax);
 
 ok = size(LMK_,1)==1 && all(LMK_==LMK) && areequal(A,A_,1e-10,'abs');
 
diff --git a/tests/ffteuler_singlewigner.m b/tests/fftso3_singlewigner.m
similarity index 91%
rename from tests/ffteuler_singlewigner.m
rename to tests/fftso3_singlewigner.m
index 6752cf2a..2a8d7d0d 100644
--- a/tests/ffteuler_singlewigner.m
+++ b/tests/fftso3_singlewigner.m
@@ -13,7 +13,7 @@
       
       f = @(a,b,c) A*wignerd([L M K],a,b,c);
       
-      [LMK,A_] = ffteuler(f,Lmax);
+      [LMK,A_] = fftso3(f,Lmax);
       
       idx = idx + 1;
       ok(idx) = all(LMK==[L M K]) && areequal(A_,A,1e-10,'abs');
diff --git a/tests/ffteuler_value.m b/tests/fftso3_value.m
similarity index 93%
rename from tests/ffteuler_value.m
rename to tests/fftso3_value.m
index 14d8e801..bf352365 100644
--- a/tests/ffteuler_value.m
+++ b/tests/fftso3_value.m
@@ -5,7 +5,7 @@
 
 % Calculate FFT
 Lmax = 12;
-[LMK,fLMK] = ffteuler(f,Lmax);
+[LMK,fLMK] = fftso3(f,Lmax);
 
 % Pick random orientation
 rng(252);

From 1ce1049decb98b1ab0460f1de5f7a8895695de70 Mon Sep 17 00:00:00 2001
From: Stefan Stoll 
Date: Fri, 5 Mar 2021 21:55:28 -0800
Subject: [PATCH 5/8] fftso3: return LMK in lexicographic order

---
 easyspin/fftso3.m      | 2 +-
 tests/fftso3_general.m | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/easyspin/fftso3.m b/easyspin/fftso3.m
index 6fd9006b..8c496975 100644
--- a/easyspin/fftso3.m
+++ b/easyspin/fftso3.m
@@ -113,7 +113,7 @@
   LMK = [];
   vals = [];
   for L = 0:Lmax
-    [idxM,idxK,vals_] = find(c{L+1});
+    [idxK,idxM,vals_] = find(c{L+1}.'); % transpose to assure lexicographic order of L,M,K
     nNonzero = numel(vals_);
     if nNonzero>0
       LMK = [LMK; repmat(L,nNonzero,1) idxM-L-1 idxK-L-1];
diff --git a/tests/fftso3_general.m b/tests/fftso3_general.m
index dfc486fd..4d21e4d8 100644
--- a/tests/fftso3_general.m
+++ b/tests/fftso3_general.m
@@ -9,8 +9,8 @@
 LMK = zeros(N,3);
 vals = zeros(N,1);
 for L = 0:Lmax
-  for K = -L:L
-    for M = -L:L
+  for M = -L:L
+    for K = -L:L
       idx = idx + 1;
       LMK(idx,:) = [L M K];
       vals(idx) = rand + 1i*rand;

From 637b1276d923b940cded7e20a191078482fc098a Mon Sep 17 00:00:00 2001
From: Stefan Stoll 
Date: Fri, 12 Mar 2021 00:12:43 -0800
Subject: [PATCH 6/8] fftso3: add Gauss-Legendre grid

---
 easyspin/fftso3.m      | 127 +++++++++++++++++++++++++++++------------
 tests/fftso3_general.m |   4 +-
 2 files changed, 92 insertions(+), 39 deletions(-)

diff --git a/easyspin/fftso3.m b/easyspin/fftso3.m
index 8c496975..683808c4 100644
--- a/easyspin/fftso3.m
+++ b/easyspin/fftso3.m
@@ -3,6 +3,7 @@
 %   c = fftso3(fcn,Lmax)
 %   [LMK,vals] = fftso3(fcn,Lmax)
 %   [c,LMK,vals] = fftso3(fcn,Lmax)
+%   ___ = fftso3(fcn,Lmax,mode)
 %
 % Takes the function fcn(alpha,beta,gamma) defined over the Euler angles
 % 0<=alpha<2*pi, 0<=beta<=pi, 0<=gamma<2*pi and calculates the coefficients
@@ -17,6 +18,8 @@
 %         in units of radians; it should be vectorized and accept array
 %         inputs for alpha and gamma
 %   Lmax  maximum L for expansion
+%   mode  beta grid mode, either 'GL' for Gauss-Legendre (Lmax+1 points)(default)
+%         or 'DH' for Driscoll-Healy (2*Lmax+2 points)
 %
 % Output:
 %   c     (Lmax+1)-element cell array of Fourier coefficients, L = 0 .. Lmax
@@ -27,63 +30,83 @@
 %   vals  values of the non-zero coefficients
 
 % References:
-% 1) P. J. Kostelec, D. N. Rockmore
+% 1) D. K. Maslen, D. N. Rockmore
+%    Generalized FFT's - A survey of some recent results
+%    Groups and Computation II
+%    DIMACS Series in Discrete Mathematics and Theoretical Computer
+%    Science, vol.28
+%    L. Finkelstein, W. M. Kantor (editors)
+%    American Mathematical Society, 1997
+%    https://doi.org/10.1007/s00041-008-9013-5
+% 2) P. J. Kostelec, D. N. Rockmore
 %    FFTs on the Rotation Group
 %    J.Fourier.Anal.Appl. 2008, 14, 145/179
 %    https://doi.org/10.1007/s00041-008-9013-5
-% 2) D.-M. Lux, C. Wülker, G. S. Chirikjian,
-%    Parallelization of the FFT on SO(3)
-%    arXiv:1808.00896 [cs.DC]
-%    https://arxiv.org/abs/1808.00896
+% 3) Z. Khalid, S. Durrani, R. A. Kennedy, Y. Wiaux, J. D. McEwen
+%    Gauss-Legendre Sampling on the Rotation Group
+%    IEEE Signal Process. Lett. 2016, 23, 207-211
+%    https://doi.org/10.1109/LSP.2015.2503295
 
-function varargout = fftso3(fcn,Lmax)
+function varargout = fftso3(fcn,Lmax,betamode)
 
 if nargin==0
   help(mfilename);
 end
 
-if nargin~=2
+if nargin~=2 && nargin~=3
   error('Two inputs are required (fcn, Lmax).');
 end
 
-B = Lmax+1; % bandwidth (0 <=L < B)
+if nargin<3
+  betamode = 'GL';
+end
+
+if numel(Lmax)~=1 || mod(Lmax,1) || Lmax<0
+  error('Second input (Lmax) must be a non-negative integer.');
+end
+
+B = Lmax + 1; % bandwidth (0 <=L < B)
 
-% Set up sampling grid over (alpha,beta,gamma)
-k = 0:2*B-2;
-alpha = 2*pi*k/(2*B-1); % radians
-gamma = 2*pi*k/(2*B-1); % radians
-[alpha,gamma] = ndgrid(alpha,gamma); % 2D array for function evaluations
+% Set up sampling grid over (alpha,gamma)
+alpha = 2*pi*(0:2*B-2)/(2*B-1); % radians
+gamma = 2*pi*(0:2*B-2)/(2*B-1); % radians
 
-k = 0:2*B-1;
-beta = pi*(2*k+1)/(4*B); % radians
+% Get knots and weights for integration over beta
+switch betamode
+  case 'DH' % Driscoll-Healy
+    beta = pi*((0:2*B-1)+1/2)/(2*B); % radians
+    nbeta = numel(beta);
+    % Calculate beta weights
+    q = 0:B-1; % summation index
+    w = zeros(1,nbeta);
+    for b = 1:nbeta
+      w(b) = 2/B * sin(beta(b)) * sum(sin((2*q+1)*beta(b))./(2*q+1));
+    end
+  case 'GL' % Gauss-Legendre
+    [z,w] = gauleg(-1,1,B);
+    beta = acos(z);
+end
 nbeta = numel(beta);
 
 % Calculate 2D inverse FFT along alpha and gamma, one for each beta
+[alpha,gamma] = ndgrid(alpha,gamma); % 2D grid for function evaluations
 S = cell(1,nbeta);
-for j = 1:nbeta
-  f = fcn(alpha,beta(j),gamma); % evaluate function over (alpha,gamma) grid
-  S_ = (2*B)^2*ifft2(f); % remove the prefactor included by ifft2()
-  S_ = fftshift(S_); % reorder to M1,M2 = -B, -B+1, ..., 0, ..., B-1
-  S{j} = S_;
-end
-
-% Calculate beta weights
-q = 0:B-1; % summation index
-wB = zeros(1,nbeta);
-for j = 1:nbeta
-  wB(j) = 2*pi/B^2 * sin(beta(j)) * sum(sin((2*q+1)*beta(j))./(2*q+1));
+for b = 1:nbeta
+  f = fcn(alpha,beta(b),gamma); % evaluate function over (alpha,gamma) grid
+  S_ = (2*pi)^2*ifft2(f);
+  S{b} = fftshift(S_); % reorder to M,K = -(B-1), ..., 0, ..., B-1
 end
 
 % Calculate Wigner d-function values for all L and beta, with matrix exponential
 d = cell(Lmax+1,nbeta);
 for L = 0:Lmax
   v = sqrt((1:2*L).*(2*L:-1:1))/2i; % off-diagonals of Jy matrix
-  Jy = diag(v,-1) - diag(v,1); % Jy matrix (M1,M2 = -Lmax,-Lmax+1,...,Lmax)
-  [U,LL] = eig(Jy);
-  LL = diag(LL).'; % equal to -L:L
-  for j = 1:nbeta
-    e = exp(-1i*beta(j)*LL);
-    d{L+1,j} = (U .* e) * U'; % equivalent to U*diag(e)*U'
+  Jy = diag(v,-1) - diag(v,1); % Jy matrix (M,K = -Lmax,-Lmax+1,...,+Lmax)
+  [U,mj] = eig(Jy);
+  mj = diag(mj).'; % equal to -L:L
+  for b = 1:nbeta
+    e = exp(-1i*beta(b)*mj);
+    d{L+1,b} = (U .* e) * U'; % equivalent to U*diag(e)*U', but slightly faster
   end
 end
 
@@ -91,11 +114,11 @@
 c = cell(Lmax+1,1);
 for L = 0:Lmax
   c_ = zeros(2*L+1);
-  idx = (Lmax+1)+(-L:L); % to access the range M,K = -L:L
-  for j = 1:nbeta % run over all beta values
-    c_ = c_ + wB(j)*d{L+1,j}.*S{j}(idx,idx);
+  idx = (Lmax+1)+(-L:L); % to access the range M,K = -L:L in S
+  for b = 1:nbeta % run over all beta values
+    c_ = c_ + w(b)*d{L+1,b}.*S{b}(idx,idx);
   end
-  c{L+1} = (2*L+1)/(8*pi*B)*c_;
+  c{L+1} = (2*L+1)/(8*pi^2)*c_;
 end
 
 % Remove numerical noise
@@ -130,3 +153,33 @@
 end
 
 end
+
+% Calculate Gauss-Legendre grid points and weights over inteval [x1,x2]
+% Press et al, Numerical Recipes in C, 2nd ed., p.152
+function [x,w] = gauleg(x1,x2,n)
+
+m = (n+1)/2;
+xm = (x2+x1)/2;
+xl = (x2-x1)/2;
+for i = 1:m
+  z = cos(pi*(i-0.25)/(n+0.5));
+  z1 = inf;
+  while abs(z-z1)>eps
+    p1 = 1;
+    p2 = 0;
+    for j = 1:n
+      p3 = p2;
+      p2 = p1;
+      p1 = ((2*j-1)*z*p2 - (j-1)*p3)/j;
+    end
+    pp = n*(z*p1-p2)/(z^2-1);
+    z1 = z;
+    z = z1 - p1/pp;
+  end
+  x(i) = xm - xl*z;
+  x(n+1-i)= xm + xl*z;
+  w(i) = 2*xl/((1-z^2)*pp^2);
+  w(n+1-i) = w(i);
+end
+
+end
diff --git a/tests/fftso3_general.m b/tests/fftso3_general.m
index 4d21e4d8..cac3be8a 100644
--- a/tests/fftso3_general.m
+++ b/tests/fftso3_general.m
@@ -17,7 +17,7 @@
     end
   end
 end
-f = @(a,b,c)fcn(LMK,vals,a,b,c);
+f = @(a,b,c)wignersum(LMK,vals,a,b,c);
 
 % Calculate Fourier transform
 [LMK_,vals_] = fftso3(f,Lmax);
@@ -26,7 +26,7 @@
 
 end
 
-function y = fcn(LMK,v,a,b,c)
+function y = wignersum(LMK,v,a,b,c)
 y = 0;
 for q = 1:numel(v)
   y = y + v(q)*wignerd(LMK(q,:),a,b,c);

From b9eef8b815e4f43a82d8dd7b37dd0b95a36670dd Mon Sep 17 00:00:00 2001
From: Stefan Stoll 
Date: Sat, 20 Mar 2021 22:51:00 -0700
Subject: [PATCH 7/8] chili_eqpopvec2: integrate fftso3 with chili_eqpopvec

---
 easyspin/private/chili_eqpopvec.m        |  10 +-
 easyspin/private/chili_eqpopvec2.m       | 158 +++++++++++++++++++++++
 easyspin/private/chili_startingvector2.m |  88 +++++++++++++
 3 files changed, 252 insertions(+), 4 deletions(-)
 create mode 100644 easyspin/private/chili_eqpopvec2.m
 create mode 100644 easyspin/private/chili_startingvector2.m

diff --git a/easyspin/private/chili_eqpopvec.m b/easyspin/private/chili_eqpopvec.m
index 1ccff317..71d88e92 100644
--- a/easyspin/private/chili_eqpopvec.m
+++ b/easyspin/private/chili_eqpopvec.m
@@ -38,7 +38,7 @@
   Mp = Potential.M;
   Kp = Potential.K;
   lambda = Potential.lambda;
-else
+elseif isnumeric(Potential)
   if size(Potential,2)~=4
     error('Potential must have four columns (L, M, K, lambda)');
   end
@@ -48,6 +48,8 @@
   lambda = Potential(:,4);
 end
 
+% Process potential
+%--------------------------------------------------------------------------
 % Assure that potential contains only terms with nonnegative K, and nonnegative M
 % for K=0. The other terms needed to render the potential real-valued are
 % implicitly supplemented in the subfunction U(a,b,c).
@@ -92,8 +94,7 @@
   if numel(idx0)~=1
     error('Exactly one orientational basis function with L=M=K=0 is allowed.');
   end
-  sqrtPeq = zeros(nOriBasis,1);
-  sqrtPeq(idx0) = 1;
+  sqrtPeq = sparse(idx0,1,1,nOriBasis,1);
   return
 end
 
@@ -187,10 +188,11 @@
     sqrtPeq(b) = Int;
   end
   
-  sqrtPeq = sparse(sqrtPeq);
   
 end
 
+sqrtPeq = sparse(sqrtPeq);
+
   % General orientational potential function (real-valued)
   % (assumes nonnegative K, and nonnegative M for K=0, and real lambda for
   % M=K=0; other terms are implicitly supplemented)
diff --git a/easyspin/private/chili_eqpopvec2.m b/easyspin/private/chili_eqpopvec2.m
new file mode 100644
index 00000000..737493be
--- /dev/null
+++ b/easyspin/private/chili_eqpopvec2.m
@@ -0,0 +1,158 @@
+% Calculates the vector representation of sqrt(Peq) in a LMK or LMKjK basis
+% under an orientational potential.
+%
+% Inputs:
+%   basis      orientational basis (structure with
+%                fields L, M, K, and possibly jK)
+%   Potential  orientational potential parameters (structure with fields
+%                L, M, K, and lambda; or n x 4 array)
+%   Opt        structure with options
+%     .useSelectionRules
+%     .PeqTolerances
+%
+% Outputs:
+%   sqrtPeq    vector representation of sqrt(Peq)
+%   nIntegrals number of evaluated 1D/2D/3D integrals
+
+function [sqrtPeq,nIntegrals] = chili_eqpopvec2(basis,Potential,Opt)
+
+% Parse options
+if nargin<3
+  Opt = struct;
+end
+if ~isfield(Opt,'useSelectionRules')
+  Opt.useSelectionRules = true;
+end
+if ~isfield(Opt,'PeqTolerances') || isempty(Opt.PeqTolerances)
+  Opt.PeqTolerances = [1e-6 1e-4 1e-4];
+end
+useSelectionRules = Opt.useSelectionRules;
+PeqTolerances = Opt.PeqTolerances;
+PeqIntThreshold = PeqTolerances(1);
+PeqIntAbsTol = PeqTolerances(2);
+PeqIntRelTol = PeqTolerances(3);
+
+% Parse potential
+if isstruct(Potential)
+  Lp = Potential.L;
+  Mp = Potential.M;
+  Kp = Potential.K;
+  lambda = Potential.lambda;
+elseif isnumeric(Potential)
+  if size(Potential,2)~=4
+    error('Potential must have four columns (L, M, K, lambda)');
+  end
+  Lp = Potential(:,1);
+  Mp = Potential(:,2);
+  Kp = Potential(:,3);
+  lambda = Potential(:,4);
+end
+
+% Process potential
+%--------------------------------------------------------------------------
+% Assure that potential contains only terms with nonnegative K, and nonnegative M
+% for K=0. The other terms needed to render the potential real-valued are
+% implicitly supplemented in the subfunction U(a,b,c).
+if any(Kp<0)
+  error('Only potential terms with nonnegative values of K are allowed. Terms with negative K required to render the potential real-valued are supplemented automatically.');
+end
+if any(Mp(Kp==0)<0)
+  error('For potential terms with K=0, M must be nonnegative. Terms with negative M required to render the potential real-valued are supplemented automatically.');
+end
+zeroMK = Mp==0 & Kp==0;
+if any(~isreal(lambda(zeroMK)))
+  error('Potential coefficients for M=K=0 must be real-valued.');
+end
+
+% Remove zero entries and constant offset
+rmv = lambda==0;
+rmv = rmv | (Lp==0 & Mp==0 & Kp==0);
+if any(rmv)
+  lambda(rmv) = [];
+  Lp(rmv) = [];
+  Mp(rmv) = [];
+  Kp(rmv) = [];
+end
+
+% Counter for the number of numerical 1D, 2D, and 3D integrals evaluated.
+nIntegrals = [0 0 0];
+
+jKbasis = isfield(basis,'jK') && ~isempty(basis.jK) && any(basis.jK);
+
+% Parse basis
+L = basis.L;
+M = basis.M;
+K = basis.K;
+if jKbasis
+  jK = basis.jK;
+end
+nOriBasis = numel(L);
+
+% Handle special case of no potential
+if ~any(lambda)
+  idx0 = find(L==0 & M==0 & K==0);
+  if numel(idx0)~=1
+    error('Exactly one orientational basis function with L=M=K=0 is allowed.');
+  end
+  sqrtPeq = sparse(idx0,1,1,nOriBasis,1);
+  return
+end
+
+% Abbreviations for 1D, 2D, and 3D integrals
+int_b = @(f) integral(f,0,pi,'AbsTol',PeqIntAbsTol,'RelTol',PeqIntRelTol);
+int_ab = @(f) integral2(f,0,2*pi,0,pi,'AbsTol',PeqIntAbsTol,'RelTol',PeqIntRelTol);
+int_bc = @(f) integral2(f,0,pi,0,2*pi,'AbsTol',PeqIntAbsTol,'RelTol',PeqIntRelTol);
+int_abc = @(f) integral3(f,0,2*pi,0,pi,0,2*pi,'AbsTol',PeqIntAbsTol,'RelTol',PeqIntRelTol);
+
+% Calculate partition sum Z for Peq = exp(-U)/Z
+zeroMp = all(Mp==0);
+zeroKp = all(Kp==0);
+if zeroMp && zeroKp
+  Z = (2*pi)^2 * int_b(@(b) exp(-U(0,b,0)).*sin(b));
+elseif zeroMp && ~zeroKp
+  Z = (2*pi) * int_bc(@(b,c) exp(-U(0,b,c)).*sin(b));
+elseif ~zeroMp && zeroKp
+  Z = (2*pi) * int_ab(@(a,b) exp(-U(a,b,0)).*sin(b));
+else
+  Z = int_abc(@(a,b,c) exp(-U(a,b,c)).*sin(b));
+end
+sqrtZ = sqrt(Z);
+
+% Calculate all Wigner coefficients up to Lmax
+c = fftso3(@(a,b,c)exp(-U(a,b,c)/2)/sqrtZ,max(L),'DH');
+
+% Rescale
+for iL = 1:numel(c)
+  c{iL} = c{iL}/sqrt((2*(iL-1)+1)/(8*pi^2));
+end
+
+% Copy coefficients into vector
+sqrtPeq = zeros(nOriBasis,1);
+idxM = M+L+1;
+idxL = K+L+1;
+for b = 1:nOriBasis
+  L_ = L(b);
+  sqrtPeq(b) = c{L_+1}(idxM(b),idxL(b));
+end
+
+% Remove imaginary part and convert to sparse form
+sqrtPeq = real(sqrtPeq);
+sqrtPeq(abs(sqrtPeq)0
+    if maxKp==0 && K~=0, continue; end
+    % Numerical integration if value is different from previous one
+    if all([L K]==cacheLK)
+      value = cacheval;
+    else
+      fun  = @(z)orifun(z,L,K,lambda);
+      if useOldIntegrator
+        value = quadl(fun,0,1,absTol); %#ok
+      else
+        value = integral(fun,0,1,'AbsTol',absTol);
+      end
+      value = value * sqrt((2*L+1)/prod(L-K+1:L+K));
+      if K==0
+        value = value/sqrt(2);
+      end
+      cacheval = value;
+      cacheLK = [L K];
+    end
+    stvec(b) = value;
+  end
+  
+end
+
+stvec = stvec/norm(stvec);
+stvec = sparse(stvec);
+
+return
+
+%===============================================================================
+% Integrand of the orientational integral over z in the starting vector.
+% (see Schneider 1989, p.19)
+%   based on the integral over gamma (valid for even K and Kp==2)
+%     \int_0^{2\pi} cos(K gamma) exp(B cos(2 gamma)) d gamma = 2 \pi besseli(K/2,B)
+function val = orifun(z,L,K,lambda)
+
+p = numel(lambda);
+
+% Potential terms with K == 0
+A = 0;
+if lambda(1), A = A + lambda(1)/2*plegendre(2,0,z); end
+if p>=3 && lambda(3), A = A + lambda(3)/2*plegendre(4,0,z); end
+
+% Potential terms with K == 2
+B = 0;
+if p>=2 && lambda(2), B = B + lambda(2)/(2*sqrt( 6))*plegendre(2,2,z); end
+if p==4 && lambda(4), B = B + lambda(4)/(6*sqrt(10))*plegendre(4,2,z); end
+
+val = plegendre(L,K,z).*exp(A).*besseli(K/2,B)*(2*pi);
+
+return

From 2a2a2fd283309671fb0e92a58126ac7b11aa1989 Mon Sep 17 00:00:00 2001
From: Stefan Stoll 
Date: Mon, 22 Mar 2021 15:38:55 -0700
Subject: [PATCH 8/8] small edits

---
 easyspin/chili.m                              |  6 +-
 ...hili_eqpopvec2.m => chili_eqpopvec_fast.m} | 10 +--
 easyspin/private/chili_startingvector2.m      | 88 -------------------
 easyspin/private/ksymmetrizer.m               |  4 +-
 4 files changed, 11 insertions(+), 97 deletions(-)
 rename easyspin/private/{chili_eqpopvec2.m => chili_eqpopvec_fast.m} (93%)
 delete mode 100644 easyspin/private/chili_startingvector2.m

diff --git a/easyspin/chili.m b/easyspin/chili.m
index 682c35d8..7970b3cd 100644
--- a/easyspin/chili.m
+++ b/easyspin/chili.m
@@ -1018,7 +1018,11 @@
     % Calculate sqrt(Peq) vector
     Opt_.useSelectionRules = Opt.useStartvecSelectionRules;
     Opt_.PeqTolerances = Opt.PeqTol;
-    [sqrtPeq,nInt] = chili_eqpopvec(Basis,Potential,Opt_);
+    if Opt.useLMKbasis
+      [sqrtPeq,nInt] = chili_eqpopvec_fast(Basis,Potential,Opt_);
+    else
+      [sqrtPeq,nInt] = chili_eqpopvec(Basis,Potential,Opt_);
+    end
     % Set up in full product basis, then prune
     StartVector = kron(sqrtPeq,SdetOp(:)/norm(SdetOp(:)));
     StartVector = StartVector(keep);
diff --git a/easyspin/private/chili_eqpopvec2.m b/easyspin/private/chili_eqpopvec_fast.m
similarity index 93%
rename from easyspin/private/chili_eqpopvec2.m
rename to easyspin/private/chili_eqpopvec_fast.m
index 737493be..7188b883 100644
--- a/easyspin/private/chili_eqpopvec2.m
+++ b/easyspin/private/chili_eqpopvec_fast.m
@@ -7,14 +7,14 @@
 %   Potential  orientational potential parameters (structure with fields
 %                L, M, K, and lambda; or n x 4 array)
 %   Opt        structure with options
-%     .useSelectionRules
 %     .PeqTolerances
 %
 % Outputs:
 %   sqrtPeq    vector representation of sqrt(Peq)
-%   nIntegrals number of evaluated 1D/2D/3D integrals
 
-function [sqrtPeq,nIntegrals] = chili_eqpopvec2(basis,Potential,Opt)
+function [sqrtPeq,nIntegrals] = chili_eqpopvec_fast(basis,Potential,Opt)
+
+nIntegrals = [];
 
 % Parse options
 if nargin<3
@@ -26,7 +26,6 @@
 if ~isfield(Opt,'PeqTolerances') || isempty(Opt.PeqTolerances)
   Opt.PeqTolerances = [1e-6 1e-4 1e-4];
 end
-useSelectionRules = Opt.useSelectionRules;
 PeqTolerances = Opt.PeqTolerances;
 PeqIntThreshold = PeqTolerances(1);
 PeqIntAbsTol = PeqTolerances(2);
@@ -74,9 +73,6 @@
   Kp(rmv) = [];
 end
 
-% Counter for the number of numerical 1D, 2D, and 3D integrals evaluated.
-nIntegrals = [0 0 0];
-
 jKbasis = isfield(basis,'jK') && ~isempty(basis.jK) && any(basis.jK);
 
 % Parse basis
diff --git a/easyspin/private/chili_startingvector2.m b/easyspin/private/chili_startingvector2.m
deleted file mode 100644
index 33d47a50..00000000
--- a/easyspin/private/chili_startingvector2.m
+++ /dev/null
@@ -1,88 +0,0 @@
-function stvec = chili_startingvector(Basis,lambda)
-
-if ~isempty(lambda) && numel(lambda)~=4
-  error('This functions works only for potentials with M=0, L=2,4, and K=0,2.');
-end
-
-isPotential = any(lambda);
-maxKp = 2;
-
-isNuc1 = isfield(Basis,'pI1');
-isNuc2 = isfield(Basis,'pI2');
-
-nBasis = numel(Basis.L);
-stvec = zeros(nBasis,1);
-  
-absTol = 1e-8; % for numerical integration
-useOldIntegrator = verLessThan('matlab','7.14'); % us quadl() instead of integral()
-
-cacheval = NaN;
-cacheLK = [NaN NaN];
-
-for b = 1:nBasis
-  L = Basis.L(b);
-  K = Basis.K(b);
-  
-  if mod(L,2) || mod(K,2), continue; end
-  if Basis.M(b)~=0, continue; end
-  if Basis.jK(b)~=1, continue; end
-  if abs(Basis.pS(b))~=1, continue; end
-  if isNuc1 && Basis.pI1(b)~=0, continue; end
-  if isNuc2 && Basis.pI2(b)~=0, continue; end
-  
-  if ~isPotential
-    % Zero potential: non-zero value only for L==K==0
-    if L~=0 || K~=0, continue; end
-    stvec(b) = 1;
-  else
-    % Zero for axial potential (max Kp == 0) and K>0
-    if maxKp==0 && K~=0, continue; end
-    % Numerical integration if value is different from previous one
-    if all([L K]==cacheLK)
-      value = cacheval;
-    else
-      fun  = @(z)orifun(z,L,K,lambda);
-      if useOldIntegrator
-        value = quadl(fun,0,1,absTol); %#ok
-      else
-        value = integral(fun,0,1,'AbsTol',absTol);
-      end
-      value = value * sqrt((2*L+1)/prod(L-K+1:L+K));
-      if K==0
-        value = value/sqrt(2);
-      end
-      cacheval = value;
-      cacheLK = [L K];
-    end
-    stvec(b) = value;
-  end
-  
-end
-
-stvec = stvec/norm(stvec);
-stvec = sparse(stvec);
-
-return
-
-%===============================================================================
-% Integrand of the orientational integral over z in the starting vector.
-% (see Schneider 1989, p.19)
-%   based on the integral over gamma (valid for even K and Kp==2)
-%     \int_0^{2\pi} cos(K gamma) exp(B cos(2 gamma)) d gamma = 2 \pi besseli(K/2,B)
-function val = orifun(z,L,K,lambda)
-
-p = numel(lambda);
-
-% Potential terms with K == 0
-A = 0;
-if lambda(1), A = A + lambda(1)/2*plegendre(2,0,z); end
-if p>=3 && lambda(3), A = A + lambda(3)/2*plegendre(4,0,z); end
-
-% Potential terms with K == 2
-B = 0;
-if p>=2 && lambda(2), B = B + lambda(2)/(2*sqrt( 6))*plegendre(2,2,z); end
-if p==4 && lambda(4), B = B + lambda(4)/(6*sqrt(10))*plegendre(4,2,z); end
-
-val = plegendre(L,K,z).*exp(A).*besseli(K/2,B)*(2*pi);
-
-return
diff --git a/easyspin/private/ksymmetrizer.m b/easyspin/private/ksymmetrizer.m
index 1fe910cd..dc6bd2db 100644
--- a/easyspin/private/ksymmetrizer.m
+++ b/easyspin/private/ksymmetrizer.m
@@ -11,7 +11,9 @@
 
 DebugMode = true;
 
-if isfield(basis,'jK')
+jKbasis = isfield(basis,'jK') && ~isempty(basis.jK) && any(basis.jK);
+
+if jKbasis
   error('This function expects an LMK basis, without jK.');
 end