src/BWRRelax/relaxNMR.cc
// --------- Spin-Pair with Spin-Pair Interaction Functions -------------
Spin-Pair with Spin-Pair Interaction Functions
Spin - Pair with Spin - Pair Interaction Functions
void Rijklds(super_op& LOp, const sys_dynamic& sys, gen_op& Ho, double* w,
         matrix& xi1s, matrix& xi2s, space_T* A1, space_T* A2,
         spin_T* T1, spin_T* T2, double* taus, double chi, int type, int level)
   {
int het = sys.heteronuclear();                  
matrix theta = sys.thetas();			
matrix phi = sys.phis(); 			
double alphaij, betaij;			
double alphakl, betakl;			
   double xi1, xi2, xi1xi2;			
   coord EA1, EA2;				
   double c1s[5];				
   double c2s[5];				
   gen_op *T1s;					
   gen_op *T2s;					
   T1s = new gen_op[5];
   T2s = new gen_op[5];
   double w0=0,w1=0,w2=0,wi=0,wj=0;		
   int ns = sys.spins();			
   int hs = sys.HS();				
   int m;					
   int ij = 0;					
   int kl = 0;					
   int i = 0;
   for(i=0; i<ns-1; i++)			
     for(int j=i+1; j<ns; j++)
       {
       xi1 = Re(xi1s.get(i,j));			
alphaij = Re(phi.get(i,j));		
betaij = Re(theta.get(i,j));		
EA1.xyz(alphaij, betaij, 0.0);
         Jcoeffs(c1s, EA1, chi);		
       for(m=-2; m<3; m++)			
         {					
         T1s[m+2] = gen_op(T1[ij].component(2,m));
         T1s[m+2].Op_base(Ho);
         }
       kl = 0;					
       for (int k=0; k<ns-1; k++)		
         for (int l=k+1; l<ns; l++)
           {
           if((ij == kl) && (type >= 0))	
             {
             xi1xi2 = xi1*xi1;			
             if(abs(level) <2)
               {
               w0 = 0.0;			
               w1 = sys.gamma(i)/GAMMA1H; 	
               w1 *= sys.Omega()*1.0e6;
               w2 = 2.0*w1;
                }
             Rmumuds(LOp,T1s,T1s,w,hs,taus,
                c1s,c1s,xi1xi2,w0,w1,w2,level,1,het);
             }
           else if((ij != kl) && (type <= 0))	
             {
             xi2 = Re(xi2s.get(k,l));		
             xi1xi2 = xi1*xi2;			
alphakl = Re(phi.get(k,l));	
betakl = Re(theta.get(k,l));	
EA2.xyz(alphakl, betakl, 0.0);
             Jcoeffs(c2s, EA2, chi);		
             for(m=-2; m<3; m++)		
               {				
               T2s[m+2] = gen_op(T2[kl].
                            component(2,m));
               T2s[m+2].Op_base(Ho);
               }
             if(abs(level) <2)
               {
               wi = sys.gamma(i)/GAMMA1H;	
               wi *= sys.Omega()*1.0e6; 	
               wj = sys.gamma(j)/GAMMA1H;
               wj *= sys.Omega()*1.0e6; 	
               w0 = wi-wj;
               w1 *= wi;
               w2 = wi+wj;
               }
             Rmumuds(LOp,T1s,T2s,w,hs,taus,
               c1s,c2s,xi1xi2,w0,w1,w2,level,0,het);
             }
           kl++;				
           }
       ij++;					
       }
   gen_op Op;					
   if(T1s)
     {
     for(int ii=0; ii<5; ii++)
       T1s[ii] = Op;
     T1s = NULL;
     }
   if(T2s)
     {
     for(int jj=0; jj<5; jj++)
       T2s[jj] = Op;
     T2s = NULL;
     }
   return;
   if(A1 != NULL) i = 0;			
   if(A2 != NULL) i = 0;			
   }
void Rijklds ( super_op & LOp , const sys_dynamic & sys , gen_op & Ho , double * w , matrix & xi1s , matrix & xi2s , space_T * A1 , space_T * A2 , spin_T * T1 , spin_T * T2 , double * taus , double chi , int type , int level ) { int het = sys . heteronuclear ( ) ; matrix theta = sys . thetas ( ) ; matrix phi = sys . phis ( ) ; double alphaij , betaij ; double alphakl , betakl ; double xi1 , xi2 , xi1xi2 ; coord EA1 , EA2 ; double c1s [ 5 ] ; double c2s [ 5 ] ; gen_op * T1s ; gen_op * T2s ; T1s = new gen_op [ 5 ] ; T2s = new gen_op [ 5 ] ; double w0 = 0 , w1 = 0 , w2 = 0 , wi = 0 , wj = 0 ; int ns = sys . spins ( ) ; int hs = sys . HS ( ) ; int m ; int ij = 0 ; int kl = 0 ; int i = 0 ; for ( i = 0 ; i < ns - 1 ; i ++ ) for ( int j = i + 1 ; j < ns ; j ++ ) { xi1 = Re ( xi1s . get ( i , j ) ) ; alphaij = Re ( phi . get ( i , j ) ) ; betaij = Re ( theta . get ( i , j ) ) ; EA1 . xyz ( alphaij , betaij , 0.0 ) ; Jcoeffs ( c1s , EA1 , chi ) ; for ( m = -2 ; m < 3 ; m ++ ) { T1s [ m + 2 ] = gen_op ( T1 [ ij ] . component ( 2 , m ) ) ; T1s [ m + 2 ] . Op_base ( Ho ) ; } kl = 0 ; for ( int k = 0 ; k < ns - 1 ; k ++ ) for ( int l = k + 1 ; l < ns ; l ++ ) { if ( ( ij == kl ) && ( type >= 0 ) ) { xi1xi2 = xi1 * xi1 ; if ( abs ( level ) < 2 ) { w0 = 0.0 ; w1 = sys . gamma ( i ) / GAMMA1H ; w1 *= sys . Omega ( ) * 1.0e6 ; w2 = 2.0 * w1 ; } Rmumuds ( LOp , T1s , T1s , w , hs , taus , c1s , c1s , xi1xi2 , w0 , w1 , w2 , level , 1 , het ) ; } else if ( ( ij != kl ) && ( type <= 0 ) ) { xi2 = Re ( xi2s . get ( k , l ) ) ; xi1xi2 = xi1 * xi2 ; alphakl = Re ( phi . get ( k , l ) ) ; betakl = Re ( theta . get ( k , l ) ) ; EA2 . xyz ( alphakl , betakl , 0.0 ) ; Jcoeffs ( c2s , EA2 , chi ) ; for ( m = -2 ; m < 3 ; m ++ ) { T2s [ m + 2 ] = gen_op ( T2 [ kl ] . component ( 2 , m ) ) ; T2s [ m + 2 ] . Op_base ( Ho ) ; } if ( abs ( level ) < 2 ) { wi = sys . gamma ( i ) / GAMMA1H ; wi *= sys . Omega ( ) * 1.0e6 ; wj = sys . gamma ( j ) / GAMMA1H ; wj *= sys . Omega ( ) * 1.0e6 ; w0 = wi - wj ; w1 *= wi ; w2 = wi + wj ; } Rmumuds ( LOp , T1s , T2s , w , hs , taus , c1s , c2s , xi1xi2 , w0 , w1 , w2 , level , 0 , het ) ; } kl ++ ; } ij ++ ; } gen_op Op ; if ( T1s ) { for ( int ii = 0 ; ii < 5 ; ii ++ ) T1s [ ii ] = Op ; T1s = NULL ; } if ( T2s ) { for ( int jj = 0 ; jj < 5 ; jj ++ ) T2s [ jj ] = Op ; T2s = NULL ; } return ; if ( A1 != NULL ) i = 0 ; if ( A2 != NULL ) i = 0 ; }
Spin-Pair with Spin-Pair Interaction Functions
Spin - Pair with Spin - Pair Interaction Functions
// Input                LOp   : Superoperator // 			sys   : Dynamic spin system //			Ho    : Static Hamiltonian //			w     : Transition frequencies //			xi1s  : Interaction constant, mu1 //			xi2s  : Interaction constant, mu2 //			A1    : Spatial tensor, mu1 //			A2    : Spatial tensor, mu2 //			T1    : Spin tensor, mu1 //			T2    : Spin tensor, mu2 //			taus  : Sytem 5 effective correlation times //			chi   : Spin system chi value // 			type  : Relaxation type to compute //				   0 = Auto & Cross Correlation //				   + = Auto Correlation Only //				   - = Cross Correlation Only //			level : Relaxation treatment level /*	             Two Spin-Pair Rank 2 Mechanisms

                      --- ---             --- --- --- ---
                      \   \               \   \   \   \
               LOp += /   /    R        = /   /   /   /   R
  	              --- ---   mu1,mu2   --- --- --- ---  ij,kl
                      mu1 mu2              i   j   k   l              */ // sosi added 3/11/96 // Flag if system is heteronuclear //int rank=2; // sosi: this is still dipole specific!! // Get dipole theta values (radians)  // Get dipole phi values (radians) // Two Euler angles for ij dipole // Two Euler angles for kl dipole // For specific interaction constants // For specific Euler angles // Set up 5 coefficients interaction 1 // Set up 5 coefficients interaction 2 // Compiler didn't like gen_op T1s[5] // These are for spin tensor components // Used for J's levels 0 & 1 // Total number of spins // Total system Hilbert space //   int ls = hs*hs;				// Total system Liouville space // z-component of ang. momentum // Sum over all dipolar pairs // and build relaxation matrix // Sum spin pairs i&j: dipole ij // Dipolar interaction constant i&j // Alpha Euler angle: dipole ij // Beta Euler angle: dipole kl //       Jcoeffs(c1s, alphaij, betaij, 0.0, chi);	// Set the 5 coefficients for ij // Set the 5 coefficients for ij // Put spin tensor for i,j into a // vector of operators in basis of Ho //         T1s[m+2] = gen_op(T1[ij].component(2,m), //                          sys.get_basis()); // dipole kl count to zero // Sum spin pairs k&l: dipole kl // Auto-correlation term // xi(ij)*xi(ij) // Need only zero, single, & double // quantum transition frequencies // Cross-correlation term // Dipolar interaction constant k&l // sosi - switched above to xi2s from xi1s 1/31/96 // xi(ij)*xi(kl) // Get the alpha Euler angle kl  // Get the beta Euler angle kl //             Jcoeffs(c2s, alphakl, betakl, 	// Set the 5 coefficients for kl //				   0.0,  chi);  // where gammakl = 0 dipolar // Set the 5 coefficients for kl // Put spin tensor for k,l into a // vector of operators in basis of Ho //               T2s[m+2] = gen_op(T2[kl]. //                            component(2,m), //                              sys.get_basis()); // Need only zero, single and double // quantum transition frequencies // !! Don't know which to use here!! // Increment second dipole // Increment first dipole // Psuedo destruction of operators // Compiler likes this used // Compiler likes this used