The above spectra were made with a general &`#'915; inversion-recovery program. In this instance the input spin system was a single ⁷Li. The quadrupolar Hamiltonian relaxes the system back to equilibrium and, for some correlation times, produces an observable dynamic frequency shift. Indeed, it almost exactly reproduces the published results of Price, Ge, and Hwang, J. Magn. Reson., 98, 134 (1992). I believe the small discrepancies are due to the inclusion of relaxation during both the delay and acquisition steps in the GAMMA simulation.

Note that this result is obtained from a generic GAMMA IR program. It doesn't use any specialized tensors, symmetry groups, etc. Its just throws together the normal pulse-delay steps for IR and was fed a spin system containing the parameters mentioned in the article. The input system (a single Li) had quadrupolar tensor assigned to its lone spin and relaxation effects were included during each step of the sequence.

This type of simulation is straight forward as there are no complex phase cycles nor worries about quadrature detection on the T1/F1 axis. Here I'll use ideal pulses (you can use square, shaped, composite, pulse trains, or any combination thereof). This particular program will run interactively and output a stack spectrum directly to the screen using Gnuplot.

#include <gamma.h>

 main(int argc, char* argv[])
   {
//			Set Up Simulation Variables

  int qn = 1;						// Query parameter
  sys_dynamic sys;					// A spin system
  String filein;					// Input filename
  filein = sys.ask_read(argc, argv, qn++);		// Ask/read in the system
  WBRExch WBR;						// Relaxation controls
  WBR.read(filein, sys);				// Read relaxation controls
  String IsoD = sys.symbol(0);				// Use 1st spin isotope type
  String IsoP = sys.symbol(0);				// Use 1st spin isotope type
  if(sys.heteronuclear())				// Ask for other possible
    {							// isotopes if heteronuclear
    query_parameter(argc, argv, qn++,			// system was input
               "\n\t180 Pulse Channel? ", IsoP);
    query_parameter(argc, argv, qn++,
        "\n\t90 Pulse, Detect Channel? ", IsoD);
    }
  double SW;
  query_parameter(argc, argv, qn++,			// system was input
           "\n\tPlot Spectral Width? ", SW);
  double FL = SW/2.0;					// Low Plot limit (Hz)
  double FH = -FL; 					// Upper Plot limit (Hz)
  int ntaus;						// Number of tau increments
  double tauinc;					// Tau time increment (sec)
  query_parameter(argc, argv, qn++,			// Ask for how many blocks
             "\n\tNumber of tau Values? ", ntaus);
  query_parameter(argc, argv, qn++,			// Ask for tau increment
    "\n\tTau Incrementation Time(s)? ", tauinc);
  int npts = 4096;

//			Set Up Mathematical Variables

  gen_op U90y = Iypuls_U(sys,IsoD,90.0);		// Propagator 90y pulse ideal
  gen_op H = Ho(sys);					// Isotropic NMR Hamiltonian
  gen_op D = Fm(sys, IsoD);				// Detection operator
  super_op L = complexi*Hsuper(H);			// L0 = -i*[Ho, ] (rad/sec)
  L += WBR.REX(sys,H);                  	        // Add in relaxation
  gen_op sigma0 = sigma_eq(sys);			// System at equilibrium
  acquire1D ACQ(D, L, sigma0);				// Set up acquisition
  row_vector data;					// A 1D data block
  matrix mx, Fmx(ntaus,npts);				// Discrete, Frequency data

//			Apply the Pulse Sequence

  gen_op sigma = Iypuls(sys,sigma0,IsoP,180.0);		// Apply 180y pulse ideal
  gen_op sigmap;					// System after 90 pulse(s)
  for(int i=0; i<ntaus; i++)				// Loop over tau values
    {
    sigmap = evolve(sigma,U90y);			// Apply 90y pulse ideal
    mx = ACQ.table(sigmap);				// Transitions table (no lb)
    data = ACQ.F(mx,npts,FL,FH);			// Get block at this tau
    Fmx.put_block(i, 0, data);				// Store block in array
    sigma = evolve(sigma-sigma0,L,tauinc) + sigma0;	// Evolve to next tau
    }

//			Output the Plot Using Gnuplot

   double time = double(ntaus-1)*tauinc;
   GP_stack(String("IR.asc"), Fmx, 0,time,FL,FH);
   ofstream gnuload("IR.gnu");				// File of gnuplot commands
   gnuload << "set data style line\n";			// Set plot to use lines
   gnuload << "set parametric\n";			// Set parametric for stack
   gnuload << "set xlabel \"W (Hz)\"\n";		// Set X axis label
   gnuload << "set ylabel \"time (sec)\"\n";		// Set Y axis label
   gnuload << "set title\"Inversion Recovery\"\n";	// Set plot title
   gnuload << "splot \"IR.asc\"\n";			// Stack 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 \"IR.gnu\"\n");			// Plot to screen
   cout << "\n\n";					// Keep the screen nice
  }

The program takes a spin system file name as input. Here is an example of such a file. You'll have to add in quadrupolar and shielding tensors as well as exchange rates if you want to include such effects.

SysName   (2) : AMX       		- 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 M
Iso(2)    (2) : 1H			- Spin Isotope Type X
Omega	  (1) : 750.00  		- Spect. Freq. in MHz (1H based)
v(0)      (1) : 500.0			- Chemical Shift (Hz)
v(1)      (1) : -500.0       		- Chemical Shift (Hz)
v(2)      (1) : 250.00       		- Chemical Shift (Hz)
J(0,1)    (1) : 10.00			- Coupling Constant
J(0,2)    (1) : 12.0			- Coupling Constant
J(1,2)    (1) : 7.0			- Coupling Constant
Coord(0)  (3) : (0.0,0.0,0.0)		- Coordinate of spin A
Coord(1)  (3) : (0.0,0.0,1.5)		- Coordinate of spin M
Coord(2)  (3) : (1.0,1.0,1.5)		- Coordinate of spin X
Taus      (3) : (1.50, 1.50, 1.00)	- Correlation times (nsec)

Rlevel    (0) : 4               - Relaxation Computation Level
Rtype     (0) : 0               - Relaxation Computation Type
RDD       (0) : 1               - Dipolar Relaxation Flag
RDDdfs    (0) : 1               - Dipolar DFS Flag
RCC       (0) : 1               - CSA Relaxation Flag
RCCdfs    (0) : 1               - CSA DFS Flag
RQQ       (0) : 1               - Quad Relaxation Flag
RQQdfs    (0) : 1               - Quad DFS Flag
RDQ       (0) : 1               - Dip-Quad Relaxation Flag
RDQdfs    (0) : 1               - Quad DFS Flag
RDC       (0) : 1               - Dip-CSA Relaxation Flag
RDCdfs    (0) : 1               - Dip-CSA DFS Flag
RQC       (0) : 1               - Quad-CSA Relaxation Flag
RQCdfs    (0) : 1               - Quad-CSA DFS Flag