src/BWRRelax/relaxanalyze.cc
// ______________________________________________________________________
// *********** RELAXATION MATRIX WHOLE ELEMENT STRING FUNCTIONS *********
// ______________________________________________________________________
RELAXATION MATRIX WHOLE ELEMENT STRING FUNCTIONS
RELAXATION MATRIX WHOLE ELEMENT STRING FUNCTIONS
void Rijkl_el(const sys_dynamic& sys, gen_op& Ho,
            int rank, spin_T* T1, spin_T* T2, std::string& Mlabel,
         std::string* line1, std::string* line2, std::string* line3, int* signs,
              int& Jterms, int& Lterms,
                              int a, int aa, int b, int bb,
                               int DFS, int Windex, int Sindex)
   {
std::string Llbl = "L";
std::string Jlbl = "J";
   gen_op *T1s;					
   gen_op *T2s;					
   T1s = new gen_op[2*rank+1];
   T2s = new gen_op[2*rank+1];
   int ns = sys.spins();			
   int hs = sys.HS();				
   int m;					
   std::string slabs;
   std::string *Wlbls;				
   Wlbls = new std::string[hs*hs];
   W_labels(Wlbls, sys, Ho, Windex);		
   int ndip = sys.dipoles();			
   int ndpairs=ndip*ndip;			
   int ij,kl;
   int maxt = 2*(2*rank+1)*(hs+1);		
   complex *strsJ;			
   int *trnsJ;			
   int *contJ;				
   std::string *Jlbs;			
   complex *strsL;			
   int *trnsL;			
   int *contL;				
   std::string *Llbs;			
   std::string *spns;			
   strsJ = new complex[ndpairs*maxt];
   trnsJ = new int[ndpairs*maxt];	
   contJ = new int[ndpairs];		
   Jlbs  = new std::string[ndpairs];		
   strsL = new complex[ndpairs*maxt];	
   trnsL = new int[ndpairs*maxt];		
   contL = new int[ndpairs];		
   Llbs  = new std::string[ndpairs];		
   spns  = new std::string[ndpairs];		
   int Jterm0=0, Lterm0=0;			
   int ijkl = 0;				
   ij=0;
   std::string *Slbls;				
   Slbls = new std::string[ns];
   Spin_labels(Slbls, sys, Sindex);		
   for(int i=0; i<ns-1; i++)			
     {
     for(int j=i+1; j<ns; j++)
       {
       for(m=-rank; m<=rank; m++)		
         {					
         T1s[m+rank] = T1[ij].component(rank,m);
         T1s[m+rank].Op_base(Ho);
         }
       kl = 0;					
       for(int k=0; k<ns-1; k++)		
         {
         for(int l=k+1; l<ns; l++)
           {
           slabs = Slbls[i] + Slbls[j]
                 + Slbls[k] + Slbls[l];
           if(ij == kl) 			
             {					
             Rel_12(hs,rank,T1s,T1s, 		
                    Jterms,strsJ,trnsJ, 	
                    Lterms,strsL,trnsL,
                              Wlbls,a,aa,b,bb);
             }
           else if(ij != kl)			
             {
             for(m=-rank; m<=rank; m++)		
               {				
               T2s[m+rank] = T2[kl].component(rank,m);
               T2s[m+rank].Op_base(Ho);
               }				
             Rel_12(hs,rank,T1s,T2s, 		
                    Jterms,strsJ,trnsJ, 	
                    Lterms,strsL,trnsL,
                              Wlbls,a,aa,b,bb);
             }
           contJ[ijkl] = Jterms-Jterm0;		
           Jlbs[ijkl] = Jlbl;			
           contL[ijkl] = Lterms-Lterm0;		
           Llbs[ijkl] = Llbl;			
           spns[ijkl] = slabs;			
           kl++;				
           ijkl++;				
           Jterm0 = Jterms;			
           Lterm0 = Lterms;			
           }
         }
       ij++;					
       }
     }
   if(Jterms && DFS >= 0)
     Rel(0, Jterms, ndpairs, strsJ, trnsJ,	
           Wlbls, contJ, spns, Jlbs, Mlabel, 	
                   line1, line2, line3, signs);
   else
     Jterms = 0;
   if(Lterms && DFS <= 0)
     Rel(Jterms, Lterms, ndpairs, strsL, trnsL,	
           Wlbls, contL, spns, Llbs, Mlabel, 	
                   line1, line2, line3, signs);
   else
     Lterms = 0;
   Rel_clean(T1s, T2s, rank);
   delete [] Wlbls;
   delete [] strsJ;
   delete [] trnsJ;
   delete [] contJ;	
   delete [] Jlbs;	
   delete [] strsL;	
   delete [] trnsL;	
   delete [] contL;	
   delete [] Llbs;		
   delete [] spns;
   delete [] Slbls;
   return;
   }
void Rijkl_el ( const sys_dynamic & sys , gen_op & Ho , int rank , spin_T * T1 , spin_T * T2 , std :: string & Mlabel , std :: string * line1 , std :: string * line2 , std :: string * line3 , int * signs , int & Jterms , int & Lterms , int a , int aa , int b , int bb , int DFS , int Windex , int Sindex ) { std :: string Llbl = " " ; std :: string Jlbl = " " ; gen_op * T1s ; gen_op * T2s ; T1s = new gen_op [ 2 * rank + 1 ] ; T2s = new gen_op [ 2 * rank + 1 ] ; int ns = sys . spins ( ) ; int hs = sys . HS ( ) ; int m ; std :: string slabs ; std :: string * Wlbls ; Wlbls = new std :: string [ hs * hs ] ; W_labels ( Wlbls , sys , Ho , Windex ) ; int ndip = sys . dipoles ( ) ; int ndpairs = ndip * ndip ; int ij , kl ; int maxt = 2 * ( 2 * rank + 1 ) * ( hs + 1 ) ; complex * strsJ ; int * trnsJ ; int * contJ ; std :: string * Jlbs ; complex * strsL ; int * trnsL ; int * contL ; std :: string * Llbs ; std :: string * spns ; strsJ = new complex [ ndpairs * maxt ] ; trnsJ = new int [ ndpairs * maxt ] ; contJ = new int [ ndpairs ] ; Jlbs = new std :: string [ ndpairs ] ; strsL = new complex [ ndpairs * maxt ] ; trnsL = new int [ ndpairs * maxt ] ; contL = new int [ ndpairs ] ; Llbs = new std :: string [ ndpairs ] ; spns = new std :: string [ ndpairs ] ; int Jterm0 = 0 , Lterm0 = 0 ; int ijkl = 0 ; ij = 0 ; std :: string * Slbls ; Slbls = new std :: string [ ns ] ; Spin_labels ( Slbls , sys , Sindex ) ; for ( int i = 0 ; i < ns - 1 ; i ++ ) { for ( int j = i + 1 ; j < ns ; j ++ ) { for ( m = - rank ; m <= rank ; m ++ ) { T1s [ m + rank ] = T1 [ ij ] . component ( rank , m ) ; T1s [ m + rank ] . Op_base ( Ho ) ; } kl = 0 ; for ( int k = 0 ; k < ns - 1 ; k ++ ) { for ( int l = k + 1 ; l < ns ; l ++ ) { slabs = Slbls [ i ] + Slbls [ j ] + Slbls [ k ] + Slbls [ l ] ; if ( ij == kl ) { Rel_12 ( hs , rank , T1s , T1s , Jterms , strsJ , trnsJ , Lterms , strsL , trnsL , Wlbls , a , aa , b , bb ) ; } else if ( ij != kl ) { for ( m = - rank ; m <= rank ; m ++ ) { T2s [ m + rank ] = T2 [ kl ] . component ( rank , m ) ; T2s [ m + rank ] . Op_base ( Ho ) ; } Rel_12 ( hs , rank , T1s , T2s , Jterms , strsJ , trnsJ , Lterms , strsL , trnsL , Wlbls , a , aa , b , bb ) ; } contJ [ ijkl ] = Jterms - Jterm0 ; Jlbs [ ijkl ] = Jlbl ; contL [ ijkl ] = Lterms - Lterm0 ; Llbs [ ijkl ] = Llbl ; spns [ ijkl ] = slabs ; kl ++ ; ijkl ++ ; Jterm0 = Jterms ; Lterm0 = Lterms ; } } ij ++ ; } } if ( Jterms && DFS >= 0 ) Rel ( 0 , Jterms , ndpairs , strsJ , trnsJ , Wlbls , contJ , spns , Jlbs , Mlabel , line1 , line2 , line3 , signs ) ; else Jterms = 0 ; if ( Lterms && DFS <= 0 ) Rel ( Jterms , Lterms , ndpairs , strsL , trnsL , Wlbls , contL , spns , Llbs , Mlabel , line1 , line2 , line3 , signs ) ; else Lterms = 0 ; Rel_clean ( T1s , T2s , rank ) ; delete [ ] Wlbls ; delete [ ] strsJ ; delete [ ] trnsJ ; delete [ ] contJ ; delete [ ] Jlbs ; delete [ ] strsL ; delete [ ] trnsL ; delete [ ] contL ; delete [ ] Llbs ; delete [ ] spns ; delete [ ] Slbls ; return ; }
______________________________________________________________________
 RELAXATION MATRIX WHOLE ELEMENT STRING FUNCTIONS
______________________________________________________________________
______________________________________________________________________ RELAXATION MATRIX WHOLE ELEMENT STRING FUNCTIONS ______________________________________________________________________
// Input 		sys   : Dynamic spin system //			Ho    : Static Hamiltonian //			rank  : Rank of the interactions //			T1    : Spin tensor, mu1 //			T2    : Spin tensor, mu2 //			Mlabel: Label for two mechanisms //			line1 : strings for the upper line of component //			line2 : strings for the middle line of component //			line3 : strings for the lower line of component //			signs : Flags for the signs of the components //			a, aa : 1st transition indices //			b, bb : 2nd transition indices //			DFS   : Flag for dynamic frequency shifts //					> 0  J Terms Only				     //					  0  Both J and L Terms (Default) //					< 0  L Terms Only //			Windex: Flag for frequency indexing //                                 >1 = Always use spin indices //                                  0 = Use spin indices if needed //                                 -1 = Never use spin indices (default) //                                 -2 = Assume extreme narrowing //                                 -3 = Assume extreme broadening //			Sindex: Flag for spin indexing //                                  0 = Use only if needed (default) //                                  1 = Always use spin indices //                                 <1 = Never use spin indices // Output		void  : // Note			      : This routine assumes that each of //				the two interactions involve 2 spins // Note			      : The above implies the rank is 2 //                                --- --- --- --- //                    12          \   \   \   \    12 //             <a,a'|R  |b,b'> =  /   /   /   /   R //	                          --- --- --- ---  ij,kl //                                 i   j   k   l // Compiler didn't like gen_op T1s[5] // These are for spin tensor components // Total number of spins // Total system Hilbert space // z-component of ang. momentum // Array of transition labels // Fill up the transition labels // Total number of dipoles // Total number of dipole pairs // Max. contributions per dip. pair // Array for interaction strengths (J's) // Array for transition indices (J's) // Number of contributors per dipole pair // Array for J labels (J, can use L,K, etc.) // Array for interaction strengths (L's) // Array for transition indices (L's) // Number of contributors per dipole pair // Array for L labels (L, can use L,Q, etc.) // Array for spin labels (ijkl, ijjk, etc.) // Array for transition indices (J's) // Number of contributors per dipole pair // Array for J labels (J, can use L,K, etc.) // Array for interaction strengths (L's) // Array for transition indices (L's) // Number of contributors per dipole pair // Array for L labels (L, can use L,Q, etc.) // Array for spin labels (ijkl, ijjk, etc.) // Total terms contributing to <a,a'|R|b,b'> // Counter for dipole-dipole interactions // strings for spin labels // Fill up spin labels // Sum spin pairs i&j: dipole ij // Put spin tensor for i,j into a // vector of operators in basis of Ho // dipole kl count to zero // Sum spin pairs k&l: dipole kl // Auto-correlation term //                      12 // Get terms for <a,a'|R    |b,b'> //			ijkl // Cross-correlation term // Put spin tensor for k,l into a // vector of operators in basis of Ho //                      12 // Get terms for <a,a'|R    |b,b'> //			ijkl // Number of terms this dipole pair // Store J label for these terms // Number of terms this dipole pair // Store L label for these terms // Store spin label for these terms // Increment second dipole // Increment total dipole pair count // Update the starting J terms counter // Update the starting L terms counter // Increment first dipole // Put J(w) terms into string arrays // line[1-3], and track signs // Add L(w) terms to string arrays // line[1-3], track signs also