propagator_exp()

template<typename T >

void cntr::propagator_exp	(	int	tstp,
		cntr::function< T > &	U,
		cntr::function< T > &	H,
		double	dt,
		int	order,
		int	kt,
		bool	fixHam = false
	)

Propagator for time-dependent free Hamiltonian

Purpose

Calculate the free propagator at time tstp from time dependent free Hamiltonian using high-order commutator-free exponential time-propagation,

see https://doi.org/10.1016/j.jcp.2011.04.006 for the description. Currently implemented versions are the second order using one exponential CF2:1 (order=2) and fourth order using two exponentials CF4:2 (order=4), see also article for more details.

Parameters

tstp	timestep at which propagator is calculated
U	propagator as an output of the function
H	Time dependent Hamiltonian
dt	time step
order	Order of approximation for commutator-free exponential, currently implemented orders = 2,4
kt	Order of integrator used for extrapolation and interpolation
fixHam	Hamiltonian is known for all times and no extrapolation is needed for the predictor/corrector

Definition at line 745 of file cntr_equilibrium_impl.hpp.

References cntr::function< T >::get_value(), interpolate_CF2(), interpolate_CF4(), cntr::function< T >::nt_, cntr::function< T >::set_value(), and cntr::function< T >::size1_.

                                                                                                                {
  int nt=U.nt_;
  int size=U.size1_;
  cdmatrix prop(size,size);

  assert(tstp<=nt);
  assert(order==2 || order==4);
  assert(size==H.size1_);

  if(tstp==-1 || tstp==0){
    prop.setIdentity();
    U.set_value(tstp,prop);
  }else{
    if(order==2){
        cdmatrix arg(size,size);
        // Get H(t+dt/2)-> Extrapolate and interpolate
        interpolate_CF2(tstp,H,arg,kt,fixHam);
        arg=arg*std::complex<double>(0.0,-1.0)*dt;
        U.get_value(tstp-1,prop);
        prop=arg.exp()*prop;
        U.set_value(tstp,prop);
    }else if(order==4){
       cdmatrix H1(size,size),H2(size,size);
       cdmatrix arg1(size,size),arg2(size,size);
       // Get H(t+dt*c1) and H(t+dt*c2) -> Extrapolate and interpolate
       interpolate_CF4(tstp,H,H1,H2,kt,fixHam);
       double a1=(3.0-2.0*sqrt(3.0))/12.0;
       double a2=(3.0+2.0*sqrt(3.0))/12.0;
       arg1=std::complex<double>(0.0,-1.0)*dt*(a1*H1+a2*H2);
       arg2=std::complex<double>(0.0,-1.0)*dt*(a2*H1+a1*H2);
       U.get_value(tstp-1,prop);
       prop=arg1.exp()*arg2.exp()*prop;
       U.set_value(tstp,prop);
    }
  }
}

Here is the call graph for this function: