src/constraints/ZconstraintForceManager.cpp

#include <cmath>
#include "constraints/ZconstraintForceManager.hpp"
#include "integrators/Integrator.hpp"
#include "utils/simError.h"
#include "utils/OOPSEConstant.hpp"
#include "utils/StringUtils.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()){
        usingZconsGap_ = true;
        zconsGap_ = simParam->getZconsGap();
    }else {
        usingZconsGap_ = false;
        zconsGap_ = 0.0;
    }

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

    //set zconsUsingSMD
    if (simParam->haveZConsUsingSMD()){
        usingSMD_ = simParam->getZconsUsingSMD();
    }else {
        usingSMD_ =false;
    }
    
    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 targetTemp;
    if (simParam->haveTargetTemp()) {
        targetTemp = simParam->getTargetTemp();
    } else {
        targetTemp = 298.0;
    }
    double zforceConstant = OOPSEConstant::kb * targetTemp / (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 = stamp[i]->getZpos();
        } else {
            param.zTargetPos = getZTargetPos(zmolIndex);
        }

        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

    // creat zconsWriter  
    fzOut = new ZConsWriter(info_, zconsOutput_.c_str());   

    if (!fzOut){
      sprintf(painCave.errMsg, "Fail to create ZConsWriter\n");
      painCave.isFatal = 1;
      simError();
    }

}

ZconstraintForceManager::~ZconstraintForceManager(){

  if (fzOut){
    delete fzOut;
  }

}

void ZconstraintForceManager::update(){
    fixedZMols_.clear();
    movingZMols_.clear();
    unzconsMols_.clear();

    for (std::map<int, ZconstraintParam>::iterator 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 = 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;
    for (std::list<ZconstraintMol>::iterator i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        zmolSet.insert(i->mol->getGlobalIndex());
    }    

    for (std::list<ZconstraintMol>::iterator i = fixedZMols_.begin(); i !=  fixedZMols_.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(mol);
        }
    }

}

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

void ZconstraintForceManager::init(){

  //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();

  currZconsTime_ = currSnapshot_->getTime();
}

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

    if (checkZConsState()){
        zeroVelocity();    
        calcTotalMassMovingZMols();
    }  

    //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() >= currZconsTime_){
        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(fixedZMols_);
        currZconsTime_ += zconsTime_;
    }
}

void ZconstraintForceManager::zeroVelocity(){

    Vector3d comVel;
    Vector3d vel;
    std::list<ZconstraintMol>::iterator i;
    Molecule* mol;
    StuntDouble* integrableObject;
    Molecule::IntegrableObjectIterator 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 = integrableObject->getVel();  
            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 = 0.0;
    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 / (totMassMovingZMols_ + totMassUnconsMols_);

    //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] -= vzMovingMols;
            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] -= vzMovingMols;
            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::IntegrableObjectIterator ii;

    for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
        mol = i->mol;
        i->fz = 0.0;
        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,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, totalFZ);
            integrableObject->addFrc(force);
        }
    }

}


void ZconstraintForceManager::doHarmonic(){
    double totalFZ;
    Vector3d force(0.0);
    Vector3d com;
    double totalFZ_local = 0;
    std::list<ZconstraintMol>::iterator i;
    StuntDouble* integrableObject;
    Molecule::IntegrableObjectIterator ii;
    Molecule* mol;
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
        mol = i->mol;
        com = mol->getCom();   
        double resPos = usingSMD_? i->cantPos : i->param.zTargetPos;
        double diff = com[whichDirection] - resPos; 
        double harmonicU = 0.5 * i->param.kz * diff * diff;
        currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
        double harmonicF = -i->param.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, totalFZ);
            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_) {
            if (usingZconsGap_) {
                i->endFixingTime = infiniteTime;
            }
            j = i++;
            newMovingZMols.push_back(*j);
            fixedZMols_.erase(j);
            
            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_) {
            if (usingZconsGap_) {
                i->endFixingTime = currSnapshot_->getTime() + zconsFixingTime_;
            }
            //this moving zconstraint molecule is about to fixed
            //moved this molecule to
            j = i++;
            newFixedZMols.push_back(*j);
            movingZMols_.erase(j);
            changed_local = 1;            
        }else {
            ++i;
        }
    }     

    //merge the lists
    fixedZMols_.insert(fixedZMols_.end(), newFixedZMols.begin(), newFixedZMols.end());
    movingZMols_.insert(movingZMols_.end(), newMovingZMols.begin(), newMovingZMols.end());

    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_ = totMassMovingZMols_local;
#endif

}

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

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

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

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