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), 
                                      indexOfZConsMols(NULL)
{

  //get properties from SimInfo
  GenericData* data;
  ZConsParaData* zConsParaData;
  DoubleData* sampleTime;
  DoubleData* tolerance;
  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);

  //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()].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[i].getGlobalIndex() == (*parameters)[i].zconsIndex){
                 molecules[i].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 i = 0; i < nMols; i++)
               if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){
                 molecules[i].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);
           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()];
  indexOfZConsMols = new int [zconsMols.size()];

  if(!fz || !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();
  
#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_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
#endif

  checkZConsState();

  //  
  fzOut = new ZConsWriter(zconsOutput.c_str());   
  
  if(!fzOut){
    sprintf( painCave.errMsg,
             "Memory allocation failure in class Zconstraint\n");
    painCave.isFatal = 1;
    simError();
  }
  
}

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

#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());

    }
  }
    
  //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(indexOfZConsMols)
    delete[] indexOfZConsMols;
    
  if (zconsMols.size() > 0){
    fz = new double[zconsMols.size()];
    indexOfZConsMols =  new int[zconsMols.size()];
    
    if(!fz || !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;
    indexOfZConsMols = NULL;
  }
  
}

#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;
}

/** 
 * Description:
 *  Reset the z coordinates
 */
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();
  
  T::integrate();

}
 

/**
 *
 *
 *
 *
 */

 
template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){

  T::calcForce(calcPot, calcStress);

  if (checkZConsState())
    zeroOutVel();

  //do zconstraint force; 
  if (haveFixedZMols())
    this->doZconstraintForce();
  
  //use harmonical poteintial to move the molecules to the specified positions
  if (haveMovingZMols())
    //this->doHarmonic();
   
  fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
      
}
 
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 massOfUncons_local;
  double massOfCurMol;
  double COM[3];
  
  totalMass_local = 0;
  totalMass = 0;
  totalMZ_local = 0;
  totalMZ = 0;
  massOfUncons_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];
    
    if(isZConstraintMol(&molecules[i]) == -1){
    
      massOfUncons_local += massOfCurMol;
    }  
    
  }
  
  
#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);  
  MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#else
  totalMass = totalMass_local;
  totalMZ = totalMZ_local;
  totalMassOfUncons = massOfUncons_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];
  
  //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);      
      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);
      }
          
    }
        
  }
  
  // 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[i]->getNAtoms(); 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[i]->getNAtoms(); 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);
        }
          
    }

  }

}

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

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

  int nMovingZMols_local;
  int nMovingZMols;

  //constrain the molecules which do not reach the specified positions  

   zsys = calcZSys();
   cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;   
    
  //Zero Out the force of z-contrained molecules    
  totalFZ_local = 0;

  //calculate the total z-contrained force of fixed z-contrained molecules
  cout << "Fixed Molecules" << endl;
  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 << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;

    }
          
  } 

  //calculate the number of atoms of moving z-constrained molecules
  nMovingZMols_local = 0;
  for(int i = 0; zconsMols.size(); i++){
    if(states[i] == zcsMoving)
        nMovingZMols_local += massOfZConsMols[i];
  }
#ifdef IS_MPI
  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#else
  totalFZ = totalFZ_local;
  nMovingZMols = nMovingZMols_local;
#endif

  force[0]= 0;
  force[1]= 0;
  force[2]= 0;
  force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);

  //modify the velocites 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++)          
       unconsAtoms[j]->addFrc(force);
     
  }      

 //modify the velocities 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++)          
       movingZAtoms[j]->addFrc(force);
     }
  }

  // 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;
        zconsAtoms[j]->addFrc(force);
      }
                
    }
        
  } 

}

template<typename T> bool ZConstraint<T>::checkZConsState(){
  double COM[3];
  double diff;
  
  bool changed;
  
  changed = false;
  
  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 = true;
    }
    else if ( diff > zconsTol && states[i] == zcsFixed){
      states[i] = zcsMoving;
        changed = true;   
    }
  
  }

  return changed;
}

template<typename T> bool ZConstraint<T>::haveFixedZMols(){
  for(int i = 0; i < zconsMols.size(); i++)
    if (states[i] == zcsFixed)
      return true;

  return false;
}


/**
 *
 */
template<typename T> bool ZConstraint<T>::haveMovingZMols(){
  for(int i = 0; i < zconsMols.size(); i++)
    if (states[i] == zcsMoving)
      return true;

  return false;
  
}

/**
  * 
  *
  */

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

  cout << "Moving Molecules" << endl;   
  for(int i = 0; i < zconsMols.size(); i++) {

    if (states[i] == zcsMoving){
      zconsMols[i]->getCOM(COM):
      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
                
                diff = COM[whichDirection] -zPos[i];
                
      harmonicU = 0.5 * kz[i] * diff * diff;  
                info->ltPot += harmonicU;

                force[whichDirection] = - kz[i] * diff / zconsMols[i]->getNAtoms();
                
      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();     

       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
         movingZAtoms[j]->addFrc(force);
    }

  }

}
Revision:	682
Committed:	Tue Aug 12 17:51:33 2003 UTC (20 years, 11 months ago) by tim
File size:	19068 byte(s)
Log Message:	added harmonical potential to z-constraint method
#	User	Rev	Content
1	tim	658	#include "Integrator.hpp"
2			#include "simError.h"
3	tim	676	#include <cmath>
4	tim	658	template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff)
5	tim	676	: T(theInfo, the_ff), fz(NULL),
6			indexOfZConsMols(NULL)
7	tim	658	{
8
9			//get properties from SimInfo
10			GenericData* data;
11	tim	682	ZConsParaData* zConsParaData;
12	tim	658	DoubleData* sampleTime;
13	tim	682	DoubleData* tolerance;
14	tim	658	StringData* filename;
15	tim	682	double COM[3];
16	tim	658
17	tim	682	//by default, the direction of constraint is z
18			// 0 --> x
19			// 1 --> y
20			// 2 --> z
21			whichDirection = 2;
22	tim	658
23	tim	682	//estimate the force constant of harmonical potential
24			double Kb = 1.986E-3 ; //in kcal/K
25
26			double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2;
27			zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox);
28	tim	676
29			//retrieve sample time of z-contraint
30	tim	682	data = info->getProperty(ZCONSTIME_ID);
31	tim	658
32			if(!data) {
33
34			sprintf( painCave.errMsg,
35			"ZConstraint error: If you use an ZConstraint\n"
36			" , you must set sample time.\n");
37			painCave.isFatal = 1;
38			simError();
39			}
40			else{
41
42			sampleTime = dynamic_cast<DoubleData*>(data);
43
44			if(!sampleTime){
45
46			sprintf( painCave.errMsg,
47			"ZConstraint error: Can not get property from SimInfo\n");
48			painCave.isFatal = 1;
49			simError();
50
51			}
52			else{
53			this->zconsTime = sampleTime->getData();
54			}
55
56			}
57
58	tim	676	//retrieve output filename of z force
59	tim	682	data = info->getProperty(ZCONSFILENAME_ID);
60	tim	658	if(!data) {
61
62
63			sprintf( painCave.errMsg,
64			"ZConstraint error: If you use an ZConstraint\n"
65			" , you must set output filename of z-force.\n");
66			painCave.isFatal = 1;
67			simError();
68
69			}
70			else{
71
72			filename = dynamic_cast<StringData*>(data);
73
74			if(!filename){
75
76			sprintf( painCave.errMsg,
77			"ZConstraint error: Can not get property from SimInfo\n");
78			painCave.isFatal = 1;
79			simError();
80
81			}
82			else{
83			this->zconsOutput = filename->getData();
84			}
85
86
87			}
88	tim	682
89			//retrieve tolerance for z-constraint molecuels
90			data = info->getProperty(ZCONSTOL_ID);
91	tim	658
92	tim	682	if(!data) {
93
94			sprintf( painCave.errMsg,
95			"ZConstraint error: can not get tolerance \n");
96			painCave.isFatal = 1;
97			simError();
98			}
99			else{
100
101			tolerance = dynamic_cast<DoubleData*>(data);
102
103			if(!tolerance){
104
105			sprintf( painCave.errMsg,
106			"ZConstraint error: Can not get property from SimInfo\n");
107			painCave.isFatal = 1;
108			simError();
109
110			}
111			else{
112			this->zconsTol = tolerance->getData();
113			}
114
115			}
116
117			//retrieve index of z-constraint molecules
118			data = info->getProperty(ZCONSPARADATA_ID);
119			if(!data) {
120
121			sprintf( painCave.errMsg,
122			"ZConstraint error: If you use an ZConstraint\n"
123			" , you must set index of z-constraint molecules.\n");
124			painCave.isFatal = 1;
125			simError();
126			}
127			else{
128
129			zConsParaData = dynamic_cast<ZConsParaData*>(data);
130
131			if(!zConsParaData){
132
133			sprintf( painCave.errMsg,
134			"ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
135			painCave.isFatal = 1;
136			simError();
137
138			}
139			else{
140
141			parameters = zConsParaData->getData();
142
143			//check the range of zconsIndex
144			//and the minimum value of index is the first one (we already sorted the data)
145			//the maximum value of index is the last one
146
147			int maxIndex;
148			int minIndex;
149			int totalNumMol;
150
151			minIndex = (*parameters)[0].zconsIndex;
152			if(minIndex < 0){
153			sprintf( painCave.errMsg,
154			"ZConstraint error: index is out of range\n");
155			painCave.isFatal = 1;
156			simError();
157			}
158
159			maxIndex = (*parameters)[parameters->size()].zconsIndex;
160
161			#ifndef IS_MPI
162			totalNumMol = nMols;
163			#else
164			totalNumMol = mpiSim->getTotNmol();
165			#endif
166
167			if(maxIndex > totalNumMol - 1){
168			sprintf( painCave.errMsg,
169			"ZConstraint error: index is out of range\n");
170			painCave.isFatal = 1;
171			simError();
172			}
173
174			//if user does not specify the zpos for the zconstraint molecule
175			//its initial z coordinate will be used as default
176			for(int i = 0; i < parameters->size(); i++){
177
178			if(!(*parameters)[i].havingZPos){
179
180			#ifndef IS_MPI
181			for(int j = 0; j < nMols; j++){
182			if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){
183			molecules[i].getCOM(COM);
184			break;
185			}
186			}
187			#else
188			//query which processor current zconstraint molecule belongs to
189			int *MolToProcMap;
190			int whichNode;
191			double initZPos;
192			MolToProcMap = mpiSim->getMolToProcMap();
193			whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
194
195			//broadcast the zpos of current z-contraint molecule
196			//the node which contain this
197
198			if (worldRank == whichNode ){
199
200			for(int i = 0; i < nMols; i++)
201			if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){
202			molecules[i].getCOM(COM);
203			break;
204			}
205
206			}
207
208			MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);
209			#endif
210
211			(*parameters)[i].zPos = COM[whichDirection];
212
213			sprintf( painCave.errMsg,
214			"ZConstraint warningr: Does not specify zpos for z-constraint molecule "
215			"initial z coornidate will be used \n");
216			painCave.isFatal = 0;
217			simError();
218
219			}
220			}
221
222			}//end if (!zConsParaData)
223			}//end if (!data)
224
225			//
226	tim	658	#ifdef IS_MPI
227			update();
228			#else
229			int searchResult;
230	tim	676
231	tim	658	for(int i = 0; i < nMols; i++){
232
233			searchResult = isZConstraintMol(&molecules[i]);
234
235			if(searchResult > -1){
236
237			zconsMols.push_back(&molecules[i]);
238			massOfZConsMols.push_back(molecules[i].getTotalMass());
239	tim	682
240			zPos.push_back((*parameters)[searchResult].zPos);
241			kz.push_back((parameters)[searchResult]. kRatio zForceConst);
242	tim	676
243	tim	658	molecules[i].getCOM(COM);
244			}
245			else
246			{
247
248			unconsMols.push_back(&molecules[i]);
249			massOfUnconsMols.push_back(molecules[i].getTotalMass());
250
251			}
252			}
253
254			fz = new double[zconsMols.size()];
255			indexOfZConsMols = new int [zconsMols.size()];
256
257			if(!fz \|\| !indexOfZConsMols){
258			sprintf( painCave.errMsg,
259			"Memory allocation failure in class Zconstraint\n");
260			painCave.isFatal = 1;
261			simError();
262			}
263
264			for(int i = 0; i < zconsMols.size(); i++)
265			indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
266
267			#endif
268	tim	676
269			//get total number of unconstrained atoms
270			int nUnconsAtoms_local;
271			nUnconsAtoms_local = 0;
272			for(int i = 0; i < unconsMols.size(); i++)
273			nUnconsAtoms_local += unconsMols[i]->getNAtoms();
274
275			#ifndef IS_MPI
276			totNumOfUnconsAtoms = nUnconsAtoms_local;
277			#else
278			MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
279			#endif
280
281	tim	682	checkZConsState();
282	tim	676
283	tim	682	//
284	tim	658	fzOut = new ZConsWriter(zconsOutput.c_str());
285
286			if(!fzOut){
287			sprintf( painCave.errMsg,
288			"Memory allocation failure in class Zconstraint\n");
289			painCave.isFatal = 1;
290			simError();
291			}
292
293			}
294
295			template<typename T> ZConstraint<T>::~ZConstraint()
296			{
297			if(fz)
298			delete[] fz;
299
300			if(indexOfZConsMols)
301			delete[] indexOfZConsMols;
302
303			if(fzOut)
304			delete fzOut;
305			}
306
307			#ifdef IS_MPI
308			template<typename T> void ZConstraint<T>::update()
309			{
310			double COM[3];
311			int index;
312
313			zconsMols.clear();
314			massOfZConsMols.clear();
315	tim	682	zPos.clear();
316			kz.clear();
317	tim	658
318			unconsMols.clear();
319			massOfUnconsMols.clear();
320
321
322			//creat zconsMol and unconsMol lists
323			for(int i = 0; i < nMols; i++){
324
325			index = isZConstraintMol(&molecules[i]);
326
327			if(index > -1){
328
329			zconsMols.push_back(&molecules[i]);
330	tim	682	zPos.push_back((*parameters)[index].zPos);
331			kz.push_back((parameters)[index].kRatio zForceConst);
332
333	tim	658	massOfZConsMols.push_back(molecules[i].getTotalMass());
334
335			molecules[i].getCOM(COM);
336			}
337			else
338			{
339
340			unconsMols.push_back(&molecules[i]);
341			massOfUnconsMols.push_back(molecules[i].getTotalMass());
342
343			}
344			}
345
346			//The reason to declare fz and indexOfZconsMols as pointer to array is
347			// that we want to make the MPI communication simple
348			if(fz)
349			delete[] fz;
350
351			if(indexOfZConsMols)
352			delete[] indexOfZConsMols;
353
354			if (zconsMols.size() > 0){
355			fz = new double[zconsMols.size()];
356			indexOfZConsMols = new int[zconsMols.size()];
357
358			if(!fz \|\| !indexOfZConsMols){
359			sprintf( painCave.errMsg,
360			"Memory allocation failure in class Zconstraint\n");
361			painCave.isFatal = 1;
362			simError();
363			}
364
365			for(int i = 0; i < zconsMols.size(); i++){
366			indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
367			}
368
369			}
370			else{
371			fz = NULL;
372			indexOfZConsMols = NULL;
373			}
374
375			}
376
377			#endif
378
379			/** Function Name: isZConstraintMol
380			** Parameter
381			** Molecule* mol
382			** Return value:
383			** -1, if the molecule is not z-constraint molecule,
384			** other non-negative values, its index in indexOfAllZConsMols vector
385			*/
386
387			template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
388			{
389			int index;
390			int low;
391			int high;
392			int mid;
393
394			index = mol->getGlobalIndex();
395
396			low = 0;
397	tim	682	high = parameters->size() - 1;
398	tim	658
399			//Binary Search (we have sorted the array)
400			while(low <= high){
401			mid = (low + high) /2;
402	tim	682	if ((*parameters)[mid].zconsIndex == index)
403	tim	658	return mid;
404	tim	682	else if ((*parameters)[mid].zconsIndex > index )
405	tim	658	high = mid -1;
406			else
407			low = mid + 1;
408			}
409
410			return -1;
411			}
412
413	tim	676	/**
414			* Description:
415			* Reset the z coordinates
416	tim	658	*/
417	tim	676	template<typename T> void ZConstraint<T>::integrate(){
418	tim	658
419	tim	676	//zero out the velocities of center of mass of unconstrained molecules
420			//and the velocities of center of mass of every single z-constrained molecueles
421			zeroOutVel();
422
423			T::integrate();
424
425	tim	658	}
426	tim	676
427	tim	658
428	tim	676	/**
429			*
430			*
431			*
432			*
433	tim	658	*/
434
435	tim	676
436			template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
437
438			T::calcForce(calcPot, calcStress);
439
440			if (checkZConsState())
441			zeroOutVel();
442
443			//do zconstraint force;
444			if (haveFixedZMols())
445			this->doZconstraintForce();
446
447			//use harmonical poteintial to move the molecules to the specified positions
448			if (haveMovingZMols())
449			//this->doHarmonic();
450
451			fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
452
453			}
454
455			template<typename T> double ZConstraint<T>::calcZSys()
456	tim	658	{
457	tim	676	//calculate reference z coordinate for z-constraint molecules
458			double totalMass_local;
459			double totalMass;
460			double totalMZ_local;
461			double totalMZ;
462			double massOfUncons_local;
463			double massOfCurMol;
464			double COM[3];
465
466			totalMass_local = 0;
467			totalMass = 0;
468			totalMZ_local = 0;
469			totalMZ = 0;
470			massOfUncons_local = 0;
471
472
473			for(int i = 0; i < nMols; i++){
474			massOfCurMol = molecules[i].getTotalMass();
475			molecules[i].getCOM(COM);
476
477			totalMass_local += massOfCurMol;
478			totalMZ_local += massOfCurMol * COM[whichDirection];
479
480			if(isZConstraintMol(&molecules[i]) == -1){
481
482			massOfUncons_local += massOfCurMol;
483			}
484
485			}
486
487
488			#ifdef IS_MPI
489			MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
490			MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
491			MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
492			#else
493			totalMass = totalMass_local;
494			totalMZ = totalMZ_local;
495			totalMassOfUncons = massOfUncons_local;
496			#endif
497	mmeineke	671
498	tim	658	double zsys;
499	tim	676	zsys = totalMZ / totalMass;
500	tim	658
501	tim	676	return zsys;
502			}
503
504			/**
505			*
506			*/
507			template<typename T> void ZConstraint<T>::thermalize( void ){
508
509			T::thermalize();
510			zeroOutVel();
511			}
512
513			/**
514			*
515			*
516			*
517			*/
518
519			template<typename T> void ZConstraint<T>::zeroOutVel(){
520
521			Atom** fixedZAtoms;
522			double COMvel[3];
523			double vel[3];
524	tim	658
525	tim	676	//zero out the velocities of center of mass of fixed z-constrained molecules
526
527	tim	658	for(int i = 0; i < zconsMols.size(); i++){
528	tim	676
529			if (states[i] == zcsFixed){
530
531			zconsMols[i]->getCOMvel(COMvel);
532			fixedZAtoms = zconsMols[i]->getMyAtoms();
533
534			for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
535			fixedZAtoms[j]->getVel(vel);
536			vel[whichDirection] -= COMvel[whichDirection];
537			fixedZAtoms[j]->setVel(vel);
538			}
539
540			}
541
542	tim	658	}
543	tim	676
544			// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
545			double MVzOfMovingMols_local;
546			double MVzOfMovingMols;
547			double totalMassOfMovingZMols_local;
548			double totalMassOfMovingZMols;
549
550			MVzOfMovingMols_local = 0;
551			totalMassOfMovingZMols_local = 0;
552	tim	658
553	tim	676	for(int i =0; i < unconsMols.size(); i++){
554			unconsMols[i]->getCOMvel(COMvel);
555			MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
556			}
557
558			for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
559
560			if (states[i] == zcsMoving){
561			zconsMols[i]->getCOMvel(COMvel);
562			MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
563			totalMassOfMovingZMols_local += massOfZConsMols[i];
564			}
565
566			}
567
568			#ifndef IS_MPI
569			MVzOfMovingMols = MVzOfMovingMols_local;
570			totalMassOfMovingZMols = totalMassOfMovingZMols_local;
571	tim	658	#else
572	tim	676	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
573			MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
574	tim	658	#endif
575
576	tim	676	double vzOfMovingMols;
577			vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
578
579			//modify the velocites of unconstrained molecules
580			Atom** unconsAtoms;
581	tim	658	for(int i = 0; i < unconsMols.size(); i++){
582	tim	676
583			unconsAtoms = unconsMols[i]->getMyAtoms();
584			for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
585			unconsAtoms[j]->getVel(vel);
586			vel[whichDirection] -= vzOfMovingMols;
587			unconsAtoms[j]->setVel(vel);
588			}
589
590			}
591
592			//modify the velocities of moving z-constrained molecuels
593			Atom** movingZAtoms;
594			for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
595
596			if (states[i] ==zcsMoving){
597
598			movingZAtoms = zconsMols[i]->getMyAtoms();
599			for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
600			movingZAtoms[j]->getVel(vel);
601			vel[whichDirection] -= vzOfMovingMols;
602			movingZAtoms[j]->setVel(vel);
603			}
604
605			}
606
607	tim	658	}
608	tim	676
609			}
610
611			template<typename T> void ZConstraint<T>::doZconstraintForce(){
612
613			Atom** zconsAtoms;
614			double totalFZ;
615			double totalFZ_local;
616			double COMvel[3];
617			double COM[3];
618			double force[3];
619			double zsys;
620
621			int nMovingZMols_local;
622			int nMovingZMols;
623
624			//constrain the molecules which do not reach the specified positions
625
626			zsys = calcZSys();
627			cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;
628
629			//Zero Out the force of z-contrained molecules
630			totalFZ_local = 0;
631
632			//calculate the total z-contrained force of fixed z-contrained molecules
633	tim	682	cout << "Fixed Molecules" << endl;
634	tim	676	for(int i = 0; i < zconsMols.size(); i++){
635
636			if (states[i] == zcsFixed){
637
638			zconsMols[i]->getCOM(COM);
639			zconsAtoms = zconsMols[i]->getMyAtoms();
640
641			fz[i] = 0;
642			for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
643			zconsAtoms[j]->getFrc(force);
644			fz[i] += force[whichDirection];
645			}
646			totalFZ_local += fz[i];
647
648			cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
649
650			}
651
652			}
653
654			//calculate the number of atoms of moving z-constrained molecules
655			nMovingZMols_local = 0;
656			for(int i = 0; zconsMols.size(); i++){
657			if(states[i] == zcsMoving)
658			nMovingZMols_local += massOfZConsMols[i];
659			}
660	tim	658	#ifdef IS_MPI
661	tim	676	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
662			MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
663	tim	658	#else
664	tim	676	totalFZ = totalFZ_local;
665			nMovingZMols = nMovingZMols_local;
666			#endif
667
668			force[0]= 0;
669			force[1]= 0;
670			force[2]= 0;
671			force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);
672
673			//modify the velocites of unconstrained molecules
674			for(int i = 0; i < unconsMols.size(); i++){
675
676			Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
677
678			for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)
679			unconsAtoms[j]->addFrc(force);
680
681			}
682
683			//modify the velocities of moving z-constrained molecules
684			for(int i = 0; i < zconsMols.size(); i++) {
685			if (states[i] == zcsMoving){
686
687			Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
688
689			for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)
690			movingZAtoms[j]->addFrc(force);
691			}
692			}
693
694			// apply negative to fixed z-constrained molecues;
695			force[0]= 0;
696			force[1]= 0;
697			force[2]= 0;
698
699			for(int i = 0; i < zconsMols.size(); i++){
700
701			if (states[i] == zcsFixed){
702
703			int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
704			zconsAtoms = zconsMols[i]->getMyAtoms();
705
706			for(int j =0; j < nAtomOfCurZConsMol; j++) {
707			force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
708			zconsAtoms[j]->addFrc(force);
709			}
710
711			}
712
713			}
714
715			}
716
717			template<typename T> bool ZConstraint<T>::checkZConsState(){
718			double COM[3];
719			double diff;
720	tim	658
721	tim	676	bool changed;
722
723			changed = false;
724
725			for(int i =0; i < zconsMols.size(); i++){
726
727	tim	658	zconsMols[i]->getCOM(COM);
728	tim	682	diff = fabs(COM[whichDirection] - zPos[i]);
729			if ( diff <= zconsTol && states[i] == zcsMoving){
730	tim	676	states[i] = zcsFixed;
731			changed = true;
732			}
733	tim	682	else if ( diff > zconsTol && states[i] == zcsFixed){
734	tim	676	states[i] = zcsMoving;
735			changed = true;
736			}
737
738	tim	658	}
739
740	tim	676	return changed;
741	tim	658	}
742	tim	676
743			template<typename T> bool ZConstraint<T>::haveFixedZMols(){
744			for(int i = 0; i < zconsMols.size(); i++)
745			if (states[i] == zcsFixed)
746			return true;
747
748			return false;
749			}
750
751
752			/**
753			*
754			*/
755			template<typename T> bool ZConstraint<T>::haveMovingZMols(){
756			for(int i = 0; i < zconsMols.size(); i++)
757			if (states[i] == zcsMoving)
758			return true;
759
760			return false;
761
762	tim	682	}
763
764			/**
765			*
766			*
767			*/
768
769			template<typename T> void ZConstraint<T>::doHarmonic(){
770			double force[3];
771			double harmonicU;
772			double COM[3];
773			double diff;
774
775			force[0] = 0;
776			force[1] = 0;
777			force[2] = 0;
778
779			cout << "Moving Molecules" << endl;
780			for(int i = 0; i < zconsMols.size(); i++) {
781
782			if (states[i] == zcsMoving){
783			zconsMols[i]->getCOM(COM):
784			cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
785
786			diff = COM[whichDirection] -zPos[i];
787
788			harmonicU = 0.5 * kz[i] * diff * diff;
789			info->ltPot += harmonicU;
790
791			force[whichDirection] = - kz[i] * diff / zconsMols[i]->getNAtoms();
792
793			Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
794
795			for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)
796			movingZAtoms[j]->addFrc(force);
797			}
798
799			}
800
801			}