traject.cc

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

//   Traject.cc	-- Stochastic simulation of QSD trajectories.
//     
//   Copyright (C) 1995  Todd Brun and Ruediger Schack
//   
//   This program is free software; you can redistribute it and/or modify
//   it under the terms of the GNU General Public License as published by
//   the Free Software Foundation; either version 2 of the License, or
//   (at your option) any later version.
//   
//   This program is distributed in the hope that it will be useful,
//   but WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//   GNU General Public License for more details.
//   
//   You should have received a copy of the GNU General Public License
//   along with this program; if not, write to the Free Software
//   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//
//   ----------------------------------------------------------------------
//   If you improve the code or make additions to it, or if you have
//   comments or suggestions, please contact us:
//
//   Dr. Todd Brun			        Tel    [phone]
//   Department of Physics                      FAX    [phone]
//   Queen Mary and Westfield College           email  [e-mail]
//   Mile End Road, London E1 4NS, UK
//
//   Dr. Ruediger Schack                        Tel    [phone]
//   Department of Mathematics                  FAX    [phone]
//   Royal Holloway, University of London       email  [e-mail]
//   Egham, Surrey TW20 0EX, UK
/////////////////////////////////////////////////////////////////////////////

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <fstream.h>

#include "Complex.h"
#include "State.h"
#include "Operator.h"
#include "ACG.h"
#include "Normal.h"
#include "Traject.h"

#define MAXSTP 10000
#define TINY 1.0e-30
#define SAFETY 0.9
#define PGROW -0.2
#define PSHRNK -0.25
#define ERRCON 1.89e-4

static const char rcsid[] = "$Id: Traject.cc,v 3.5 1996/11/22 10:34:37 mrigo Exp $";

Trajectory::Trajectory()	// Constructor
{
  error("A trajectory must be constructed with initial data!");
}

AdaptiveStep::AdaptiveStep(const State& psi, const Operator& theH, int theNL,
	const Operator* theL, double theEpsilon)
//
// Constructor for AdaptiveStep; this contains all the information needed
// to integrate the QSD equation for a single adaptive timestep.  The
// Hamiltonian and Lindblad operators are stored as private data; the
// State psi is used to initialize the temporary States used in the
// Runge-Kutta integration; nL is the number of Lindblad operators, and
// is also stored as private data; and the integration is performed up to
// an accuracy epsilon.
//
{
  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;
  }
  epsilon = theEpsilon;
  stochasticFlag=0;
  numdtsUsed = 0;
  dtNumList = 0;
  listDim = 0;
  listPtr = 0;

  ak2 = ak3 = ak4 = ak5 = psi;		// initialize temporaries
  ytemp = yerr = y = yout = dydt = psi;
  temp1 = newsum = psi;
}

AdaptiveStep::~AdaptiveStep()	// AdaptiveStep destructor
{
  if( listDim != 0 )
#ifndef NON_GNU_DELETE
    delete[] dtNumList;		// deallocate memory
#else
    delete[listDim] dtNumList;	// deallocate memory
#endif
}

AdaptiveJump::~AdaptiveJump()	// AdaptiveJump destructor
{
  if( nL != 0 )
#ifndef NON_GNU_DELETE
    delete[] lindVal;		// deallocate memory
#else
    delete[nL] lindVal;		// deallocate memory
#endif
  if( outFile != 0 ) fclose(outFile);
}

AdaptiveOrthoJump::~AdaptiveOrthoJump()     // AdaptiveOrthoJump destructor
{
  if( nL != 0 )
#ifndef NON_GNU_DELETE
    delete[] lindVal;		// deallocate memory
#else
    delete[nL] lindVal;		// deallocate memory
#endif
  if( outFile != 0 ) fclose(outFile);
}

Trajectory::Trajectory(const State& thePsiIni, double thedt,
  IntegrationStep& theStepper, ComplexRandom* theRand, double theT0)
//
// Constructor for trajectory class.  This stores all the
// information needed to calculate a QSD trajectory.  The initial state,
// size of the output step, and number of dts per output step are stored
// as member data, and the size of dt is calculated; the seed for the
// random number generator is stored, as is the initial time; and the
// IntegrationStep contains all the information about the Hamiltonian and
// Lindblad operators, as well as the integration algorithm.
//
{
  psi = thePsiIni;
  dt = thedt;
  rndm = theRand;
  t0 = theT0;
  stepper = &theStepper;
}

IntegrationStep::~IntegrationStep()	// destructor
{
#ifndef NON_GNU_DELETE
  delete [] L;				// dispose of allocated memory
  delete [] Ldag;
  delete [] dxi;
#else // ---- NON GNU CODE ----
  delete [nL] L;			// dispose of allocated memory
  delete [nL] Ldag;
  delete [nL] dxi;
#endif // ---- NON GNU CODE ----
}

void AdaptiveStep::operator()(State& psi, double t, double dt,
	double& dtlast, ComplexRandom* rndm)
//
// This member function calculates a single QSD integration step.  It
// invokes the member function odeint (adapted from the Numerical Recipes
// routine) to advance the deterministic part of the equation in time,
// and uses the supplied function rndm to advance the stochastic
// part (with noise dxi).  The state psi is evolved from t to t+dt; dtlast
// contains the actual physical stepsize used by the adaptive stepsize
// routine.  (In case of a constant stepsize stepper, dtlast=dt.)
//
{
  int nok, nbad;

  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;
  odeint(psi,t,t+dt,epsilon,dtlast,0.0,nok,nbad,dxi);
	// generate deterministic evolution (and also stochastic term)
  psi += newsum;				// add stochastic term
  psi.normalize();
}

void AdaptiveStochStep::operator()(State& psi, double t, double dt,
	double& dtlast, ComplexRandom* rndm)
//
// This member function calculates a single QSD integration step.  It
// invokes the member function odeint (adapted from the Numerical Recipes
// routine) to advance the deterministic part of the equation in time,
// and uses the supplied function rndm to advance the stochastic
// part (with noise dxi).  The state psi is evolved from t to t+dt; dtlast
// contains the actual physical stepsize used by the adaptive stepsize
// routine.  (In case of a constant stepsize stepper, dtlast=dt.)
//
{
  int nok, nbad;

  newsum = psi;					// initialize temp
  odeint(psi,t,t+dt,epsilon,dtlast,0.0,nok,nbad,rndm);
	// generate deterministic evolution (and also stochastic)
}

void AdaptiveJump::operator()(State& psi, double t, double dt,
	double& dtlast, ComplexRandom* rndm)
//
// This member function calculates a single Quantum Jumps integration step.
// It invokes the member function odeint (adapted from the Numerical Recipes
// routine) to advance the deterministic part of the equation in time,
// and uses the supplied function rndm to advance the stochastic
// part (with noise dxi).  The state psi is evolved from t to t+dt; dtlast
// contains the actual physical stepsize used by the adaptive stepsize
// routine.  (In case of a constant stepsize stepper, dtlast=dt.)
//
{
  int nok, nbad;

  odeint(psi,t,t+dt,epsilon,dtlast,0.0,nok,nbad,rndm);
	// generate deterministic evolution (and also stochastic)
}

void AdaptiveOrthoJump::operator()(State& psi, double t, double dt,
	double& dtlast, ComplexRandom* rndm)
//
// This member function calculates a single QSD Quantum Jumps integration step.
// QSD Quantum Jumps as the same deterministic part as QSD but the stochastic
// part correspond to jumps.
// It invokes the member function odeint (adapted from the Numerical Recipes
// routine) to advance the deterministic part of the equation in time,
// and uses the supplied function rndm to advance the stochastic
// part (with noise dxi).  The state psi is evolved from t to t+dt; dtlast
// contains the actual physical stepsize used by the adaptive stepsize
// routine.  (In case of a constant stepsize stepper, dtlast=dt.)
//
{
  int nok, nbad;

  odeint(psi,t,t+dt,epsilon,dtlast,0.0,nok,nbad,rndm);
	// generate deterministic evolution (and also stochastic)
}

State Trajectory::getState()
// Returns the value of psi stored in the trajectory.
{
  return psi;
}

void Trajectory::plotExp_obsolete( int nX, const Operator* X, FILE** fp, int* pipe,
	int dtsPerStep, int nOfSteps, int move,
	double delta, int width, double moveEps)
//
// This member function calculates a trajectory consisting of nOfSteps
// output steps, and at each step stores the expectation values of
// the nX operators X into files.  The array `pipe'
// instructs plotExp which expectation values it should print to the
// standard output.  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.  After completing the
// trajectory, the initial state psi is reset to the final state.
//
{
  double t=t0;				// time
  double dtlast=dt;			// keep track of numerical timestep
  double xpipe[4];			// output to standard output
  if( nX > 0 ) {
    for( int k=0; k<4; k++ )
      if( pipe[k] < 1 || pipe[k] > 4*nX )
        cerr << "Warning: illegal argument 'pipe' in Trajectory::plotExp." 
          << endl;
  }

  State psi1 = psi;			// temporary state

  for( int 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;
    fprintf( fp[i], "%lG %lG %lG %lG %lG\n", 
	    t, real(expec), imag(expec), real(expec2), imag(expec2) ); // file
    for( int k=0; k<4; k++ ) {
      if( pipe[k] == 4*i+1 ) xpipe[k] = real(expec);          // standard i/o
      else if( pipe[k] == 4*i+2 ) xpipe[k] = imag(expec);
      else if( pipe[k] == 4*i+3 ) xpipe[k] = real(expec2);
      else if( pipe[k] == 4*i+4 ) xpipe[k] = imag(expec2);
    }
  }
  if( nX > 0 )
    printf("%lG %lG %lG %lG %lG %d %d\n",
      t,xpipe[0],xpipe[1],xpipe[2],xpipe[3],psi.size(),(*stepper).getNumdts());

  for (int n=1; n<=nOfSteps; n++) {
    for (int n1=0; n1<dtsPerStep; n1++) {
      (*stepper)(psi,t,dt,dtlast,rndm);
      t += dt;
      for(int k=0; k<move; k++) {		// use moving basis & cutoff
        psi.adjustCutoff(k,delta,width);
        psi.recenter(k,moveEps);
      }
    }
    for( int 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;
      fprintf( fp[i], "%lG %lG %lG %lG %lG\n", 
	      t, real(expec),imag(expec), real(expec2),imag(expec2) ); // file
      for( int k=0; k<4; k++ ) {
	if( pipe[k] == 4*i+1 ) xpipe[k] = real(expec);          // standard i/o
	else if( pipe[k] == 4*i+2 ) xpipe[k] = imag(expec);
	else if( pipe[k] == 4*i+3 ) xpipe[k] = real(expec2);
	else if( pipe[k] == 4*i+4 ) xpipe[k] = imag(expec2);
      }
    }
    if( nX > 0 )
      printf("%lG %lG %lG %lG %lG %d %d\n",
        t,xpipe[0],xpipe[1],xpipe[2],xpipe[3],psi.size(),(*stepper).getNumdts());
  }
//
// Save state to output file
//
//  ofstream os("saved_state.dat", ios::out );
//  if (!os) {
//    cerr << "Can't open output file:  saved_state.dat" << endl;
//    exit(1);
//  }
//  os << psi;
}


void Trajectory::plotExp( int nX, const Operator* X, char** fname, int* pipe,
	int dtsPerStep, int nOfSteps, int move,
	double delta, int width, double moveEps, char* savedState)
//
// This member function calculates a trajectory consisting of nOfSteps
// output steps, and at each step stores the expectation values of
// the nX operators X into files.  The array `pipe'
// instructs plotExp which expectation values it should print to the
// standard output.  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.  After completing the
// trajectory, the initial state psi is reset to the final state.
// Unless the string variable `savedState' is equal to the Null pointer 
// (which is the default), the final state is saved to the file `savedState'.
{
  double t=t0;				// time
  double dtlast=dt;			// keep track of numerical timestep
  double xpipe[4];			// output to standard output
  FILE** fp;
  fp = new FILE*[nX];                   // nX files are used

  if( nX > 0 ) 
    for ( int i=0; i<nX; i++) fp[i] = fopen(fname[i],"w"); 

  if( nX > 0 ) {
    for( int k=0; k<4; k++ )
      if( pipe[k] < 1 || pipe[k] > 4*nX )
        cerr << "Warning: illegal argument 'pipe' in Trajectory::plotExp." 
          << endl;
  }

  State psi1 = psi;			// temporary state

  for( int 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;
    fprintf( fp[i], "%lG %lG %lG %lG %lG\n", 
	    t, real(expec), imag(expec), real(expec2), imag(expec2) ); // file
    for( int k=0; k<4; k++ ) {
      if( pipe[k] == 4*i+1 ) xpipe[k] = real(expec);          // standard i/o
      else if( pipe[k] == 4*i+2 ) xpipe[k] = imag(expec);
      else if( pipe[k] == 4*i+3 ) xpipe[k] = real(expec2);
      else if( pipe[k] == 4*i+4 ) xpipe[k] = imag(expec2);
    }
  }
  if( nX > 0 )
    printf("%lG %lG %lG %lG %lG %d %d\n",
      t,xpipe[0],xpipe[1],xpipe[2],xpipe[3],psi.size(),(*stepper).getNumdts());

  for (int n=1; n<=nOfSteps; n++) {
    for (int n1=0; n1<dtsPerStep; n1++) {
      (*stepper)(psi,t,dt,dtlast,rndm);
      t += dt;
      for(int k=0; k<move; k++) {		// use moving basis & cutoff
        psi.adjustCutoff(k,delta,width);
        psi.recenter(k,moveEps);
      }
    }
    for( int 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;
      fprintf( fp[i], "%lG %lG %lG %lG %lG\n", 
	      t, real(expec),imag(expec), real(expec2),imag(expec2) ); // file
      for( int k=0; k<4; k++ ) {
	if( pipe[k] == 4*i+1 ) xpipe[k] = real(expec);          // standard i/o
	else if( pipe[k] == 4*i+2 ) xpipe[k] = imag(expec);
	else if( pipe[k] == 4*i+3 ) xpipe[k] = real(expec2);
	else if( pipe[k] == 4*i+4 ) xpipe[k] = imag(expec2);
      }
    }
    if( nX > 0 )
      printf("%lG %lG %lG %lG %lG %d %d\n",
        t,xpipe[0],xpipe[1],xpipe[2],xpipe[3],psi.size(),(*stepper).getNumdts());
  }

  if( nX > 0 ) {
    for ( int i=0; i<nX; i++) fclose(fp[i]); 
    delete fp;
  }
  // Save state to output file
  if (savedState) {
    ofstream os(savedState, ios::out );
    if (!os) {
      cerr << "Can't open output file: " << savedState << endl;
      exit(1);
    }
    os << psi;
  }
}

void Trajectory::sumExp( int nX, const Operator* X, char** fname,
	int dtsPerStep, int nOfSteps, 
        int nTrajectory, int nTrajSave, int ReadFile,