src/BWRRelax/relaxNMR.cc
// ____________________________________________________________________________
// B            RELAXATION SUPEROPERATORS VIA DOUBLE COMMUTATORS
// ____________________________________________________________________________
B            RELAXATION SUPEROPERATORS VIA DOUBLE COMMUTATORS
B RELAXATION SUPEROPERATORS VIA DOUBLE COMMUTATORS
void R_4s(super_op& LOp, int rank, gen_op* T1s, gen_op* T2s, matrix& J12)
					       LOp    = LOp   + R4
					          out	   in     12
	 Note				T1s, T2s, & LOp assumed in proper bases
                   --- ---
                   \   \       m    m     -m
	     R4    = /   /   (-1)  [ T , [ T   ,  ] ] J  (w )
	       1,2   --- ---          1     2,p        12  p
		       p   m
 where p sums over transitions, m sums over angular momentum components,
 and J values include xi1 and xi2
   {
   int hs = T1s[0].dim();			 Get Hilbert space size
   matrix mx(hs, hs, 0.0);			 Null Hilbert space matrix
   basis bs = T1s[0].get_basis();		 Current basis of spin tensors
   gen_op nullOp;				 Set up a null operator
   int k=0, m=0;				 Angular momentum indices
   double Jwp;
   complex Taaa;
   gen_op *T2ps;				 Components for specific transitions
   T2ps = new gen_op[2*rank+1];			     (2*l+1 of these for each p)
   for(int a=0; a<hs; a++)			 Sum over transitions a-aa: p
     for(int aa=0; aa<hs; aa++)
       {
       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
         Taaa = T2s[k].get(a,aa);		                       m
         if(Re(Taaa) || Im(Taaa))		           T2ps[k] == T
           {					                       2,p
           mx.put(Taaa, a, aa);			 
           T2ps[k] = gen_op(mx,bs);		 These are either sparse or null arrays
           mx.put(0.0, a, aa);			 so we insure that null components are
           }					 recognized (and not used in computation).
         }
       Jwp = Re(J12.get(a,aa));			 Get spectral density for this transition
       R_CC_0(T1s, T2ps, LOp, rank, Jwp);	 Now, perform sum over m, and add to LOp
       }
   return;
   }
void R_4s ( super_op & LOp , int rank , gen_op * T1s , gen_op * T2s , matrix & J12 ) { int hs = T1s [ 0 ] . dim ( ) ; matrix mx ( hs , hs , 0.0 ) ; basis bs = T1s [ 0 ] . get_basis ( ) ; gen_op nullOp ; int k = 0 , m = 0 ; double Jwp ; complex Taaa ; gen_op * T2ps ; T2ps = new gen_op [ 2 * rank + 1 ] ; for ( int a = 0 ; a < hs ; a ++ ) for ( int aa = 0 ; aa < hs ; aa ++ ) { for ( m = - rank , k = 0 ; m <= rank ; m ++ , k ++ ) { T2ps [ k ] = nullOp ; Taaa = T2s [ k ] . get ( a , aa ) ; if ( Re ( Taaa ) || Im ( Taaa ) ) { mx . put ( Taaa , a , aa ) ; T2ps [ k ] = gen_op ( mx , bs ) ; mx . put ( 0.0 , a , aa ) ; } } Jwp = Re ( J12 . get ( a , aa ) ) ; R_CC_0 ( T1s , T2ps , LOp , rank , Jwp ) ; } return ; }
____________________________________________________________________________
B            RELAXATION SUPEROPERATORS VIA DOUBLE COMMUTATORS
____________________________________________________________________________
____________________________________________________________________________ B RELAXATION SUPEROPERATORS VIA DOUBLE COMMUTATORS ____________________________________________________________________________
// Input		LOp   : Current relaxation superoperator //			rank  : Rank of the interactions involved // 			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 4 relax. superop for interactions //			        1 & 2 added to input superoperator LOp // //				       LOp    = LOp   + R4 //				          out	   in     12 // // Note				T1s, T2s, & LOp assumed in proper bases //                   --- --- //                   \   \       m    m     -m //	     R4    = /   /   (-1)  [ T , [ T   ,  ] ] J  (w ) //	       1,2   --- ---          1     2,p        12  p //		       p   m // where p sums over transitions, m sums over angular momentum components, // and J values include xi1 and xi2 // Get Hilbert space size // Null Hilbert space matrix // Current basis of spin tensors // Set up a null operator // Angular momentum indices // Components for specific transitions //     (2*l+1 of these for each p) // Sum over transitions a-aa: 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 null components are // recognized (and not used in computation). // Get spectral density for this transition // Now, perform sum over m, and add to LOp