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, 7 months ago) by tim
File size:	18107 byte(s)
Log Message:	ZConstraintForceManager in progress
#	Content
1	#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	}