GAMMA has the ability to treat both mutual and non-mutual exchange. This is done in two ways. First, it supplies functions which return exchange superoperators for a given spin system. These are then used for system evolution steps under exchange (and relaxation) effects. Second, GAMMA supplies extended spin systems - systems that are made of interacting spin systems. Each individual system has defined spin and spin-spin interactions and the extended system has interactions between spins in differing systems.

Mutual Exchange

The above results were from a GAMMA mutual exchange simulation using an input ABC spin system. It was used for comparison to the article by Gerhard Binsch

"A Unified Theory of Exchange Effects on Nuclear Magnetic Resonance Line Shapes", J. Am. Chem. Soc., 91, 1304-1309, (1969)

The figure contains three spectra placed side by side for comparison. The middle spectrum (K=5/sec) nicely reproduces Figure 1 of the Binsch article. Here is the source code. It is a typical GAMMA 1D NMR simulation program in Liouville space as it uses an exchange superoperator. It runs interactively and plots to the screen using gnuplot. The code isn't totally general in that it detects protons, pulses all spins, and plots only between 0-100 Hz. It also adds in an additional linewidth that Binsch used. You can make modifications to generalize it as needed. Note that it must be run three times with differing input rates to get the three spectra that made the composite figure above. I made the code lines that have to do with mutual exchange in red...

#include <gamma.h>

 int main(int argc, char* argv[], int argn)
   {
   sys_dynamic sys;                             // Declare a sys dynamic
   sys.ask_read(argc, argv, 1);                 // Ask/Read in sys
   gen_op H = Ho(sys);                          // Isotropic Hamiltonian
   gen_op D = Fm(sys, "1H");                    // Detect F-, protons
   gen_op sigma0 = sigma_eq(sys);               // System at equilibrium
   gen_op sigma = Iypuls(sys,sigma0,90.0);      // Apply 90y pulse ideal
   super_op L = complexi*Hsuper(H);             // Liouvillian
   L += Kex(sys, H.get_basis());                // Add in exchange
   acquire1D ACQ(D, L, 1.e-3);                  // Set up acquisition
   TTable1D TT = ACQ.table(sigma);              // Transitions table (no lb)
   cout << "\n\n\t\t\t\tSystem Transitions\n"   // Look at the transitions
        << TT;
   TT.broaden(1.25,0);                          // Add .8sec T2 (1.25/sec R2)
   row_vector data = TT.F(4096,0,100.0);        // Frequency spectrum
   GP_1D("gp.gnu", data, 0, 0,100.0);           // * Output data to gp.gnu, gnuplot
   ofstream gnuload("gnu.dat");                 // * File of gnuplot commands
   gnuload << "set data style line\n";          // * Command to plot using lines
   gnuload << "set xlabel \"v(Hz)\"\n";         // * Command to set X axis label
   gnuload << "set title\"Binsch Test\"\n";     // * Command to set plot title
   gnuload << "plot \"gp.gnu\"\n";              // * Command to plot gp.gnu
   gnuload << "pause -1 \' To Exit \n";         // * Command to wait
   gnuload << "exit\n";                         // * Command to exit
   gnuload.close();                             // * Close gnuplot command file
   system("gnuplot \"gnu.dat\"\n");             // * Plot to screen
   cout << "\n\n";                              // Keep the screen nice
   }

Code lines that deal with plotting have a * at the start of their comments. You may replace these with other output types or even simpler lines that omit plot labels. The program takes a spin system file name as input. Here is one for the Binsch article

SysName    (2) : ABC_Binsch              - Name of the Spin System
NSpins     (0) : 3                       - Number of Spins in the System
Iso(0)     (2) : 1H                      - Spin Isotope Type A
Iso(1)     (2) : 1H                      - Spin Isotope Type B
Iso(2)     (2) : 1H                      - Spin Isotope Type C
Omega      (1) : 500.00                  - Spect. Freq. in MHz (1H based)
v(0)       (1) : 19.6                    - Chemical Shift (Hz)
v(1)       (1) : 22.9                    - Chemical Shift (Hz)
v(2)       (1) : 52.7                    - Chemical Shift (Hz)
J(0,1)     (1) : 12.3                    - Coupling Constant
J(0,2)     (1) : 7.3                     - Coupling Constant
J(1,2)     (1) : -17.7                   - Coupling Constant
MutExch(0) (2) : (0,1)                   - Mutual Exchange process
Kex(0)     (1) : 0.5                     - Mutual Exchange Rate (1/Sec)

In the above system, only one exchange process is defined: spin 0 <--> spin 1 with an exchange rate of 0.5/sec. You may define more spins in the process and any number of exchange processes. For example, MutExch(1) with a value (0,2,1) would set a second exchange process with 0 -> 2 -> 1 -> 0 and one would set an associated exchange rate with parameter Kex(1).

Note that in this program the only things specific to exchange were 1.) the declaration of exchange processes in the spin system and 2.) addition of the exchange superoperator to the system Liouvillian. There is nothing preventing you from adding in relaxation, using shaped pulses, performing a multi-dimensional NMR simulation, etc. etc. The rest of the GAMMA machinery is still at your disposal.

Non-Mutual Exchange

Instead, one just follows the same programming ideas but uses an extended spin system, multi_sys. Each multi_sys spin system contains any number of independent spin systems that may be involved in non-mutual exchange with each other. Each system may be assigned a relative population and GAMMA will automatically increase the dimension (spin Hilbert space) of the problem to be the sum of the individual system Hilbert space (as opposed to the product of all spin Hilbert spaces).

For example, compare this non-mutual exchange program with the mutual exchange program shown previously on this page.

#include <gamma.h>

int main()
  {
  multi_sys sys;                                // Declare a spin system
  sys.read("ATPADP.msys");                      // Read it in
  gen_op H = Ho(sys);                           // Isotropic Hamiltonian
  gen_op D = Fm(sys, "31P");                    // Detection F-, 31P
  gen_op sigma0 = sigma_eq(sys);                // System at equilibrium
  gen_op sigmaP = Iypuls(sys,sigma0,90.0);      // Apply a hard 90y pulse
  super_op L = -complexi*Lo(sys);               // Commutation H superop
  L += Xnm(sys);                                // Add in non-mutual exchange
  acquire1D ACQ(D,L,1.e-3);                     // Set up acquisition
  TTable1D TT= ACQ.table(sigmaP);               // Transitions table (no lb)
   cout << "\n\n\t\t\t\tSystem Transitions\n"   // Look at the transitions
        << TT;
  TT.broaden(8.0,1);                            // Add in 8Hz lwhh
  double p_low  = 0.0, p_high = 23.0;           // Plot limits in PPM
  double f_low  = 40.3*p_low;                   // Low plot limit in Hz
  double f_high = 40.3*p_high;                  // Hight plot limit in Hz
  string fname = "ATP.asc";                     // Specrum data file
  row_vector data = TT.F(4096,f_low,f_high);    // Frequency acquisition
  GP_1D(fname, data, 0, p_low, p_high);         // *Write in gnuplot  (ASCII)
  ofstream gnuload("ATP.gnu");                  // *File of gnuplot commands
  gnuload << "set data style lines\n";          // *Set plot to use lines
  gnuload << "set noparametric\n";              // *Set parametric for stack
  gnuload << "set xlabel \"w (PPM)\"\n";        // *Set X axis label
  gnuload << "set title\"ATP <-> ADP + P\"\n";  // *Set plot title
  gnuload << "plot \"" << fname << "\"\n";      // *Plot in gnuplot
  gnuload << "pause -1 \' To Exit \n";		// *Pause before exit
  gnuload << "exit\n";                          // *Exit this mess
  gnuload.close();                              // *Close gnuplot command file
  system("gnuplot \"ATP.gnu\"\n");              // *Plot to screen
  cout << "\n\n";                               // Keep the screen nice
  }

Other than some cosmetic changes, this simple program is nearly identical to the mutual exchange program (which is identical to a typical 1D NMR simulation program in GAMMA). I have highlighted the main code differences in red. Note that the system has been switched to type multi_sys and the exchange superoperator to Xnm for non-mutual exchange. Each component of the multi_sys is itself a spin system that may be undergoing relaxation and/or involved in mutual exchange!

The ASCII file(s) that define these systems are slightly more complicated, but only slightly. For example, here is a spin system that might be used in the above program.

MSysName   (2) : ATP/ADP                - Name of the multi system
NComp      (0) : 3                      - Number of components in multi system
Popul(0)   (1) : 0.8                    - Population of component 0
Popul(1)   (1) : 1.0                    - Population of component 1
Popul(2)   (1) : 1.0                    - Population of component 3
Exch(0)    (2) : (0<=>1+2)              - Non-mutual exchange scheme
Smap(0,0)  (2) : (0)0(1)0               - Cmp 0 0th spin <-> Cmp 1 0th spin
Smap(0,1)  (2) : (0)1(1)1               - Cmp 0 1st spin <-> Cmp 1 1st spin
Smap(0,2)  (2) : (0)2(2)0               - Cmp 0 2nd spin <-> Cmp 2 0th spin
Kex_nm(0)  (1) : 50.0                   - Exchange rate (1/sec)

                                  ATP

[0]SysName   (2) : ATP                  - Name of the Spin System
[0]NSpins    (0) : 3                    - Number of Spins in the System
[0]Iso(0)    (2) : 31P                  - Spin 0 Isotope Type
[0]Iso(1)    (2) : 31P                  - Spin 1 Isotope Type
[0]Iso(2)    (2) : 31P                  - Spin 2 Isotope Type
[0]Omega     (1) : 100.00               - Spect. Freq. in MHz (1H based)
[0]PPM(0)    (1) : 11.0                 - Chemical Shift (ppm)
[0]PPM(1)    (1) : 19.4                 - Chemical Shift (ppm)
[0]PPM(2)    (1) : 4.8                  - Chemical Shift (ppm)
[0]J(0,1)    (1) : 15.0                 - Coupling Constant
[0]J(0,2)    (1) : 0.0                  - Coupling Constant
[0]J(1,2)    (1) : 15.0                 - Coupling Constant

                                ADP

[1]SysName   (2) : ADP                  - Name of the Spin System
[1]NSpins    (0) : 2                    - Number of Spins in the System
[1]Iso(0)    (2) : 31P                  - Spin 0 Isotope Type (alpha)
[1]Iso(1)    (2) : 31P                  - Spin 1 Isotope Type (beta)
[1]Omega     (1) : 100.00               - Spect. Freq. in MHz (1H based)
[1]PPM(0)    (1) : 10.4                 - Chemical Shift (ppm)
[1]PPM(1)    (1) : 4.4                  - Chemical Shift (ppm)
[1]J(0,1)    (1) : 18.0                 - Coupling Constant

                            Free Phospate

[2]SysName   (2) : free_P               - Name of the Spin System
[2]NSpins    (0) : 1                    - Number of Spins in the System
[2]Iso(0)    (2) : 31P                  - Spin 0 Isotope Type
[2]Omega          (1) : 100.00          - Spect. Freq. in MHz (1H based)
[2]PPM(0)    (1) : 2.9                  - Chemical Shift (ppm)

As you can see, the multisystem is defined with specified spin systems, their relative populations, and their relative amounts. I have colored those parameters in red. Then each component system has it' own (usual) parameters except that they are prefaced with [#] to indicate which system they belong to. In the file shown the first (colored purple) has its parameters prefixed with [0], the second (colored blue) prefixed with [1], and so on. This indexing scheme allows users to define as many spin systems as needed in a single file. Alternatively, one may put each component system in its own external file (see class multi_sys documentation).

The spin system data presented above is roughly taken from the article

K.V. Vasavada, J.I. Kaplan, and B.D. Nageswara Rao, J. Magn. Reson., 41, 467-482, (1980)

to simulate the exhange between ATP, ADP, and free phosphate.

Below are the results from the above program and input file for three different exchange rates. Note that it must be run three times with differing input rates to get the three spectra that made the composite figure below. Also, this does not accurately reflect the experimental exchange spectra... why? Because we have not included any relaxation in our simulation. ATP and ADP are both large enough that their motion in a liquid will produce 31P relaxation via shift anisotropy that broadens their resonances, at least it more that that from free phosphate.

By adding the CSA relaxation superoperator to the above program (1 line of code) and defining shift anisotropy tensors and correlation times in our spin systems one might do a lot better. For example, that is how I generated the spectra in the following figure, although this was done in a single program run with Adobe FrameMaker output followed by some quick cosmetic work.

Here are some other exchange programs you might want to give a try. The icons show what they produce.

A-B Mutual Exchange	A-B Non-Mutual Exchange