state.cc

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

  cout << "data pointer " << data << ".\n";
  cout << "My Pointer " << myPointer << ".\n";
  cout << "Dimensions:\n";
  if (nFreedoms > 1 ) {
    for (i=0; i<nFreedoms; i++)
      cout << "  " << nDims[i] << ".\n";
    cout << "Sizes:\n";
    for (i=0; i<nFreedoms; i++)
      cout << "  " << sizes[i] << ".\n";
    cout << "Skips:\n";
    for (i=0; i<nFreedoms; i++)
      cout << "  " << nSkips[i] << ".\n";
    cout << "Partial dimensions:\n";
    for (i=0; i<nFreedoms; i++)
      cout << "  " << partDims[i] << ".\n";
    cout << "Centers of coordinates:\n";
    for (i=0; i<nFreedoms; i++)
      cout << "  " << coords[i] << ".\n";
  }
}

// Friend I/O functions

ostream& operator<<( ostream& s, const State& a )
//
// Output a State to an output stream
//
{
  int i;

  s << a.nFreedoms << endl;
  s << a.totalDim << endl;
  s << a.maxSize << endl;
  s << a.mySize << endl;
  s << a.nSkip << endl;
  s << a.coord << endl;
  s << a.myType << endl;
  if( a.nFreedoms > 1 )
    for(i=0; i<a.nFreedoms; i++) {
      s << a.nDims[i] << endl;
      s << a.sizes[i] << endl;
      s << a.nSkips[i] << endl;
      s << a.partDims[i] << endl;
      s << a.freedomTypes[i] << endl;
      s << a.coords[i] << endl;
    }
  for (i=0; i<a.maxSize; i++)
    s << a.data[i] << endl;
  return s;
}

istream& operator>>( istream& s, State& a )
//
// Input a State from an input stream
//
{
  int i;

  a.free();			// free up State
  s >> a.nFreedoms;
  s >> a.totalDim;
  s >> a.maxSize;
  s >> a.mySize;
  s >> a.nSkip;
  s >> a.coord;
  s >> a.myType;
  if( a.totalDim > 0 ) {
    a.data = new Complex[a.totalDim];
    a.myPointer = a.data;
  }
  if( a.nFreedoms > 1 ) {
    a.nDims = new int[a.nFreedoms];
    a.sizes = new int[a.nFreedoms];
    a.nSkips = new int[a.nFreedoms];
    a.partDims = new int[a.nFreedoms];
    a.freedomTypes = new FreedomType[a.nFreedoms];
    a.coords = new Complex[a.nFreedoms];
    for(i=0; i<a.nFreedoms; i++) {
      s >> a.nDims[i];
      s >> a.sizes[i];
      s >> a.nSkips[i];
      s >> a.partDims[i];
      s >> a.freedomTypes[i];
      s >> a.coords[i];
    }
  }
  for (i=0; i<a.maxSize; i++)
    s >> a.data[i];
  return s;
}

// apply a PrimaryOperator to the State

void State::apply(PrimaryOperator& theOp, int hc, int theFreedom, FreedomType theType, double t)
{
  int i,j;
  if( nFreedoms == 1 ) {		// 1 freedom case
#ifndef OPTIMIZE_QSD
    if( (theType != myType) && (theType != ALL) )
      error("Incompatible operator type in State::apply!");
#endif
    theOp.applyTo(*this,hc,t);
  }
  else {				// multi-freedom case
#ifndef OPTIMIZE_QSD
    if( (theFreedom < 0) || (theFreedom >= nFreedoms) )
      error("Illegal degree of freedom requested in State::apply!");
    if( (theType != freedomTypes[theFreedom]) && (theType != ALL) )
      error("Incompatible operator type in State::apply!");
#endif
    mySize = sizes[theFreedom];		// set up 1D state for PrimaryOperator
    nSkip = nSkips[theFreedom];
    coord = coords[theFreedom];
    for (i=0; i<nSkip; i++) {				// loop over inner and
    for (j=0; j<maxSize; j += partDims[theFreedom]) {	// outer indices
      myPointer = data + i + j;
      if( i+j > maxSize ) error("Gone off end of array!");
      theOp.applyTo(*this,hc,t);
    } }
    mySize = maxSize;
    nSkip = 1;
    coord = coords[0];
    myPointer = data;
  }
}

State operator*(const Complex& a, const State& b) {	// multiplication
  State result(b);
  return (result *= a);
}

State operator*(const State& b, const Complex& a) {	// multiplication
  State result(b);
  return (result *= a);
}

State operator*(double a, const State& b) {		// multiplication
  State result(b);
  return (result *= a);
}

State operator*(const State& b, double a) {		// multiplication
  State result(b);
  return (result *= a);
}

State operator*(ImaginaryUnit a, const State& b) {	// multiplication
  State result(b);
  return (result *= a);
}

State operator*(const State& b, ImaginaryUnit a) {	// multiplication
  State result(b);
  return (result *= a);
}

State operator+(const State& a, const State& b) {	// addition
  State result(a);
  return (result += b);
}

State operator-(const State& a, const State& b) {	// subtraction
  State result(a);
  return (result -= b);
}

State operator+(const State& a) {			// unary +
  State result(a);
  return result;
}

State operator-(const State& a) {			// unary -
  State result(a);
  for (int i=0; i<result.maxSize; i++) result.data[i] = - result.data[i];
  return result;
}

int State::size() { return mySize; }	// return state size

int State::getSize(int theFreedom)
{
#ifndef OPTIMIZE_QSD
  if( (theFreedom < 0) || (theFreedom >= nFreedoms) )
    error("Illegal freedom passed to State::getSize!");
#endif
  return sizes[theFreedom];
}

Complex State::centerCoords() { return coord; }	// return center of coordinates

Complex State::getCoords(int theFreedom)
{
#ifndef OPTIMIZE_QSD
  if( (theFreedom<0) || (theFreedom >= nFreedoms) )
    error("Illegal freedom requested in State::getCoords!");
#endif
  if( nFreedoms > 1 ) return coords[theFreedom];
  return coord;
}

void State::setCoords(Complex& alpha, int theFreedom)
//
// Reset the value of the phase space coordinate for a freedom.
{
#ifndef OPTIMIZE_QSD
  if( (theFreedom<0) || (theFreedom >= nFreedoms) )
    error("Illegal freedom requested in State::setCoords!");
#endif
  if( nFreedoms == 1 )
    coord = alpha;
  else {
    coords[theFreedom] = alpha;
    if( theFreedom == 0 ) coord = alpha;
  }
}

void State::displaceCoords(Complex& alpha, int theFreedom)
//
// Reset the value of the phase space coordinate for a freedom.
{
#ifndef OPTIMIZE_QSD
  if( (theFreedom<0) || (theFreedom >= nFreedoms) )
    error("Illegal freedom requested in State::setCoords!");
#endif
  if( nFreedoms == 1 )
    coord += alpha;
  else {
    coords[theFreedom] += alpha;
    if( theFreedom == 0 ) coord += alpha;
  }
}

// Move the center of coordinates to a new position in phase space.

void State::moveToCoords(const Complex& alpha, int theFreedom, double shiftAccuracy)
{
#ifndef OPTIMIZE_QSD
  if( (theFreedom < 0) || (theFreedom >= nFreedoms) )
    error("Illegal degree of freedom passed in State::moveToCoords!");
  if( nFreedoms == 1 ) {
    if( myType != FIELD )
      error("Can only change basis of FIELD freedoms in State::moveToCoords!");
  }
  else {
    if( freedomTypes[theFreedom] != FIELD )
      error("Can only change basis of FIELD freedoms in State::moveToCoords!");
  }
#endif
  Complex displacement;
  displacement = getCoords(theFreedom) - alpha;
  moveCoords(displacement,theFreedom,shiftAccuracy);
}

// Move the center of coordinates to match those of another state.

void State::centerOn(State& psi, double shiftAccuracy)
{
#ifndef OPTIMIZE_QSD
  if( nFreedoms != psi.nFreedoms )
    error("Mismatched number of degrees of freedom in State::centerOn!");
#endif
  Complex alpha;
  if( nFreedoms == 1 ) {
    if( (myType != FIELD) || (psi.myType != FIELD) )
      error("Can't change basis of non-FIELD type in State::centerOn!");
    alpha = centerCoords() - psi.centerCoords();
    moveCoords(alpha,0,shiftAccuracy);
  }
  else {
    for (int i=0; i<nFreedoms; i++) {
      if( freedomTypes[i] != psi.freedomTypes[i] )
        error("Mismatched types in State::centerOn!");
      if( freedomTypes[i] == FIELD ) {
        alpha = getCoords(i) - psi.getCoords(i);
        moveCoords(alpha,i);
      }
    }
  }
}

// Recenter the coordinates at the current mean position in phase space.

void State::recenter(int theFreedom, double shiftAccuracy)
{
#ifndef OPTIMIZE_QSD
  if( (theFreedom >= nFreedoms) || (theFreedom < 0) )
    error("Illegal degree of freedom in State::recenter!");
  if( nFreedoms == 1 ) {
    if (myType != FIELD) error("Only type FIELD can be recentered!");
  }
  else {
    if (freedomTypes[theFreedom] != FIELD)
      error("Only type FIELD can be recentered!");
  }
#endif
  int i,j,k,l;
  State dummy=(*this);
  Complex magnitude, alpha;
  magnitude = (*this)*(*this);
  if( magnitude == 0 ) error("Can't recenter zero state!");
  if( nFreedoms == 1 ) {		// 1 freedom case
    if( mySize > intSqrtMax ) setIntSqrtMax(mySize);
    for (k=1; k<mySize; k++)	// apply local annihilation op.
      dummy.myPointer[k-1] = intSqrtList[k]*dummy.myPointer[k];
    dummy.myPointer[mySize-1] = 0;
  }
  else {				// multi-freedom case
    mySize = sizes[theFreedom];
    nSkip = nSkips[theFreedom];
    if( mySize > intSqrtMax ) setIntSqrtMax(mySize);
    for (i=0; i<nSkip; i++)		// apply local annihilation op
    for (j=0; j<maxSize; j+=partDims[theFreedom]) {	// loop over other
      dummy.myPointer = dummy.data + i + j;		// indices
      for (k=1,l=nSkip; k<mySize; k++,l+=nSkip)
        dummy.myPointer[l-nSkip] = intSqrtList[k]*dummy.myPointer[l];
      dummy[mySize-1] = 0;
    }
    mySize = maxSize;
    nSkip = 1;
    dummy.myPointer = dummy.data;
  }
  alpha = - ((*this)*dummy)/magnitude;	// get relative coords in phase space
  if( abs(alpha) > shiftAccuracy )
    moveCoords(alpha,theFreedom, shiftAccuracy);	// shift basis
}

// Move the center of coordinates by a displacement alpha in
// phase space.

void State::moveCoords(const Complex& alpha, int theFreedom, double shiftAccuracy)
{
#ifndef OPTIMIZE_QSD
  if( (theFreedom >= nFreedoms) || (theFreedom < 0) )
    error("Illegal freedom passed to State::moveCoords!");
#endif
  Complex mag=(*this)*(*this);
  if( mag == 0 ) {				// move zero state
    if( nFreedoms == 1 ) coord -= alpha;
    if( nFreedoms > 1 ) {
      coords[theFreedom] -= alpha;
      coord = coords[0];
    }
    return;
  }
  int numshifts;				// Break basis shift into
  numshifts = int(abs(alpha)/shiftAccuracy)+1;	// sequence of small shifts
  Complex delta=alpha/numshifts;		// according to accuracy
  Complex deltacc=conj(delta);			// requirement.
  Complex theCoord, phase=0;
  if (nFreedoms == 1) {				// 1 freedom case
#ifndef OPTIMIZE_QSD
    if( myType != FIELD) error("Can only move coordinates for FIELD type.");
#endif
    theCoord = coord;
    if( mySize > intSqrtMax) setIntSqrtMax(mySize);
    for (int i=0; i<numshifts; i++) {
      moveStep(delta,deltacc);			// 1 small shift
      phase += imag(delta*conj(theCoord));	// sum up phase changes
      theCoord -= delta;
    }
    coord = theCoord;
  }
  else {					// multi-freedom case
#ifndef OPTIMIZE_QSD
    if( freedomTypes[theFreedom] != FIELD )
      error("Can only move coordinates for FIELD type.");
#endif
    mySize = sizes[theFreedom];
    nSkip = nSkips[theFreedom];
    theCoord = coords[theFreedom];
    if( mySize > intSqrtMax) setIntSqrtMax(mySize);
    for (int i=0; i<numshifts; i++) {
      for (int j=0; j<nSkip; j++)		// loop over other indices
      for (int k=0; k<maxSize; k+=partDims[theFreedom]) {
        myPointer = data + j + k;
        moveStep(delta,deltacc);		// 1 small shift
      }
      phase += imag(delta*conj(theCoord));	// sum up phase changes
      theCoord -= delta;
    }
    coords[theFreedom] = theCoord;
    if( theFreedom == 0 ) coord = coords[0];
    myPointer = data;
    mySize = maxSize;
    nSkip = 1;
  }
  phase.timesI();
  (*this) *= exp(phase);			// remove phase change
  normalize();					// normalize state
}

// Shift the center of coordinates by a single infinitesimal step

void State::moveStep(Complex& delta, Complex& deltacc)
{
  int i,j;
  Complex psizero,psinth,temp,previous;
  psizero = (*this)[0] - deltacc*(*this)[1];
  psinth = (*this)[mySize-1] + delta*intSqrtList[mySize-1]*(*this)[mySize-2];
  previous = (*this)[0];
  for (i=1,j=nSkip; i<(mySize-1); i++,j+=nSkip) {	// infinitesimal
    temp = myPointer[j];				// displacement
    myPointer[j] += delta*intSqrtList[i]*previous;
    myPointer[j] -= deltacc*intSqrtList[i+1]*myPointer[j+nSkip];
    previous = temp;
  }
  (*this)[0] = psizero;
  (*this)[mySize-1] = psinth;
}

void State::error(const char* message) const	// error handling
{
  cerr << message << "\n";
  exit(1);
}