anchor.diff:  		do t_counter = 1,size(time_dependence,2)
   			propagators(:,:,t_counter,p_counter) = 0
   			do op_counter = 1,size(time_dependence,1)
  -				propagators(:,:,t_counter,p_counter) = &
  -					propagators(:,:,t_counter,p_counter) &
  -					+ p(op_counter,p_counter) &
  -					* time_dependence(op_counter,t_counter) &
  -					* operators(:,:,op_counter)
  +				propagators(:,:,t_counter,p_counter) =&
  +                    propagators(:,:,t_counter,p_counter) +&
  +                    p(op_counter,p_counter) *&
  +                    time_dependence(op_counter,t_counter) *&
  +                    operators(:,:,op_counter)
   			end do
   			call complexexp(propagators(:,:,t_counter,p_counter)) ! this passes a HERMITIAN argument, and gets back the exponent of the 2*pi*i* that argument
   			!call exact(propagators(:,:,t_counter,p_counter)) ! this passes a HERMITIAN argument, and gets back the exponent
   			! for now, i'm not actually propagating, just exponentiating
   			if (t_counter.gt.1) then
  -				propagators(:,:,t_counter,p_counter) = &
  -					matmul(&
  -					propagators(:,:,t_counter,p_counter),&
  -					propagators(:,:,t_counter-1,p_counter))
  +				propagators(:,:,t_counter,p_counter) = matmul( propagators(:,:,t_counter,p_counter), propagators(:,:,t_counter-1,p_counter))
   			endif
   		enddo
   		!write(*,*) 'done propagating parameter ',p_counter,' out of ',size(p,2)
  
positive.code: 	t_steps = size(propagators,3)
  	do p_counter = 1,size(propagators,4)
  		if (size(C,3).gt.1) then
  			call sandwich_internal(&
  				conjg(transpose(C(:,:,p_counter))),&
  				!propagators(:,:,size(propagators,3),p_counter),&
  				conjg(transpose(propagators(:,:,size(propagators,3),p_counter))),&
  				end_sandwich)
  		else
  			call sandwich_internal(&
  				conjg(transpose(C(:,:,1))),&
  				!propagators(:,:,size(propagators,3),p_counter),&
  				conjg(transpose(propagators(:,:,size(propagators,3),p_counter))),&
  				end_sandwich)
  		endif
  		do t_counter = 1,t_steps
  			call sandwich_internal(&
  				H_c,&
  				!propagators(:,:,t_counter,p_counter),&
  				conjg(transpose(propagators(:,:,t_counter,p_counter))),&
  				tempmat)
  			if (size(rho,3).gt.1) then
  				call commutator_internal(tempmat,rho(:,:,1,p_counter))
  			else
  				call commutator_internal(tempmat,rho(:,:,1,1))
  			endif
  			tempmat = matmul(end_sandwich,tempmat)
  			call trace_internal(tempmat,&
  				gradient(t_counter,p_counter))
  		enddo
  	enddo
  	gradient = gradient * (0,1)
  
negative.code: integer p_counter,t_counter
  double complex tempmat(size(propagators,1),size(propagators,2)) 
  write(*,*) '------------------------------------------------------------'
  write (*,*) 'sandwiching rho of size ',size(rho,1),size(rho,2),size(rho,3),size(rho,4),&
      ' with propagator of size ',size(propagators,1),size(propagators,2),&
      size(propagators,3),size(propagators,4)
  if (allocated(propagators)) then
  	write (*,*) 'propagators is allocated for size ',size(propagators)
  else
  