src/LSLib/LSprop.cc
// Input		eLt	: Exponential Liouvillian to relax
//				  the density matrix
//			sigmaeq : Density matrix at equilibrium
//				  (or steady state density matrix)
//			Ho      : Operator in Hilbert space
// Output		LOp     : Superoperator which fullfills
//
//	   eLt * {sigma-sigmaeq} + sigmaeq  = LOp * sigma
//
// Note				: LOp constructed in the current basis of
//				  sigmaeq & in the default Liouville basis
// Note				: sigmaeq is copied so trace is unchanged
Input          eLt     : Exponential Liouvillian to relax
the density matrix
sigmaeq : Density matrix at equilibrium
(or steady state density matrix)
Ho      : Operator in Hilbert space
Output         LOp     : Superoperator which fullfills

eLt * {sigma-sigmaeq} + sigmaeq  = LOp * sigma

Note                           : LOp constructed in the current basis of
sigmaeq & in the default Liouville basis
Note                           : sigmaeq is copied so trace is unchanged
Input eLt : Exponential Liouvillian to relax the density matrix sigmaeq : Density matrix at equilibrium ( or steady state density matrix ) Ho : Operator in Hilbert space Output LOp : Superoperator which fullfills eLt * { sigma - sigmaeq } + sigmaeq = LOp * sigma Note : LOp constructed in the current basis of sigmaeq & in the default Liouville basis Note : sigmaeq is copied so trace is unchanged
super_op R_prop(super_op& eLt, gen_op sigmaeq)
  {
  eLt.set_HBR();				
  super_op LOp(eLt);				
  eLt.LOp_base(sigmaeq);			
  set_trace(sigmaeq, 1.0);			
  int hs = sigmaeq.size();			
  int a,aa,b,bb,g,gg,aaa=0,bbb=0,ggg=0;
  complex Rel;
  int nd = sigmaeq.in_EBR();			
  for(b=0; b<hs; b++)				
    for(bb=0; bb<hs; bb++)
      {
      if(b==bb) 				
        {					
        aaa = 0;		
        for(a=0; a<hs; a++)
          for(aa=0; aa<hs; aa++)
            {
            Rel = LOp.get(aaa,bbb);
            if(!nd || a==aa)			
              Rel += sigmaeq.get(a,aa);		
            ggg=0;
            for(g=0; g<hs; g++)
              for(gg=0; gg<hs; gg++)
                {
                Rel -= eLt.get(aaa,ggg)*sigmaeq.get(g,gg);
                ggg++;
                }
            LOp.put(aaa,bbb,Rel);
            aaa++;
            }
        }
      bbb++;
      }
  return LOp;
  }
super_op R_prop ( super_op & eLt , gen_op sigmaeq ) { eLt . set_HBR ( ) ; super_op LOp ( eLt ) ; eLt . LOp_base ( sigmaeq ) ; set_trace ( sigmaeq , 1.0 ) ; int hs = sigmaeq . size ( ) ; int a , aa , b , bb , g , gg , aaa = 0 , bbb = 0 , ggg = 0 ; complex Rel ; int nd = sigmaeq . in_EBR ( ) ; for ( b = 0 ; b < hs ; b ++ ) for ( bb = 0 ; bb < hs ; bb ++ ) { if ( b == bb ) { aaa = 0 ; for ( a = 0 ; a < hs ; a ++ ) for ( aa = 0 ; aa < hs ; aa ++ ) { Rel = LOp . get ( aaa , bbb ) ; if ( ! nd || a == aa ) Rel += sigmaeq . get ( a , aa ) ; ggg = 0 ; for ( g = 0 ; g < hs ; g ++ ) for ( gg = 0 ; gg < hs ; gg ++ ) { Rel -= eLt . get ( aaa , ggg ) * sigmaeq . get ( g , gg ) ; ggg ++ ; } LOp . put ( aaa , bbb , Rel ) ; aaa ++ ; } } bbb ++ ; } return LOp ; }
Input          eLt     : Exponential Liouvillian to relax
the density matrix
sigmaeq : Density matrix at equilibrium
(or steady state density matrix)
Ho      : Operator in Hilbert space
Output         LOp     : Superoperator which fullfills
Input eLt : Exponential Liouvillian to relax the density matrix sigmaeq : Density matrix at equilibrium ( or steady state density matrix ) Ho : Operator in Hilbert space Output LOp : Superoperator which fullfills
// Put eLt in default Liouville basis // Copy the original superoperator // Put sigma into eLt Hilbert basis // Set the trace to 1 // Get the Hilbert space size // Test if sigmaeq is diagonal // *should add a check for i-matrix so prop is i-matrix returned //  if(mx.stored_type() != i_matrix_type) // *actually need to test (1-eLt)*sigmaeq and see if it is diagonal, not sigmaeq! // should set nd=0; // For nd=1, diagonal  // <aaa|LOp|bbb> = <aaa|eLt|bbb> if |bb> != |b> // so do not make any modifications // If non-diagonal sigmaeq or |a> = |aa> must // add <aaa|sigmaeq|1> to each <aaa|LOp|bbb>