src/BWRRelax/relaxNMR.cc
// ------------------ Level 3 via Double Commutators --------------------
//                    (Applies Secular Approximation)
Level 3 via Double Commutators
(Applies Secular Approximation)
Level 3 via Double Commutators ( Applies Secular Approximation )
void R_3s(super_op& LOp, double* w, int rank, gen_op* T1s, gen_op* T2s, matrix& J12)
					       LOp    = LOp   + R3
					          out	   in     12
	 Note				T1s, T2s, and LOp assumed in proper bases
            --- --- ---
            \   \   \       m                 m       -m
    R3    = /   /   /   (-1)  delta        [ T   , [ T    ,  ] ] J  (w  )
	1,2   --- --- ---            w , -w     1,p     2,p'        12  p'
	       p   p'  m              p    p'
 where p and p' sum over transitions, m sums over angular momentum components,
 and J values include xi1 and xi2.  The secular approximation restricts the
 sum to contain only terms where the transition frequencies for p & p' cancel.
   {
   double cutoff = 1.e-9;			 Secular approximation cutoff
   int hs = T1s[0].dim();			 Get Hilbert space size
   matrix mx1(hs, hs, 0.0);			 Two Zeroed Hilbert space matrices
   matrix mx2(hs, hs, 0.0);
   basis bs = T1s[0].get_basis();		 Current working basis
   gen_op nullOp;				 A NULL operator
   double wp, wpp;				 Transition frequencies for p & p'
   double Jwp;					 Spectral density, transition p (& p')
   complex Ta_aa, Tb_bb;
   gen_op *T1ps, *T2ps;				 Components of the two spin tensors
   T1ps = new gen_op[2*rank+1]; 		 (2*l+1 of these for each p)
   T2ps = new gen_op[2*rank+1];			 (2*l+1 of these for each p')
   int m=0, k=0;				 Angular momentum indicies
   for(int a=0; a<hs; a++)			 Sum over transitions a-aa: p
     for(int aa=0; aa<hs; aa++)
       {
       wp = w[a] - w[aa];			 Transition frequency for a-aa: p
       Jwp = Re(J12.get(a,aa));			 Get spectral density for transition p
       for(m=-rank, k=0; m<=rank; m++, k++)	 Perform a pre-summation over m
         {					 and set up the spin tensor operator
         T1ps[k] = nullOp;			 components for this transition
         Ta_aa = T1s[k].get(a,aa);		                       m
         if(Re(Ta_aa) || Im(Ta_aa))		           T1ps[k] == T
           {					                       1,p
           mx1.put(Ta_aa, a, aa);		 
           T1ps[k] = gen_op(mx1,bs);		 These are either sparse or null arrays
           mx1.put(0.0, a, aa);			 so we insure that null components are
           }					 recognized (and not used in computation).
         }
       for(int b=0; b<hs; b++)			 Sum over transitions b-bb: p'
         for(int bb=0; bb<hs; bb++)
           {
           wpp = w[b] - w[bb];			 Transition b-bb: p' frequency
	   if(fabs(wp+wpp) < cutoff)		 Apply the secular approximation: p = -p'
             {
             for(m=-rank,k=0; m<=rank; m++,k++)	 Perform a pre-summation over m
               {				 and set up the spin tensor operator
               T2ps[k] = nullOp;		 components for this transition
               Tb_bb = T2s[k].get(b,bb);	                       m
               if(Re(Tb_bb) || Im(Tb_bb))	           T2ps[k] == T
                 {				                       2,p'
                 mx2.put(Tb_bb, b, bb);		
                 T2ps[k] = gen_op(mx2,bs);	 These are either sparse or null arrays
                 mx2.put(0.0, b, bb);		 so we insure that the null components are
                 }				 recognized (and not used in computaiton).
               }
             R_CC_0(T1ps, T2ps, LOp, rank, Jwp); Perform sum over m, add LOp    to LOp
             }								      pp'
           }					 Next transition p'
       }					 Next transition p
   return;
   }
void R_3s ( super_op & LOp , double * w , int rank , gen_op * T1s , gen_op * T2s , matrix & J12 ) { double cutoff = 1.e-9 ; int hs = T1s [ 0 ] . dim ( ) ; matrix mx1 ( hs , hs , 0.0 ) ; matrix mx2 ( hs , hs , 0.0 ) ; basis bs = T1s [ 0 ] . get_basis ( ) ; gen_op nullOp ; double wp , wpp ; double Jwp ; complex Ta_aa , Tb_bb ; gen_op * T1ps , * T2ps ; T1ps = new gen_op [ 2 * rank + 1 ] ; T2ps = new gen_op [ 2 * rank + 1 ] ; int m = 0 , k = 0 ; for ( int a = 0 ; a < hs ; a ++ ) for ( int aa = 0 ; aa < hs ; aa ++ ) { wp = w [ a ] - w [ aa ] ; Jwp = Re ( J12 . get ( a , aa ) ) ; for ( m = - rank , k = 0 ; m <= rank ; m ++ , k ++ ) { T1ps [ k ] = nullOp ; Ta_aa = T1s [ k ] . get ( a , aa ) ; if ( Re ( Ta_aa ) || Im ( Ta_aa ) ) { mx1 . put ( Ta_aa , a , aa ) ; T1ps [ k ] = gen_op ( mx1 , bs ) ; mx1 . put ( 0.0 , a , aa ) ; } } for ( int b = 0 ; b < hs ; b ++ ) for ( int bb = 0 ; bb < hs ; bb ++ ) { wpp = w [ b ] - w [ bb ] ; if ( fabs ( wp + wpp ) < cutoff ) { for ( m = - rank , k = 0 ; m <= rank ; m ++ , k ++ ) { T2ps [ k ] = nullOp ; Tb_bb = T2s [ k ] . get ( b , bb ) ; if ( Re ( Tb_bb ) || Im ( Tb_bb ) ) { mx2 . put ( Tb_bb , b , bb ) ; T2ps [ k ] = gen_op ( mx2 , bs ) ; mx2 . put ( 0.0 , b , bb ) ; } } R_CC_0 ( T1ps , T2ps , LOp , rank , Jwp ) ; } } } return ; }
Level 3 via Double Commutators
(Applies Secular Approximation)
Level 3 via Double Commutators ( Applies Secular Approximation )
// Input		LOp   : Current relaxation superoperator //			w     : Vector of system energies //				in rad/sec in the LAB frame //			rank  : Rank of the two interactions // 			T1s   : Spin tensor components, 1st interaction // 			T2s   : Spin tensor components, 2nd interaction //			J12   : Matrix of reduced spectral density  //				function values for 2 interactions //				times the two interaction constants // Return		void  : Level 3 relaxation superop for interactions //			        1 & 2 added to input superoperator LOp // //				       LOp    = LOp   + R3 //				          out	   in     12 // // Note				T1s, T2s, and LOp assumed in proper bases //            --- --- --- //            \   \   \       m                 m       -m //    R3    = /   /   /   (-1)  delta        [ T   , [ T    ,  ] ] J  (w  ) //	1,2   --- --- ---            w , -w     1,p     2,p'        12  p' //	       p   p'  m              p    p' // where p and p' sum over transitions, m sums over angular momentum components, // and J values include xi1 and xi2.  The secular approximation restricts the // sum to contain only terms where the transition frequencies for p & p' cancel. // Secular approximation cutoff // Get Hilbert space size // Two Zeroed Hilbert space matrices // Current working basis // A NULL operator // Transition frequencies for p & p' // Spectral density, transition p (& p') // Components of the two spin tensors // (2*l+1 of these for each p) // (2*l+1 of these for each p') // Angular momentum indicies // Sum over transitions a-aa: p // Transition frequency for a-aa: p // Get spectral density for transition p // Perform a pre-summation over m // and set up the spin tensor operator // components for this transition //                       m //           T1ps[k] == T //                       1,p //  // These are either sparse or null arrays // so we insure that null components are // recognized (and not used in computation). // Sum over transitions b-bb: p' // Transition b-bb: p' frequency // Apply the secular approximation: p = -p' // Perform a pre-summation over m // and set up the spin tensor operator // components for this transition //                       m //           T2ps[k] == T //                       2,p' // // These are either sparse or null arrays // so we insure that the null components are // recognized (and not used in computaiton). // Perform sum over m, add LOp    to LOp //			      pp' // Next transition p' // Next transition p