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, 5 months ago) by tim
File size:	18858 byte(s)
Log Message:	more work in zconstraint
#	Content
1	#include <cmath>
2	#include "constraints/ZconstraintForceManager.hpp"
3	#include "integrators/Integrator.hpp"
4	#include "utils/simError.h"
5	#include "utils/OOPSEConstant.hpp"
6
7	namespace oopse {
8	ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info) {
9	currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
10	Globals* simParam = info_->getSimParams();
11
12	if (simParam->haveDt()){
13	dt_ = simParam->getDt();
14	} else {
15	sprintf(painCave.errMsg,
16	"Integrator Error: dt is not set\n");
17	painCave.isFatal = 1;
18	simError();
19	}
20
21	if (simParam->haveZconstraintTime()){
22	zconsTime_ = simParam->getZconsTime();
23	}
24	else{
25	sprintf(painCave.errMsg,
26	"ZConstraint error: If you use a ZConstraint,\n"
27	"\tyou must set zconsTime.\n");
28	painCave.isFatal = 1;
29	simError();
30	}
31
32	if (simParam->haveZconsTol()){
33	zconsTol_ = simParam->getZconsTol();
34	}
35	else{
36	zconsTol_ = 0.01;
37	sprintf(painCave.errMsg,
38	"ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
39	"\tOOPSE will use a default value of %f.\n"
40	"\tTo set the tolerance, use the zconsTol variable.\n",
41	zconsTol_);
42	painCave.isFatal = 0;
43	simError();
44	}
45
46	//set zcons gap
47	if (simParam->haveZConsGap()){
48	zconsGap_ = simParam->getZconsGap();
49	}else {
50	zconsGap_ = false;
51	}
52
53	//set zcons fixtime
54	if (simParam->haveZConsFixTime()){
55	zconsFixingTime_ = simParam->getZconsFixtime();
56	} else {
57	zconsFixingTime_ = infiniteTime;
58	}
59
60	//set zconsUsingSMD
61	if (simParam->haveZConsUsingSMD()){
62	usingSMD_ = simParam->getZconsUsingSMD();
63	}else {
64	usingSMD_ =false;
65	}
66
67	std::string zconsOutput = getPrefix(info_->getFinalConfigFileName()) + ".fz";
68
69	//estimate the force constant of harmonical potential
70	Mat3x3d hmat = currSnapshot_->getHmat();
71	double halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2;
72	double 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	}