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), infiniteTime(1e31) {
    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:	1926
Committed:	Tue Jan 11 21:03:34 2005 UTC (19 years, 6 months ago) by tim
File size:	18842 byte(s)
Log Message:	change const static double data member to const double
#	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	#include "utils/StringUtils.hpp"
7	namespace oopse {
8	ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info), infiniteTime(1e31) {
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	usingZconsGap_ = true;
49	zconsGap_ = simParam->getZconsGap();
50	}else {
51	usingZconsGap_ = false;
52	zconsGap_ = 0.0;
53	}
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	zconsOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".fz";
70
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	param.zTargetPos = stamp[i]->getZpos();
91	} else {
92	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	fzOut = new ZConsWriter(info_, zconsOutput_.c_str());
124
125	if (!fzOut){
126	sprintf(painCave.errMsg, "Fail to create ZConsWriter\n");
127	painCave.isFatal = 1;
128	simError();
129	}
130
131	}
132
133	ZconstraintForceManager::~ZconstraintForceManager(){
134
135	if (fzOut){
136	delete fzOut;
137	}
138
139	}
140
141	void ZconstraintForceManager::update(){
142	fixedZMols_.clear();
143	movingZMols_.clear();
144	unzconsMols_.clear();
145
146	for (std::map<int, ZconstraintParam>::iterator i = allZMolIndices_.begin(); i != allZMolIndices_.end(); ++i) {
147	#ifdef IS_MPI
148	if (info_->getMolToProc(i->first) == worldRank) {
149	#endif
150	ZconstraintMol zmol;
151	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	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	for (std::list<ZconstraintMol>::iterator i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
172	zmolSet.insert(i->mol->getGlobalIndex());
173	}
174
175	for (std::list<ZconstraintMol>::iterator i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
176	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	if (zmolSet.find(mol->getGlobalIndex()) == zmolSet.end()) {
183	unzconsMols_.push_back(mol);
184	}
185	}
186
187	}
188
189	bool ZconstraintForceManager::isZMol(Molecule* mol){
190	return allZMolIndices_.find(mol->getGlobalIndex()) == allZMolIndices_.end() ? false : true;
191	}
192
193	void ZconstraintForceManager::init(){
194
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	currZconsTime_ = currSnapshot_->getTime();
200	}
201
202	void ZconstraintForceManager::calcForces(bool needPotential, bool needStress){
203	ForceManager::calcForces(needPotential, needStress);
204
205	if (usingZconsGap_){
206	updateZPos();
207	}
208
209	if (checkZConsState()){
210	zeroVelocity();
211	calcTotalMassMovingZMols();
212	}
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	if (currSnapshot_->getTime() >= currZconsTime_){
226	std::list<ZconstraintMol>::iterator i;
227	Vector3d com;
228	for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
229	com = i->mol->getCom();
230	i->zpos = com[whichDirection];
231	}
232
233	fzOut->writeFZ(fixedZMols_);
234	currZconsTime_ += zconsTime_;
235	}
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	Molecule::IntegrableObjectIterator ii;
246
247	//zero out the velocities of center of mass of fixed z-constrained molecules
248	for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
249	mol = i->mol;
250	comVel = mol->getComVel();
251	for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
252	integrableObject = mol->nextIntegrableObject(ii)) {
253	vel = integrableObject->getVel();
254	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	double pzMovingMols_local = 0.0;
262	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	comVel = mol->getComVel();
274	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	double vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_);
285
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	vel[whichDirection] -= vzMovingMols;
294	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	vel[whichDirection] -= vzMovingMols;
306	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	Molecule::IntegrableObjectIterator ii;
331
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	i->fz += force[whichDirection];
340	}
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	force[whichDirection] = -getZFOfMovingMols(mol,totalFZ);
370	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	force[whichDirection] = -getZFOfMovingMols(mol, totalFZ);
382	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	StuntDouble* integrableObject;
396	Molecule::IntegrableObjectIterator ii;
397	Molecule* mol;
398	for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
399	mol = i->mol;
400	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	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	force[whichDirection] = getHFOfUnconsMols(mol, totalFZ);
431	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	j = i++;
454	newMovingZMols.push_back(*j);
455	fixedZMols_.erase(j);
456
457	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	//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	changed_local = 1;
477	}else {
478	++i;
479	}
480	}
481
482	//merge the lists
483	fixedZMols_.insert(fixedZMols_.end(), newFixedZMols.begin(), newFixedZMols.end());
484	movingZMols_.insert(movingZMols_.end(), newMovingZMols.begin(), newMovingZMols.end());
485
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	totMassMovingZMols_ = totMassMovingZMols_local;
540	#endif
541
542	}
543
544	double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){
545	return totalForce * sd->getMass() / mol->getMass();
546	}
547
548	double ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, double totalForce){
549	return totalForce * mol->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
550	}
551
552	double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){
553	return totalForce * sd->getMass() / mol->getMass();
554	}
555
556	double ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, double totalForce){
557	return totalForce * mol->getMass() / totMassUnconsMols_;
558	}
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	assert(mol);
580	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	}