OOPSE/libmdtools/Integrator.cpp

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

#ifdef IS_MPI
#include "mpiSimulation.hpp"
#include <unistd.h>
#endif //is_mpi

#ifdef PROFILE
#include "mdProfile.hpp"
#endif // profile

#include "Integrator.hpp"
#include "simError.h"


template<typename T> Integrator<T>::Integrator(SimInfo* theInfo,
                                               ForceFields* the_ff){
  info = theInfo;
  myFF = the_ff;
  isFirst = 1;

  molecules = info->molecules;
  nMols = info->n_mol;

  // give a little love back to the SimInfo object

  if (info->the_integrator != NULL){
    delete info->the_integrator;
  }

  nAtoms = info->n_atoms;

  // check for constraints

  constrainedA = NULL;
  constrainedB = NULL;
  constrainedDsqr = NULL;
  moving = NULL;
  moved = NULL;
  oldPos = NULL;

  nConstrained = 0;

  checkConstraints();
}

template<typename T> Integrator<T>::~Integrator(){
  if (nConstrained){
    delete[] constrainedA;
    delete[] constrainedB;
    delete[] constrainedDsqr;
    delete[] moving;
    delete[] moved;
    delete[] oldPos;
  }
}

template<typename T> void Integrator<T>::checkConstraints(void){
  isConstrained = 0;

  Constraint* temp_con;
  Constraint* dummy_plug;
  temp_con = new Constraint[info->n_SRI];
  nConstrained = 0;
  int constrained = 0;

  SRI** theArray;
  for (int i = 0; i < nMols; i++){

          theArray = (SRI * *) molecules[i].getMyBonds();
    for (int j = 0; j < molecules[i].getNBonds(); j++){
      constrained = theArray[j]->is_constrained();

      if (constrained){
        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a(dummy_plug->get_a());
        temp_con[nConstrained].set_b(dummy_plug->get_b());
        temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());

        nConstrained++;
        constrained = 0;
      }
    }

    theArray = (SRI * *) molecules[i].getMyBends();
    for (int j = 0; j < molecules[i].getNBends(); j++){
      constrained = theArray[j]->is_constrained();

      if (constrained){
        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a(dummy_plug->get_a());
        temp_con[nConstrained].set_b(dummy_plug->get_b());
        temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());

        nConstrained++;
        constrained = 0;
      }
    }

    theArray = (SRI * *) molecules[i].getMyTorsions();
    for (int j = 0; j < molecules[i].getNTorsions(); j++){
      constrained = theArray[j]->is_constrained();

      if (constrained){
        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a(dummy_plug->get_a());
        temp_con[nConstrained].set_b(dummy_plug->get_b());
        temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());

        nConstrained++;
        constrained = 0;
      }
    }
  }


  if (nConstrained > 0){
    isConstrained = 1;

    if (constrainedA != NULL)
      delete[] constrainedA;
    if (constrainedB != NULL)
      delete[] constrainedB;
    if (constrainedDsqr != NULL)
      delete[] constrainedDsqr;

    constrainedA = new int[nConstrained];
    constrainedB = new int[nConstrained];
    constrainedDsqr = new double[nConstrained];

    for (int i = 0; i < nConstrained; i++){
      constrainedA[i] = temp_con[i].get_a();
      constrainedB[i] = temp_con[i].get_b();
      constrainedDsqr[i] = temp_con[i].get_dsqr();
    }


    // save oldAtoms to check for lode balancing later on.

    oldAtoms = nAtoms;

    moving = new int[nAtoms];
    moved = new int[nAtoms];

    oldPos = new double[nAtoms * 3];
  }

  delete[] temp_con;
}


template<typename T> void Integrator<T>::integrate(void){

  double runTime = info->run_time;
  double sampleTime = info->sampleTime;
  double statusTime = info->statusTime;
  double thermalTime = info->thermalTime;
  double resetTime = info->resetTime;


  double currSample;
  double currThermal;
  double currStatus;
  double currReset;

  int calcPot, calcStress;

  tStats = new Thermo(info);
  statOut = new StatWriter(info);
  dumpOut = new DumpWriter(info);

  atoms = info->atoms;

  dt = info->dt;
  dt2 = 0.5 * dt;

  readyCheck();

  // initialize the forces before the first step

  calcForce(1, 1);

  //temp test
  tStats->getPotential();
  
  if (nConstrained){
    preMove();
    constrainA();
    calcForce(1, 1);
    constrainB();
  }
  
  if (info->setTemp){
    thermalize();
  }

  calcPot     = 0;
  calcStress  = 0;
  currSample  = sampleTime + info->getTime();
  currThermal = thermalTime+ info->getTime();
  currStatus  = statusTime + info->getTime();
  currReset   = resetTime  + info->getTime();

  dumpOut->writeDump(info->getTime());
  statOut->writeStat(info->getTime());


#ifdef IS_MPI
  strcpy(checkPointMsg, "The integrator is ready to go.");
  MPIcheckPoint();
#endif // is_mpi

  while (info->getTime() < runTime){
    if ((info->getTime() + dt) >= currStatus){
      calcPot = 1;
      calcStress = 1;
    }

#ifdef PROFILE
    startProfile( pro1 );
#endif
    
    integrateStep(calcPot, calcStress);

#ifdef PROFILE
    endProfile( pro1 );

    startProfile( pro2 );
#endif // profile

    info->incrTime(dt);

    if (info->setTemp){
      if (info->getTime() >= currThermal){
        thermalize();
        currThermal += thermalTime;
      }
    }

    if (info->getTime() >= currSample){
      dumpOut->writeDump(info->getTime());
      currSample += sampleTime;
    }

    if (info->getTime() >= currStatus){
      statOut->writeStat(info->getTime());
      calcPot = 0;
      calcStress = 0;
      currStatus += statusTime;
    }

    if (info->resetIntegrator){
      if (info->getTime() >= currReset){
        this->resetIntegrator();
        currReset += resetTime;
      }
    }
    
#ifdef PROFILE
    endProfile( pro2 );
#endif //profile

#ifdef IS_MPI
    strcpy(checkPointMsg, "successfully took a time step.");
    MPIcheckPoint();
#endif // is_mpi
  }

  delete dumpOut;
  delete statOut;
}

template<typename T> void Integrator<T>::integrateStep(int calcPot,
                                                       int calcStress){
  // Position full step, and velocity half step

#ifdef PROFILE
  startProfile(pro3);
#endif //profile

  preMove();

#ifdef PROFILE
  endProfile(pro3);

  startProfile(pro4);
#endif // profile

  moveA();

#ifdef PROFILE
  endProfile(pro4);
  
  startProfile(pro5);
#endif//profile


#ifdef IS_MPI
  strcpy(checkPointMsg, "Succesful moveA\n");
  MPIcheckPoint();
#endif // is_mpi


  // calc forces

  calcForce(calcPot, calcStress);

#ifdef IS_MPI
  strcpy(checkPointMsg, "Succesful doForces\n");
  MPIcheckPoint();
#endif // is_mpi

#ifdef PROFILE
  endProfile( pro5 );

  startProfile( pro6 );
#endif //profile

  // finish the velocity  half step

  moveB();

#ifdef PROFILE
  endProfile(pro6);
#endif // profile

#ifdef IS_MPI
  strcpy(checkPointMsg, "Succesful moveB\n");
  MPIcheckPoint();
#endif // is_mpi
}


template<typename T> void Integrator<T>::moveA(void){
  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double vel[3], pos[3], frc[3];
  double mass;

  for (i = 0; i < nAtoms; i++){
    atoms[i]->getVel(vel);
    atoms[i]->getPos(pos);
    atoms[i]->getFrc(frc);

    mass = atoms[i]->getMass();

    for (j = 0; j < 3; j++){
      // velocity half step
      vel[j] += (dt2 * frc[j] / mass) * eConvert;
      // position whole step
      pos[j] += dt * vel[j];
    }

    atoms[i]->setVel(vel);
    atoms[i]->setPos(pos);

    if (atoms[i]->isDirectional()){
      dAtom = (DirectionalAtom *) atoms[i];

      // get and convert the torque to body frame

      dAtom->getTrq(Tb);
      dAtom->lab2Body(Tb);

      // get the angular momentum, and propagate a half step

      dAtom->getJ(ji);

      for (j = 0; j < 3; j++)
        ji[j] += (dt2 * Tb[j]) * eConvert;

      this->rotationPropagation( dAtom, ji );

      dAtom->setJ(ji);
    }
  }

  if (nConstrained){
    constrainA();
  }
}


template<typename T> void Integrator<T>::moveB(void){
  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double vel[3], frc[3];
  double mass;

  for (i = 0; i < nAtoms; i++){
    atoms[i]->getVel(vel);
    atoms[i]->getFrc(frc);

    mass = atoms[i]->getMass();

    // velocity half step
    for (j = 0; j < 3; j++)
      vel[j] += (dt2 * frc[j] / mass) * eConvert;

    atoms[i]->setVel(vel);

    if (atoms[i]->isDirectional()){
      dAtom = (DirectionalAtom *) atoms[i];

      // get and convert the torque to body frame

      dAtom->getTrq(Tb);
      dAtom->lab2Body(Tb);

      // get the angular momentum, and propagate a half step

      dAtom->getJ(ji);

      for (j = 0; j < 3; j++)
        ji[j] += (dt2 * Tb[j]) * eConvert;


      dAtom->setJ(ji);
    }
  }

  if (nConstrained){
    constrainB();
  }
}

template<typename T> void Integrator<T>::preMove(void){
  int i, j;
  double pos[3];

  if (nConstrained){
    for (i = 0; i < nAtoms; i++){
      atoms[i]->getPos(pos);

      for (j = 0; j < 3; j++){
        oldPos[3 * i + j] = pos[j];
      }
    }
  }
}

template<typename T> void Integrator<T>::constrainA(){
  int i, j;
  int done;
  double posA[3], posB[3];
  double velA[3], velB[3];
  double pab[3];
  double rab[3];
  int a, b, ax, ay, az, bx, by, bz;
  double rma, rmb;
  double dx, dy, dz;
  double rpab;
  double rabsq, pabsq, rpabsq;
  double diffsq;
  double gab;
  int iteration;

  for (i = 0; i < nAtoms; i++){
    moving[i] = 0;
    moved[i] = 1;
  }

  iteration = 0;
  done = 0;
  while (!done && (iteration < maxIteration)){
    done = 1;
    for (i = 0; i < nConstrained; i++){
      a = constrainedA[i];
      b = constrainedB[i];

      ax = (a * 3) + 0;
      ay = (a * 3) + 1;
      az = (a * 3) + 2;

      bx = (b * 3) + 0;
      by = (b * 3) + 1;
      bz = (b * 3) + 2;

      if (moved[a] || moved[b]){
        atoms[a]->getPos(posA);
        atoms[b]->getPos(posB);

        for (j = 0; j < 3; j++)
          pab[j] = posA[j] - posB[j];

        //periodic boundary condition

        info->wrapVector(pab);

        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];

        rabsq = constrainedDsqr[i];
        diffsq = rabsq - pabsq;

        // the original rattle code from alan tidesley
        if (fabs(diffsq) > (tol * rabsq * 2)){
          rab[0] = oldPos[ax] - oldPos[bx];
          rab[1] = oldPos[ay] - oldPos[by];
          rab[2] = oldPos[az] - oldPos[bz];

          info->wrapVector(rab);

          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];

          rpabsq = rpab * rpab;


          if (rpabsq < (rabsq * -diffsq)){
#ifdef IS_MPI
            a = atoms[a]->getGlobalIndex();
            b = atoms[b]->getGlobalIndex();
#endif //is_mpi
            sprintf(painCave.errMsg,
                    "Constraint failure in constrainA at atom %d and %d.\n", a,
                    b);
            painCave.isFatal = 1;
            simError();
          }

          rma = 1.0 / atoms[a]->getMass();
          rmb = 1.0 / atoms[b]->getMass();

          gab = diffsq / (2.0 * (rma + rmb) * rpab);

          dx = rab[0] * gab;
          dy = rab[1] * gab;
          dz = rab[2] * gab;

          posA[0] += rma * dx;
          posA[1] += rma * dy;
          posA[2] += rma * dz;

          atoms[a]->setPos(posA);

          posB[0] -= rmb * dx;
          posB[1] -= rmb * dy;
          posB[2] -= rmb * dz;

          atoms[b]->setPos(posB);

          dx = dx / dt;
          dy = dy / dt;
          dz = dz / dt;

          atoms[a]->getVel(velA);

          velA[0] += rma * dx;
          velA[1] += rma * dy;
          velA[2] += rma * dz;

          atoms[a]->setVel(velA);

          atoms[b]->getVel(velB);

          velB[0] -= rmb * dx;
          velB[1] -= rmb * dy;
          velB[2] -= rmb * dz;

          atoms[b]->setVel(velB);

          moving[a] = 1;
          moving[b] = 1;
          done = 0;
        }
      }
    }

    for (i = 0; i < nAtoms; i++){
      moved[i] = moving[i];
      moving[i] = 0;
    }

    iteration++;
  }

  if (!done){
    sprintf(painCave.errMsg,
            "Constraint failure in constrainA, too many iterations: %d\n",
            iteration);
    painCave.isFatal = 1;
    simError();
  }

}

template<typename T> void Integrator<T>::constrainB(void){
  int i, j;
  int done;
  double posA[3], posB[3];
  double velA[3], velB[3];
  double vxab, vyab, vzab;
  double rab[3];
  int a, b, ax, ay, az, bx, by, bz;
  double rma, rmb;
  double dx, dy, dz;
  double rvab;
  double gab;
  int iteration;

  for (i = 0; i < nAtoms; i++){
    moving[i] = 0;
    moved[i] = 1;
  }

  done = 0;
  iteration = 0;
  while (!done && (iteration < maxIteration)){
    done = 1;

    for (i = 0; i < nConstrained; i++){
      a = constrainedA[i];
      b = constrainedB[i];

      ax = (a * 3) + 0;
      ay = (a * 3) + 1;
      az = (a * 3) + 2;

      bx = (b * 3) + 0;
      by = (b * 3) + 1;
      bz = (b * 3) + 2;

      if (moved[a] || moved[b]){
        atoms[a]->getVel(velA);
        atoms[b]->getVel(velB);

        vxab = velA[0] - velB[0];
        vyab = velA[1] - velB[1];
        vzab = velA[2] - velB[2];

        atoms[a]->getPos(posA);
        atoms[b]->getPos(posB);

        for (j = 0; j < 3; j++)
          rab[j] = posA[j] - posB[j];

        info->wrapVector(rab);

        rma = 1.0 / atoms[a]->getMass();
        rmb = 1.0 / atoms[b]->getMass();

        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;

        gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);

        if (fabs(gab) > tol){
          dx = rab[0] * gab;
          dy = rab[1] * gab;
          dz = rab[2] * gab;

          velA[0] += rma * dx;
          velA[1] += rma * dy;
          velA[2] += rma * dz;

          atoms[a]->setVel(velA);

          velB[0] -= rmb * dx;
          velB[1] -= rmb * dy;
          velB[2] -= rmb * dz;

          atoms[b]->setVel(velB);

          moving[a] = 1;
          moving[b] = 1;
          done = 0;
        }
      }
    }

    for (i = 0; i < nAtoms; i++){
      moved[i] = moving[i];
      moving[i] = 0;
    }

    iteration++;
  }

  if (!done){
    sprintf(painCave.errMsg,
            "Constraint failure in constrainB, too many iterations: %d\n",
            iteration);
    painCave.isFatal = 1;
    simError();
  }
}

template<typename T> void Integrator<T>::rotationPropagation
( DirectionalAtom* dAtom, double ji[3] ){

  double angle;
  double A[3][3], I[3][3];

  // use the angular velocities to propagate the rotation matrix a
  // full time step

  dAtom->getA(A);
  dAtom->getI(I);

  // rotate about the x-axis
  angle = dt2 * ji[0] / I[0][0];
  this->rotate( 1, 2, angle, ji, A );

  // rotate about the y-axis
  angle = dt2 * ji[1] / I[1][1];
  this->rotate( 2, 0, angle, ji, A );

  // rotate about the z-axis
  angle = dt * ji[2] / I[2][2];
  this->rotate( 0, 1, angle, ji, A);

  // rotate about the y-axis
  angle = dt2 * ji[1] / I[1][1];
  this->rotate( 2, 0, angle, ji, A );

  // rotate about the x-axis
  angle = dt2 * ji[0] / I[0][0];
  this->rotate( 1, 2, angle, ji, A );

  dAtom->setA( A  );
}

template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
                                                double angle, double ji[3],
                                                double A[3][3]){
  int i, j, k;
  double sinAngle;
  double cosAngle;
  double angleSqr;
  double angleSqrOver4;
  double top, bottom;
  double rot[3][3];
  double tempA[3][3];
  double tempJ[3];

  // initialize the tempA

  for (i = 0; i < 3; i++){
    for (j = 0; j < 3; j++){
      tempA[j][i] = A[i][j];
    }
  }

  // initialize the tempJ

  for (i = 0; i < 3; i++)
    tempJ[i] = ji[i];

  // initalize rot as a unit matrix

  rot[0][0] = 1.0;
  rot[0][1] = 0.0;
  rot[0][2] = 0.0;

  rot[1][0] = 0.0;
  rot[1][1] = 1.0;
  rot[1][2] = 0.0;

  rot[2][0] = 0.0;
  rot[2][1] = 0.0;
  rot[2][2] = 1.0;

  // use a small angle aproximation for sin and cosine

  angleSqr = angle * angle;
  angleSqrOver4 = angleSqr / 4.0;
  top = 1.0 - angleSqrOver4;
  bottom = 1.0 + angleSqrOver4;

  cosAngle = top / bottom;
  sinAngle = angle / bottom;

  rot[axes1][axes1] = cosAngle;
  rot[axes2][axes2] = cosAngle;

  rot[axes1][axes2] = sinAngle;
  rot[axes2][axes1] = -sinAngle;

  // rotate the momentum acoording to: ji[] = rot[][] * ji[]

  for (i = 0; i < 3; i++){
    ji[i] = 0.0;
    for (k = 0; k < 3; k++){
      ji[i] += rot[i][k] * tempJ[k];
    }
  }

  // rotate the Rotation matrix acording to:
  //            A[][] = A[][] * transpose(rot[][])


  // NOte for as yet unknown reason, we are performing the
  // calculation as:
  //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])

  for (i = 0; i < 3; i++){
    for (j = 0; j < 3; j++){
      A[j][i] = 0.0;
      for (k = 0; k < 3; k++){
        A[j][i] += tempA[i][k] * rot[j][k];
      }
    }
  }
}

template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
  myFF->doForces(calcPot, calcStress);
}

template<typename T> void Integrator<T>::thermalize(){
  tStats->velocitize();
}

template<typename T> double Integrator<T>::getConservedQuantity(void){
  return tStats->getTotalE();
}
template<typename T> string Integrator<T>::getAdditionalParameters(void){
  //By default, return a null string
  //The reason we use string instead of char* is that if we use char*, we will
  //return a pointer point to local variable which might cause problem
  return string();
}
Revision:	1057
Committed:	Tue Feb 17 19:23:44 2004 UTC (20 years, 4 months ago) by tim
File size:	17369 byte(s)
Log Message:	adding function shakeF in order to remove the constraint force along bond direction
#	User	Rev	Content
1	mmeineke	558	#include <iostream>
2	gezelter	829	#include <stdlib.h>
3			#include <math.h>
4	mmeineke	558
5			#ifdef IS_MPI
6			#include "mpiSimulation.hpp"
7			#include <unistd.h>
8			#endif //is_mpi
9
10	chuckv	892	#ifdef PROFILE
11			#include "mdProfile.hpp"
12			#endif // profile
13
14	mmeineke	558	#include "Integrator.hpp"
15			#include "simError.h"
16
17
18	tim	725	template<typename T> Integrator<T>::Integrator(SimInfo* theInfo,
19			ForceFields* the_ff){
20	mmeineke	558	info = theInfo;
21			myFF = the_ff;
22			isFirst = 1;
23
24			molecules = info->molecules;
25			nMols = info->n_mol;
26
27			// give a little love back to the SimInfo object
28	tim	725
29			if (info->the_integrator != NULL){
30			delete info->the_integrator;
31			}
32	tim	837
33	mmeineke	558	nAtoms = info->n_atoms;
34
35			// check for constraints
36	tim	725
37			constrainedA = NULL;
38			constrainedB = NULL;
39	mmeineke	558	constrainedDsqr = NULL;
40	tim	725	moving = NULL;
41			moved = NULL;
42			oldPos = NULL;
43
44	mmeineke	558	nConstrained = 0;
45
46			checkConstraints();
47			}
48
49	tim	725	template<typename T> Integrator<T>::~Integrator(){
50			if (nConstrained){
51	mmeineke	558	delete[] constrainedA;
52			delete[] constrainedB;
53			delete[] constrainedDsqr;
54			delete[] moving;
55			delete[] moved;
56	mmeineke	561	delete[] oldPos;
57	mmeineke	558	}
58			}
59
60	tim	725	template<typename T> void Integrator<T>::checkConstraints(void){
61	mmeineke	558	isConstrained = 0;
62
63	tim	725	Constraint* temp_con;
64			Constraint* dummy_plug;
65	mmeineke	558	temp_con = new Constraint[info->n_SRI];
66			nConstrained = 0;
67			int constrained = 0;
68	tim	725
69	mmeineke	558	SRI** theArray;
70	tim	725	for (int i = 0; i < nMols; i++){
71	tim	1057
72			theArray = (SRI * *) molecules[i].getMyBonds();
73	tim	725	for (int j = 0; j < molecules[i].getNBonds(); j++){
74	mmeineke	558	constrained = theArray[j]->is_constrained();
75	mmeineke	594
76	tim	725	if (constrained){
77			dummy_plug = theArray[j]->get_constraint();
78			temp_con[nConstrained].set_a(dummy_plug->get_a());
79			temp_con[nConstrained].set_b(dummy_plug->get_b());
80			temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());
81	mmeineke	594
82	tim	725	nConstrained++;
83			constrained = 0;
84			}
85	mmeineke	558	}
86
87	tim	725	theArray = (SRI * *) molecules[i].getMyBends();
88			for (int j = 0; j < molecules[i].getNBends(); j++){
89	mmeineke	558	constrained = theArray[j]->is_constrained();
90	tim	725
91			if (constrained){
92			dummy_plug = theArray[j]->get_constraint();
93			temp_con[nConstrained].set_a(dummy_plug->get_a());
94			temp_con[nConstrained].set_b(dummy_plug->get_b());
95			temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());
96
97			nConstrained++;
98			constrained = 0;
99	mmeineke	558	}
100			}
101
102	tim	725	theArray = (SRI * *) molecules[i].getMyTorsions();
103			for (int j = 0; j < molecules[i].getNTorsions(); j++){
104	mmeineke	558	constrained = theArray[j]->is_constrained();
105	tim	725
106			if (constrained){
107			dummy_plug = theArray[j]->get_constraint();
108			temp_con[nConstrained].set_a(dummy_plug->get_a());
109			temp_con[nConstrained].set_b(dummy_plug->get_b());
110			temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());
111
112			nConstrained++;
113			constrained = 0;
114	mmeineke	558	}
115			}
116			}
117
118	tim	1057
119	tim	725	if (nConstrained > 0){
120	mmeineke	558	isConstrained = 1;
121
122	tim	725	if (constrainedA != NULL)
123			delete[] constrainedA;
124			if (constrainedB != NULL)
125			delete[] constrainedB;
126			if (constrainedDsqr != NULL)
127			delete[] constrainedDsqr;
128	mmeineke	558
129	tim	725	constrainedA = new int[nConstrained];
130			constrainedB = new int[nConstrained];
131	mmeineke	558	constrainedDsqr = new double[nConstrained];
132	tim	725
133			for (int i = 0; i < nConstrained; i++){
134	mmeineke	558	constrainedA[i] = temp_con[i].get_a();
135			constrainedB[i] = temp_con[i].get_b();
136			constrainedDsqr[i] = temp_con[i].get_dsqr();
137			}
138
139	tim	725
140	chrisfen	999	// save oldAtoms to check for lode balancing later on.
141	tim	725
142	mmeineke	558	oldAtoms = nAtoms;
143	tim	725
144	mmeineke	558	moving = new int[nAtoms];
145	tim	725	moved = new int[nAtoms];
146	mmeineke	558
147	tim	725	oldPos = new double[nAtoms * 3];
148	mmeineke	558	}
149	tim	725
150	mmeineke	558	delete[] temp_con;
151			}
152
153
154	tim	725	template<typename T> void Integrator<T>::integrate(void){
155	mmeineke	558
156	tim	725	double runTime = info->run_time;
157			double sampleTime = info->sampleTime;
158			double statusTime = info->statusTime;
159	mmeineke	558	double thermalTime = info->thermalTime;
160	mmeineke	746	double resetTime = info->resetTime;
161	mmeineke	558
162	mmeineke	746
163	mmeineke	558	double currSample;
164			double currThermal;
165			double currStatus;
166	mmeineke	746	double currReset;
167	tim	837
168	mmeineke	558	int calcPot, calcStress;
169
170	tim	725	tStats = new Thermo(info);
171			statOut = new StatWriter(info);
172			dumpOut = new DumpWriter(info);
173	mmeineke	558
174	mmeineke	561	atoms = info->atoms;
175
176			dt = info->dt;
177	mmeineke	558	dt2 = 0.5 * dt;
178
179	mmeineke	784	readyCheck();
180
181	mmeineke	558	// initialize the forces before the first step
182
183	tim	677	calcForce(1, 1);
184	tim	781
185	tim	1035	//temp test
186			tStats->getPotential();
187
188	tim	781	if (nConstrained){
189			preMove();
190			constrainA();
191	tim	837	calcForce(1, 1);
192	tim	781	constrainB();
193			}
194	mmeineke	843
195	tim	725	if (info->setTemp){
196	tim	677	thermalize();
197	mmeineke	558	}
198	tim	725
199	mmeineke	558	calcPot = 0;
200			calcStress = 0;
201	mmeineke	711	currSample = sampleTime + info->getTime();
202			currThermal = thermalTime+ info->getTime();
203			currStatus = statusTime + info->getTime();
204	mmeineke	746	currReset = resetTime + info->getTime();
205	mmeineke	558
206	tim	725	dumpOut->writeDump(info->getTime());
207			statOut->writeStat(info->getTime());
208	mmeineke	559
209
210			#ifdef IS_MPI
211	tim	725	strcpy(checkPointMsg, "The integrator is ready to go.");
212	mmeineke	559	MPIcheckPoint();
213			#endif // is_mpi
214
215	tim	725	while (info->getTime() < runTime){
216			if ((info->getTime() + dt) >= currStatus){
217	mmeineke	558	calcPot = 1;
218			calcStress = 1;
219			}
220	mmeineke	561
221	chuckv	892	#ifdef PROFILE
222			startProfile( pro1 );
223			#endif
224
225	tim	725	integrateStep(calcPot, calcStress);
226
227	chuckv	892	#ifdef PROFILE
228			endProfile( pro1 );
229
230			startProfile( pro2 );
231			#endif // profile
232
233	mmeineke	643	info->incrTime(dt);
234	mmeineke	558
235	tim	725	if (info->setTemp){
236			if (info->getTime() >= currThermal){
237			thermalize();
238			currThermal += thermalTime;
239	mmeineke	558	}
240			}
241
242	tim	725	if (info->getTime() >= currSample){
243			dumpOut->writeDump(info->getTime());
244	mmeineke	558	currSample += sampleTime;
245			}
246
247	tim	725	if (info->getTime() >= currStatus){
248	tim	837	statOut->writeStat(info->getTime());
249			calcPot = 0;
250	mmeineke	558	calcStress = 0;
251			currStatus += statusTime;
252	tim	837	}
253	mmeineke	559
254	mmeineke	746	if (info->resetIntegrator){
255			if (info->getTime() >= currReset){
256			this->resetIntegrator();
257			currReset += resetTime;
258			}
259			}
260	chuckv	892
261			#ifdef PROFILE
262			endProfile( pro2 );
263			#endif //profile
264	mmeineke	746
265	mmeineke	559	#ifdef IS_MPI
266	tim	725	strcpy(checkPointMsg, "successfully took a time step.");
267	mmeineke	559	MPIcheckPoint();
268			#endif // is_mpi
269	mmeineke	558	}
270
271	mmeineke	561	delete dumpOut;
272			delete statOut;
273	mmeineke	558	}
274
275	tim	725	template<typename T> void Integrator<T>::integrateStep(int calcPot,
276			int calcStress){
277	mmeineke	558	// Position full step, and velocity half step
278	chuckv	892
279			#ifdef PROFILE
280			startProfile(pro3);
281			#endif //profile
282
283	tim	725	preMove();
284	mmeineke	558
285	chuckv	892	#ifdef PROFILE
286			endProfile(pro3);
287
288			startProfile(pro4);
289			#endif // profile
290
291	mmeineke	558	moveA();
292
293	chuckv	892	#ifdef PROFILE
294			endProfile(pro4);
295
296			startProfile(pro5);
297			#endif//profile
298	tim	725
299
300	mmeineke	614	#ifdef IS_MPI
301	tim	725	strcpy(checkPointMsg, "Succesful moveA\n");
302	mmeineke	614	MPIcheckPoint();
303			#endif // is_mpi
304
305	tim	725
306	mmeineke	558	// calc forces
307
308	tim	725	calcForce(calcPot, calcStress);
309	mmeineke	558
310	mmeineke	614	#ifdef IS_MPI
311	tim	725	strcpy(checkPointMsg, "Succesful doForces\n");
312	mmeineke	614	MPIcheckPoint();
313			#endif // is_mpi
314
315	chuckv	892	#ifdef PROFILE
316			endProfile( pro5 );
317	tim	725
318	chuckv	892	startProfile( pro6 );
319			#endif //profile
320
321	mmeineke	558	// finish the velocity half step
322	tim	725
323	mmeineke	558	moveB();
324	tim	725
325	chuckv	892	#ifdef PROFILE
326			endProfile(pro6);
327			#endif // profile
328	tim	725
329	mmeineke	614	#ifdef IS_MPI
330	tim	725	strcpy(checkPointMsg, "Succesful moveB\n");
331	mmeineke	614	MPIcheckPoint();
332			#endif // is_mpi
333	mmeineke	558	}
334
335
336	tim	725	template<typename T> void Integrator<T>::moveA(void){
337	gezelter	600	int i, j;
338	mmeineke	558	DirectionalAtom* dAtom;
339	gezelter	600	double Tb[3], ji[3];
340			double vel[3], pos[3], frc[3];
341			double mass;
342	mmeineke	558
343	tim	725	for (i = 0; i < nAtoms; i++){
344			atoms[i]->getVel(vel);
345			atoms[i]->getPos(pos);
346			atoms[i]->getFrc(frc);
347	mmeineke	567
348	gezelter	600	mass = atoms[i]->getMass();
349	mmeineke	594
350	tim	725	for (j = 0; j < 3; j++){
351	gezelter	600	// velocity half step
352	tim	725	vel[j] += (dt2 * frc[j] / mass) * eConvert;
353	gezelter	600	// position whole step
354			pos[j] += dt * vel[j];
355			}
356	mmeineke	594
357	tim	725	atoms[i]->setVel(vel);
358			atoms[i]->setPos(pos);
359	gezelter	600
360	tim	725	if (atoms[i]->isDirectional()){
361			dAtom = (DirectionalAtom *) atoms[i];
362	mmeineke	558
363			// get and convert the torque to body frame
364	mmeineke	597
365	tim	725	dAtom->getTrq(Tb);
366			dAtom->lab2Body(Tb);
367
368	mmeineke	558	// get the angular momentum, and propagate a half step
369	gezelter	600
370	tim	725	dAtom->getJ(ji);
371	gezelter	600
372	tim	725	for (j = 0; j < 3; j++)
373	gezelter	600	ji[j] += (dt2 * Tb[j]) * eConvert;
374	tim	725
375	mmeineke	778	this->rotationPropagation( dAtom, ji );
376	gezelter	600
377	tim	725	dAtom->setJ(ji);
378			}
379	mmeineke	558	}
380	mmeineke	768
381			if (nConstrained){
382			constrainA();
383			}
384	mmeineke	558	}
385
386
387	tim	725	template<typename T> void Integrator<T>::moveB(void){
388	gezelter	600	int i, j;
389	mmeineke	558	DirectionalAtom* dAtom;
390	gezelter	600	double Tb[3], ji[3];
391			double vel[3], frc[3];
392			double mass;
393	mmeineke	558
394	tim	725	for (i = 0; i < nAtoms; i++){
395			atoms[i]->getVel(vel);
396			atoms[i]->getFrc(frc);
397	mmeineke	558
398	gezelter	600	mass = atoms[i]->getMass();
399
400	mmeineke	558	// velocity half step
401	tim	725	for (j = 0; j < 3; j++)
402			vel[j] += (dt2 * frc[j] / mass) * eConvert;
403	gezelter	600
404	tim	725	atoms[i]->setVel(vel);
405	mmeineke	597
406	tim	725	if (atoms[i]->isDirectional()){
407			dAtom = (DirectionalAtom *) atoms[i];
408
409	tim	837	// get and convert the torque to body frame
410	gezelter	600
411	tim	725	dAtom->getTrq(Tb);
412			dAtom->lab2Body(Tb);
413	gezelter	600
414			// get the angular momentum, and propagate a half step
415
416	tim	725	dAtom->getJ(ji);
417	gezelter	600
418	tim	725	for (j = 0; j < 3; j++)
419	gezelter	600	ji[j] += (dt2 * Tb[j]) * eConvert;
420	mmeineke	597
421	tim	725
422			dAtom->setJ(ji);
423	mmeineke	558	}
424			}
425	mmeineke	768
426			if (nConstrained){
427			constrainB();
428			}
429	mmeineke	558	}
430
431	tim	725	template<typename T> void Integrator<T>::preMove(void){
432	gezelter	600	int i, j;
433			double pos[3];
434	mmeineke	558
435	tim	725	if (nConstrained){
436			for (i = 0; i < nAtoms; i++){
437			atoms[i]->getPos(pos);
438	mmeineke	561
439	tim	725	for (j = 0; j < 3; j++){
440			oldPos[3 * i + j] = pos[j];
441	gezelter	600	}
442			}
443	tim	725	}
444	gezelter	600	}
445
446	tim	645	template<typename T> void Integrator<T>::constrainA(){
447	mmeineke	787	int i, j;
448	mmeineke	558	int done;
449	gezelter	600	double posA[3], posB[3];
450			double velA[3], velB[3];
451	mmeineke	572	double pab[3];
452			double rab[3];
453	mmeineke	563	int a, b, ax, ay, az, bx, by, bz;
454	mmeineke	558	double rma, rmb;
455			double dx, dy, dz;
456	mmeineke	561	double rpab;
457	mmeineke	558	double rabsq, pabsq, rpabsq;
458			double diffsq;
459			double gab;
460			int iteration;
461
462	tim	725	for (i = 0; i < nAtoms; i++){
463	mmeineke	558	moving[i] = 0;
464	tim	725	moved[i] = 1;
465	mmeineke	558	}
466	mmeineke	567
467	mmeineke	558	iteration = 0;
468			done = 0;
469	tim	725	while (!done && (iteration < maxIteration)){
470	mmeineke	558	done = 1;
471	tim	725	for (i = 0; i < nConstrained; i++){
472	mmeineke	558	a = constrainedA[i];
473			b = constrainedB[i];
474	mmeineke	563
475	tim	725	ax = (a * 3) + 0;
476			ay = (a * 3) + 1;
477			az = (a * 3) + 2;
478	mmeineke	563
479	tim	725	bx = (b * 3) + 0;
480			by = (b * 3) + 1;
481			bz = (b * 3) + 2;
482
483			if (moved[a] \|\| moved[b]){
484			atoms[a]->getPos(posA);
485			atoms[b]->getPos(posB);
486
487			for (j = 0; j < 3; j++)
488	gezelter	600	pab[j] = posA[j] - posB[j];
489	mmeineke	567
490	tim	725	//periodic boundary condition
491	mmeineke	567
492	tim	725	info->wrapVector(pab);
493	mmeineke	572
494	tim	725	pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
495	mmeineke	558
496	tim	725	rabsq = constrainedDsqr[i];
497			diffsq = rabsq - pabsq;
498	mmeineke	567
499	tim	725	// the original rattle code from alan tidesley
500			if (fabs(diffsq) > (tol * rabsq * 2)){
501			rab[0] = oldPos[ax] - oldPos[bx];
502			rab[1] = oldPos[ay] - oldPos[by];
503			rab[2] = oldPos[az] - oldPos[bz];
504	mmeineke	558
505	tim	725	info->wrapVector(rab);
506	mmeineke	567
507	tim	725	rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
508	mmeineke	558
509	tim	725	rpabsq = rpab * rpab;
510	mmeineke	558
511	mmeineke	563
512	tim	725	if (rpabsq < (rabsq * -diffsq)){
513	mmeineke	558	#ifdef IS_MPI
514	tim	725	a = atoms[a]->getGlobalIndex();
515			b = atoms[b]->getGlobalIndex();
516	mmeineke	558	#endif //is_mpi
517	tim	725	sprintf(painCave.errMsg,
518			"Constraint failure in constrainA at atom %d and %d.\n", a,
519			b);
520			painCave.isFatal = 1;
521			simError();
522			}
523	mmeineke	558
524	tim	725	rma = 1.0 / atoms[a]->getMass();
525			rmb = 1.0 / atoms[b]->getMass();
526	mmeineke	567
527	tim	725	gab = diffsq / (2.0 * (rma + rmb) * rpab);
528	mmeineke	567
529	mmeineke	572	dx = rab[0] * gab;
530			dy = rab[1] * gab;
531			dz = rab[2] * gab;
532	mmeineke	558
533	tim	725	posA[0] += rma * dx;
534			posA[1] += rma * dy;
535			posA[2] += rma * dz;
536	mmeineke	558
537	tim	725	atoms[a]->setPos(posA);
538	mmeineke	558
539	tim	725	posB[0] -= rmb * dx;
540			posB[1] -= rmb * dy;
541			posB[2] -= rmb * dz;
542	gezelter	600
543	tim	725	atoms[b]->setPos(posB);
544	gezelter	600
545	mmeineke	558	dx = dx / dt;
546			dy = dy / dt;
547			dz = dz / dt;
548
549	tim	725	atoms[a]->getVel(velA);
550	mmeineke	558
551	tim	725	velA[0] += rma * dx;
552			velA[1] += rma * dy;
553			velA[2] += rma * dz;
554	mmeineke	558
555	tim	725	atoms[a]->setVel(velA);
556	gezelter	600
557	tim	725	atoms[b]->getVel(velB);
558	gezelter	600
559	tim	725	velB[0] -= rmb * dx;
560			velB[1] -= rmb * dy;
561			velB[2] -= rmb * dz;
562	gezelter	600
563	tim	725	atoms[b]->setVel(velB);
564	gezelter	600
565	tim	725	moving[a] = 1;
566			moving[b] = 1;
567			done = 0;
568			}
569	mmeineke	558	}
570			}
571	tim	725
572			for (i = 0; i < nAtoms; i++){
573	mmeineke	558	moved[i] = moving[i];
574			moving[i] = 0;
575			}
576
577			iteration++;
578			}
579
580	tim	725	if (!done){
581			sprintf(painCave.errMsg,
582			"Constraint failure in constrainA, too many iterations: %d\n",
583			iteration);
584	mmeineke	558	painCave.isFatal = 1;
585			simError();
586			}
587	mmeineke	768
588	mmeineke	558	}
589
590	tim	725	template<typename T> void Integrator<T>::constrainB(void){
591	mmeineke	787	int i, j;
592	mmeineke	558	int done;
593	gezelter	600	double posA[3], posB[3];
594			double velA[3], velB[3];
595	mmeineke	558	double vxab, vyab, vzab;
596	mmeineke	572	double rab[3];
597	mmeineke	563	int a, b, ax, ay, az, bx, by, bz;
598	mmeineke	558	double rma, rmb;
599			double dx, dy, dz;
600	mmeineke	787	double rvab;
601	mmeineke	558	double gab;
602			int iteration;
603
604	tim	725	for (i = 0; i < nAtoms; i++){
605	mmeineke	558	moving[i] = 0;
606			moved[i] = 1;
607			}
608
609			done = 0;
610	mmeineke	561	iteration = 0;
611	tim	725	while (!done && (iteration < maxIteration)){
612	mmeineke	567	done = 1;
613
614	tim	725	for (i = 0; i < nConstrained; i++){
615	mmeineke	558	a = constrainedA[i];
616			b = constrainedB[i];
617
618	tim	725	ax = (a * 3) + 0;
619			ay = (a * 3) + 1;
620			az = (a * 3) + 2;
621	mmeineke	563
622	tim	725	bx = (b * 3) + 0;
623			by = (b * 3) + 1;
624			bz = (b * 3) + 2;
625	mmeineke	563
626	tim	725	if (moved[a] \|\| moved[b]){
627			atoms[a]->getVel(velA);
628			atoms[b]->getVel(velB);
629	mmeineke	558
630	tim	725	vxab = velA[0] - velB[0];
631			vyab = velA[1] - velB[1];
632			vzab = velA[2] - velB[2];
633	gezelter	600
634	tim	725	atoms[a]->getPos(posA);
635			atoms[b]->getPos(posB);
636	gezelter	600
637	tim	725	for (j = 0; j < 3; j++)
638	gezelter	600	rab[j] = posA[j] - posB[j];
639	mmeineke	558
640	tim	725	info->wrapVector(rab);
641	mmeineke	558
642	tim	725	rma = 1.0 / atoms[a]->getMass();
643			rmb = 1.0 / atoms[b]->getMass();
644	mmeineke	558
645	tim	725	rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
646	gezelter	600
647	tim	725	gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);
648	gezelter	600
649	tim	725	if (fabs(gab) > tol){
650			dx = rab[0] * gab;
651			dy = rab[1] * gab;
652			dz = rab[2] * gab;
653
654			velA[0] += rma * dx;
655			velA[1] += rma * dy;
656			velA[2] += rma * dz;
657
658			atoms[a]->setVel(velA);
659
660			velB[0] -= rmb * dx;
661			velB[1] -= rmb * dy;
662			velB[2] -= rmb * dz;
663
664			atoms[b]->setVel(velB);
665
666			moving[a] = 1;
667			moving[b] = 1;
668			done = 0;
669			}
670	mmeineke	558	}
671			}
672
673	tim	725	for (i = 0; i < nAtoms; i++){
674	mmeineke	558	moved[i] = moving[i];
675			moving[i] = 0;
676			}
677	tim	725
678	mmeineke	558	iteration++;
679			}
680
681	tim	725	if (!done){
682			sprintf(painCave.errMsg,
683			"Constraint failure in constrainB, too many iterations: %d\n",
684			iteration);
685	mmeineke	558	painCave.isFatal = 1;
686			simError();
687	tim	725	}
688	mmeineke	558	}
689
690	mmeineke	778	template<typename T> void Integrator<T>::rotationPropagation
691			( DirectionalAtom* dAtom, double ji[3] ){
692
693			double angle;
694			double A[3][3], I[3][3];
695
696			// use the angular velocities to propagate the rotation matrix a
697			// full time step
698
699			dAtom->getA(A);
700			dAtom->getI(I);
701	tim	837
702			// rotate about the x-axis
703	mmeineke	778	angle = dt2 * ji[0] / I[0][0];
704	tim	837	this->rotate( 1, 2, angle, ji, A );
705
706	mmeineke	778	// rotate about the y-axis
707			angle = dt2 * ji[1] / I[1][1];
708			this->rotate( 2, 0, angle, ji, A );
709	tim	837
710	mmeineke	778	// rotate about the z-axis
711			angle = dt * ji[2] / I[2][2];
712			this->rotate( 0, 1, angle, ji, A);
713	tim	837
714	mmeineke	778	// rotate about the y-axis
715			angle = dt2 * ji[1] / I[1][1];
716			this->rotate( 2, 0, angle, ji, A );
717	tim	837
718	mmeineke	778	// rotate about the x-axis
719			angle = dt2 * ji[0] / I[0][0];
720			this->rotate( 1, 2, angle, ji, A );
721	tim	837
722			dAtom->setA( A );
723	mmeineke	778	}
724
725	tim	725	template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
726			double angle, double ji[3],
727			double A[3][3]){
728			int i, j, k;
729	mmeineke	558	double sinAngle;
730			double cosAngle;
731			double angleSqr;
732			double angleSqrOver4;
733			double top, bottom;
734			double rot[3][3];
735			double tempA[3][3];
736			double tempJ[3];
737
738			// initialize the tempA
739
740	tim	725	for (i = 0; i < 3; i++){
741			for (j = 0; j < 3; j++){
742	gezelter	600	tempA[j][i] = A[i][j];
743	mmeineke	558	}
744			}
745
746			// initialize the tempJ
747
748	tim	725	for (i = 0; i < 3; i++)
749			tempJ[i] = ji[i];
750
751	mmeineke	558	// initalize rot as a unit matrix
752
753			rot[0][0] = 1.0;
754			rot[0][1] = 0.0;
755			rot[0][2] = 0.0;
756
757			rot[1][0] = 0.0;
758			rot[1][1] = 1.0;
759			rot[1][2] = 0.0;
760	tim	725
761	mmeineke	558	rot[2][0] = 0.0;
762			rot[2][1] = 0.0;
763			rot[2][2] = 1.0;
764	tim	725
765	mmeineke	558	// use a small angle aproximation for sin and cosine
766
767	tim	725	angleSqr = angle * angle;
768	mmeineke	558	angleSqrOver4 = angleSqr / 4.0;
769			top = 1.0 - angleSqrOver4;
770			bottom = 1.0 + angleSqrOver4;
771
772			cosAngle = top / bottom;
773			sinAngle = angle / bottom;
774
775			rot[axes1][axes1] = cosAngle;
776			rot[axes2][axes2] = cosAngle;
777
778			rot[axes1][axes2] = sinAngle;
779			rot[axes2][axes1] = -sinAngle;
780	tim	725
781	mmeineke	558	// rotate the momentum acoording to: ji[] = rot[][] * ji[]
782	tim	725
783			for (i = 0; i < 3; i++){
784	mmeineke	558	ji[i] = 0.0;
785	tim	725	for (k = 0; k < 3; k++){
786	mmeineke	558	ji[i] += rot[i][k] * tempJ[k];
787			}
788			}
789
790	tim	837	// rotate the Rotation matrix acording to:
791	mmeineke	558	// A[][] = A[][] * transpose(rot[][])
792
793
794	mmeineke	561	// NOte for as yet unknown reason, we are performing the
795	mmeineke	558	// calculation as:
796			// transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
797
798	tim	725	for (i = 0; i < 3; i++){
799			for (j = 0; j < 3; j++){
800	gezelter	600	A[j][i] = 0.0;
801	tim	725	for (k = 0; k < 3; k++){
802			A[j][i] += tempA[i][k] * rot[j][k];
803	mmeineke	558	}
804			}
805			}
806			}
807	tim	677
808	tim	725	template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
809			myFF->doForces(calcPot, calcStress);
810	tim	677	}
811
812			template<typename T> void Integrator<T>::thermalize(){
813	tim	725	tStats->velocitize();
814	tim	677	}
815	tim	763
816			template<typename T> double Integrator<T>::getConservedQuantity(void){
817			return tStats->getTotalE();
818	mmeineke	768	}
819	tim	837	template<typename T> string Integrator<T>::getAdditionalParameters(void){
820			//By default, return a null string
821			//The reason we use string instead of char* is that if we use char*, we will
822			//return a pointer point to local variable which might cause problem
823			return string();
824			}