StollLab / EasySpin

MATLAB toolbox for Electron Paramagnetic Resonance (EPR) spectroscopy
http://easyspin.org
MIT License
48 stars 25 forks source link

Please document the parameter `Opt.FuzzLevel` #279

Closed micb25 closed 1 year ago

micb25 commented 1 year ago

Description According to this post, the yet undocumented Opt.FuzzLevel parameter was introduced in EasySpin 6.0.0-dev.22 to better handle situations with degenerate energy levels:

Starting with 6.0.0-dev.22, the handling of this is now improved. EasySpin now adds a tiny bit of numerical noise to non-zero Hamiltonian matrix elements if it suspects degeneracies.You can control the noise level with the undocumented option Opt.FuzzLevel. Its default is set to 1e-10.

I am fully aware that in most cases, i.e., in a normal EPR simulation, the impact of Opt.FuzzLevel is totally negligible. However, I think that knowing about the existence and behavior of this parameter is crucial for EasySpin users. Currently, I am working on single-crystal EPR difference spectra of trinuclear systems with tiny g shifts. In some of these simulations, I found a difference in the spectra depending on either Opt.FuzzLevel = 1e-10 (default value) or Opt.FuzzLevel = 0 (manually switched off). In general, it might get harder to spot errors and reproduce simulations by adding randomness.

For me, it makes more sense to fix degenerate energy levels by adjusting simulation parameters such as J or g, because this leads to more consistent results. Thus, in cases of degenerate energy levels it might be helpful to inform the user about the issue and let the user decide how to fix it: either by the user himself/herself (adjusting simulation parameters) or by the program (fuzzing). In my opinion, the fuzzing should be switched off by default (Opt.FuzzLevel = 0).

Suggestions

Documentation

EasySpin

stestoll commented 1 year ago

Good points and good suggestions, thanks!

Can you please post a minimal script that illustrates a case where the difference between simulations with Opt.FuzzLevel=1e-10 and Opt.FuzzLevel=0 are non-neglibigle? We can then use this as a test case to improve upon the current implementation.

micb25 commented 1 year ago

Thanks! I will try to post a working minimal script by next week.

We had an internal discussion yesterday, in which we discussed another solution to the issue. It would be great if the fuzzing could be somehow controlled and thus be more consistent, i.e., by applying the same randomness multiple times. As a programmer, I personally would solve this issue by a seed value which generates exactly the same randomness for a given seed. In this way, the user can either provide a custom seed value (which always generates the same fuzzing) or not. For the latter case, EasySpin could just generate a random seed value for the fuzzing every time.

For my simulations, I calculate four single-crystal EPR difference spectra and add them all together. Unfortunately, with a small random variation in each of the eight Hamiltonians, it can happen that the sum spectra are slightly different. Hence, it would be beneficial, if the same fuzzing could be applied every single time.

stestoll commented 1 year ago

You can seed the raandom number generator yourself before the simulation, using rng(). This will give reproducible results. By design, EasySpin does not set seeds internally - exactly for this reason.

micb25 commented 1 year ago

First of all, thank you very much for the tip with the seed for rng() which solved another issue for us!

Can you please post a minimal script that illustrates a case where the difference between simulations with Opt.FuzzLevel=1e-10 and Opt.FuzzLevel=0 are non-neglibigle? We can then use this as a test case to improve upon the current implementation.

I attached a minimal script for a trinuclear system with a degenerate ground state (see below). The simulations in the first column use only different values for Opt.FuzzLevel to lift the degenerate states. In the second column, the degenerate states are lifted by a tiny variation in one of the three J values and the Opt.FuzzLevel is also varied.

For the default value Opt.FuzzLevel = 1e-10 (first row) there are small changes in the simulated intensity for both cases. For our simulations, however, we add several of these difference spectra together which in sum basically lead to different results depending on Opt.FuzzLevel. As you can imagine, by not using a fixed seed for the random number generator the issue gets much worse in terms of reproducibility.

Screenshot Screenshot_20230125_091243

Script

clc, clear, clf

Exp.CrystalOrientation = [ 0 0 0];
Exp.Range = [240 380];
Exp.nPoints = 10000;
Exp.mwFreq = 9.4;
Exp.Temperature = 20;

FuzzLevels = [1e-10 1e-12 1e-14 1e-16 0];
rng(65725345); % "fix" randomness

lwpp = [10 10];
gvec1 = [ 4.1105 0 0];
gvec2 = [ 2.0161 0 0];
gvec3 = [ 6.2049 0 0];
dvec = [0 0 0];
Jval = 8.4242e+06;
gvala = [2.016399328628010   2.048602157921936   2.201603656381605];
gvalb = [2.016400671371990   2.048597842078065   2.201604968993395];

Sys.S = [1/2 1/2 1/2];
Sys.gFrame = [ gvec1 ; gvec2; gvec3 ];
Sys.dvec = [dvec; dvec; dvec];
Sys.lwpp = lwpp;

tiledlayout(numel(FuzzLevels), 2);
for j=0:1
    Sys.J = [ Jval Jval Jval+j ];
    for f=1:numel(FuzzLevels)
        Opt.FuzzLevel = FuzzLevels(f);
        nexttile((f-1)*2+j+1); hold on
        title(sprintf('FuzzLevel = %.0e, Jdiff = %i', Opt.FuzzLevel, j));
        for N=1:30
            Sys.g = [gvala; gvala; gvala]; [x,y1] = pepper(Sys, Exp, Opt);
            Sys.g = [gvalb; gvalb; gvalb]; [x,y2] = pepper(Sys, Exp, Opt);
            plot(x, y2-y1);
        end
    end
end
stestoll commented 1 year ago

Thanks for the example script.

This system with Exp.CrystalOrientation=[0,0,0] consists of three equivalent spins, i.e. the Hamiltonian doesn't change under any permutation of the three spins. If the orientation is tilted away from [0,0,0], then the spins are no longer equivalent (due to different g tensor orientations), and the problem disappears.

With the [0,0,0] orientation, the energy levels fall into two separate groups of four, with different transformation properties under permutation of the thee spins: under the permutation group S3, a group of 4 has A1 symmetry and the other group of 4 as E symmetry. The problem is due to the degeneracies within the E manifold. Blocking out the Hamiltonian into A1 and E parts (which would be easily doable) won't directly solve the problem of the degeneracies.

We will have to dig in a bit more to figure out how to solve this problem in a rigorous and general way.

As for now, your strategy of breaking the permutation symmetry at the parameter level seems to be a successful workaround.

stestoll commented 1 year ago

Here's a minimal script for this system with degeneracies:

clc, clear

g1Frame = [0 0 0];
g2Frame = [2*pi/3 0 0];
g3Frame = [4*pi/3 0 0];
J = 1e5;
g = [2.02, 2.05, 2.20];

Sys.S = [1/2 1/2 1/2];
Sys.g = [g; g; g];
Sys.gFrame = [g1Frame; g2Frame; g3Frame];
Sys.J = [J J J];

B0 = [0;0;350];
H = ham(Sys,B0);
energies = eig(H)/1e3
stestoll commented 1 year ago

We'll fix documentation of Opt.FuzzLevel with this issue.

For the symmetry determination, I opened a new issue #280.