src/constraints/ZconstraintForceManager.cpp

#include <cmath>
#include "constraints/ZconstraintForceManager.hpp"
#include "integrators/Integrator.hpp"
#include "utils/simError.h"
#include "utils/OOPSEConstant.hpp"

namespace oopse {
ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info) {
    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
    Globals* simParam = info_->getSimParams();

    if (simParam->haveDt()){
        dt_ = simParam->getDt();
    } else {
        sprintf(painCave.errMsg,
                "Integrator Error: dt is not set\n");
        painCave.isFatal = 1;
        simError();    
    }

    if (simParam->haveZconstraintTime()){
        zconsTime_ = simParam->getZconsTime();
    }
    else{
        sprintf(painCave.errMsg,
        "ZConstraint error: If you use a ZConstraint,\n"
        "\tyou must set zconsTime.\n");
        painCave.isFatal = 1;
        simError();
    }

    if (simParam->haveZconsTol()){
        zconsTol_ = simParam->getZconsTol();
    }
    else{
        zconsTol_ = 0.01;
        sprintf(painCave.errMsg,
            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
            "\tOOPSE will use a default value of %f.\n"
            "\tTo set the tolerance, use the zconsTol variable.\n",
            zconsTol_);
        painCave.isFatal = 0;
        simError();      
    }

    //set zcons gap
    if (simParam->haveZConsGap()){
        zconsGap_ = simParam->getZconsGap();
    }else {
        zconsGap_ = false;
    }

    //set zcons fixtime
    if (simParam->haveZConsFixTime()){
        zconsFixingTime_ = simParam->getZconsFixtime();
    } else {
        zconsFixingTime_ = infiniteTime;
    }

    //set zconsUsingSMD
    if (simParam->haveZConsUsingSMD()){
        usingSMD_ = simParam->getZconsUsingSMD();
    }else {
        usingSMD_ =false;
    }
    
    std::string zconsOutput = getPrefix(info_->getFinalConfigFileName()) + ".fz";

    //estimate the force constant of harmonical potential
    Mat3x3d hmat = currSnapshot_->getHmat();
    double halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2;
    double zforceConstant = OOPSEConstant::kb * info->target_temp / (halfOfLargestBox * halfOfLargestBox);
         
    int nZConstraints = simParam->getNzConstraints();
    ZconStamp** stamp = simParam->getZconStamp();
    //
    for (int i = 0; i < nZConstraints; i++){

        ZconstraintParam param;
        int zmolIndex = stamp[i]->getMolIndex();
        if (stamp[i]->haveZpos()) {
            param.zTargetPos = getZTargetPos(zmolIndex);
        } else {
            param.zTargetPos = stamp[i]->getZpos();
        }

        param.kz = zforceConstant * stamp[i]->getKratio();

        if (stamp[i]->haveCantVel()) {
            param.cantVel = stamp[i]->getCantVel();
        } else {
            param.cantVel = 0.0;
        }

        allZMolIndices_.insert(std::make_pair(zmolIndex, param));
    }

    //create fixedMols_, movingMols_ and unconsMols lists 
    update();
    
    //calculate masss of unconstraint molecules in the whole system (never change during the simulation)
    double totMassUnconsMols_local = 0.0;    
    std::vector<Molecule*>::iterator j;
    for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
        totMassUnconsMols_local += (*j)->getMass();
    }    
#ifndef IS_MPI
    totMassUnconsMols_ = totMassUnconsMols_local;
#else
    MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_DOUBLE,
        MPI_SUM, MPI_COMM_WORLD);  
#endif


}

ZconstraintForceManager::~ZConstraint(){

  if (fzOut){
    delete fzOut;
  }

}

void ZconstraintForceManager::update(){
    fixedZMols_.clear();
    movingZMols_.clear();
    unzconsMols_.clear();
    
    std::map<int, ZconstraintParam>::iterator i;
    for (i = allZMolIndices_.begin(); i != allZMolIndices_.end(); ++i) {
#ifdef IS_MPI
        if (info_->getMolToProc(i->first) == worldRank) {
#endif
            ZconstraintMol zmol;
            zmol->mol = info_->getMoleculeByGlobalIndex(i->first);
            assert(zmol->mol);
            zmol->param = info_->getMoleculeByGlobalIndex(i->second);
            zmol->cantPos = zmol->param->zTargetPos; /**@todo fixed me when zmol migrate, it is incorrect*/
            Vector3d com = zmol->mol->getCom();
            double diff = fabs(zmol->param->zTargetPos - com[whichDirection]);
            if (diff < zconsTol_) {
                fixedZMols_.push_back(zmol);
            } else {
                movingZMols_.push_back(zmol);            
            }

#ifdef IS_MPI
        }
#endif
    }

    calcTotalMassMovingZMols();

    std::set<int> zmolSet;
    std::list<ZconstraintMol>::iterator i;
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        zmolSet.insert(i->mol->getGlobalIndex());
    }    

    for ( i = fixedZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        zmolSet.insert(i->mol->getGlobalIndex());
    }

    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    for(mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
        if (zmolSet->find(mol->getGlobalIndex()) == zmolSet.end()) {
            unzconsMols_.push_back(zmolSet);
        }
    }

}

bool ZconstraintForceManager::isZMol(Molecule* mol){
    return allZMolIndices.find(mol->getGlobalIndex()) == allZMolIndices_.end() ? false : true;
}

void ZconstraintForceManager::integrate(){
  // creat zconsWriter  
  fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);   

  if (!fzOut){
    sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
    painCave.isFatal = 1;
    simError();
  }

  //zero out the velocities of center of mass of unconstrained molecules 
  //and the velocities of center of mass of every single z-constrained molecueles
  zeroVelocity();

  curZconsTime = zconsTime + info->getTime();

  T::integrate();
}

void ZconstraintForceManager::calcForce(bool needPotential, bool needStress){
    ForceManager::calcForce(needPotential, needStress);
    
    if (usingZconsGap_){
        updateZPos();
    }

    if (checkZConsState()){
        zeroVelocity();    
        forcePolicy->update();
    }  

    //do zconstraint force; 
    if (haveFixedZMols()){
        doZconstraintForce();
    }

    //use external force to move the molecules to the specified positions
    if (haveMovingZMols()){
        doHarmonic();
    }

    //write out forces and current positions of z-constraint molecules    
    if (currSnapshot_->getTime() >= curZconsTime_){
        std::list<ZconstraintMol>::iterator i;
        Vector3d com;
        for(i = fixedZMols_.begin(); i != fixedZMols_.end() ++i) {
            com = i->mol->getCom();
            i->zpos = com[whichDirection];
        }
        
        fzOut->writeFZ(info->getTime(), fixedZMols_);
        curZconsTime_ += zconsTime_;
    }
}

void ZconstraintForceManager::zeroVelocity(){

    Vector3d comVel;
    Vector3d vel;
    std::list<ZconstraintMol>::iterator i;
    Molecule* mol;
    StuntDouble* integrableObject;
    Molecule::IntegrableObjectIteraot ii;

    //zero out the velocities of center of mass of fixed z-constrained molecules
    for(i = fixedZMols_.begin(); i != fixedZMols_.end() ++i) {
        mol = i->mol;
        comVel = mol->getComVel();
        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {
            vel[whichDirection] -= comVel[whichDirection];
            integrableObject->setVel(vel);
        }
    }

    // calculate the vz of center of mass of moving molecules(include unconstrained molecules 
    // and moving z-constrained molecules)  
    double pzMovingMols_local;
    double pzMovingMols;
    
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        mol = i->mol;        
        comVel = mol->getComVel();
        pzMovingMols_local +=  mol->getMass() * comVel[whichDirection];   
    }

    std::vector<Molecule*>::iterator j;
    for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
        mol =*j;
        comVel = mol->getCoMVel();
        pzMovingMols_local += mol->getMass() * comVel[whichDirection];
    }
    
#ifndef IS_MPI
    pzMovingMols = pzMovingMols_local;
#else
    MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_DOUBLE,
        MPI_SUM, MPI_COMM_WORLD);
#endif

    double vzMovingMols = pzMovingMols / totMassMovingMols_;

    //modify the velocities of moving z-constrained molecuels
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        mol = i->mol;
        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {

            vel = integrableObject->getVel();
            vel[whichDirection] -= vzOfMovingMols;
            integrableObject->setVel(vel); 
        }
    }

    //modify the velocites of unconstrained molecules  
    for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
        mol =*j;
        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {

            vel = integrableObject->getVel();
            vel[whichDirection] -= vzOfMovingMols;
            integrableObject->setVel(vel); 
        }
    }
    
}


void ZconstraintForceManager::doZconstraintForce(){
  Atom** zconsAtoms;
  double totalFZ; 
  double totalFZ_local;
  Vector3d com;
  Vector3d force(0.0);

  //constrain the molecules which do not reach the specified positions  

  //Zero Out the force of z-contrained molecules    
  totalFZ_local = 0;


    //calculate the total z-contrained force of fixed z-contrained molecules
    std::list<ZconstraintMol>::iterator i;
    Molecule* mol;
    StuntDouble* integrableObject;
    Molecule::IntegrableObjectIteraot ii;

    for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
        mol = i->mol;
        com = mol->getCom();        

        i->fz = 0.0;

        totalFZ_local += harmonicF;

        //adjust force
        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {

            force = integrableObject->getFrc();    
            i->fz += force[whichDirection] 
        }
         totalFZ_local += i->fz;
    }

  //calculate total z-constraint force
#ifdef IS_MPI
  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#else
  totalFZ = totalFZ_local;
#endif


     // apply negative to fixed z-constrained molecues;
    for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
        mol = i->mol;
        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {

            force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz);
            integrableObject->addFrc(force);
        }
    }

  //modify the forces of moving z-constrained molecules
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        mol = i->mol;
        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {

            force[whichDirection] = -getZFOfMovingMols(mol, integrableObject, totalFZ);
            integrableObject->addFrc(force);
        }
    }

  //modify the forces of unconstrained molecules
    std::vector<Molecule*>::iterator j;
    for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
        mol =*j;
        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {

            force[whichDirection] = -getZFOfMovingMols(mol, integrableObject, totalFZ);
            integrableObject->addFrc(force);
        }
    }

}


void ZconstraintForceManager::doHarmonic(){
    double harmonicU;
    double harmonicF;
    double diff;
    double totalFZ_local;
    double totalFZ;

    Vector3d force(0.0);
    Vector3d com;

    double totalFZ_local = 0;

    std::list<ZconstraintMol>::iterator i;
    Molecule* mol;
     StuntDouble* integrableObject;
    Molecule::IntegrableObjectIteraot ii;
    
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        mol = i->mol;
        com = mol->getCom();        
        diff = com[whichDirection] - resPos; 
        harmonicU = 0.5 * kz * diff * diff;
        currSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
        harmonicF = -kz * diff;
        totalFZ_local += harmonicF;

        //adjust force
        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {

            force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF);
            integrableObject->addFrc(force);            
        }
    }

#ifndef IS_MPI
  totalFZ = totalFZ_local;
#else
  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
#endif

    //modify the forces of unconstrained molecules
    std::vector<Molecule*>::iterator j;
    for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
        mol = *j;
        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; 
            integrableObject = mol->nextIntegrableObject(ii)) {

            force[whichDirection] = getHFOfUnconsMols(mol, harmonicF);
            integrableObject->addFrc(force);            
        }
    }

}

bool ZconstraintForceManager::checkZConsState(){
    Vector3d com;
    double diff;
    int changed_local = 0;

    std::list<ZconstraintMol>::iterator i;
    std::list<ZconstraintMol>::iterator j;
    
    std::list<ZconstraintMol> newMovingZMols;
    for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end();) {
        com = i->mol->getCom();
        diff = fabs(com[whichDirection] - i->param.zTargetPos);
        if (diff > zconsTol_) {
            //this fixed zconstraint molecule is about to moving
            //moved this molecule to 
            j = i++;
            newMovingZMols.push_back(fixedZMols_.erase(j));
            if (usingZconsGap_) {
                i->endFixingTime = infiniteTime;
            }
            changed_local = 1;            
        }else {
            ++i;
        }
    }  

    std::list<ZconstraintMol> newFixedZMols;
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end();) {
        com = i->mol->getCom();
        diff = fabs(com[whichDirection] - i->param.zTargetPos);
        if (diff <= zconsTol_) {
            //this moving zconstraint molecule is about to fixed
            //moved this molecule to 
            j = i++;
            newFixedZMols.push_back(movingZMols_.erase(j));
            if (usingZconsGap_) {
                i->endFixingTime = currSnapshot_->getTime() + zconsFixingTime_;
            }
            changed_local = 1;            
        }else {
            ++i;
        }
    }     

    //merge the lists
    fixedZMols_.merge(newFixedZMols);
    movingZMols_.merge(newFixedZMols);

    int changed;
#ifndef IS_MPI
    changed = changed_local; 
#else
    MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#endif

    return (changed > 0);
}

bool ZconstraintForceManager::haveFixedZMols(){
  int havingFixed;
  int havingFixed_local = fixedZMols_.empty() ? 0 : 1;

#ifndef IS_MPI
  havingFixed = havingFixed_local;
#else
  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
                MPI_COMM_WORLD);
#endif

  return havingFixed > 0;
}


bool ZconstraintForceManager::haveMovingZMols(){
  int havingMoving_local;
  int havingMoving;

  havingMoving_local = movingZMols_.empty()? 0 : 1;

#ifndef IS_MPI
  havingMoving = havingMoving_local;
#else
  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
                MPI_COMM_WORLD);
#endif

  return havingMoving > 0;
}

void ZconstraintForceManager::calcTotalMassMovingZMols(){

    double totMassMovingZMols_local = 0.0;
    std::list<ZconstraintMol>::iterator i;
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        totMassMovingZMols_local += i->mol->getMass();
    }
    
#ifdef IS_MPI
    MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_DOUBLE,
                MPI_SUM, MPI_COMM_WORLD);
#else
    totMassMovingZMols_ = massOfMovingZAtoms_local;
#endif

}

double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){
  return totalForce * sd->getMass() / mol->getMass();
}

double ZconstraintForceManager::getZFOfMovingMols(StuntDouble*sd, double totalForce){
  return totalForce * sd->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
}

double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){
  return totalForce * sd->getMass() / mol->getMass();
}

double ZconstraintForceManager::getHFOfUnconsMols(StuntDouble*sd, double totalForce){
  return totalForce * sd->getMass() / totMassUnconsMols_;
}

void ZconstraintForceManager::updateZPos(){
    double curTime = currSnapshot_->getTime();
    std::list<ZconstraintMol>::iterator i;
    for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
        i->param.zTargetPos += zconsGap_;     
    }  
}

void ZconstraintForceManager::updateCantPos(){
    std::list<ZconstraintMol>::iterator i;
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        i->cantPos += i->param.cantVel * dt_;
    }
}

double ZconstraintForceManager::getZTargetPos(int index){
    double zTargetPos;
#ifndef IS_MPI    
    Molecule* mol = info_->getMoleculeByGlobalIndex(index);
    assert(!mol);
    Vector3d com = mol->getCom();
    zTargetPos = com[whichDirection];
#else
    int whicProc = info_->getMolToProc(index);
    MPI_Bcast(&zTargetPos, 1, MPI_DOUBLE, whicProc, MPI_COMM_WORLD);
#endif
    return zTargetPos;
}

}
Revision:	1910
Committed:	Fri Jan 7 21:50:13 2005 UTC (19 years, 8 months ago) by tim
File size:	18107 byte(s)
Log Message:	ZConstraintForceManager in progress
#	User	Rev	Content
1	tim	1910	#include <cmath>
2			#include "constraints/ZconstraintForceManager.hpp"
3			#include "integrators/Integrator.hpp"
4			#include "utils/simError.h"
5			#include "utils/OOPSEConstant.hpp"
6
7			namespace oopse {
8			ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info) {
9			currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
10			Globals* simParam = info_->getSimParams();
11
12			if (simParam->haveDt()){
13			dt_ = simParam->getDt();
14			} else {
15			sprintf(painCave.errMsg,
16			"Integrator Error: dt is not set\n");
17			painCave.isFatal = 1;
18			simError();
19			}
20
21			if (simParam->haveZconstraintTime()){
22			zconsTime_ = simParam->getZconsTime();
23			}
24			else{
25			sprintf(painCave.errMsg,
26			"ZConstraint error: If you use a ZConstraint,\n"
27			"\tyou must set zconsTime.\n");
28			painCave.isFatal = 1;
29			simError();
30			}
31
32			if (simParam->haveZconsTol()){
33			zconsTol_ = simParam->getZconsTol();
34			}
35			else{
36			zconsTol_ = 0.01;
37			sprintf(painCave.errMsg,
38			"ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
39			"\tOOPSE will use a default value of %f.\n"
40			"\tTo set the tolerance, use the zconsTol variable.\n",
41			zconsTol_);
42			painCave.isFatal = 0;
43			simError();
44			}
45
46			//set zcons gap
47			if (simParam->haveZConsGap()){
48			zconsGap_ = simParam->getZconsGap();
49			}else {
50			zconsGap_ = false;
51			}
52
53			//set zcons fixtime
54			if (simParam->haveZConsFixTime()){
55			zconsFixingTime_ = simParam->getZconsFixtime();
56			} else {
57			zconsFixingTime_ = infiniteTime;
58			}
59
60			//set zconsUsingSMD
61			if (simParam->haveZConsUsingSMD()){
62			usingSMD_ = simParam->getZconsUsingSMD();
63			}else {
64			usingSMD_ =false;
65			}
66
67			std::string zconsOutput = getPrefix(info_->getFinalConfigFileName()) + ".fz";
68
69			//estimate the force constant of harmonical potential
70			Mat3x3d hmat = currSnapshot_->getHmat();
71			double halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2;
72			double zforceConstant = OOPSEConstant::kb * info->target_temp / (halfOfLargestBox * halfOfLargestBox);
73
74			int nZConstraints = simParam->getNzConstraints();
75			ZconStamp** stamp = simParam->getZconStamp();
76			//
77			for (int i = 0; i < nZConstraints; i++){
78
79			ZconstraintParam param;
80			int zmolIndex = stamp[i]->getMolIndex();
81			if (stamp[i]->haveZpos()) {
82			param.zTargetPos = getZTargetPos(zmolIndex);
83			} else {
84			param.zTargetPos = stamp[i]->getZpos();
85			}
86
87			param.kz = zforceConstant * stamp[i]->getKratio();
88
89			if (stamp[i]->haveCantVel()) {
90			param.cantVel = stamp[i]->getCantVel();
91			} else {
92			param.cantVel = 0.0;
93			}
94
95			allZMolIndices_.insert(std::make_pair(zmolIndex, param));
96			}
97
98			//create fixedMols_, movingMols_ and unconsMols lists
99			update();
100
101			//calculate masss of unconstraint molecules in the whole system (never change during the simulation)
102			double totMassUnconsMols_local = 0.0;
103			std::vector<Molecule*>::iterator j;
104			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
105			totMassUnconsMols_local += (*j)->getMass();
106			}
107			#ifndef IS_MPI
108			totMassUnconsMols_ = totMassUnconsMols_local;
109			#else
110			MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_DOUBLE,
111			MPI_SUM, MPI_COMM_WORLD);
112			#endif
113
114
115			}
116
117			ZconstraintForceManager::~ZConstraint(){
118
119			if (fzOut){
120			delete fzOut;
121			}
122
123			}
124
125			void ZconstraintForceManager::update(){
126			fixedZMols_.clear();
127			movingZMols_.clear();
128			unzconsMols_.clear();
129
130			std::map<int, ZconstraintParam>::iterator i;
131			for (i = allZMolIndices_.begin(); i != allZMolIndices_.end(); ++i) {
132			#ifdef IS_MPI
133			if (info_->getMolToProc(i->first) == worldRank) {
134			#endif
135			ZconstraintMol zmol;
136			zmol->mol = info_->getMoleculeByGlobalIndex(i->first);
137			assert(zmol->mol);
138			zmol->param = info_->getMoleculeByGlobalIndex(i->second);
139			zmol->cantPos = zmol->param->zTargetPos; /*@todo fixed me when zmol migrate, it is incorrect/
140			Vector3d com = zmol->mol->getCom();
141			double diff = fabs(zmol->param->zTargetPos - com[whichDirection]);
142			if (diff < zconsTol_) {
143			fixedZMols_.push_back(zmol);
144			} else {
145			movingZMols_.push_back(zmol);
146			}
147
148			#ifdef IS_MPI
149			}
150			#endif
151			}
152
153			calcTotalMassMovingZMols();
154
155			std::set<int> zmolSet;
156			std::list<ZconstraintMol>::iterator i;
157			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
158			zmolSet.insert(i->mol->getGlobalIndex());
159			}
160
161			for ( i = fixedZMols_.begin(); i != movingZMols_.end(); ++i) {
162			zmolSet.insert(i->mol->getGlobalIndex());
163			}
164
165			SimInfo::MoleculeIterator mi;
166			Molecule* mol;
167			for(mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
168			if (zmolSet->find(mol->getGlobalIndex()) == zmolSet.end()) {
169			unzconsMols_.push_back(zmolSet);
170			}
171			}
172
173			}
174
175			bool ZconstraintForceManager::isZMol(Molecule* mol){
176			return allZMolIndices.find(mol->getGlobalIndex()) == allZMolIndices_.end() ? false : true;
177			}
178
179			void ZconstraintForceManager::integrate(){
180			// creat zconsWriter
181			fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);
182
183			if (!fzOut){
184			sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
185			painCave.isFatal = 1;
186			simError();
187			}
188
189			//zero out the velocities of center of mass of unconstrained molecules
190			//and the velocities of center of mass of every single z-constrained molecueles
191			zeroVelocity();
192
193			curZconsTime = zconsTime + info->getTime();
194
195			T::integrate();
196			}
197
198			void ZconstraintForceManager::calcForce(bool needPotential, bool needStress){
199			ForceManager::calcForce(needPotential, needStress);
200
201			if (usingZconsGap_){
202			updateZPos();
203			}
204
205			if (checkZConsState()){
206			zeroVelocity();
207			forcePolicy->update();
208			}
209
210			//do zconstraint force;
211			if (haveFixedZMols()){
212			doZconstraintForce();
213			}
214
215			//use external force to move the molecules to the specified positions
216			if (haveMovingZMols()){
217			doHarmonic();
218			}
219
220			//write out forces and current positions of z-constraint molecules
221			if (currSnapshot_->getTime() >= curZconsTime_){
222			std::list<ZconstraintMol>::iterator i;
223			Vector3d com;
224			for(i = fixedZMols_.begin(); i != fixedZMols_.end() ++i) {
225			com = i->mol->getCom();
226			i->zpos = com[whichDirection];
227			}
228
229			fzOut->writeFZ(info->getTime(), fixedZMols_);
230			curZconsTime_ += zconsTime_;
231			}
232			}
233
234			void ZconstraintForceManager::zeroVelocity(){
235
236			Vector3d comVel;
237			Vector3d vel;
238			std::list<ZconstraintMol>::iterator i;
239			Molecule* mol;
240			StuntDouble* integrableObject;
241			Molecule::IntegrableObjectIteraot ii;
242
243			//zero out the velocities of center of mass of fixed z-constrained molecules
244			for(i = fixedZMols_.begin(); i != fixedZMols_.end() ++i) {
245			mol = i->mol;
246			comVel = mol->getComVel();
247			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
248			integrableObject = mol->nextIntegrableObject(ii)) {
249			vel[whichDirection] -= comVel[whichDirection];
250			integrableObject->setVel(vel);
251			}
252			}
253
254			// calculate the vz of center of mass of moving molecules(include unconstrained molecules
255			// and moving z-constrained molecules)
256			double pzMovingMols_local;
257			double pzMovingMols;
258
259			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
260			mol = i->mol;
261			comVel = mol->getComVel();
262			pzMovingMols_local += mol->getMass() * comVel[whichDirection];
263			}
264
265			std::vector<Molecule*>::iterator j;
266			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
267			mol =*j;
268			comVel = mol->getCoMVel();
269			pzMovingMols_local += mol->getMass() * comVel[whichDirection];
270			}
271
272			#ifndef IS_MPI
273			pzMovingMols = pzMovingMols_local;
274			#else
275			MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_DOUBLE,
276			MPI_SUM, MPI_COMM_WORLD);
277			#endif
278
279			double vzMovingMols = pzMovingMols / totMassMovingMols_;
280
281			//modify the velocities of moving z-constrained molecuels
282			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
283			mol = i->mol;
284			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
285			integrableObject = mol->nextIntegrableObject(ii)) {
286
287			vel = integrableObject->getVel();
288			vel[whichDirection] -= vzOfMovingMols;
289			integrableObject->setVel(vel);
290			}
291			}
292
293			//modify the velocites of unconstrained molecules
294			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
295			mol =*j;
296			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
297			integrableObject = mol->nextIntegrableObject(ii)) {
298
299			vel = integrableObject->getVel();
300			vel[whichDirection] -= vzOfMovingMols;
301			integrableObject->setVel(vel);
302			}
303			}
304
305			}
306
307
308			void ZconstraintForceManager::doZconstraintForce(){
309			Atom** zconsAtoms;
310			double totalFZ;
311			double totalFZ_local;
312			Vector3d com;
313			Vector3d force(0.0);
314
315			//constrain the molecules which do not reach the specified positions
316
317			//Zero Out the force of z-contrained molecules
318			totalFZ_local = 0;
319
320
321			//calculate the total z-contrained force of fixed z-contrained molecules
322			std::list<ZconstraintMol>::iterator i;
323			Molecule* mol;
324			StuntDouble* integrableObject;
325			Molecule::IntegrableObjectIteraot ii;
326
327			for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
328			mol = i->mol;
329			com = mol->getCom();
330
331			i->fz = 0.0;
332
333			totalFZ_local += harmonicF;
334
335			//adjust force
336			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
337			integrableObject = mol->nextIntegrableObject(ii)) {
338
339			force = integrableObject->getFrc();
340			i->fz += force[whichDirection]
341			}
342			totalFZ_local += i->fz;
343			}
344
345			//calculate total z-constraint force
346			#ifdef IS_MPI
347			MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
348			#else
349			totalFZ = totalFZ_local;
350			#endif
351
352
353			// apply negative to fixed z-constrained molecues;
354			for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
355			mol = i->mol;
356			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
357			integrableObject = mol->nextIntegrableObject(ii)) {
358
359			force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz);
360			integrableObject->addFrc(force);
361			}
362			}
363
364			//modify the forces of moving z-constrained molecules
365			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
366			mol = i->mol;
367			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
368			integrableObject = mol->nextIntegrableObject(ii)) {
369
370			force[whichDirection] = -getZFOfMovingMols(mol, integrableObject, totalFZ);
371			integrableObject->addFrc(force);
372			}
373			}
374
375			//modify the forces of unconstrained molecules
376			std::vector<Molecule*>::iterator j;
377			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
378			mol =*j;
379			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
380			integrableObject = mol->nextIntegrableObject(ii)) {
381
382			force[whichDirection] = -getZFOfMovingMols(mol, integrableObject, totalFZ);
383			integrableObject->addFrc(force);
384			}
385			}
386
387			}
388
389
390			void ZconstraintForceManager::doHarmonic(){
391			double harmonicU;
392			double harmonicF;
393			double diff;
394			double totalFZ_local;
395			double totalFZ;
396
397			Vector3d force(0.0);
398			Vector3d com;
399
400			double totalFZ_local = 0;
401
402			std::list<ZconstraintMol>::iterator i;
403			Molecule* mol;
404			StuntDouble* integrableObject;
405			Molecule::IntegrableObjectIteraot ii;
406
407			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
408			mol = i->mol;
409			com = mol->getCom();
410			diff = com[whichDirection] - resPos;
411			harmonicU = 0.5 * kz * diff * diff;
412			currSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
413			harmonicF = -kz * diff;
414			totalFZ_local += harmonicF;
415
416			//adjust force
417			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
418			integrableObject = mol->nextIntegrableObject(ii)) {
419
420			force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF);
421			integrableObject->addFrc(force);
422			}
423			}
424
425			#ifndef IS_MPI
426			totalFZ = totalFZ_local;
427			#else
428			MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
429			#endif
430
431			//modify the forces of unconstrained molecules
432			std::vector<Molecule*>::iterator j;
433			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
434			mol = *j;
435			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
436			integrableObject = mol->nextIntegrableObject(ii)) {
437
438			force[whichDirection] = getHFOfUnconsMols(mol, harmonicF);
439			integrableObject->addFrc(force);
440			}
441			}
442
443			}
444
445			bool ZconstraintForceManager::checkZConsState(){
446			Vector3d com;
447			double diff;
448			int changed_local = 0;
449
450			std::list<ZconstraintMol>::iterator i;
451			std::list<ZconstraintMol>::iterator j;
452
453			std::list<ZconstraintMol> newMovingZMols;
454			for ( i = fixedZMols_.begin(); i != fixedZMols_.end();) {
455			com = i->mol->getCom();
456			diff = fabs(com[whichDirection] - i->param.zTargetPos);
457			if (diff > zconsTol_) {
458			//this fixed zconstraint molecule is about to moving
459			//moved this molecule to
460			j = i++;
461			newMovingZMols.push_back(fixedZMols_.erase(j));
462			if (usingZconsGap_) {
463			i->endFixingTime = infiniteTime;
464			}
465			changed_local = 1;
466			}else {
467			++i;
468			}
469			}
470
471			std::list<ZconstraintMol> newFixedZMols;
472			for ( i = movingZMols_.begin(); i != movingZMols_.end();) {
473			com = i->mol->getCom();
474			diff = fabs(com[whichDirection] - i->param.zTargetPos);
475			if (diff <= zconsTol_) {
476			//this moving zconstraint molecule is about to fixed
477			//moved this molecule to
478			j = i++;
479			newFixedZMols.push_back(movingZMols_.erase(j));
480			if (usingZconsGap_) {
481			i->endFixingTime = currSnapshot_->getTime() + zconsFixingTime_;
482			}
483			changed_local = 1;
484			}else {
485			++i;
486			}
487			}
488
489			//merge the lists
490			fixedZMols_.merge(newFixedZMols);
491			movingZMols_.merge(newFixedZMols);
492
493			int changed;
494			#ifndef IS_MPI
495			changed = changed_local;
496			#else
497			MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
498			#endif
499
500			return (changed > 0);
501			}
502
503			bool ZconstraintForceManager::haveFixedZMols(){
504			int havingFixed;
505			int havingFixed_local = fixedZMols_.empty() ? 0 : 1;
506
507			#ifndef IS_MPI
508			havingFixed = havingFixed_local;
509			#else
510			MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
511			MPI_COMM_WORLD);
512			#endif
513
514			return havingFixed > 0;
515			}
516
517
518			bool ZconstraintForceManager::haveMovingZMols(){
519			int havingMoving_local;
520			int havingMoving;
521
522			havingMoving_local = movingZMols_.empty()? 0 : 1;
523
524			#ifndef IS_MPI
525			havingMoving = havingMoving_local;
526			#else
527			MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
528			MPI_COMM_WORLD);
529			#endif
530
531			return havingMoving > 0;
532			}
533
534			void ZconstraintForceManager::calcTotalMassMovingZMols(){
535
536			double totMassMovingZMols_local = 0.0;
537			std::list<ZconstraintMol>::iterator i;
538			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
539			totMassMovingZMols_local += i->mol->getMass();
540			}
541
542			#ifdef IS_MPI
543			MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_DOUBLE,
544			MPI_SUM, MPI_COMM_WORLD);
545			#else
546			totMassMovingZMols_ = massOfMovingZAtoms_local;
547			#endif
548
549			}
550
551			double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){
552			return totalForce * sd->getMass() / mol->getMass();
553			}
554
555			double ZconstraintForceManager::getZFOfMovingMols(StuntDouble*sd, double totalForce){
556			return totalForce * sd->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
557			}
558
559			double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){
560			return totalForce * sd->getMass() / mol->getMass();
561			}
562
563			double ZconstraintForceManager::getHFOfUnconsMols(StuntDouble*sd, double totalForce){
564			return totalForce * sd->getMass() / totMassUnconsMols_;
565			}
566
567			void ZconstraintForceManager::updateZPos(){
568			double curTime = currSnapshot_->getTime();
569			std::list<ZconstraintMol>::iterator i;
570			for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
571			i->param.zTargetPos += zconsGap_;
572			}
573			}
574
575			void ZconstraintForceManager::updateCantPos(){
576			std::list<ZconstraintMol>::iterator i;
577			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
578			i->cantPos += i->param.cantVel * dt_;
579			}
580			}
581
582			double ZconstraintForceManager::getZTargetPos(int index){
583			double zTargetPos;
584			#ifndef IS_MPI
585			Molecule* mol = info_->getMoleculeByGlobalIndex(index);
586			assert(!mol);
587			Vector3d com = mol->getCom();
588			zTargetPos = com[whichDirection];
589			#else
590			int whicProc = info_->getMolToProc(index);
591			MPI_Bcast(&zTargetPos, 1, MPI_DOUBLE, whicProc, MPI_COMM_WORLD);
592			#endif
593			return zTargetPos;
594			}
595
596			}