src/BWRRelax/relaxRF.cc
// --------------- Spin-Pair with Spin Functions -------------------
Spin-Pair with Spin Functions
Spin - Pair with Spin Functions
void Rrfijkds(super_op& LOp, const sys_dynamic& sys, gen_op& Heff, double* w,
         double Wrflab, 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 rank=2;
matrix theta = sys.thetas();			
matrix phi = sys.phis(); 			
double alphaij, betaij;				
   double xi1, xi2, xi1xi2=0;			
   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;					
   double cutoff = 1e-12;
   double mWrf=0;
   double J[3];
   matrix* J12 = NULL;
   if(abs(level) > 1)                           
     {
     J12 = new matrix[5];
     Heff.eigvals(w);                           
     }
   int ij = 0;					
   for(int 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(Heff);
         }
       for(int k=0; k<ns; k++)			
         {
         xi2 = Re(xi2s.get(k,k));		
         xi1xi2 = xi1*xi2;			
         if(fabs(xi1xi2) > cutoff)		
           {
           EA2 = (A2[k]).PASys_EA();		
           Jcoeffs(c2s, EA2, chi);		
           for(m=-2; m<3; m++)			
             {					
             T2s[m+2] = gen_op(T2[k].component(2,m));
             T2s[m+2].Op_base(Heff);
             }
           if(abs(level) > 1)
             for(m=-rank; m<=rank; m++)		
               {				
               mWrf = -1.0*double(m)*Wrflab;
               J12[m+2]=J_red_shft(w,mWrf,
 		             hs,taus,c1s,c1s,1);
               J12[m+2] *= complex(xi1xi2);	
               }
           else
             {
             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;
             }
           Rrfmumu(LOp,T1s,T2s,J12,J,w,rank,level,0);
           }
         } 					
       }					
   return;
   if(A1==NULL) type=0;				
   }
void Rrfijkds ( super_op & LOp , const sys_dynamic & sys , gen_op & Heff , double * w , double Wrflab , 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 rank = 2 ; matrix theta = sys . thetas ( ) ; matrix phi = sys . phis ( ) ; double alphaij , betaij ; double xi1 , xi2 , xi1xi2 = 0 ; 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 ; double cutoff = 1e-12 ; double mWrf = 0 ; double J [ 3 ] ; matrix * J12 = NULL ; if ( abs ( level ) > 1 ) { J12 = new matrix [ 5 ] ; Heff . eigvals ( w ) ; } int ij = 0 ; for ( int 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 ( Heff ) ; } for ( int k = 0 ; k < ns ; k ++ ) { xi2 = Re ( xi2s . get ( k , k ) ) ; xi1xi2 = xi1 * xi2 ; if ( fabs ( xi1xi2 ) > cutoff ) { EA2 = ( A2 [ k ] ) . PASys_EA ( ) ; Jcoeffs ( c2s , EA2 , chi ) ; for ( m = -2 ; m < 3 ; m ++ ) { T2s [ m + 2 ] = gen_op ( T2 [ k ] . component ( 2 , m ) ) ; T2s [ m + 2 ] . Op_base ( Heff ) ; } if ( abs ( level ) > 1 ) for ( m = - rank ; m <= rank ; m ++ ) { mWrf = -1.0 * double ( m ) * Wrflab ; J12 [ m + 2 ] = J_red_shft ( w , mWrf , hs , taus , c1s , c1s , 1 ) ; J12 [ m + 2 ] *= complex ( xi1xi2 ) ; } else { 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 ; } Rrfmumu ( LOp , T1s , T2s , J12 , J , w , rank , level , 0 ) ; } } } return ; if ( A1 == NULL ) type = 0 ; }
Spin-Pair with Spin Functions
Spin - Pair with Spin Functions
// Input                LOp   : Superoperator // 			sys   : Dynamic spin system //			Heff  : Effective Hamiltonian //			w     : Transition frequencies //			Wrflab: Lab frame rf-field frequency //			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 //	       Mechanism 1 = Spin-Pair; Mechanism 2 = Spin /*                    --- ---             --- --- ---
                      \   \               \   \   \
               LOp += /   /    R        = /   /   /   R
  	              --- ---   mu1,mu2   --- --- ---  ij,k
                      mu1 mu2              i   j   k
*/ // sosi: these three variables are still dipole specific!! // Get dipole theta values (radians)  // Get dipole phi values (radians) // Two Euler angles for ij 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 // Needed for higher level computations // Frequencies w in field rotating frame // Sum over all dipolar pairs // 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*/ // Put spin tensor for i,j into a // vector of operators in basis of Heff //         T1s[m+2] = gen_op((T1[ij].component(2,m)).matrix()); //         T1s[m+2] = gen_op(T1[ij].component(2,m), //                          sys.get_basis()); // Sum over spins k // Get spin k (mu2) interaction constant // xi(mu1,i)*xi(mu2,j) // Only add non-trivial contributions // Get spin k (mu2) space tensor Euler angles // Set the 5 coefficients for k // Put spin tensor for k into a // vector of operators in basis of Heff //             T2s[m+2] = gen_op((T2[k].component(2,m)).matrix()); //             T2s[m+2] = gen_op(T2[k].component(2,m), //                                 sys.get_basis()); // Set up all the reduced spectral // density functions J(w-m*Wrf)  // sosi 3-10-96 why isnt this Q and c1s with c2s?? // Scale by the interaction constants // Need only zero, single, & double // quantum transition frequencies // !! Don't know which to use here!! // Increment spin k (mu2) // Increment pair ij (mu1) // Compiler likes these used