Exact calculation of dipolar interactions considering g-tensors and frames

[akboudalis]

Currently, calculation of dipolar interactions with diptensor does not consider: -the anisotropy of the g-tensors -the relative orientations of the g-frames -the orientation of the interspin vector within the molecular frame (needed for larger than 2-spin systems).

Adding such a feature (potentially by augmenting diptensor) would be very useful for the exact calculation of dipolar interactions, at least within the point-dipole model.

Thanks in advance!

[stestoll]

It should be straightforward to extend diptensor to accommodate this, without disrupting the existing input argument syntax. Here's a possible suggestion. r is the distance vector between the two spins, in units of nm:

% returning 3x3 dipolar coupling tensor
ee = diptensor(g1,g2,r);    % using 3x3 g matrices
ee = diptensor({g1,gFrame1},{g2,gFrame2},r);  % using g principal values and Euler angles

% returning principal values and Euler angles of dipolar coupling tensor
[ee,eeFrame] = diptensor(...);

[akboudalis]

Yes, I believe this would be the basis of a straightforward syntax.

A couple of points:

1. From the documentation of the current version of diptensor, I see the syntax for r is [a; b; c], in nm dimensions. I would see a utility in separating the orientational part from the length, i.e. using an input considering 3 Euler angles + 1 length variable. The reason is that it is often sufficient to start working with idealized symmetries (much more easily defined by Euler angles), leaving the length as a fitting variable for an experimental spectrum. Of course, the Euler angles can also be fitted, if the data permit it.

2. As for the format of the dipolar coupling tensor, in my experience, when calculated exactly, it may need to be decomposed in a symmetric, antisymmetric and a traceless part, just like the exchange tensor (e.g. for different orientations of dissimilar spins). In these cases it cannot be expressed in three principal values + three Euler angles, and all nine elements are needed. In practice, and for more than two spins, I find it practical to calculate the 3x3 tensor in the molecular frame and add it to the other 3x3 exchange terms (also brought to the molecular frame): E.g ee12 = J12 + antisymm12 + ... + dip12 ee13 = J13 + antisymm13 + ... + dip13 ee23 = J23 + antisymm23 + ... + dip23 Sys.ee = [ee12; ee13; ee23] all being 3x3 matrices in the molecular frame. In other words, I would rather see eeFrame as input and a 3x3 matrix as an output.

I hope the above makes sense.

[stestoll]

1. I think this is a special application that can be addressed easily with one extra line of code before calling diptensor.
2. You are right, the output tensor is not necessarily symmetric. So diptensor should provide only the 3x3 matrix as output. I will add an example in the documentation on how to decompose the 3x3 matrix into its isotropic, symmetric and antisymmetric components.
3. I don't see why eeFrame would be needed as input. The dipolar tensor is calculated in this function from the inter-spin vector (defined in the molecular frame), so it doesn't make sense to provide information about the tensor orientation.

[stestoll]

Implemented in 0905ccb . Let me know what you think.

[akboudalis]

Regarding point 2, this should be just fine.

Points 1 and 3 are actually connected since my point was to define the interspin vector orientation by Euler angles. I suppose you are suggesting that if I want to define a set of Euler angles eul = [a, b, c] Then I use

rotmZYZ = eul2rotm(eul,'ZYZ') % Rotation matrix 
rij = [0;0;1]'*rotmZYZ % Unit vector, starting from parallel to the z-axis
rij = rij * dist % Same, but in nm

Would that be compatible with the current diptensor input format?

[stestoll]

The easiest is to take the polar angles that describe the direction from site 1 to site2, represented in the molecular frame

theta = 60*pi/180;  % angle down from the molecular z axis, in radians
phi = 45*pi/180;   % angle in molecular xy plane, in radians
r = 0.8;  % distance, in nm
n12 = ang2vec(phi,theta);  % unit vector from 1 to 2, represented in mol frame
rvec = n12*r;
T = diptensor(g1,g2,rvec)

If you have the Euler angles that relate the site 1 frame (molecular frame) to the site 2 frame, then phi and theta correspond to the first and second Euler angle. You can either use the above code snippet with these two angles or use erot with all three Euler angles (the third one has no effect):

angles = [45 60 10]*pi/180;  % phi, theta, chi in radians
R_1to2 = erot(angles);
n12 = R_1to2(3,:).';

[stestoll]

I am closing this issue as completed. Happy to reopen if needed!