Inconsistency of PCS fields calculated from different susceptibility tensors

[patroklus97]

Dear SpinDynamics community,

I was investigating a paramagnetic Eu(III) coordination complex (C3 symmetry) and wanted to calculate the PCS field and values using Spinach. To provide a appropriate spin density, g- and susceptibility tensor, I did a CASSCF calculation with Orca (details below) using the 7 4f orbitals as active space. Since S>1, ZFS occurs and Orca outputs an updated geometry in the D-Tensor axis frame (which should save me the later adaption of the tensors via the Euler angles). I redid the CAS calculation using this geometry.

CAS setup:

! DKH DKH-DEF2-TZVPP RIJCOSX AUTOAUX SlowConv TIGHTSCF KDIIS MoRead

%moinp "Eu_pre_CAS.gbw"

%pal
nprocs 32
end

%maxcore 6000

%basis
NewGTO Eu "SARC-DKH-TZVPP" end
end

%rel
picturechange true
FiniteNuc true
end

%casscf
nel 6
norb 7
mult 7,5,3,1
nroots 7,140,588,490
bweight = 0.00571429,0.114286,0.48,0.4
...


The obtained g-factors have the following values: 0.809722 0.953096 0.977340

The susceptibility tensor in cm³*K/mol at 273 K yielded:
chi_eu =
[ 1.409858 -0.004367 0.020120
-0.004367 1.414043 0.001456
0.020120 0.001456 1.796160 ]

It was then converted to A³ via the following line:
chi_eu=(1.660539067*chi_eu)/273;

PCS_chi_orca.png (220.78 KiB) Viewed 127144 times

Using these values, Spinach calculated a PCS field, which seems to be reasonable (oriented along the C3 axis) and agrees with literature cases of comparable systems. Experimental results from Bleany Analysis are relatively close to the respective PCS values. The only thing that is off is the discrepancy of the theoretical and experimental axial susceptibilities. mat2axrh (using chi_ax = zz-(xx+yy)/2) yields 2.35e-33 m³. The experimental value is -1.80E-32 m³ (agreeing with literature values). So first of all, it seems to be smaller by the order of 10, and second the sign is opposite. According to my understanding, the sign should be negative with a negative torus in the PCS field. Does someone have an idea, where this difference comes from? Am I doing somethin wrong with the ordering of the tensor values? That's actually something, that confused me even more: when I used Haeberlen notation (zz- iso >= xx - iso >= yy - iso) instead of Mehring notation (xx <= yy <= zz) for a series of Ln(III) complexes (Nd, Sm, Eu, Yb) I always got the opposite sign of what I expected and was known from literature.

That's basically my main problem. Something that I then wanted to try was to compare those results to PCS fields that I would get from approximating chi via the built-in g2chi-function based on the g-tensor and spin. This resulted in a PCS field, which appears to be rotated by 90° and opposite sign. chi_ax is again positive (1.99e-32 m³), but at least the value itself now comes way closer to the experimental one.

PCS_chi_Spinach.png (170.75 KiB) Viewed 127144 times

I then even considered involving low-lying excited SO states. Using single_aniso, I was able to obtain a pseudospin Hamiltonian g-Tensor taking into account the first three degenerate SO states after the ground state (198.390 cm-1, 346.182 cm-1, 353.886 cm-1):
1.66609838 1.64718265 1.05881485

In this case, the orientation seems to be right, whereas the chi_ax value (9.85e-32 m³) seems to be a bit large. The value again is positive, but in this case this would be expected looking at the PCS cone (basically the previous one rotated by 90°, the torus is positive; I somehow could'nt upload more than 3 files).

In conclusion, I can't seem to reproduce the PCS cone obtained from the Orca susceptibility tensor and moreover I can't make sense of the accompanying chi_ax values.


I hope I described the encountered problems sufficiently. Any ideas or insights to mistakes I might have made would be appreciated. Thank you in advance!

I attached my Spinach code.

Kind regards

Adrian

[kuprov]

A couple of questions:

1. The solver uses periodic boundary conditions, and it's necessary to add sufficient margins to prevent the PCS field from wrapping around. Try increasing the pad_size parameter in ocparse function and see if it makes any difference.

2. Is your probability density normalised (see line 21 in dft_density.m example)? How many unpaired electrons do you have? An appropriate multiplier needs to be introduced.

Best wishes,
Ilya.

[patroklus97]

Hello Ilya, thanks for your response!

1. Unfortunately increasing the pad_size did not change anything for chi_ax (or otherwise really) (I tried 3 and 4).

2. Neither did the normalization of the density (is it really necessary though in case I exported it via orca_plot as .spindens.3d? since I couldn't find this step in the porphyrin_example_3.m and I think ocparse already includes this step).

The Eu(III) system has 6 unpaired electrons. Where would this multiplier need to be introduced other than for chi=g2chi(g,273,6/2)?

A lot of steps would of course be more convenient for me if I calculated the parameters via DFT (using eprnmr). As far as I'm concerned the automatic extraction in Spinach does not work for Orca CAS output files, right? Maybe I could try to mimic the DFT output...

I still have the feeling that I'm doing something wrong either with the provided form of chi, the conversion or the method of determining chi_ax...
Is there maybe an example-file demonstrating the extraction of a susceptibility tensor generated with Orca or Molcas using single_aniso?

Best wishes,

Adrian.

[kuprov]

I'll take a better look when I get to a computer - probably something simple like XYZ axis order, angstroms vs nm, etc.

For a lanthanide, the point susceptibility centre should be a decent approximation. Try plotting the output of point PCS calculation against the distributed one - is there any obvious numerical factor in the slope? If the two methods broadly agree, then it's indeed an incorrect chi tensor. It could be that there are conformational averaging effects present and a single geometry is not a good approximation.

[patroklus97]

Thanks, that would be great! And definitely a relief if it was something that simple. I attached the spin density file (as Eu_III_spindens.m, needs to be renamed to Eu_III.spindens.3d)

So I assume in this case with distributed you mean the output from the Kuprov equation and not the multipole model (lpcs)? In this case the plot looks like this:

PCS_point_vs_distr.png (32.31 KiB) Viewed 127119 times

But why would it be an incorrect chi tensor if the two methods agree and the point model is a decent approximation for lanthanides?

Regarding conformational averaging: I actually did some meta dynamics calculations for this complex and it doesn't seem to have a lot of strikingly different low lying conformers (apart from the inversion of helicity (picture on the right)):

conformers.png (153.97 KiB) Viewed 127119 times

[kuprov]

Alas, the perfect agreement with the point model only leaves one possibility - for some reason, there is a problem with coordinates and / or chi: either orientation, or amplitude, or both. Essentially, that agreement plot means that two completely mathematically dissimilar solvers have returned identical answers, so in all likelihood both solvers are doing it right.

If there's still no match to experiment or common sense, this can only mean that the input variables have a problem. Can you make sure that chi and XYZ are in the same reference frame, that there are no stray minus signs anywhere, that definitions of axiality and rhombicity match (they can sometimes differ by a factor of two) and not refer to similar-sounding "anisotropy and asymmetry" or "span and skew", etc.

Essentially, the agreement between point model and distributed solver gives Spinach mathematics a clean bill of health, this is not a program bug. It's something with the physics of the system or the definitions of the parameters that go in.

[patroklus97]

Ok, so I redid the setup and there actually seemed to be a problem with the reference frames. This at least solved the problem of different orientations for the Orca chi and Spinach chi. But some main issues still remain. To demonstrate it more clearly, I provided the additional example of the same complex but with Nd(III).

Nd.png (327.79 KiB) Viewed 127107 times

In the case of Nd(III) I get the same PCS field both for the chi tensor provided by Orca and the one generated by Spinach (g2chi) based on g values. To calculate the axial and rhombic parts of the tensors, I use the following definitions, where the eigenvalues are sorted in Mehring order (xx < yy < zz; implemented in the MatLab function mat2axrh function):

chi_ax = zz-(xx+yy)/2
chi_rh = (xx-yy)

This way I get chi_ax values that still differ by a factor of 10. The one from the g2chi tensor actually comes really close to the experimental one.

In the Eu(III) case, I get PCS cones with different signs and the shape of the torus differs slightly. Again, the g2chi value comes closer to the experimental one (apart from the sign).

Eu.png (283.34 KiB) Viewed 127107 times

Now I have three major questions:

1. why are the values from the Orca chi still off by a factor of 10? I checked the unit conversion from cm³*K/mol to A³ multiple times. (odd thing is, the produced PCS values are actually better and the ones from the g2chi tensor are actually too big)

2. what are the limits of the g2chi method? is it a good approximation in these cases? Since it's based on the Curie law I assume it would completely neglect spin orbit coupling?

3. I was wondering about the differences of the signs of the chi_ax values. I realized that chi_ax can never be negative if the Mehring order is used. So can I assume most experimental values are referring to the Haeberlen notation? If I wanted to use this one, would I have to supply a traceless tensor? And would the definition of chi_ax change to 3*zz/2?

Best wishes,

Adrian

[kuprov]

Sorry about the slow reply - teaching season just started here!

The g2chi approximation is indeed quite crude - there are many other mechanisms that can contribute to d2E/dB2, not just spin of course. Electronic structure theory packages would pick those mechanisms up, and that is part of the difference. That said, much would depend on zero-field splitting, and that's of course a parameter that even relativistic electronic structure theory would struggle with.

Conventions is something that we are trying to steer clear from! Too many of them in the literature - this is why Spinach always insists on getting the full chi matrix. On any particular paper, my guess would be as good as yours.

So yes, the summary I think is: blame the chi. As far as I can see from the above, Spinach performs nominally and diligently uses the susceptibility that it is given.

[patroklus97]

No problem, thanks for taking the time to answer!

Thanks for the insights! I think I have now found the source of discrepancy regarding the magnitude of experimental and calculated chi_ax values. While browsing some Kernel utility functions I stumbled across cgsppm2ang, which made me realize that the chi values of quantum chemistry packages are apparently quoted in the form of the irrational susceptibility and need to be multiplied by 4*pi. This seems to account for the missing factor of around 10...

As a side note/question: I was actually a little confused regarding the conversion from cgs-ppm to A³ in the cgs2ppm function. Wouldn't cm³ convert to 1e24 A³ and not 1e18 (as mentioned in the calculation line)?

Furthermore, I think I found a small error in tsm2param: Since in line 43 D was lastly defined as diag(D), the matrix notation for the indices in lines 51 and 52 doesn't work anymore.

For chi I now decided to use a convention from one of your papers (traceless part of chi, chi_ax=3*zz/2, |xx|<|yy|<|zz|), which seems to work quite well and feels like a good representation of those quantities. I double checked my values and found that my signs should actually be correct and the disagreement of theory and experiment rather lies in the fitting procedure of the experimental side and the reference values that were initially used.

[kuprov]

Ah, but remember ppm - that's another factor of 1e-6.

Well spotted for the tsm2param function! Will fix now and push to SVN. I was just staring at it and thinking "how could I miss this?" I blame teaching.