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 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 = 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

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

}

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(){


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