Pseudocontact shift calculation from spin density

[ppikulova]

Dear Prof. Kuprov,

I have been working with Spinach's paramagnetic NMR functionalities and have some suggestions about the function ocparse.m used to import spin density from a .3d Orca output and a theoretical question about the PCS calculation.

First, a technicality regarding the handling of box extents. I have a .3d file with 80 grid points along each dimension, the X range being [-7 A, 14 A]. After the density is padded with zeroes, the following line of code calculates new extents of the box in Angstrom:

Code: Select all

ext=(2*pad_factor+1)*ext;

Say I have pad_factor=1, this means adding 80 grid points on each side. Physically, this must correspond to extending the box by 21 Angstrom to each side, so the new X range should be [-28 A, 35 A]. The code above gives (2*1+1)*[-7 A, 14 A] = [-21 A, 42 A]. With the incorrect box extents, subsequent calculation of PCS with kpcs.m gives nonsensical results which change wildly with pad_factor. The original code can work only if the cube is distributed symmetrically around 0, which is accidentally the case for the Cu porphyrin example (it confused me why the same thing is working for the example .3d file but not for mine).

For my purposes, I modified the code so that the new box extents are calculated as:

Code: Select all

ext=[ext(1)-abs(ext(2)-ext(1))*pad_factor ext(2)+abs(ext(2)-ext(1))*pad_factor...
ext(3)-abs(ext(4)-ext(3))*pad_factor ext(4)+abs(ext(4)-ext(3))*pad_factor...
ext(5)-abs(ext(6)-ext(5))*pad_factor ext(6)+abs(ext(6)-ext(5))*pad_factor];

which can also handle spin density cubes not centered around origin of the Cartesian system. I suggest that it would be a good idea to update the official version of ocparse.m in a similar way to make it more generally usable.

I also believe that the following line of code in ocparse.m

Code: Select all

density=permute(density,[2 3 1]);

should actually be

Code: Select all

density=permute(density,[3 2 1]);

given how data is stored in the .3d file with X as the outer loop and Z as the inner loop. I found this half empirically, using the former results in the following deformed and wrongly placed spin density (sorry about the crazy connectivity, I used a high threshold for bond lengths):

spd_231.png (42.27 KiB) Viewed 714 times

while using the latter works fine (positive spd located in a dxy orbital at the metal atom with delocalization to coordinated Cl- ligands):

spd_321.png (39.02 KiB) Viewed 714 times

So it seems that the original code incorrectly switches the X and Y directions. (Again, I believe that this incidentally does not produce an error for the Cu porphyrin example because the model is axially symmetric around Z).

With these issues fixed, the procedure outlined in the Cu porphyrin example finally gives reasonable pseudocontact shifts which converge with increasing pad_factor for my system, so I now look forward to playing with it more

Last, I have a question about the physics. I am curious about the link between the paramagnetic centre probability distribution the PDE operates with and spin density output from QC calculations. I noticed ocparse.m takes the absolute value of the spin density. In your opinion, would it have meaning to use the original spin density with both signs for the PCS calculation? Granted, any beta spin density is a few orders of magnitude smaller than alpha spin density, so this effect should be relatively small. Still, I would be really grateful for any insight about this last point so that I can understand the calculation better

Best regards
Petra Pikulová

[kuprov]

OK, looking at the extents matter... you are right of course, thanks for spotting that bug! Indeed, the example file has symmetric extents and therefore works by pure accident. I'll update the extent handling.

[kuprov]

Also correct for the XYZ order - that would teach me not to use highly symmetric examples for testing...

[kuprov]

Finally on physics. The partial differential equation that is being solved was derived for a classical probability density of a classical magnetic susceptibility centre. In quantum mechanical calculations, we assume that one or more unpaired electrons can be treated in that way. So the density does strictly need to be positive because of its physical meaning.

[ppikulova]

Hello, interesting.

Of cource, I am influenced by quantum chemical calculations in thinking about the negative spin density. Though I feel a little bit excused by the fact that you draw the parallel yourself in one of the papers.

The system is a doublet, which makes the situation not so hard to imagine. I know the point dipole equation is derived by Boltzmann averaging of the electron magnetic moment over the electron spin states (and averaging over all molecular orientations) and the PCS field results from what is left over. My understanding of the physical sense of the PDE is that instead of one point dipole localized at the paramagnetic centre, we now have a sort of point dipole at every point in space where the probability of finding the unpaired electron is finite (and we suppose that the susceptibility is the same at every point in space). Is this accurate enough?

My thinking was that if QC gives negative spin density delocalized to some region of the molecule and this is incidentally close to the probe nucleus, artificially turning it into positive spin density might skew the resulting PCS somewhat. But I repeat that I wouldn't even expect it to be a large effect, as the negative spin density is small. It was not a very practically intended question, only driven by curiosity.

Thanks for your consideration of this non-expert question

Petra P.

[kuprov]

Yes, the PDE essentially treats the problem as multiple tiny susceptibility centres at each voxel - a distributed susceptibility density. To make the problem tractable, we assume the same tensor everywhere, and only the amplitude is variable. Of course, the amplitude has to stay positive.