traject.cc

使用量子轨道方法计算量子主方程的C++库

// interval.  derivs is the routine for calculating the right-hand side of
// the QSD equation, while rkqs is the name of the stepper routine
// to be used.
//
// Adapted from Numerical Recipes routine ODEINT.C
//
{
  int nstp;
  double t,hnext,hdid,h;

  t = t1;
  h = h1; if( h > (t2-t1) ) h = (t2-t1);
  nok = nbad = 0;
  y = ystart;
  for (nstp=1;nstp<=MAXSTP;nstp++) {
	if ((t+h-t2)*(t+h-t1) > 0.0) h=t2-t;
	rkqs(t,h,eps,hdid,hnext);		// deterministic part
	if (hdid == h) ++nok; else ++nbad;
	y.normalize();				// renormalize
	if ((t-t2)*(t2-t1) >= 0.0) {
		ystart = y;
		h1 = hnext;
		if( listPtr < listDim ) {	// if keeping track...
		  dtNumList[listPtr] = nok+nbad;  // ...store number of steps
		  listPtr++;
		}
		return;
	}
	if (fabs(hnext) <= hmin) {
	  cerr << "Step size too small in odeint" << endl;
	  exit(1);
	}
	h=hnext;
  }
  cerr << "Too many steps in routine odeint" << endl;
  exit(1);
}

void AdaptiveStochStep::odeint(State& ystart, double t1, double t2, double eps, double& h1,
	double hmin, int& nok, int& nbad, ComplexRandom* rndm) 
//
// Runge-Kutta driver with adaptive stepsize control  Integrate starting values
// ystart from t1 to t2 with accuracy eps.  h1 should be set as a guessed
// first stepsize, hmin as the minimum allowed stepsize (can be zero).  On
// output nok and nbad are the number of good and bad (but retried and fixed)
// steps taken, and ystart is replaced by values at the end of the integration
// interval.  derivs is the routine for calculating the right-hand side of
// the QSD equation, while rkqs is the name of the stepper routine
// to be used.
//
// Also, generate one stochastic timestep per deterministic timestep;
//
// Adapted from Numerical Recipes routine ODEINT.C
//
{
  int nstp;
  double t,hnext,hdid,h;
  double sqrtdt;

  t = t1;
  h = h1; if( h > (t2-t1) ) h = (t2-t1);
  nok = nbad = 0;
  y = ystart;
  for (nstp=1;nstp<=MAXSTP;nstp++) {
    if ((t+h-t2)*(t+h-t1) > 0.0) h=t2-t;
    if( rndm != 0 ) {				// any noise?
      stochasticFlag = 1;
	  // set flag for derivs to generate stoch. part and store in newsum
      for (int i=0; i<nL; i++)			// for each Lindblad...
        dxi[i] = (*rndm)();			// ...independent noise
    }
    newsum = 0;
    rkqs(t,h,eps,hdid,hnext);			// deterministic part
    if (hdid == h) ++nok; else ++nbad;
    if( rndm != 0 ) {
      sqrtdt = sqrt(hdid);
      newsum *= sqrtdt;
      y += newsum;
    }
    y.normalize();				// renormalize
    if ((t-t2)*(t2-t1) >= 0.0) {
      ystart = y;
      h1 = hnext;
      if( listPtr < listDim ) {			// if keeping track...
        dtNumList[listPtr] = nok+nbad;		// ...store number of steps
        listPtr++;
      }
      return;
    }
    if (fabs(hnext) <= hmin) {
      cerr << "Step size too small in odeint" << endl;
      exit(1);
    }
    h=hnext;
  }
  cerr << "Too many steps in routine odeint" << endl;
  exit(1);
}


void AdaptiveJump::odeint(State& ystart, double t1, double t2, double eps, double& h1,
	double hmin, int& nok, int& nbad, ComplexRandom* rndm) 
//
// Runge-Kutta driver with adaptive stepsize control  Integrate starting values
// ystart from t1 to t2 with accuracy eps.  h1 should be set as a guessed
// first stepsize, hmin as the minimum allowed stepsize (can be zero).  On
// output nok and nbad are the number of good and bad (but retried and fixed)
// steps taken, and ystart is replaced by values at the end of the integration
// interval.  derivs is the routine for calculating the right-hand side of
// the QSD equation, while rkqs is the name of the stepper routine
// to be used.
//
// Also, generate one stochastic timestep per deterministic timestep;
//
// Adapted from Numerical Recipes routine ODEINT.C
//
{
  int nstp,count;
  double t,hnext,hdid,h;
  double prob,prob2,total;

  t = t1;
  h = h1; if( h > (t2-t1) ) h = (t2-t1);
  nok = nbad = 0;
  y = ystart;
  for (nstp=1;nstp<=MAXSTP;nstp++) {
    if ((t+h-t2)*(t+h-t1) > 0.0) h=t2-t;
    if( rndm != 0 ) {				// any noise?
      stochasticFlag = 1;			// have derivs calc <Ldag L>
    }
    newsum = y;                                 // keep the actual value of psi
    rkqs(t,h,eps,hdid,hnext);			// deterministic part
    if (hdid == h) ++nok; else ++nbad;
    if( rndm != 0 ) {				// any noise?
      total = 0;
      for(count=0; count<nL; count++)
	total += lindVal[count];		// sum up <Ldag L>
// Two different choice for the test-jump
// they should give the same results
      prob = 1.0 - real(y*y);
//      prob = total*hdid;
//
      if( real((*rndm)()) < prob ) {		// make a jump?
	y = newsum;
        prob = real((*rndm)());
	for(count=0; prob>0; count++) {
          prob2 = lindVal[count]/total;
          if( prob < prob2 ) {			// jump !
            y *= L[count];			// L|psy>
            if( outFile != 0 )			// output jump time to file
              fprintf(outFile,"%lf  %d\n",t,count);
          }
          prob -= prob2;
        }
	if (t2-t > hnext) hnext = t2-t;  // After a jump the step h
	                                 // is re-initialized 
      }
    }
    y.normalize();				// renormalize
    if ((t-t2)*(t2-t1) >= 0.0) {
      ystart = y;
      h1 = hnext;
      if( listPtr < listDim ) {			// if keeping track...
        dtNumList[listPtr] = nok+nbad;		// ...store number of steps
        listPtr++;
      }
      return;
    }
    if (fabs(hnext) <= hmin) {
      cerr << "Step size too small in odeint" << endl;
      exit(1);
    }
    h=hnext;
  }
  cerr << "Too many steps in routine odeint" << endl;
  exit(1);
}

void AdaptiveOrthoJump::odeint(State& ystart, double t1, double t2, double eps, double& h1,
	double hmin, int& nok, int& nbad, ComplexRandom* rndm) 
//
// Runge-Kutta driver with adaptive stepsize control  Integrate starting values
// ystart from t1 to t2 with accuracy eps.  h1 should be set as a guessed
// first stepsize, hmin as the minimum allowed stepsize (can be zero).  On
// output nok and nbad are the number of good and bad (but retried and fixed)
// steps taken, and ystart is replaced by values at the end of the integration
// interval.  derivs is the routine for calculating the right-hand side of
// the orthogonal jump equation, while rkqs is the name of the stepper routine
// to be used.
//
// Also, generate one stochastic timestep per deterministic timestep;
//
// Adapted from Numerical Recipes routine ODEINT.C
//
{
  int nstp,count;
  double t,hnext,hdid,h;
  double prob,prob2,total;
  Complex expect;

  t = t1;
  h = h1; if( h > (t2-t1) ) h = (t2-t1);
  nok = nbad = 0;
  y = ystart;
  for (nstp=1;nstp<=MAXSTP;nstp++) {
    if ((t+h-t2)*(t+h-t1) > 0.0) h=t2-t;
    if( rndm != 0 ) {				// any noise?
      stochasticFlag = 1;			// have derivs calc <Ldag L>
    }
    newsum = y;
    rkqs(t,h,eps,hdid,hnext);			// deterministic part
    if (hdid == h) ++nok; else ++nbad;
    if( rndm != 0 ) {				// any noise?
      total = 0;
      for(count=0; count<nL; count++)
	total += lindVal[count];		// sum up <Ldag L>-<Ldag><L>
      prob = total*hdid;
      if( real((*rndm)()) < prob ) {		// make a jump?
	y = newsum;
        prob = real((*rndm)());
        for(count=0; prob>0; count++) {
          prob2 = lindVal[count]/total;
          if( prob < prob2 ) {		        // jump!
            y *= L[count];	                // L |psi>
	    expect = newsum*y;                  // <L>
	    newsum *= expect;                   // <L>|psi>
	    y -= newsum;                        // L|psi> - <L>|psi>
            if( outFile != 0 )			// output jump time to file
              fprintf(outFile,"%lf  %d\n",t,count);
	  }
          prob -= prob2;
        }
	if (t2-t > hnext) hnext = t2-t;  // After a jump the step h
	                                 // is re-initialized 
      }
    }
    y.normalize();				// renormalize
    if ((t-t2)*(t2-t1) >= 0.0) {
      ystart = y;
      h1 = hnext;
      if( listPtr < listDim ) {			// if keeping track...
        dtNumList[listPtr] = nok+nbad;		// ...store number of steps
        listPtr++;
      }
      return;
    }
    if (fabs(hnext) <= hmin) {
      cerr << "Step size too small in odeint" << endl;
      exit(1);
    }
    h=hnext;
  }
  cerr << "Too many steps in routine odeint" << endl;
  exit(1);
}

void AdaptiveStep::dtListSet( int theDim )
//
// Sets up a list to store the number of deterministic steps per
// stochastic timestep
//
{
#ifndef OPTIMIZE_QSD
  if( theDim < 1 )
    error("Negative dimension specified in AdaptiveStep::dtListSet!");
#endif
  listDim = theDim;
  dtNumList = new int[theDim];
}

void AdaptiveStep::dtListRead()
//
// Prints out the number of deterministic steps per stochastic timestep
//
{
  for( int i=0; i<listPtr; i++ )
    cout << dtNumList[i] << endl;
}

int AdaptiveStep::dtListElem( int theElem )
//
// Prints out the number of deterministic steps in a particular
// stochastic timestep
//
{
#ifndef OPTIMIZE_QSD
  if( (theElem < 0) || (theElem >= listPtr) )
    error("Illegal list element in AdaptiveStep::dtListElem!");
#endif
  return dtNumList[theElem];
}

void AdaptiveStep::dtListClear()
//
// Clear the list of numbers of deterministic steps per stochastic timestep
//
{
  if( listDim != 0 )
#ifndef NON_GNU_DELETE
    delete[] dtNumList;			// deallocate memory
#else
    delete[listDim] dtNumList;		// deallocate memory
#endif
  listDim = 0;
  listPtr = 0;
  dtNumList = 0;
}

void AdaptiveStep::dtListReset()
//
// Reset the array of numbers without deallocating it
//
{
  listPtr = 0;
}

void IntegrationStep::derivs(double t, State& psi, State& dpsi)
//
// Computes the change of the state psi in a timestep dt due to the
// deterministic part of the QSD equation.
//
// Also, when stochasticFlag is set, computes the sum for the stochastic part.
//
{
  Complex expect;
  double sum=0;
  Complex stochSum=0;

  dpsi = psi;
  dpsi *= H(t);
  dpsi *= M_IM;
  for (int count=0; count<nL; count++) {
    temp1 = psi;
    temp1 *= L[count];
    expect = psi*temp1;
    if( stochasticFlag != 0 ) {
      temp0 = temp1;
      temp0 *= dxi[count];
      newsum += temp0;
      stochSum -= expect*dxi[count];
    }
    sum -= norm(expect);
    temp0 = temp1;
    temp0 *= conj(expect);
    dpsi += temp0;
    temp1 *= Ldag[count];
    temp1 *= -0.5;
    dpsi += temp1;
  }
  sum *= 0.5;
  temp0 = psi;
  temp0 *= sum;
  dpsi += temp0;
  if( stochasticFlag != 0 ) {
    temp0 = psi;
    temp0 *= stochSum;
    newsum += temp0;
  }
  stochasticFlag = 0;
}

void AdaptiveJump::derivs(double t, State& psi, State& dpsi)
//
// Computes the change of the state psi in a timestep dt due to the
// deterministic part of the quantum jumps evolution.
//
// Also, when stochasticFlag is set, computes the sum for the stochastic part.
//
{
  dpsi = psi;
  dpsi *= H(t);
  dpsi *= M_IM;                                     // -iH|psi>
  for (int count=0; count<nL; count++) {
    temp1 = psi;
    temp1 *= L[count];                              // L|psi>
    if (stochasticFlag == 1)
      lindVal[count] = real(temp1*temp1);           // <Ldag L>
    temp1 *= Ldag[count];                           // Ldag L|psi>
    temp1 *= -0.5;
    dpsi += temp1;                                  // -0.5 Ldag L|psi>  
  }
  stochasticFlag = 0;
}

void AdaptiveOrthoJump::derivs(double t, State& psi, State& dpsi)
//
// Computes the change of the state psi in a timestep dt due to the
// deterministic part of the orthogonal jump evolution.
//
// Also, when stochasticFlag is set, computes the sum for the stochastic part.
//
{
  Complex expect;

  dpsi = psi;
  dpsi *= H(t);
  dpsi *= M_IM;                 // dpsi = -iH|psi>
  for (int count=0; count<nL; count++) {
    temp1 = psi;
    temp1 *= L[count];          // L|psi>
    expect = psi*temp1;         // <L>=<psi|L|psi>
    if (stochasticFlag == 1)
      lindVal[count] = real(temp1*temp1)-norm(expect);  // <LdagL>-<Ldag><L>
    temp0 = temp1;              
    temp0 *= conj(expect);      // <Ldag>L|psi>
    dpsi += temp0;              // dpsi += <Ldag> L |psi>
    temp1 *= Ldag[count];       // Ldag L |psi>
    temp1 *= -0.5;
    dpsi += temp1;              // dpsi += -0.5 Ldag L |psi> 
  }
  stochasticFlag = 0;
}

void Trajectory::error(char* message)	// print error message and exit
{
  cerr << message << endl;
  exit(1);
}

void IntegrationStep::error(char* message)	// print error message
{						// and exit
  cerr << message << endl;
  exit(1);
}

#undef MAXSTP
#undef TINY
#undef SAFETY
#undef PGROW
#undef PSHRNK
#undef ERRCON