src/BWRRelax/relaxRF.cc
// ----------------- Level 3 by Double Commutators ----------------------
//                  (Applies Secular Approximation)
Level 3 by Double Commutators
(Applies Secular Approximation)
Level 3 by Double Commutators ( Applies Secular Approximation )
void Rrf_3s(super_op& LOp, double* w, int rank, gen_op* T1s, gen_op* T2s, matrix* J12)
					LOp    = LOp   + LOp
					   out	    in      12
	 Note				T1s, T2s, & LOp assumed in proper bases
   {
 sosi UNDER CONSTRUCTION: THIS  DOESN'T WORK YET
   int hs = T1s[0].dim();			 Get Hilbert space size
   matrix mx1(hs, hs, 0.0);			 Two zero matrices of the
   matrix mx2(hs, hs, 0.0);			 Hilbert space size
   basis bs = T1s[0].get_basis();		 Current Basis
   gen_op nullOp;				 A NULL operator
   double delwaa, delwbb;
   int k=0, kk=0;
   double Jwp;
   complex Ta_aa, Tb_bb;
   gen_op *T1ps;
   T1ps = new gen_op[2*rank+1];
   gen_op *T2ps;
   T2ps = new gen_op[2*rank+1];
double Js[5];
   for(int a=0; a<hs; a++)			 Sum over transitions a-aa: p
     for(int aa=0; aa<hs; aa++)
       {
       delwaa = w[a] - w[aa];			 Transition a-aa: p frequency
       Jwp = Re(J12.get(a,aa));
       k=0;
       for(int m=-rank; m<=rank; m++)		 Sum over all tensor components
         {					 and generate the T1ps eigenoperators
         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);
           }
Js[m+2] = Re(J12[m+2].get(a,aa));
         k++;
         }
       for (int b=0; b<hs; b++)			 Sum over transitions b-bb: p'
         for (int bb=0; bb<hs; bb++)
           {
           delwbb = w[b] - w[bb];		 Transition b-bb: p' frequency
	   if(fabs(delwaa+delwbb) < 1.e-9)	 Secular approximation: p = -p'
             {
             kk=0;
             for (int mm=-rank; mm<=rank; mm++)	 Sum over all tensor components
               {				 and generate the T2ps eigenoperators
               T2ps[kk] = nullOp;
               Tb_bb = T2s[kk].get(b,bb);
               if(Re(Tb_bb) || Im(Tb_bb))
                 {
                 mx2.put(Tb_bb, b, bb);
                 T2ps[kk] = gen_op(mx2,bs);
                 mx2.put(0.0, b, bb);
                 }
               kk++;
               }
 R_CC_0(T1ps, T2ps, LOp, rank, Js);
               R_CC_0(T1ps, T2ps, LOp, rank, Jwp);
 the above step should be reviewed.
             }
           }
       }
   return;
   if(LOp.dim() && J12==NULL) hs=0;		 Compiler likes these used 
   }
void Rrf_3s ( super_op & LOp , double * w , int rank , gen_op * T1s , gen_op * T2s , matrix * J12 ) { 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 delwaa , delwbb ; int k = 0 , kk = 0 ; complex Ta_aa , Tb_bb ; gen_op * T1ps ; T1ps = new gen_op [ 2 * rank + 1 ] ; gen_op * T2ps ; T2ps = new gen_op [ 2 * rank + 1 ] ; double Js [ 5 ] ; for ( int a = 0 ; a < hs ; a ++ ) for ( int aa = 0 ; aa < hs ; aa ++ ) { delwaa = w [ a ] - w [ aa ] ; k = 0 ; for ( int m = - rank ; m <= rank ; m ++ ) { 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 ) ; } Js [ m + 2 ] = Re ( J12 [ m + 2 ] . get ( a , aa ) ) ; k ++ ; } for ( int b = 0 ; b < hs ; b ++ ) for ( int bb = 0 ; bb < hs ; bb ++ ) { delwbb = w [ b ] - w [ bb ] ; if ( fabs ( delwaa + delwbb ) < 1.e-9 ) { kk = 0 ; for ( int mm = - rank ; mm <= rank ; mm ++ ) { T2ps [ kk ] = nullOp ; Tb_bb = T2s [ kk ] . get ( b , bb ) ; if ( Re ( Tb_bb ) || Im ( Tb_bb ) ) { mx2 . put ( Tb_bb , b , bb ) ; T2ps [ kk ] = gen_op ( mx2 , bs ) ; mx2 . put ( 0.0 , b , bb ) ; } kk ++ ; } } } } return ; if ( LOp . dim ( ) && J12 == NULL ) hs = 0 ; }
Level 3 by Double Commutators
(Applies Secular Approximation)
Level 3 by Double Commutators ( Applies Secular Approximation )
// Input		LOp   : Current relaxation superoperator //			w     : Double vector of enery levels //				(rad/sec) for static Hamiltonian //				in the LAB frame (range in 10**8 Hz) //			rank  : Rank of the two interactions // 			T1s   : Spin tensor components, 1st interaction // 			T2s   : Spin tensor components, 2nd interaction //			J12   : Matrices of reduced spectral density  //				function values for 2 interactions //				times the two interaction constants //				evaluated at m based frequencies // Void			LOp   : Relaxation superoperator elements //				for interactions 1 & 2 added to input //				superoperator elements // //				LOp    = LOp   + LOp //				   out	    in      12 // // Note				T1s, T2s, & LOp assumed in proper bases // sosi UNDER CONSTRUCTION: THIS  DOESN'T WORK YET // Get Hilbert space size // Two zero matrices of the // Hilbert space size // Current Basis // A NULL operator //   double Jwp; // Sum over transitions a-aa: p // Transition a-aa: p frequency //       Jwp = Re(J12.get(a,aa)); // Sum over all tensor components // and generate the T1ps eigenoperators // Sum over transitions b-bb: p' // Transition b-bb: p' frequency // Secular approximation: p = -p' // Sum over all tensor components // and generate the T2ps eigenoperators // R_CC_0(T1ps, T2ps, LOp, rank, Js); //               R_CC_0(T1ps, T2ps, LOp, rank, Jwp); // the above step should be reviewed. // Compiler likes these used 