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;
  IndexData* index;
  DoubleData* sampleTime;
  StringData* filename; 
  
  //retrieve index of z-constraint molecules
  data = info->getProperty("zconsindex");
  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{
    index = dynamic_cast<IndexData*>(data);
    
    if(!index){

      sprintf( painCave.errMsg,
                 "ZConstraint error: Can not get property from SimInfo\n");
      painCave.isFatal = 1;
      simError();  
    
    }
    else{
           
      indexOfAllZConsMols = index->getIndexData();
      
      //the maximum value of index is the last one(we sorted the index data in SimSetup.cpp)
      int maxIndex;
        int minIndex;
      int totalNumMol;

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

#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(); 
                 
      }
      
    }
         
  }
  
  //retrieve sample time of z-contraint 
  data = info->getProperty("zconstime");
  
  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");
  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();
    }
    

  }
  
#ifdef IS_MPI
  update(); 
#else  
  int searchResult;
  double COM[3];
      
  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());  
      
      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


  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();
  
  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]);      
      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 = indexOfAllZConsMols.size() - 1;
  
  //Binary Search (we have sorted the array)  
  while(low <= high){
    mid = (low + high) /2;
    if (indexOfAllZConsMols[mid] == index)
      return mid;
    else if (indexOfAllZConsMols[mid] > 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
  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 <= ztol && states[i] == zcsMoving){
      states[i] = zcsFixed;
        changed = true;
    }
    else if ( diff > ztol && 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;
  
}
Revision:	676
Committed:	Mon Aug 11 19:40:06 2003 UTC (20 years, 11 months ago) by tim
File size:	15171 byte(s)
Log Message:	* empty log message *
#	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			IndexData* index;
12			DoubleData* sampleTime;
13			StringData* filename;
14
15	tim	676	//retrieve index of z-constraint molecules
16	tim	658	data = info->getProperty("zconsindex");
17			if(!data) {
18
19			sprintf( painCave.errMsg,
20			"ZConstraint error: If you use an ZConstraint\n"
21			" , you must set index of z-constraint molecules.\n");
22			painCave.isFatal = 1;
23			simError();
24			}
25			else{
26			index = dynamic_cast<IndexData*>(data);
27
28			if(!index){
29
30			sprintf( painCave.errMsg,
31			"ZConstraint error: Can not get property from SimInfo\n");
32			painCave.isFatal = 1;
33			simError();
34
35			}
36			else{
37	tim	660
38	tim	658	indexOfAllZConsMols = index->getIndexData();
39	tim	660
40			//the maximum value of index is the last one(we sorted the index data in SimSetup.cpp)
41			int maxIndex;
42	tim	676	int minIndex;
43	tim	660	int totalNumMol;
44	tim	676
45			minIndex = indexOfAllZConsMols[0];
46			if(minIndex < 0){
47			sprintf( painCave.errMsg,
48			"ZConstraint error: index is out of range\n");
49			painCave.isFatal = 1;
50			simError();
51			}
52
53	tim	660	maxIndex = indexOfAllZConsMols[indexOfAllZConsMols.size() - 1];
54
55			#ifndef IS_MPI
56			totalNumMol = nMols;
57			#else
58			totalNumMol = mpiSim->getTotNmol();
59			#endif
60
61			if(maxIndex > totalNumMol - 1){
62			sprintf( painCave.errMsg,
63			"ZConstraint error: index is out of range\n");
64			painCave.isFatal = 1;
65			simError();
66
67			}
68
69	tim	658	}
70
71			}
72
73	tim	676	//retrieve sample time of z-contraint
74	tim	658	data = info->getProperty("zconstime");
75
76			if(!data) {
77
78			sprintf( painCave.errMsg,
79			"ZConstraint error: If you use an ZConstraint\n"
80			" , you must set sample time.\n");
81			painCave.isFatal = 1;
82			simError();
83			}
84			else{
85
86			sampleTime = dynamic_cast<DoubleData*>(data);
87
88			if(!sampleTime){
89
90			sprintf( painCave.errMsg,
91			"ZConstraint error: Can not get property from SimInfo\n");
92			painCave.isFatal = 1;
93			simError();
94
95			}
96			else{
97			this->zconsTime = sampleTime->getData();
98			}
99
100			}
101
102
103	tim	676	//retrieve output filename of z force
104	tim	658	data = info->getProperty("zconsfilename");
105			if(!data) {
106
107
108			sprintf( painCave.errMsg,
109			"ZConstraint error: If you use an ZConstraint\n"
110			" , you must set output filename of z-force.\n");
111			painCave.isFatal = 1;
112			simError();
113
114			}
115			else{
116
117			filename = dynamic_cast<StringData*>(data);
118
119			if(!filename){
120
121			sprintf( painCave.errMsg,
122			"ZConstraint error: Can not get property from SimInfo\n");
123			painCave.isFatal = 1;
124			simError();
125
126			}
127			else{
128			this->zconsOutput = filename->getData();
129			}
130
131
132			}
133
134			#ifdef IS_MPI
135			update();
136			#else
137			int searchResult;
138	tim	676	double COM[3];
139
140	tim	658	for(int i = 0; i < nMols; i++){
141
142			searchResult = isZConstraintMol(&molecules[i]);
143
144			if(searchResult > -1){
145
146			zconsMols.push_back(&molecules[i]);
147			massOfZConsMols.push_back(molecules[i].getTotalMass());
148	tim	676
149	tim	658	molecules[i].getCOM(COM);
150			}
151			else
152			{
153
154			unconsMols.push_back(&molecules[i]);
155			massOfUnconsMols.push_back(molecules[i].getTotalMass());
156
157			}
158			}
159
160			fz = new double[zconsMols.size()];
161			indexOfZConsMols = new int [zconsMols.size()];
162
163			if(!fz \|\| !indexOfZConsMols){
164			sprintf( painCave.errMsg,
165			"Memory allocation failure in class Zconstraint\n");
166			painCave.isFatal = 1;
167			simError();
168			}
169
170			for(int i = 0; i < zconsMols.size(); i++)
171			indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
172
173			#endif
174	tim	676
175			//get total number of unconstrained atoms
176			int nUnconsAtoms_local;
177			nUnconsAtoms_local = 0;
178			for(int i = 0; i < unconsMols.size(); i++)
179			nUnconsAtoms_local += unconsMols[i]->getNAtoms();
180
181			#ifndef IS_MPI
182			totNumOfUnconsAtoms = nUnconsAtoms_local;
183			#else
184			MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
185			#endif
186
187
188	tim	658
189			fzOut = new ZConsWriter(zconsOutput.c_str());
190
191			if(!fzOut){
192			sprintf( painCave.errMsg,
193			"Memory allocation failure in class Zconstraint\n");
194			painCave.isFatal = 1;
195			simError();
196			}
197
198			}
199
200			template<typename T> ZConstraint<T>::~ZConstraint()
201			{
202			if(fz)
203			delete[] fz;
204
205			if(indexOfZConsMols)
206			delete[] indexOfZConsMols;
207
208			if(fzOut)
209			delete fzOut;
210			}
211
212			#ifdef IS_MPI
213			template<typename T> void ZConstraint<T>::update()
214			{
215			double COM[3];
216			int index;
217
218			zconsMols.clear();
219			massOfZConsMols.clear();
220
221			unconsMols.clear();
222			massOfUnconsMols.clear();
223
224
225			//creat zconsMol and unconsMol lists
226			for(int i = 0; i < nMols; i++){
227
228			index = isZConstraintMol(&molecules[i]);
229
230			if(index > -1){
231
232			zconsMols.push_back(&molecules[i]);
233			massOfZConsMols.push_back(molecules[i].getTotalMass());
234
235			molecules[i].getCOM(COM);
236			}
237			else
238			{
239
240			unconsMols.push_back(&molecules[i]);
241			massOfUnconsMols.push_back(molecules[i].getTotalMass());
242
243			}
244			}
245
246			//The reason to declare fz and indexOfZconsMols as pointer to array is
247			// that we want to make the MPI communication simple
248			if(fz)
249			delete[] fz;
250
251			if(indexOfZConsMols)
252			delete[] indexOfZConsMols;
253
254			if (zconsMols.size() > 0){
255			fz = new double[zconsMols.size()];
256			indexOfZConsMols = new int[zconsMols.size()];
257
258			if(!fz \|\| !indexOfZConsMols){
259			sprintf( painCave.errMsg,
260			"Memory allocation failure in class Zconstraint\n");
261			painCave.isFatal = 1;
262			simError();
263			}
264
265			for(int i = 0; i < zconsMols.size(); i++){
266			indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
267			}
268
269			}
270			else{
271			fz = NULL;
272			indexOfZConsMols = NULL;
273			}
274
275			}
276
277			#endif
278
279			/** Function Name: isZConstraintMol
280			** Parameter
281			** Molecule* mol
282			** Return value:
283			** -1, if the molecule is not z-constraint molecule,
284			** other non-negative values, its index in indexOfAllZConsMols vector
285			*/
286
287			template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
288			{
289			int index;
290			int low;
291			int high;
292			int mid;
293
294			index = mol->getGlobalIndex();
295
296			low = 0;
297			high = indexOfAllZConsMols.size() - 1;
298
299			//Binary Search (we have sorted the array)
300			while(low <= high){
301			mid = (low + high) /2;
302			if (indexOfAllZConsMols[mid] == index)
303			return mid;
304			else if (indexOfAllZConsMols[mid] > index )
305			high = mid -1;
306			else
307			low = mid + 1;
308			}
309
310			return -1;
311			}
312
313	tim	676	/**
314			* Description:
315			* Reset the z coordinates
316	tim	658	*/
317	tim	676	template<typename T> void ZConstraint<T>::integrate(){
318	tim	658
319	tim	676	//zero out the velocities of center of mass of unconstrained molecules
320			//and the velocities of center of mass of every single z-constrained molecueles
321			zeroOutVel();
322
323			T::integrate();
324
325	tim	658	}
326	tim	676
327	tim	658
328	tim	676	/**
329			*
330			*
331			*
332			*
333	tim	658	*/
334
335	tim	676
336			template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
337
338			T::calcForce(calcPot, calcStress);
339
340			if (checkZConsState())
341			zeroOutVel();
342
343			//do zconstraint force;
344			if (haveFixedZMols())
345			this->doZconstraintForce();
346
347			//use harmonical poteintial to move the molecules to the specified positions
348			if (haveMovingZMols())
349			//this->doHarmonic();
350
351			fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
352
353			}
354
355			template<typename T> double ZConstraint<T>::calcZSys()
356	tim	658	{
357	tim	676	//calculate reference z coordinate for z-constraint molecules
358			double totalMass_local;
359			double totalMass;
360			double totalMZ_local;
361			double totalMZ;
362			double massOfUncons_local;
363			double massOfCurMol;
364			double COM[3];
365
366			totalMass_local = 0;
367			totalMass = 0;
368			totalMZ_local = 0;
369			totalMZ = 0;
370			massOfUncons_local = 0;
371
372
373			for(int i = 0; i < nMols; i++){
374			massOfCurMol = molecules[i].getTotalMass();
375			molecules[i].getCOM(COM);
376
377			totalMass_local += massOfCurMol;
378			totalMZ_local += massOfCurMol * COM[whichDirection];
379
380			if(isZConstraintMol(&molecules[i]) == -1){
381
382			massOfUncons_local += massOfCurMol;
383			}
384
385			}
386
387
388			#ifdef IS_MPI
389			MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
390			MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
391			MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
392			#else
393			totalMass = totalMass_local;
394			totalMZ = totalMZ_local;
395			totalMassOfUncons = massOfUncons_local;
396			#endif
397	mmeineke	671
398	tim	658	double zsys;
399	tim	676	zsys = totalMZ / totalMass;
400	tim	658
401	tim	676	return zsys;
402			}
403
404			/**
405			*
406			*/
407			template<typename T> void ZConstraint<T>::thermalize( void ){
408
409			T::thermalize();
410			zeroOutVel();
411			}
412
413			/**
414			*
415			*
416			*
417			*/
418
419			template<typename T> void ZConstraint<T>::zeroOutVel(){
420
421			Atom** fixedZAtoms;
422			double COMvel[3];
423			double vel[3];
424	tim	658
425	tim	676	//zero out the velocities of center of mass of fixed z-constrained molecules
426
427	tim	658	for(int i = 0; i < zconsMols.size(); i++){
428	tim	676
429			if (states[i] == zcsFixed){
430
431			zconsMols[i]->getCOMvel(COMvel);
432			fixedZAtoms = zconsMols[i]->getMyAtoms();
433
434			for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
435			fixedZAtoms[j]->getVel(vel);
436			vel[whichDirection] -= COMvel[whichDirection];
437			fixedZAtoms[j]->setVel(vel);
438			}
439
440			}
441
442	tim	658	}
443	tim	676
444			// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
445			double MVzOfMovingMols_local;
446			double MVzOfMovingMols;
447			double totalMassOfMovingZMols_local;
448			double totalMassOfMovingZMols;
449
450			MVzOfMovingMols_local = 0;
451			totalMassOfMovingZMols_local = 0;
452	tim	658
453	tim	676	for(int i =0; i < unconsMols.size(); i++){
454			unconsMols[i]->getCOMvel(COMvel);
455			MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
456			}
457
458			for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
459
460			if (states[i] == zcsMoving){
461			zconsMols[i]->getCOMvel(COMvel);
462			MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
463			totalMassOfMovingZMols_local += massOfZConsMols[i];
464			}
465
466			}
467
468			#ifndef IS_MPI
469			MVzOfMovingMols = MVzOfMovingMols_local;
470			totalMassOfMovingZMols = totalMassOfMovingZMols_local;
471	tim	658	#else
472	tim	676	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
473			MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
474	tim	658	#endif
475
476	tim	676	double vzOfMovingMols;
477			vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
478
479			//modify the velocites of unconstrained molecules
480			Atom** unconsAtoms;
481	tim	658	for(int i = 0; i < unconsMols.size(); i++){
482	tim	676
483			unconsAtoms = unconsMols[i]->getMyAtoms();
484			for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
485			unconsAtoms[j]->getVel(vel);
486			vel[whichDirection] -= vzOfMovingMols;
487			unconsAtoms[j]->setVel(vel);
488			}
489
490			}
491
492			//modify the velocities of moving z-constrained molecuels
493			Atom** movingZAtoms;
494			for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
495
496			if (states[i] ==zcsMoving){
497
498			movingZAtoms = zconsMols[i]->getMyAtoms();
499			for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
500			movingZAtoms[j]->getVel(vel);
501			vel[whichDirection] -= vzOfMovingMols;
502			movingZAtoms[j]->setVel(vel);
503			}
504
505			}
506
507	tim	658	}
508	tim	676
509			}
510
511			template<typename T> void ZConstraint<T>::doZconstraintForce(){
512
513			Atom** zconsAtoms;
514			double totalFZ;
515			double totalFZ_local;
516			double COMvel[3];
517			double COM[3];
518			double force[3];
519			double zsys;
520
521			int nMovingZMols_local;
522			int nMovingZMols;
523
524			//constrain the molecules which do not reach the specified positions
525
526			zsys = calcZSys();
527			cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;
528
529			//Zero Out the force of z-contrained molecules
530			totalFZ_local = 0;
531
532			//calculate the total z-contrained force of fixed z-contrained molecules
533			for(int i = 0; i < zconsMols.size(); i++){
534
535			if (states[i] == zcsFixed){
536
537			zconsMols[i]->getCOM(COM);
538			zconsAtoms = zconsMols[i]->getMyAtoms();
539
540			fz[i] = 0;
541			for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
542			zconsAtoms[j]->getFrc(force);
543			fz[i] += force[whichDirection];
544			}
545			totalFZ_local += fz[i];
546
547			cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
548
549			}
550
551			}
552
553			//calculate the number of atoms of moving z-constrained molecules
554			nMovingZMols_local = 0;
555			for(int i = 0; zconsMols.size(); i++){
556			if(states[i] == zcsMoving)
557			nMovingZMols_local += massOfZConsMols[i];
558			}
559	tim	658	#ifdef IS_MPI
560	tim	676	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
561			MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
562	tim	658	#else
563	tim	676	totalFZ = totalFZ_local;
564			nMovingZMols = nMovingZMols_local;
565			#endif
566
567			force[0]= 0;
568			force[1]= 0;
569			force[2]= 0;
570			force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);
571
572			//modify the velocites of unconstrained molecules
573			for(int i = 0; i < unconsMols.size(); i++){
574
575			Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
576
577			for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)
578			unconsAtoms[j]->addFrc(force);
579
580			}
581
582			//modify the velocities of moving z-constrained molecules
583			for(int i = 0; i < zconsMols.size(); i++) {
584			if (states[i] == zcsMoving){
585
586			Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
587
588			for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)
589			movingZAtoms[j]->addFrc(force);
590			}
591			}
592
593			// apply negative to fixed z-constrained molecues;
594			force[0]= 0;
595			force[1]= 0;
596			force[2]= 0;
597
598			for(int i = 0; i < zconsMols.size(); i++){
599
600			if (states[i] == zcsFixed){
601
602			int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
603			zconsAtoms = zconsMols[i]->getMyAtoms();
604
605			for(int j =0; j < nAtomOfCurZConsMol; j++) {
606			force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
607			zconsAtoms[j]->addFrc(force);
608			}
609
610			}
611
612			}
613
614			}
615
616			template<typename T> bool ZConstraint<T>::checkZConsState(){
617			double COM[3];
618			double diff;
619	tim	658
620	tim	676	bool changed;
621
622			changed = false;
623
624			for(int i =0; i < zconsMols.size(); i++){
625
626	tim	658	zconsMols[i]->getCOM(COM);
627	tim	676	diff = fabs(COM[whichDirection] - ZPos[i]);
628			if ( diff <= ztol && states[i] == zcsMoving){
629			states[i] = zcsFixed;
630			changed = true;
631			}
632			else if ( diff > ztol && states[i] == zcsFixed){
633			states[i] = zcsMoving;
634			changed = true;
635			}
636
637	tim	658	}
638
639	tim	676	return changed;
640	tim	658	}
641	tim	676
642			template<typename T> bool ZConstraint<T>::haveFixedZMols(){
643			for(int i = 0; i < zconsMols.size(); i++)
644			if (states[i] == zcsFixed)
645			return true;
646
647			return false;
648			}
649
650
651			/**
652			*
653			*/
654			template<typename T> bool ZConstraint<T>::haveMovingZMols(){
655			for(int i = 0; i < zconsMols.size(); i++)
656			if (states[i] == zcsMoving)
657			return true;
658
659			return false;
660
661			}