src/constraints/ZconstraintForceManager.cpp

 /*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
#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(){
  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:	1930
Committed:	Wed Jan 12 22:41:40 2005 UTC (19 years, 5 months ago) by gezelter
File size:	20309 byte(s)
Log Message:	merging new_design branch into OOPSE-2.0
#	User	Rev	Content
1	gezelter	1930	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41
42			#include <cmath>
43			#include "constraints/ZconstraintForceManager.hpp"
44			#include "integrators/Integrator.hpp"
45			#include "utils/simError.h"
46			#include "utils/OOPSEConstant.hpp"
47			#include "utils/StringUtils.hpp"
48			namespace oopse {
49			ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info), infiniteTime(1e31) {
50			currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
51			Globals* simParam = info_->getSimParams();
52
53			if (simParam->haveDt()){
54			dt_ = simParam->getDt();
55			} else {
56			sprintf(painCave.errMsg,
57			"Integrator Error: dt is not set\n");
58			painCave.isFatal = 1;
59			simError();
60			}
61
62			if (simParam->haveZconstraintTime()){
63			zconsTime_ = simParam->getZconsTime();
64			}
65			else{
66			sprintf(painCave.errMsg,
67			"ZConstraint error: If you use a ZConstraint,\n"
68			"\tyou must set zconsTime.\n");
69			painCave.isFatal = 1;
70			simError();
71			}
72
73			if (simParam->haveZconsTol()){
74			zconsTol_ = simParam->getZconsTol();
75			}
76			else{
77			zconsTol_ = 0.01;
78			sprintf(painCave.errMsg,
79			"ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
80			"\tOOPSE will use a default value of %f.\n"
81			"\tTo set the tolerance, use the zconsTol variable.\n",
82			zconsTol_);
83			painCave.isFatal = 0;
84			simError();
85			}
86
87			//set zcons gap
88			if (simParam->haveZConsGap()){
89			usingZconsGap_ = true;
90			zconsGap_ = simParam->getZconsGap();
91			}else {
92			usingZconsGap_ = false;
93			zconsGap_ = 0.0;
94			}
95
96			//set zcons fixtime
97			if (simParam->haveZConsFixTime()){
98			zconsFixingTime_ = simParam->getZconsFixtime();
99			} else {
100			zconsFixingTime_ = infiniteTime;
101			}
102
103			//set zconsUsingSMD
104			if (simParam->haveZConsUsingSMD()){
105			usingSMD_ = simParam->getZconsUsingSMD();
106			}else {
107			usingSMD_ =false;
108			}
109
110			zconsOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".fz";
111
112			//estimate the force constant of harmonical potential
113			Mat3x3d hmat = currSnapshot_->getHmat();
114			double halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2;
115			double targetTemp;
116			if (simParam->haveTargetTemp()) {
117			targetTemp = simParam->getTargetTemp();
118			} else {
119			targetTemp = 298.0;
120			}
121			double zforceConstant = OOPSEConstant::kb * targetTemp / (halfOfLargestBox * halfOfLargestBox);
122
123			int nZconstraints = simParam->getNzConstraints();
124			ZconStamp** stamp = simParam->getZconStamp();
125			//
126			for (int i = 0; i < nZconstraints; i++){
127
128			ZconstraintParam param;
129			int zmolIndex = stamp[i]->getMolIndex();
130			if (stamp[i]->haveZpos()) {
131			param.zTargetPos = stamp[i]->getZpos();
132			} else {
133			param.zTargetPos = getZTargetPos(zmolIndex);
134			}
135
136			param.kz = zforceConstant * stamp[i]->getKratio();
137
138			if (stamp[i]->haveCantVel()) {
139			param.cantVel = stamp[i]->getCantVel();
140			} else {
141			param.cantVel = 0.0;
142			}
143
144			allZMolIndices_.insert(std::make_pair(zmolIndex, param));
145			}
146
147			//create fixedMols_, movingMols_ and unconsMols lists
148			update();
149
150			//calculate masss of unconstraint molecules in the whole system (never change during the simulation)
151			double totMassUnconsMols_local = 0.0;
152			std::vector<Molecule*>::iterator j;
153			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
154			totMassUnconsMols_local += (*j)->getMass();
155			}
156			#ifndef IS_MPI
157			totMassUnconsMols_ = totMassUnconsMols_local;
158			#else
159			MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_DOUBLE,
160			MPI_SUM, MPI_COMM_WORLD);
161			#endif
162
163			// creat zconsWriter
164			fzOut = new ZConsWriter(info_, zconsOutput_.c_str());
165
166			if (!fzOut){
167			sprintf(painCave.errMsg, "Fail to create ZConsWriter\n");
168			painCave.isFatal = 1;
169			simError();
170			}
171
172			}
173
174			ZconstraintForceManager::~ZconstraintForceManager(){
175
176			if (fzOut){
177			delete fzOut;
178			}
179
180			}
181
182			void ZconstraintForceManager::update(){
183			fixedZMols_.clear();
184			movingZMols_.clear();
185			unzconsMols_.clear();
186
187			for (std::map<int, ZconstraintParam>::iterator i = allZMolIndices_.begin(); i != allZMolIndices_.end(); ++i) {
188			#ifdef IS_MPI
189			if (info_->getMolToProc(i->first) == worldRank) {
190			#endif
191			ZconstraintMol zmol;
192			zmol.mol = info_->getMoleculeByGlobalIndex(i->first);
193			assert(zmol.mol);
194			zmol.param = i->second;
195			zmol.cantPos = zmol.param.zTargetPos; /*@todo fixed me when zmol migrate, it is incorrect/
196			Vector3d com = zmol.mol->getCom();
197			double diff = fabs(zmol.param.zTargetPos - com[whichDirection]);
198			if (diff < zconsTol_) {
199			fixedZMols_.push_back(zmol);
200			} else {
201			movingZMols_.push_back(zmol);
202			}
203
204			#ifdef IS_MPI
205			}
206			#endif
207			}
208
209			calcTotalMassMovingZMols();
210
211			std::set<int> zmolSet;
212			for (std::list<ZconstraintMol>::iterator i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
213			zmolSet.insert(i->mol->getGlobalIndex());
214			}
215
216			for (std::list<ZconstraintMol>::iterator i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
217			zmolSet.insert(i->mol->getGlobalIndex());
218			}
219
220			SimInfo::MoleculeIterator mi;
221			Molecule* mol;
222			for(mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
223			if (zmolSet.find(mol->getGlobalIndex()) == zmolSet.end()) {
224			unzconsMols_.push_back(mol);
225			}
226			}
227
228			}
229
230			bool ZconstraintForceManager::isZMol(Molecule* mol){
231			return allZMolIndices_.find(mol->getGlobalIndex()) == allZMolIndices_.end() ? false : true;
232			}
233
234			void ZconstraintForceManager::init(){
235
236			//zero out the velocities of center of mass of unconstrained molecules
237			//and the velocities of center of mass of every single z-constrained molecueles
238			zeroVelocity();
239
240			currZconsTime_ = currSnapshot_->getTime();
241			}
242
243			void ZconstraintForceManager::calcForces(bool needPotential, bool needStress){
244			ForceManager::calcForces(needPotential, needStress);
245
246			if (usingZconsGap_){
247			updateZPos();
248			}
249
250			if (checkZConsState()){
251			zeroVelocity();
252			calcTotalMassMovingZMols();
253			}
254
255			//do zconstraint force;
256			if (haveFixedZMols()){
257			doZconstraintForce();
258			}
259
260			//use external force to move the molecules to the specified positions
261			if (haveMovingZMols()){
262			doHarmonic();
263			}
264
265			//write out forces and current positions of z-constraint molecules
266			if (currSnapshot_->getTime() >= currZconsTime_){
267			std::list<ZconstraintMol>::iterator i;
268			Vector3d com;
269			for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
270			com = i->mol->getCom();
271			i->zpos = com[whichDirection];
272			}
273
274			fzOut->writeFZ(fixedZMols_);
275			currZconsTime_ += zconsTime_;
276			}
277			}
278
279			void ZconstraintForceManager::zeroVelocity(){
280
281			Vector3d comVel;
282			Vector3d vel;
283			std::list<ZconstraintMol>::iterator i;
284			Molecule* mol;
285			StuntDouble* integrableObject;
286			Molecule::IntegrableObjectIterator ii;
287
288			//zero out the velocities of center of mass of fixed z-constrained molecules
289			for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
290			mol = i->mol;
291			comVel = mol->getComVel();
292			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
293			integrableObject = mol->nextIntegrableObject(ii)) {
294			vel = integrableObject->getVel();
295			vel[whichDirection] -= comVel[whichDirection];
296			integrableObject->setVel(vel);
297			}
298			}
299
300			// calculate the vz of center of mass of moving molecules(include unconstrained molecules
301			// and moving z-constrained molecules)
302			double pzMovingMols_local = 0.0;
303			double pzMovingMols;
304
305			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
306			mol = i->mol;
307			comVel = mol->getComVel();
308			pzMovingMols_local += mol->getMass() * comVel[whichDirection];
309			}
310
311			std::vector<Molecule*>::iterator j;
312			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
313			mol =*j;
314			comVel = mol->getComVel();
315			pzMovingMols_local += mol->getMass() * comVel[whichDirection];
316			}
317
318			#ifndef IS_MPI
319			pzMovingMols = pzMovingMols_local;
320			#else
321			MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_DOUBLE,
322			MPI_SUM, MPI_COMM_WORLD);
323			#endif
324
325			double vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_);
326
327			//modify the velocities of moving z-constrained molecuels
328			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
329			mol = i->mol;
330			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
331			integrableObject = mol->nextIntegrableObject(ii)) {
332
333			vel = integrableObject->getVel();
334			vel[whichDirection] -= vzMovingMols;
335			integrableObject->setVel(vel);
336			}
337			}
338
339			//modify the velocites of unconstrained molecules
340			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
341			mol =*j;
342			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
343			integrableObject = mol->nextIntegrableObject(ii)) {
344
345			vel = integrableObject->getVel();
346			vel[whichDirection] -= vzMovingMols;
347			integrableObject->setVel(vel);
348			}
349			}
350
351			}
352
353
354			void ZconstraintForceManager::doZconstraintForce(){
355			double totalFZ;
356			double totalFZ_local;
357			Vector3d com;
358			Vector3d force(0.0);
359
360			//constrain the molecules which do not reach the specified positions
361
362			//Zero Out the force of z-contrained molecules
363			totalFZ_local = 0;
364
365
366			//calculate the total z-contrained force of fixed z-contrained molecules
367			std::list<ZconstraintMol>::iterator i;
368			Molecule* mol;
369			StuntDouble* integrableObject;
370			Molecule::IntegrableObjectIterator ii;
371
372			for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
373			mol = i->mol;
374			i->fz = 0.0;
375			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
376			integrableObject = mol->nextIntegrableObject(ii)) {
377
378			force = integrableObject->getFrc();
379			i->fz += force[whichDirection];
380			}
381			totalFZ_local += i->fz;
382			}
383
384			//calculate total z-constraint force
385			#ifdef IS_MPI
386			MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
387			#else
388			totalFZ = totalFZ_local;
389			#endif
390
391
392			// apply negative to fixed z-constrained molecues;
393			for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
394			mol = i->mol;
395			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
396			integrableObject = mol->nextIntegrableObject(ii)) {
397
398			force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz);
399			integrableObject->addFrc(force);
400			}
401			}
402
403			//modify the forces of moving z-constrained molecules
404			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
405			mol = i->mol;
406			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
407			integrableObject = mol->nextIntegrableObject(ii)) {
408
409			force[whichDirection] = -getZFOfMovingMols(mol,totalFZ);
410			integrableObject->addFrc(force);
411			}
412			}
413
414			//modify the forces of unconstrained molecules
415			std::vector<Molecule*>::iterator j;
416			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
417			mol =*j;
418			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
419			integrableObject = mol->nextIntegrableObject(ii)) {
420
421			force[whichDirection] = -getZFOfMovingMols(mol, totalFZ);
422			integrableObject->addFrc(force);
423			}
424			}
425
426			}
427
428
429			void ZconstraintForceManager::doHarmonic(){
430			double totalFZ;
431			Vector3d force(0.0);
432			Vector3d com;
433			double totalFZ_local = 0;
434			std::list<ZconstraintMol>::iterator i;
435			StuntDouble* integrableObject;
436			Molecule::IntegrableObjectIterator ii;
437			Molecule* mol;
438			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
439			mol = i->mol;
440			com = mol->getCom();
441			double resPos = usingSMD_? i->cantPos : i->param.zTargetPos;
442			double diff = com[whichDirection] - resPos;
443			double harmonicU = 0.5 * i->param.kz * diff * diff;
444			currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
445			double harmonicF = -i->param.kz * diff;
446			totalFZ_local += harmonicF;
447
448			//adjust force
449			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
450			integrableObject = mol->nextIntegrableObject(ii)) {
451
452			force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF);
453			integrableObject->addFrc(force);
454			}
455			}
456
457			#ifndef IS_MPI
458			totalFZ = totalFZ_local;
459			#else
460			MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
461			#endif
462
463			//modify the forces of unconstrained molecules
464			std::vector<Molecule*>::iterator j;
465			for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
466			mol = *j;
467			for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
468			integrableObject = mol->nextIntegrableObject(ii)) {
469
470			force[whichDirection] = getHFOfUnconsMols(mol, totalFZ);
471			integrableObject->addFrc(force);
472			}
473			}
474
475			}
476
477			bool ZconstraintForceManager::checkZConsState(){
478			Vector3d com;
479			double diff;
480			int changed_local = 0;
481
482			std::list<ZconstraintMol>::iterator i;
483			std::list<ZconstraintMol>::iterator j;
484
485			std::list<ZconstraintMol> newMovingZMols;
486			for ( i = fixedZMols_.begin(); i != fixedZMols_.end();) {
487			com = i->mol->getCom();
488			diff = fabs(com[whichDirection] - i->param.zTargetPos);
489			if (diff > zconsTol_) {
490			if (usingZconsGap_) {
491			i->endFixingTime = infiniteTime;
492			}
493			j = i++;
494			newMovingZMols.push_back(*j);
495			fixedZMols_.erase(j);
496
497			changed_local = 1;
498			}else {
499			++i;
500			}
501			}
502
503			std::list<ZconstraintMol> newFixedZMols;
504			for ( i = movingZMols_.begin(); i != movingZMols_.end();) {
505			com = i->mol->getCom();
506			diff = fabs(com[whichDirection] - i->param.zTargetPos);
507			if (diff <= zconsTol_) {
508			if (usingZconsGap_) {
509			i->endFixingTime = currSnapshot_->getTime() + zconsFixingTime_;
510			}
511			//this moving zconstraint molecule is about to fixed
512			//moved this molecule to
513			j = i++;
514			newFixedZMols.push_back(*j);
515			movingZMols_.erase(j);
516			changed_local = 1;
517			}else {
518			++i;
519			}
520			}
521
522			//merge the lists
523			fixedZMols_.insert(fixedZMols_.end(), newFixedZMols.begin(), newFixedZMols.end());
524			movingZMols_.insert(movingZMols_.end(), newMovingZMols.begin(), newMovingZMols.end());
525
526			int changed;
527			#ifndef IS_MPI
528			changed = changed_local;
529			#else
530			MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
531			#endif
532
533			return (changed > 0);
534			}
535
536			bool ZconstraintForceManager::haveFixedZMols(){
537			int havingFixed;
538			int havingFixed_local = fixedZMols_.empty() ? 0 : 1;
539
540			#ifndef IS_MPI
541			havingFixed = havingFixed_local;
542			#else
543			MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
544			MPI_COMM_WORLD);
545			#endif
546
547			return havingFixed > 0;
548			}
549
550
551			bool ZconstraintForceManager::haveMovingZMols(){
552			int havingMoving_local;
553			int havingMoving;
554
555			havingMoving_local = movingZMols_.empty()? 0 : 1;
556
557			#ifndef IS_MPI
558			havingMoving = havingMoving_local;
559			#else
560			MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
561			MPI_COMM_WORLD);
562			#endif
563
564			return havingMoving > 0;
565			}
566
567			void ZconstraintForceManager::calcTotalMassMovingZMols(){
568
569			double totMassMovingZMols_local = 0.0;
570			std::list<ZconstraintMol>::iterator i;
571			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
572			totMassMovingZMols_local += i->mol->getMass();
573			}
574
575			#ifdef IS_MPI
576			MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_DOUBLE,
577			MPI_SUM, MPI_COMM_WORLD);
578			#else
579			totMassMovingZMols_ = totMassMovingZMols_local;
580			#endif
581
582			}
583
584			double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){
585			return totalForce * sd->getMass() / mol->getMass();
586			}
587
588			double ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, double totalForce){
589			return totalForce * mol->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
590			}
591
592			double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){
593			return totalForce * sd->getMass() / mol->getMass();
594			}
595
596			double ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, double totalForce){
597			return totalForce * mol->getMass() / totMassUnconsMols_;
598			}
599
600			void ZconstraintForceManager::updateZPos(){
601			double curTime = currSnapshot_->getTime();
602			std::list<ZconstraintMol>::iterator i;
603			for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
604			i->param.zTargetPos += zconsGap_;
605			}
606			}
607
608			void ZconstraintForceManager::updateCantPos(){
609			std::list<ZconstraintMol>::iterator i;
610			for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
611			i->cantPos += i->param.cantVel * dt_;
612			}
613			}
614
615			double ZconstraintForceManager::getZTargetPos(int index){
616			double zTargetPos;
617			#ifndef IS_MPI
618			Molecule* mol = info_->getMoleculeByGlobalIndex(index);
619			assert(mol);
620			Vector3d com = mol->getCom();
621			zTargetPos = com[whichDirection];
622			#else
623			int whicProc = info_->getMolToProc(index);
624			MPI_Bcast(&zTargetPos, 1, MPI_DOUBLE, whicProc, MPI_COMM_WORLD);
625			#endif
626			return zTargetPos;
627			}
628
629			}