traject.cc

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

        int move, double delta, int width, double moveEps)
//
// This member function calculates nTrajectory trajectories consisting of 
// nOfSteps output steps, and at each step stores the expectation values 
// of the nX operators X.
// Every nTrajSave trajectory, the mean is computed and saved into files.
// If ReadFile=1 then the files `fname[i]' are read at the beginig of the 
// computation and their values take into account for the mean computation.  
// The variable move contains the number of freedoms for which the moving 
// basis and changing cutoff algorithms should be used;
// delta is the probability threshold for the cutoff adjustment, and
// width is the size of the "pad" for the cutoff.  
//
{
  double t=t0;				// time
  double dtlast=dt;			// keep track of numerical timestep
  State psi0 = psi;                     // initial state
  State psi1 = psi;		        // temporary state
  int nTrajEff = 0;                     // number of trajectories
  FILE** fp;
  fp = new FILE*[nX];                   // nX files are used
  int i,j,k,n;

  // Intialisation of the table OutputExp 
  double** OutputExpR;                  // table of real expectations values
  double** OutputExpI;                  // table of imag expectations values
  OutputExpR = new double*[nOfSteps+1];
  OutputExpI = new double*[nOfSteps+1];
  for ( j=0; j<=nOfSteps; j++) OutputExpR[j] = new double[2*nX];
  for ( j=0; j<=nOfSteps; j++) OutputExpI[j] = new double[2*nX];
  for ( n=0; n<=nOfSteps; n++)
    for ( i=0; i<2*nX; i++) OutputExpR[n][i] = 0.0;
  for ( n=0; n<=nOfSteps; n++)
    for ( i=0; i<2*nX; i++) OutputExpI[n][i] = 0.0;

  //
  // Read the previous results
  //
  if (ReadFile) {
    // open the files in read mode
    for ( i=0; i<nX; i++) fp[i] = fopen(fname[i],"r");
    // check for the files open
    int OpenAllFiles = 1;
    for ( i=0; i<nX; i++) if (fp[i]==NULL) OpenAllFiles=0;
    if (OpenAllFiles) {
      // number of trajectories already computed
      char dummy[100];
      for ( i=0; i<nX; i++) 
	fscanf( fp[i], "%s %d ", dummy, &nTrajEff );
      // read the preceding expectations values
      for ( n=0; n<=nOfSteps; n++)
	for ( i=0; i<nX; i++)
	  fscanf( fp[i], "%lf %lf %lf %lf %lf", 
		 &OutputExpR[n][i], &OutputExpI[n][i], 
		 &OutputExpR[n][nX+i], &OutputExpI[n][nX+i]);
      // close the files
      for ( i=0; i<nX; i++) fclose(fp[i]);
      // transform mean values  
      for ( n=0; n<=nOfSteps; n++)
	for ( i=0; i<2*nX; i++) OutputExpR[n][i] *= nTrajEff;
      for ( n=0; n<=nOfSteps; n++)
	for ( i=0; i<2*nX; i++) OutputExpI[n][i] *= nTrajEff;
      cout << " Using " << nTrajEff << " previous trajectories " << endl;
    }
  }
  //
  // Computation of the expectation values 
  //
  for( k=0; k<nTrajectory; k++ ) {
    cout << "Computing trajectory " << k+1  << endl;
    ++nTrajEff;
    t = t0;
    psi = psi0;

    n=0; 
    for( i=0; i<nX; i++ ) {		// compute expectation values for
      psi1 = psi;			// output...
      psi1 *= X[i];
      Complex expec = psi*psi1;
      psi1 *= X[i];
      Complex expec2 = psi*psi1 - expec*expec;
      OutputExpR[n][i] += real(expec);
      OutputExpI[n][i] += imag(expec);
      OutputExpR[n][nX+i] += real(expec2);
      OutputExpI[n][nX+i] += imag(expec2);
    }

    for ( n=1; n<=nOfSteps; n++) {
      for (int n1=0; n1<dtsPerStep; n1++) {
        (*stepper)(psi,t,dt,dtlast,rndm);
        t += dt;
        for( j=0; j<move; j++) {		// use moving basis & cutoff
          psi.adjustCutoff(j,delta,width);
          psi.recenter(j,moveEps);
        }
      }
      for( i=0; i<nX; i++ ) {		// compute expectation values for
        psi1 = psi;			// output...
        psi1 *= X[i];
        Complex expec = psi*psi1;
        psi1 *= X[i];
        Complex expec2 = psi*psi1 - expec*expec;
        OutputExpR[n][i] += real(expec);
        OutputExpI[n][i] += imag(expec);
        OutputExpR[n][nX+i] += real(expec2);
        OutputExpI[n][nX+i] += imag(expec2);
      }
    }
    //
    // Temporary save of results every nTrajSave 
    //
    if ((nTrajSave>0) && ((k+1) % nTrajSave)==0) {
      for ( i=0; i<nX; i++) fp[i] = fopen(fname[i],"w"); 
      t=t0;
      for ( i=0; i<nX; i++) 
        fprintf( fp[i], "Number_of_Trajectories  %d \n", nTrajEff);
      for ( n=0; n<=nOfSteps; n++) {
        for ( i=0; i<nX; i++) 
          fprintf( fp[i], "%lG %lG %lG %lG %lG\n", t, 
		  OutputExpR[n][i]/nTrajEff, OutputExpI[n][i]/nTrajEff, 
		  OutputExpR[n][nX+i]/nTrajEff, OutputExpI[n][nX+i]/nTrajEff);
	t += dtsPerStep*dt;
      } 
      for ( i=0; i<nX; i++) fclose(fp[i]); 
    }
  }

  // Computation of the mean value
  for ( n=0; n<=nOfSteps; n++)
    for ( i=0; i<2*nX; i++) OutputExpR[n][i] /= nTrajEff;
  for ( n=0; n<=nOfSteps; n++)
    for ( i=0; i<2*nX; i++) OutputExpI[n][i] /= nTrajEff;

  //  
  // Writing the results in the files
  //
  for ( i=0; i<nX; i++) fp[i] = fopen(fname[i],"w"); 
  t=t0;
  for ( i=0; i<nX; i++) 
    fprintf( fp[i], "Number_of_Trajectories  %d \n", nTrajEff);
  for ( n=0; n<=nOfSteps; n++) {
    for ( i=0; i<nX; i++) 
      fprintf( fp[i], "%lG %lG %lG %lG %lG\n", t, 
	      OutputExpR[n][i], OutputExpI[n][i], 
	      OutputExpR[n][nX+i], OutputExpI[n][nX+i]);
    t += dtsPerStep*dt;
  } 
  for ( i=0; i<nX; i++) fclose(fp[i]); 
  delete fp;

  // Destruction of the table OutputExp
  for ( j=0; j<nOfSteps; j++) delete OutputExpR[j];
  for ( j=0; j<nOfSteps; j++) delete OutputExpI[j];
  delete OutputExpR;
  delete OutputExpI;
}

Order2Step::Order2Step(const State& psi, const Operator& theH, int theNL,
	const Operator* theL)
{
  H = theH;				// store private data
  nL = theNL;
  L = new Operator[nL];
  Ldag = new Operator[nL];
  dxi = new Complex[nL];
  for (int i=0; i<nL; i++) {
    L[i] = theL[i];
    Ldag[i] = L[i].hc();
    dxi[i] = 0;
  }
  stochasticFlag=0;
  numdtsUsed=0;

  psi2=psi;				// initialize temporaries
  dpsi=psi; 
}

void Order2Step::operator()(State& psi, double t, double dt,
	double& dtlast, ComplexRandom* rndm)
//
// Second order integration of the deterministic part: 
//      see Numerical Recipes, 2nd edition, Eq. (16.1.2).
// Euler for the stochastic part.
{
  Complex lexpect;
 
   
  if( rndm != 0 ) {				// any noise?
    double sqrtdt = sqrt(dt);
    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] = sqrtdt*(*rndm)();		// ...independent noise
    }
  }
  newsum = psi;					// initialize temp
  newsum = 0;
  derivs(t,psi,dpsi);				// also calc. stoch. terms
  dpsi *= (0.5*dt);
  // dpsi = k_1/2 = h/2 f(x_n,y_n) (cf. Numerical Recipes)
  psi2 = psi;
  psi2 += dpsi;       // psi2 = y_n + k_1/2 
  derivs(t+0.5*dt,psi2,dpsi);
  dpsi *= dt;
  psi += dpsi;      // dpsi = k_2 (cf. Numerical Recipes)
  psi += newsum;    // newsum = (Stoch. part)
  psi.normalize();
  numdtsUsed++;
}

Order4Step::Order4Step(const State& psi, const Operator& theH, int theNL,
	const Operator* theL)
{
  H = theH;				// store private data
  nL = theNL;
  L = new Operator[nL];
  Ldag = new Operator[nL];
  dxi = new Complex[nL];
  for (int i=0; i<nL; i++) {
    L[i] = theL[i];
    Ldag[i] = L[i].hc();
    dxi[i] = 0;
  }
  stochasticFlag=0;
  numdtsUsed=0;

  psi0=psi; 				// initialize temporaries
  psi2=psi; 
  dpsi=psi;
}

void Order4Step::operator()(State& psi, double t, double dt,
	double& dtlast, ComplexRandom* rndm)
//
// Fourth order integration of the deterministic part: 
//      see Numerical Recipes, 2nd edition, Eq. (16.1.3).
// Euler for the stochastic part.
// psi is replaced by the new state. psi1 is temporary storage.
{
  Complex lexpect;
  
  if( rndm != 0 ) {				// any noise?
    double sqrtdt = sqrt(dt);
    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] = sqrtdt*(*rndm)();		// ...independent noise
    }
  }
  newsum = psi;					// initialize temp
  newsum = 0;
  psi0 = psi;
  derivs(t,psi,dpsi);			// Also calc. stoch. terms
  dpsi *= (dt/6);
  psi += dpsi;          // dpsi = k1/6

  dpsi *= 3;
  psi2 = psi0;
  psi2 += dpsi;          // psi2 = psi0 + k1/2
  derivs(t+0.5*dt,psi2,dpsi);
  dpsi *= (dt/3);
  psi += dpsi;           // dpsi = k2/3

  dpsi *= 1.5;
  psi2 = psi0;
  psi2 += dpsi;          // psi2 = psi0 + k2/2
  derivs(t+0.5*dt,psi2,dpsi);
  dpsi *= (dt/3);
  psi += dpsi;           // dpsi = k3/3

  dpsi *= 3;
  psi2 = psi0;
  psi2 += dpsi;           // psi2 = psi0 + k3
  derivs(t+0.5*dt,psi2,dpsi);
  dpsi *= (dt/6);
  psi += dpsi;           // psi += k4/6
  psi += newsum;         // psi += (Stoch. part)
  psi.normalize();
  numdtsUsed++;
}

int IntegrationStep::getNumdts()
//
// Return number of deterministic timesteps performed since last
// call to getNumdts
{
  int i;
  i = numdtsUsed;
  numdtsUsed = 0;
  return i;
}

void AdaptiveStep::rkck(double t, double h)
//
// Given values for a state y and its derivative dydt known at t, use the
// fifth-order Cash-Karp Runge-Kutta method to advance the solution over an
// interval h and return the incremented state as yout.  Also return an
// estimate of the local truncation error in yout using the embedded
// fourth-order method.  The user supplies the routine derivs(t,y,dydt) which
// returns the derivative dydt at t.
//
// Adapted from the Numerical Recipes routine RKCK.C
//
{
  static double a2=0.2,a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2,
	b31=3.0/8.0,b32=9.0/40.0,b41=4.0,b42 = -4.0,b43=1.2,
	b51 = -55.0/81.0, b52=-25.0/9.0, b53 = -175.0/81.0,b54=35.0/27.0,
	b61=-1631.0/11264.0,b62=35.0/256.0,b63=-115.0/7168.0,
	b64=1265.0/4096.0,b65=253.0/4096.0,c1=37888.0/11417.0,
	c3=5120.0/529.0,c4=10240.0/19481.0,c6=512.0/1771.0,
	dc5 = -554.0/1771.0;
  double dc1=-831.0/18944.0,dc3=831.0/17920.0,
	dc4=-831.0/5120.0,dc6=277.0/2048.0;

  ytemp = y;
  derivs(t,ytemp,dydt);
  dydt *= (h*b21);
  ytemp += dydt;				// First step
  derivs(t+a2*h,ytemp,ak2);			// Second step
  ytemp = y;
  dydt *= b31;
  ak2 *= (h*b32);
  ytemp += dydt;
  ytemp += ak2;
  derivs(t+a3*h,ytemp,ak3);			// Third step
  dydt *= b41;
  ak2 *= b42;
  ak3 *= (h*b43);
  ytemp = y;
  ytemp += dydt;
  ytemp += ak2;
  ytemp += ak3;
  derivs(t+a4*h,ytemp,ak4);			// Fourth step
  dydt *= b51;
  ak2 *= b52;
  ak3 *= b53;
  ak4 *= (h*b54);
  ytemp = y;
  ytemp += dydt;
  ytemp += ak2;
  ytemp += ak3;
  ytemp += ak4;
  derivs(t+a5*h,ytemp,ak5);			// Fifth step
  dydt *= b61;
  ak2 *= b62;
  ak3 *= b63;
  ak4 *= b64;
  ak5 *= (h*b65);
  ytemp = y;
  ytemp += dydt;
  ytemp += ak2;
  ytemp += ak3;
  ytemp += ak4;
  ytemp += ak5;
  derivs(t+a6*h,ytemp,ak2);			// Sixth step
//
// Accumulate increments with proper weights
//
  dydt *= c1;
  ak3 *= c3;
  ak4 *= c4;
  ak2 *= (h*c6);
  yout = y;
  yout += dydt;
  yout += ak3;
  yout += ak4;
  yout += ak2;
//
// Estimate error as difference between fourth and fifth order methods.
//
  dydt *= dc1;
  ak3 *= dc3;
  ak4 *= dc4;
  ak5 *= dc5;
  ak2 *= dc6;
  yerr = dydt;
  yerr += ak3;
  yerr += ak4;
  yerr += ak5;
  yerr += ak2;
}

void AdaptiveStep::rkqs(double& t, double htry, double eps,
	double& hdid, double& hnext)
//
// Fifth-order Runge-Kutta step with monitoring of local truncation error to
// ensure accuracy and adjust stepsize.  Input are the State y and its
// derivative dydt at the starting value of the time t.  Also input are the
// stepsize to be attempted, htry, and the required accuracy, eps. On output,
// y and t are replaced by their new values, hdid is the stepsize that was
// actually accomplished, and hnext is the estimated next stepsize.  derivs
// is the user-supplied routine that computes the right-hand side derivative.
//
// Adapted from the Numerical Recipes routine RKQS.C
//
{
  double errmax,h,tnew;

  h=htry;
  for (;;) {
	rkck(t,h);
	errmax = sqrt(real(yerr*yerr));
	errmax /= eps;
	if (errmax > 1.0) {
		h=SAFETY*h*pow(errmax,PSHRNK);
//		if (h < 0.1*h) h *= 0.1;
		tnew=t+h;
		if (tnew == t) {
		  cerr << "stepsize underflow in rkqs" << endl;
		  cerr << "errmax = " << errmax << endl;
		  cerr << "h = " << h << endl;
		  cerr << "eps = " << eps << endl;
		  exit(1);
		}
		continue;
	} else {
		if (errmax > ERRCON) hnext=SAFETY*h*pow(errmax,PGROW);
		else hnext=5.0*h;
		t += (hdid=h);
		y = yout;
		numdtsUsed++;
		break;
	}
  }
}


void AdaptiveStep::odeint(State& ystart, double t1, double t2, double eps, double& h1,
	double hmin, int& nok, int& nbad, Complex* dxi)
//
// 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