OOPSE/libmdtools/ZConstraint.cpp

#include "Integrator.hpp"
#include "simError.h"
#include <cmath>
template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) 
                                    : T(theInfo, the_ff), fz(NULL), curZPos(NULL),
                                 indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0)
{

  //get properties from SimInfo
  GenericData* data;
  ZConsParaData* zConsParaData;
  DoubleData* sampleTime;
  DoubleData* tolerance;
  StringData* policy;
  StringData* filename; 
  double COM[3];

  //by default, the direction of constraint is z
  // 0 --> x
  // 1 --> y
  // 2 --> z
  whichDirection = 2;

  //estimate the force constant of harmonical potential
  double Kb = 1.986E-3 ; //in kcal/K
  
  double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2;
  zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox);

  //creat force substraction policy
  data = info->getProperty(ZCONSFORCEPOLICY_ID);
  if(!data){
    sprintf( painCave.errMsg,
               "ZConstraint Warning: User does not set force substraction policy, "
               "PolicyByMass is used\n");
    painCave.isFatal = 0;
    simError();      

    forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
  }
  else{
    policy = dynamic_cast<StringData*>(data);
    
    if(!policy){
      sprintf( painCave.errMsg,
                 "ZConstraint Error: Convertion from GenericData to StringData failure, "
                 "PolicyByMass is used\n");
      painCave.isFatal = 0;
      simError();      

      forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
    }
    else{
      if(policy->getData() == "BYNUMBER")
        forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
      else if(policy->getData() == "BYMASS")
        forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
      else{
        sprintf( painCave.errMsg,
                  "ZConstraint Warning: unknown force substraction policy, "
                  "average force substraction policy is used\n");
        painCave.isFatal = 0;
        simError();      
      }   
    }
  }
  
  
  //retrieve sample time of z-contraint 
  data = info->getProperty(ZCONSTIME_ID);
  
  if(!data) {
      
    sprintf( painCave.errMsg,
               "ZConstraint error: If you use an ZConstraint\n"
               " , you must set sample time.\n");
    painCave.isFatal = 1;
    simError();      
  }
  else{
  
    sampleTime = dynamic_cast<DoubleData*>(data);
     
    if(!sampleTime){

      sprintf( painCave.errMsg,
                 "ZConstraint error: Can not get property from SimInfo\n");
      painCave.isFatal = 1;
      simError();  
       
    }
    else{
      this->zconsTime = sampleTime->getData();
    }

  }
  
  //retrieve output filename of z force
  data = info->getProperty(ZCONSFILENAME_ID);
  if(!data) {

      
    sprintf( painCave.errMsg,
               "ZConstraint error: If you use an ZConstraint\n"
               " , you must set output filename of z-force.\n");
    painCave.isFatal = 1;
    simError();  

  }
  else{

    filename = dynamic_cast<StringData*>(data);
    
    if(!filename){

      sprintf( painCave.errMsg,
                 "ZConstraint error: Can not get property from SimInfo\n");
      painCave.isFatal = 1;
      simError();  
        
    }
    else{
      this->zconsOutput = filename->getData();
    }
    
  }

  //retrieve tolerance for z-constraint molecuels
  data = info->getProperty(ZCONSTOL_ID);
  
  if(!data) {
      
    sprintf( painCave.errMsg,
               "ZConstraint error: can not get tolerance \n");
    painCave.isFatal = 1;
    simError();      
  }
  else{
  
    tolerance = dynamic_cast<DoubleData*>(data);
     
    if(!tolerance){

      sprintf( painCave.errMsg,
                 "ZConstraint error: Can not get property from SimInfo\n");
      painCave.isFatal = 1;
      simError();  
       
    }
    else{
      this->zconsTol = tolerance->getData();
    }

  }
  
  //retrieve index of z-constraint molecules
  data = info->getProperty(ZCONSPARADATA_ID);
  if(!data) {

    sprintf( painCave.errMsg,
               "ZConstraint error: If you use an ZConstraint\n"
               " , you must set index of z-constraint molecules.\n");
    painCave.isFatal = 1;
    simError();  
  }
  else{
  
    zConsParaData = dynamic_cast<ZConsParaData*>(data);
    
    if(!zConsParaData){

      sprintf( painCave.errMsg,
                 "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
      painCave.isFatal = 1;
      simError();  
    
    }
    else{
       
      parameters = zConsParaData->getData();

      //check the range of zconsIndex
      //and the minimum value of index is the first one (we already sorted the data)
      //the maximum value of index is the last one

      int maxIndex;
      int minIndex;
      int totalNumMol;

      minIndex = (*parameters)[0].zconsIndex;
      if(minIndex < 0){
        sprintf( painCave.errMsg,
               "ZConstraint error: index is out of range\n");
        painCave.isFatal = 1;
        simError(); 
        }

      maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;

#ifndef IS_MPI
      totalNumMol = nMols;
#else
      totalNumMol = mpiSim->getTotNmol();   
#endif      
      
      if(maxIndex > totalNumMol - 1){
        sprintf( painCave.errMsg,
               "ZConstraint error: index is out of range\n");
        painCave.isFatal = 1;
        simError();                  
      }

      //if user does not specify the zpos for the zconstraint molecule
      //its initial z coordinate  will be used as default
      for(int i = 0; i < parameters->size(); i++){

        if(!(*parameters)[i].havingZPos){
#ifndef IS_MPI
          for(int j = 0; j < nMols; j++){
            if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
              molecules[j].getCOM(COM);
        break; 
            }
          }
#else
            //query which processor current zconstraint molecule belongs to
          int *MolToProcMap;
          int whichNode;
    double initZPos;
          MolToProcMap = mpiSim->getMolToProcMap();
          whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
        
          //broadcast the zpos of current z-contraint molecule
          //the node which contain this 
           
          if (worldRank == whichNode ){
            
            for(int j = 0; j < nMols; j++)
              if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
                molecules[j].getCOM(COM);
          break; 
              }
         
          }

          MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);        
#endif
            
          (*parameters)[i].zPos = COM[whichDirection];

          sprintf( painCave.errMsg,
                     "ZConstraint warningr: Does not specify zpos for z-constraint molecule "
                     "initial z coornidate will be used \n");
           painCave.isFatal = 0;
           simError();  
    
  }
      }
      
    }//end if (!zConsParaData)
  }//end  if (!data)
             
//  
#ifdef IS_MPI
  update(); 
#else  
  int searchResult;
      
  for(int i = 0; i < nMols; i++){
    
    searchResult = isZConstraintMol(&molecules[i]); 
    
    if(searchResult > -1){
    
      zconsMols.push_back(&molecules[i]);      
      massOfZConsMols.push_back(molecules[i].getTotalMass());  

      zPos.push_back((*parameters)[searchResult].zPos);
//       cout << "index: "<< (*parameters)[searchResult].zconsIndex 
//              <<"\tzPos = " << (*parameters)[searchResult].zPos << endl; 
      kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
      
      molecules[i].getCOM(COM);
    }
    else
    {
    
      unconsMols.push_back(&molecules[i]);
      massOfUnconsMols.push_back(molecules[i].getTotalMass());

    }
  }

  fz = new double[zconsMols.size()];
  curZPos = new double[zconsMols.size()]; 
  indexOfZConsMols = new int [zconsMols.size()];

  if(!fz || !curZPos || !indexOfZConsMols){
    sprintf( painCave.errMsg,
             "Memory allocation failure in class Zconstraint\n");
    painCave.isFatal = 1;
    simError();
  }

  //determine the states of z-constraint molecules
  for(int i = 0; i < zconsMols.size(); i++){
    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();

    zconsMols[i]->getCOM(COM);
    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
      states.push_back(zcsFixed);
    else
      states.push_back(zcsMoving);
  }
  
#endif 

  //get total masss of unconstraint molecules
  double totalMassOfUncons_local;
  totalMassOfUncons_local = 0;
  
  for(int i = 0; i < unconsMols.size(); i++)
    totalMassOfUncons_local += unconsMols[i]->getTotalMass();
    
#ifndef IS_MPI
  totalMassOfUncons = totalMassOfUncons_local;
#else
  MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1,
                      MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
#endif


  //get total number of unconstrained atoms
  int nUnconsAtoms_local;
  nUnconsAtoms_local = 0;
  for(int i = 0; i < unconsMols.size(); i++)
    nUnconsAtoms_local += unconsMols[i]->getNAtoms();
    
#ifndef IS_MPI
  totNumOfUnconsAtoms = nUnconsAtoms_local;
#else
  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, 
                      MPI_INT,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();
  }

  forcePolicy->update();
}

template<typename T> ZConstraint<T>::~ZConstraint()
{
  if(fz)
    delete[] fz;

  if(curZPos)
    delete[] curZPos;
  
  if(indexOfZConsMols)
    delete[] indexOfZConsMols;
  
  if(fzOut)
    delete fzOut;
  
  if(forcePolicy)
    delete forcePolicy;
}


/**
 * 
 */

#ifdef IS_MPI
template<typename T> void ZConstraint<T>::update()
{
  double COM[3];
  int index;
  
  zconsMols.clear();
  massOfZConsMols.clear();
  zPos.clear();
  kz.clear();
  
  unconsMols.clear();
  massOfUnconsMols.clear();
 

  //creat zconsMol and unconsMol lists
  for(int i = 0; i < nMols; i++){
    
    index = isZConstraintMol(&molecules[i]); 
    
    if(index > -1){
    
      zconsMols.push_back(&molecules[i]);      
      zPos.push_back((*parameters)[index].zPos);
      kz.push_back((*parameters)[index].kRatio * zForceConst);
      
      massOfZConsMols.push_back(molecules[i].getTotalMass());  
       
      molecules[i].getCOM(COM);
    }
    else
    {
    
      unconsMols.push_back(&molecules[i]);
      massOfUnconsMols.push_back(molecules[i].getTotalMass());

    }
  }

  //determine the states of z-constraint molecules
  for(int i = 0; i < zconsMols.size(); i++){
     zconsMols[i]->getCOM(COM);
      if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
      states.push_back(zcsFixed);
     else
      states.push_back(zcsMoving);
  }

    
  //The reason to declare fz and indexOfZconsMols as pointer to array is
  // that we want to make the MPI communication simple
  if(fz)
    delete[] fz;
  
  if(curZPos)
    delete[] curZPos;
    
  if(indexOfZConsMols)
    delete[] indexOfZConsMols;
    
  if (zconsMols.size() > 0){
    fz = new double[zconsMols.size()];
   curZPos = new double[zconsMols.size()];
    indexOfZConsMols =  new int[zconsMols.size()];
    
    if(!fz || !curZPos || !indexOfZConsMols){
      sprintf( painCave.errMsg,
               "Memory allocation failure in class Zconstraint\n");
      painCave.isFatal = 1;
      simError();
    }
        
    for(int i = 0; i < zconsMols.size(); i++){
      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
    }

  }
  else{
    fz = NULL;
   curZPos = NULL;
    indexOfZConsMols = NULL;
  }
  
  //
  forcePolicy->update();
  
}

#endif

/**
 *  Function Name: isZConstraintMol
 *  Parameter
 *    Molecule* mol
 *  Return value:
 *    -1, if the molecule is not z-constraint molecule, 
 *    other non-negative values, its index in indexOfAllZConsMols vector
 */

template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
{
  int index;
  int low;
  int high;
  int mid;

  index = mol->getGlobalIndex();
  
  low = 0;
  high = parameters->size() - 1;
  
  //Binary Search (we have sorted the array)  
  while(low <= high){
    mid = (low + high) /2;
    if ((*parameters)[mid].zconsIndex == index)
      return mid;
    else if ((*parameters)[mid].zconsIndex > index )
       high = mid -1;
    else     
      low = mid + 1; 
  }
  
  return -1;
}

template<typename T> void ZConstraint<T>::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
  zeroOutVel();

  curZconsTime = zconsTime + info->getTime();
  
  T::integrate();

}
 

/**
 *
 *
 *
 *
 */ 
template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
  double zsys;
  double COM[3];
  double force[3];
  double zSysCOMVel;

  T::calcForce(calcPot, calcStress);

  if (checkZConsState()){
    
#ifdef IS_MPI
    if(worldRank == 0){
#endif
//       std::cerr << "\n"
//              << "*******************************************\n"
//              << " about to call zeroOutVel()\n"
//              << "*******************************************\n"
//              << "\n";
#ifdef IS_MPI
    }
#endif
    zeroOutVel();

#ifdef IS_MPI
    if(worldRank == 0){
#endif
//       std::cerr << "\n"
//              << "*******************************************\n"
//              << " finished zeroOutVel()\n"
//              << "*******************************************\n"
//              << "\n";
#ifdef IS_MPI
    }
#endif
    
    forcePolicy->update();
  }  
  
  zsys = calcZSys();
  zSysCOMVel = calcSysCOMVel();
#ifdef IS_MPI
  if(worldRank == 0){
#endif
//     cout << "---------------------------------------------------------------------" <<endl;
//     cout << "current time: " << info->getTime() << endl;
//     cout << "center of mass at z: " << zsys << endl;     
//     cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;

#ifdef IS_MPI
  }
#endif

  //do zconstraint force; 
  if (haveFixedZMols())
    this->doZconstraintForce();
     
  //use harmonical poteintial to move the molecules to the specified positions
  if (haveMovingZMols())
    this->doHarmonic();

  //write out forces and current positions of z-constraint molecules
  if(info->getTime() >= curZconsTime){    
   for(int i = 0; i < zconsMols.size(); i++){
      zconsMols[i]->getCOM(COM);
    curZPos[i] = COM[whichDirection];

    //if the z-constraint molecule is still moving, just record its force
    if(states[i] == zcsMoving){
         fz[i] = 0;
      Atom** movingZAtoms;
      movingZAtoms = zconsMols[i]->getMyAtoms();
      for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
           movingZAtoms[j]->getFrc(force);
           fz[i] += force[whichDirection];
      }
     }
   }
    fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos);
   curZconsTime += zconsTime;
  }

  zSysCOMVel = calcSysCOMVel();  
#ifdef IS_MPI
  if(worldRank == 0){
#endif
 //    cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
#ifdef IS_MPI
  }
#endif

}


/**
 * 
 */
 
template<typename T> double ZConstraint<T>::calcZSys()
{
  //calculate reference z coordinate for z-constraint molecules
  double totalMass_local;
  double totalMass;
  double totalMZ_local;
  double totalMZ;
  double massOfCurMol;
  double COM[3];
  
  totalMass_local = 0;
  totalMZ_local = 0;
  
  for(int i = 0; i < nMols; i++){
    massOfCurMol = molecules[i].getTotalMass(); 
    molecules[i].getCOM(COM); 
     
    totalMass_local += massOfCurMol;
    totalMZ_local += massOfCurMol * COM[whichDirection];

  }
 
  
#ifdef IS_MPI  
  MPI_Allreduce(&totalMass_local, &totalMass, 1, 
                      MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&totalMZ_local, &totalMZ, 1,
                      MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
 #else
  totalMass = totalMass_local;
  totalMZ = totalMZ_local;
 #endif  

  double zsys;
  zsys = totalMZ / totalMass;

  return zsys;
}

/**
 *
 */
template<typename T> void ZConstraint<T>::thermalize( void ){

  T::thermalize();
  zeroOutVel();
}

/**
 *
 */

template<typename T> void ZConstraint<T>::zeroOutVel(){

  Atom** fixedZAtoms;  
  double COMvel[3];
  double vel[3];
  double zSysCOMVel;

  //zero out the velocities of center of mass of fixed z-constrained molecules
  
  for(int i = 0; i < zconsMols.size(); i++){

    if (states[i] == zcsFixed){ 

     zconsMols[i]->getCOMvel(COMvel);      
    //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;

      fixedZAtoms = zconsMols[i]->getMyAtoms(); 
    
      for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
        fixedZAtoms[j]->getVel(vel);
       vel[whichDirection] -= COMvel[whichDirection];
       fixedZAtoms[j]->setVel(vel);
      }

    zconsMols[i]->getCOMvel(COMvel);
    //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
    }
  
  }

    //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; 

  zSysCOMVel = calcSysCOMVel();
#ifdef IS_MPI
  if(worldRank == 0){
#endif
//     cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
#ifdef IS_MPI
  }
#endif
       
  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
  double MVzOfMovingMols_local;
  double MVzOfMovingMols;
  double totalMassOfMovingZMols_local;
  double totalMassOfMovingZMols;
      
  MVzOfMovingMols_local = 0;
  totalMassOfMovingZMols_local = 0;

  for(int i =0; i < unconsMols.size(); i++){
    unconsMols[i]->getCOMvel(COMvel);
    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
  } 

  for(int i = 0; i < zconsMols.size(); i++){
    if (states[i] == zcsMoving){
      zconsMols[i]->getCOMvel(COMvel);
      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];   
      totalMassOfMovingZMols_local += massOfZConsMols[i];       
    }
    
  }

#ifndef IS_MPI
  MVzOfMovingMols = MVzOfMovingMols_local;
  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
#else
  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
#endif

  double vzOfMovingMols;
  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);

  //modify the velocites of unconstrained molecules  
  Atom** unconsAtoms;
  for(int i = 0; i < unconsMols.size(); i++){
  
    unconsAtoms = unconsMols[i]->getMyAtoms();
    for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
      unconsAtoms[j]->getVel(vel);
      vel[whichDirection] -= vzOfMovingMols;
      unconsAtoms[j]->setVel(vel);
    }
  
  }  

  //modify the velocities of moving z-constrained molecuels
  Atom** movingZAtoms;
  for(int i = 0; i < zconsMols.size(); i++){

    if (states[i] ==zcsMoving){
   
      movingZAtoms = zconsMols[i]->getMyAtoms();
     for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
        movingZAtoms[j]->getVel(vel);
        vel[whichDirection] -= vzOfMovingMols;
       movingZAtoms[j]->setVel(vel);
    }
    
   }

 }


  zSysCOMVel = calcSysCOMVel();
#ifdef IS_MPI
  if(worldRank == 0){
#endif
//     cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
#ifdef IS_MPI
  }
#endif

}

/**
 * 
 */

template<typename T> void ZConstraint<T>::doZconstraintForce(){

  Atom** zconsAtoms;
  double totalFZ; 
  double totalFZ_local;
  double COMvel[3];   
  double COM[3];
  double force[3];


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

  for(int i = 0; i < zconsMols.size(); i++){
    
    if (states[i] == zcsFixed){
      
      zconsMols[i]->getCOM(COM);
      zconsAtoms = zconsMols[i]->getMyAtoms();  

      fz[i] = 0;      
      for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
        zconsAtoms[j]->getFrc(force);
        fz[i] += force[whichDirection];      
      } 
      totalFZ_local += fz[i];

      //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i] 
      //      <<"\tcurrent zpos: " << COM[whichDirection] 
      //       << "\tcurrent fz: " <<fz[i] << endl;


    }
    
  } 

  //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;
  force[0]= 0;
  force[1]= 0;
  force[2]= 0;

  for(int i = 0; i < zconsMols.size(); i++){

    if (states[i] == zcsFixed){  
  
      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
      zconsAtoms = zconsMols[i]->getMyAtoms();  
    
      for(int j =0; j < nAtomOfCurZConsMol; j++) {
        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
        force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
        zconsAtoms[j]->addFrc(force);
      }
    
    }
  
  } 

//   cout << "after zero out z-constraint force on fixed z-constraint molecuels "
//        << "total force is " << calcTotalForce() << endl;

  //calculate the number of atoms of moving z-constrained molecules
  int nMovingZAtoms_local;
  int nMovingZAtoms;
  
  nMovingZAtoms_local = 0;
  for(int i = 0; i < zconsMols.size(); i++)
    if(states[i] == zcsMoving)
     nMovingZAtoms_local += zconsMols[i]->getNAtoms();
  
#ifdef IS_MPI
  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, 
                      MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#else
  nMovingZAtoms = nMovingZAtoms_local;
#endif

  force[0]= 0;
  force[1]= 0;
  force[2]= 0;

  //modify the forces of unconstrained molecules
  for(int i = 0; i < unconsMols.size(); i++){
     
     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
     
     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){          
       //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
       force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
       unconsAtoms[j]->addFrc(force);
     }
     
  }      

 //modify the forces of moving z-constrained molecules
  for(int i = 0; i < zconsMols.size(); i++) {
    if (states[i] == zcsMoving){
    
      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();     

      for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
        movingZAtoms[j]->addFrc(force);
      }
    }
  }

  //cout << "after substracting z-constraint force from moving molecuels "
  //      << "total force is " << calcTotalForce()  << endl;

}

/**
  * 
  *
  */

template<typename T> void ZConstraint<T>::doHarmonic(){
  double force[3];
  double harmonicU;
  double harmonicF;
  double COM[3];
  double diff;
  double totalFZ_local;
  double totalFZ;
  
  force[0] = 0;
  force[1] = 0;
  force[2] = 0;

  totalFZ_local = 0;

  for(int i = 0; i < zconsMols.size(); i++) {

    if (states[i] == zcsMoving){
      zconsMols[i]->getCOM(COM);
//       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] 
//         << "\tcurrent zpos: " << COM[whichDirection] << endl;

      diff = COM[whichDirection] -zPos[i];
    
      harmonicU = 0.5 * kz[i] * diff * diff;  
      info->lrPot += harmonicU;

      harmonicF =  - kz[i] * diff;
      totalFZ_local += harmonicF;

       //adjust force
    
      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();     

       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){          
        //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
         movingZAtoms[j]->addFrc(force);
       }
    }

  }

#ifndef IS_MPI
  totalFZ = totalFZ_local;
#else
  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
#endif

  force[0]= 0;
  force[1]= 0;
  force[2]= 0;

  //modify the forces of unconstrained molecules
  for(int i = 0; i < unconsMols.size(); i++){
     
     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
     
     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){          
       //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
       force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
       unconsAtoms[j]->addFrc(force);     
     }
  }   

}

/**
 * 
 */

template<typename T> bool ZConstraint<T>::checkZConsState(){
  double COM[3];
  double diff;
  
  int changed_local;
  int changed;
  
  changed_local = 0;
  
  for(int i =0; i < zconsMols.size(); i++){

    zconsMols[i]->getCOM(COM);
    diff = fabs(COM[whichDirection] - zPos[i]);  
    if (  diff <= zconsTol && states[i] == zcsMoving){
      states[i] = zcsFixed;
     changed_local = 1;
    }
    else if ( diff > zconsTol && states[i] == zcsFixed){
      states[i] = zcsMoving;
     changed_local = 1;   
    }
  
  }

#ifndef IS_MPI
  changed =changed_local; 
#else
  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#endif

  return (changed > 0);

}

template<typename T> bool ZConstraint<T>::haveFixedZMols(){

  int havingFixed_local;
  int havingFixed;

  havingFixed_local = 0;

  for(int i = 0; i < zconsMols.size(); i++)
    if (states[i] == zcsFixed){
      havingFixed_local = 1;
    break;
    }

#ifndef IS_MPI
  havingFixed = havingFixed_local;
#else
  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#endif

  return (havingFixed > 0);
}


/**
 *
 */
template<typename T> bool ZConstraint<T>::haveMovingZMols(){

  int havingMoving_local;
  int havingMoving;

  havingMoving_local = 0;

  for(int i = 0; i < zconsMols.size(); i++)
    if (states[i] == zcsMoving){
      havingMoving_local = 1;
    break;
    }

#ifndef IS_MPI
  havingMoving = havingMoving_local;
#else
  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
#endif

  return (havingMoving > 0);
  
}

/**
 *
 */

template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
{
  double MVzOfMovingMols_local;
  double MVzOfMovingMols;
  double totalMassOfMovingZMols_local;
  double totalMassOfMovingZMols;
  double COMvel[3];
      
  MVzOfMovingMols_local = 0;
  totalMassOfMovingZMols_local = 0;

  for(int i =0; i < unconsMols.size(); i++){
    unconsMols[i]->getCOMvel(COMvel);
    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
  } 

  for(int i = 0; i < zconsMols.size(); i++){

    if (states[i] == zcsMoving){
      zconsMols[i]->getCOMvel(COMvel);
      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];   
      totalMassOfMovingZMols_local += massOfZConsMols[i];       
    }
    
  }

#ifndef IS_MPI
  MVzOfMovingMols = MVzOfMovingMols_local;
  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
#else
  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
#endif

  double vzOfMovingMols;
  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);

  return vzOfMovingMols;
}

/**
 *
 */

template<typename T> double ZConstraint<T>::calcSysCOMVel()
{
  double COMvel[3];
  double tempMVz_local;
  double tempMVz;
  double massOfZCons_local;
  double massOfZCons;


 tempMVz_local = 0;

  for(int i =0 ; i < nMols; i++){
    molecules[i].getCOMvel(COMvel);
   tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
  }

  massOfZCons_local = 0;
  
  for(int i = 0; i < massOfZConsMols.size(); i++){
    massOfZCons_local += massOfZConsMols[i];
  }
#ifndef IS_MPI
  massOfZCons = massOfZCons_local;
  tempMVz = tempMVz_local;
#else
  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#endif

  return tempMVz /(totalMassOfUncons + massOfZCons);
}

/**
 *
 */

template<typename T> double ZConstraint<T>::calcTotalForce(){

  double force[3];  
  double totalForce_local;
  double totalForce;

  totalForce_local = 0;

  for(int i = 0; i < nAtoms; i++){
    atoms[i]->getFrc(force);
    totalForce_local += force[whichDirection];
  }

#ifndef IS_MPI
  totalForce = totalForce_local;
#else
  MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#endif

  return totalForce;

}

/**
 *
 */

template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
  //calculate the number of atoms of moving z-constrained molecules
  int nMovingZAtoms_local;
  int nMovingZAtoms;
  
  nMovingZAtoms_local = 0;
  for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
    if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
     nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
  
#ifdef IS_MPI
  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#else
  nMovingZAtoms = nMovingZAtoms_local;
#endif
  totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;

#ifdef IS_MPI
  if(worldRank == 0){
#endif
 //    std::cerr << "\n"
//            << "*******************************************\n"
//            << " fiished Policy by numbr()\n"
//            << "*******************************************\n"
//            << "\n";
#ifdef IS_MPI
  }
#endif
}

template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
  return totalForce / mol->getNAtoms();
}

template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, double totalForce){
  return totalForce / totNumOfMovingAtoms;
}

template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
    return totalForce / mol->getNAtoms();
}

template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, double totalForce){
  return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
}

/**
 *
 */

template<typename T> void ZConstraint<T>::PolicyByMass::update(){
  //calculate the number of atoms of moving z-constrained molecules
  double massOfMovingZAtoms_local;
  double massOfMovingZAtoms;
  
  massOfMovingZAtoms_local = 0;
  for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
    if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
     massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
  
#ifdef IS_MPI
  MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#else
  massOfMovingZAtoms = massOfMovingZAtoms_local;
#endif
  totMassOfMovingAtoms = massOfMovingZAtoms_local + zconsIntegrator->totalMassOfUncons;
}

template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
  return totalForce * atom->getMass() / mol->getTotalMass();
}

template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols( Atom* atom, double totalForce){
    return totalForce * atom->getMass() / totMassOfMovingAtoms;
}

template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
  return totalForce * atom->getMass() / mol->getTotalMass();
}

template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, double totalForce){
    return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
}

Revision:	726
Committed:	Tue Aug 26 20:37:30 2003 UTC (20 years, 10 months ago) by tim
File size:	32141 byte(s)
Log Message:	set default force substraction policy to PolicyByMass