-
+
+- EasySpin +
- Documentation +
- Publications +
- Website +
- Academy +
- Forum +
++Fast Fourier transform over the rotation group SO(3). +
+ + ++c = fftso3(fcn,Lmax) +[LMK,vals] = fftso3(fcn,Lmax) +[c,LMK,vals] = fftso3(fcn,Lmax) ++ + +
+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. fcn should be able to handle array inputs for alpha and gamma.
+
+The Fourier expansion is in terms of Wigner functions +
+ +![[eqn]](eqn/fftso31.png)
+
+Lmax is the maxmium value of L.
+
+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]](eqn/fftso32.png)
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.
+
+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)10*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 ++ + + +
+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. +
+ + +Opt structure has changed from Opt.Output to Opt.separate and has more options: 'components', 'transitions', 'sites' and 'orientations'. There is extended support for returning subspectra across all EasySpin simulation functions.Major bug fixes
diff --git a/easyspin/chili.m b/easyspin/chili.m index 820c1151..4869bca8 100644 --- a/easyspin/chili.m +++ b/easyspin/chili.m @@ -1007,7 +1007,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/fftso3.m b/easyspin/fftso3.m new file mode 100644 index 00000000..683808c4 --- /dev/null +++ b/easyspin/fftso3.m @@ -0,0 +1,185 @@ +% fftso3 Fast Fourier transform over the rotation group SO(3) +% +% 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 +% c^L_MK of its tuncated expansion in terms of Wigner D-functions, D^L_MK: +% +% fcn(alpha,beta,gamma) = sum_(L,M,K) c^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; 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 +% c{L+1} contains the (2L+1)x(2L+1)-dimensional matrix of expansion +% coefficients c^L_MK, with M and K ordered in ascending order: M,K = +% -L,-L+1,..,L. +% LMK L, M, and K values of the non-zero coefficients +% vals values of the non-zero coefficients + +% References: +% 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 +% 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,betamode) + +if nargin==0 + help(mfilename); +end + +if nargin~=2 && nargin~=3 + error('Two inputs are required (fcn, Lmax).'); +end + +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,gamma) +alpha = 2*pi*(0:2*B-2)/(2*B-1); % radians +gamma = 2*pi*(0:2*B-2)/(2*B-1); % 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 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 (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 + +% Calculate Fourier coefficients +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 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^2)*c_; +end + +% Remove numerical noise +maxcoeff = cellfun(@(x)max(abs(x),[],'all'),c); +threshold = 1e-13*max(maxcoeff); +for L = 0:Lmax + idx = abs(c{L+1})