Program for chemical exchange

[Stefan Stoll]

Hi all,

While EasySpin currently does not support the simulation of EPR spectra for chemical exchange, I have spent some time a while back to implement it. Attached the code, with the main function exchange() and several example. It works reasonably well, and might get incorporated into EasySpin at some point in the future. Meanwhile, download it from here and enjoy!

If you publish simulations based on this code, please cite dx.doi.org/10.1021/jp3104358 . That's where the code has first been used.

Stefan

[katarkon]

I have made some corrections in exchange.m function. Now it correctly works with separate linewidths for each site. Also, the site's populations are calculated automatically from 2D Sys.k matrix (by default) or defined manually using Sys.Pop field. In other cases the equal populations are assumed.
Known restrictions.
1. Only general methods (1 and 2) now supported.
2. Exchange sites should be defined manually (by Sys.A matix). Sys.xc seems to be working too, but it was not closely tested.

Small example how to use the modified function

Code: Select all

%% DOI: 10.1039/B407972H Dalton Trans., 2004, 2957-2962
% Simulation for complex 7
clear
% defining exchanging sites
%
% every one have own g-factor
Sys.g(1,:) = 2.00635;
Sys.g(2,:) = 2.0041;
% two equivalent phosphoruses and one hydrogen HFCs
Sys.Nucs = '31P,1H';
Sys.n = [2 1];
Sys.A(1,:) = mt2mhz([24.4 3.4]/10);
Sys.A(2,:) = mt2mhz([18.4 2.4]/10);
% and own lw
Sys.lwpp(1,:) = 0.3;
Sys.lwpp(2,:) = 0.1;
% defining rate constants
% 2D matrix have the form
% |  0	k<21>	k<31>   ...  k<n1> |
% | k<12> 0     k<32>   ...  k<n2> |
% | ...  ...    ...     ...   ...  |
% | k<1n> k<2n> k<3n>   ...    0   |
%
% k<ij> means rate constant of reaction
% Site<i> -> Site<j>
% in MHz (or 10^-6 s^-1)
 
% for mutual exchange k<ij>=k<ji> and simplified form may be used
% Sys.k = [k<12> k<13> ... k<1n> k<23> k<24> ... k<34> ...];
Sys.k = [0 45; 90 0]; %corresponds to 270K

%Populations are calculated automatically from 2D k matrix
%If necessary, Sys.Pop may be used instead
%Sys.Pop = [0.5 1];
Exp.nPoints = 4096;
Exp.Freq = 9.65; % in GHz

Opt.Method = 2;
exchange(Sys,Exp,Opt);

UPD. Some wrapper functions added, function exchange() was patched for correct working with nuclei with I=0.
exchange1() - works with isotope mixtures. It uses undocumented function isotopologues() and call exchange() function for each isotopologue with further accumulation. Exp.Range must be defined, separated output was not tested.
exchange2() - convolute each line with additional broadening. Sys.lb must be defined. It contains n elements for n lines in spectrum according to their appearance with field increasing. Each one is additional broadening in mT (or zero if no additional broadening for this line). Exp.Range must be defined. Separated output and isotope mixtures are not supported yet.
UPD. The function exchange() function now works with additional Gaussian broadening. Use Sys.lwg field to define it. Separate Gaussian broadening for each site is not supported.
UPD. The wrapper exchange2() was modified. If only one element in Sys.lb is defined it assumed to be a maximal line broadening. Minimal broadening assumed to be a zero. Other broadenings estimates according to fast motion regime broadening. If sign of Sys.A is positive, the broadening increase with the field increasing, if negative - in opposite order. Also, if Sys.lb contains n element for n lines, it will be used as previously.
In work. Simulation of additional broadening caused by 'thumbling effect' from g and A principal values and correlation time (like function chili).

[katarkon]

Chemical exchange simulation - getting started.
Installation.
Download the archives from the upper posts. Extract the 'chemical exchange' folder into your easyspin working directory. Replace the 'exchange.m'
file by one from second 'exchange.zip' archive. Compile files for your platform. Just type 'mex(exlinesum.c)' and 'mex(exlinesum2.c)' into matlab command window. Now you can write your own scripts in the 'chemical exchange' folder.
Usage of exchange function.
First let's define the exchanging sites.
The spin system for each site can be defined either by permutation of the nuclei of by manual definition of each site.
Definition of nuclei permutation is follow:

Code: Select all

Sys.g = 2.0023;
Sys.Nucs = '14N,14N,1H,1H,1H';
Sys.n = [1 1 3 1 1];
Sys.A = [15 5 1 2 5]; %in MHz
Sys.xc = [2 1 3 5 4];

This code defines two-site exchange, the nuclei 1 and 2 are swaped as well as 4 and 5 ones whereas nucleus 3 is the same for both sites.

Code: Select all

Sys.Nucs = '1H,1H,1H';
Sys.n = [1 1 1];
Sys.A = [15 5 1]; %in MHz
Sys.xc = [2 3 1; 3 1 2];

This code defines three-site exchange, second site is a result of swapind of the nuclei 1-2 and 2-3, third one is a result of swapind of the nuclei 1-3 and 2-1.
The other way is defining the exchange sites manually.

Code: Select all

Sys.Nucs = '1H,1H,1H';
Sys.n = [1 1 1];
Sys.g(1,:) = 2.0013;
Sys.g(2,:) = 2.0023;
Sys.g(3,:) = 2.0033;
Sys.A(1,:) = [1 2 3]; %in MHz
Sys.A(2,:) = [2 3 1]; %in MHz
Sys.A(3,:) = [3 2 1]; %in MHz

Next step is defining the line width (in mT):

Code: Select all

Sys.lwpp = 0.015;

or

Code: Select all

Sys.lw = 0.015;

If separate line width is necessary for every site, use the follow sintax

Code: Select all

Sys.lwpp(1,:) = 0.015;
Sys.lwpp(2,:) = 0.025;
Sys.lwpp(3,:) = 0.005;

At last, the exchange rates must be defined. It can be defined as 2D (N x N) matrix in which non-diagonal elements are rate constants for each site-to-site exchange reaction. The direct constants (Site 1 -> Site 2) constants are placed in the columns whereas reverse ones (Site2 -> Site 1) in the rows. Diagonal elements may have any values (i.e. zeroes):

Code: Select all

Sys.k = [0 1 1; 2 0 2; 1 1 0];

This code defines k<12>=2MHz, k<13>=1MHz, k<21>=1MHz, k<23>=1MHz, k<31>=1MHz, k<32>=2MHz where k<ij> means the rate constant for Site i -> Site j reaction.
Another way assumes that k<ij>=k<ji> for mutual exchange. The sintax is quite simple:

Code: Select all

Sys.k = [k<12> k<13> ... k<1n> k<23> k<24> ... k<34> ...];

The populations of each site is calculated by default if rate constants are defined by 2D matrix. Also they may be defined manually by defining the weight coefficient for each site:

Code: Select all

Sys.Pop=[1 2 0.5 ...];

At last, experimental parameters should be given:

Code: Select all

Exp.nPoints = 4096; %number of points
Exp.mwFreq = 9.5; %microwave frequency in GHh
Exp.Range = [337.5 340.5]; % the range in mT, calculates automatically if not given

And some options as well:

Code: Select all

Opt.Method = 2; 
Opt.reduce = 0; 
Opt.Output = 'separate';

Opt.Method can be a number from -3 to 3 and means follow:
-3 Explicit method, two-site mutual exchange, Matlab accumulation
-2 Explicit method, two-site mutual exchange, mex accumulation
-1 Liouville method, two-site mutual exchange, Matlab accumulation
0 Liouville method, two-site mutual exchange, mex accumulation
1 Liouville method, general n-site exchange, mex accumulation
2 Liouville method, general n-site exchange, Matlab accumulation
3 Bloch/McConnell, two-site mutual exchange
Opt.reduce (1 by default) - reduction of dublicates of lines. If set to zero, the reduction is neglected.
Opt.Output ('summed' by default). If set to 'separate', the each line draws as separate spectrum, else all lines are summed in single spectrum.
Mex acumulation works faster but incompatible wits separate output. Only general exchange methods (1 and 2) are compatible with separate linewidths for each site.
And finally, the function exchange can be called.

Code: Select all

exchange(Sys,Exp,Opt);

Just draw the simulated spectrum.

Code: Select all

[x,y]=exchange(Sys,Exp,Opt);

Returns spectrum as a pair of x,y values for thuther processing.

[katarkon]

The article https://doi.org/10.1016/j.molstruc.2017.07.001 came out. It illustrates the usage of modified exchange() function.