OOPSE/libmdtools/ZConstraint.cpp

#include "Integrator.hpp"
#include "simError.h"

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; 
  
  
  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();
    }
         
  }
  
  //retrive 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();
    }

  }
  
  
  //retrive 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();
    }
    

  }

 
  //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[2];
    
    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;

#ifndef IS_MPI  
  for(int i = 0; i < nMols; i++){
    
    if(isZConstraintMol(&molecules[i]) > -1 ){
      molecules[i].getCOM(COM); 
      allRefZ.push_back(COM[2] - zsys);  
    }
    
  }
#else
 
  int whichNode;
  enum CommType { RequestMolZPos, EndOfRequest} status;
  //int status;
  double zpos;
  int localIndex;
  MPI_Status ierr;
  int tag = 0;
  
  if(worldRank == 0){
    
    int globalIndexOfCurMol;
    int *MolToProcMap;
    MolToProcMap = mpiSim->getMolToProcMap();
    
    for(int i = 0; i < indexOfAllZConsMols.size(); i++){
      
      whichNode = MolToProcMap[indexOfAllZConsMols[i]];
      globalIndexOfCurMol = indexOfAllZConsMols[i];
      
      if(whichNode == 0){
        
        for(int j = 0; j < nMols; j++)
          if(molecules[j].getGlobalIndex() == globalIndexOfCurMol){
            localIndex = j;
            break;
          }
                  
        molecules[localIndex].getCOM(COM);
        allRefZ.push_back(COM[2] - zsys);  
              
      }
      else{
        status = RequestMolZPos;
        MPI_Send(&status, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
        MPI_Send(&globalIndexOfCurMol, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
        MPI_Recv(&zpos, 1, MPI_DOUBLE_PRECISION, whichNode, tag, MPI_COMM_WORLD, &ierr);
        
        allRefZ.push_back(zpos - zsys);
      
      }
              
    } //End of Request Loop
    
    //Send ending request message to slave nodes    
    status = EndOfRequest;
    for(int i =1; i < mpiSim->getNumberProcessors(); i++)
      MPI_Send(&status, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
     
  }
  else{
  
    int whichMol;
    bool done = false;

    while (!done){  
      
      MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
    
      switch (status){
          
        case RequestMolZPos : 
          
          MPI_Recv(&whichMol, 1, MPI_INT, 0, tag, MPI_COMM_WORLD,&ierr);
          
          for(int i = 0; i < nMols; i++)
            if(molecules[i].getGlobalIndex() == whichMol){
              localIndex = i;
              break;
            }
          
          molecules[localIndex].getCOM(COM);
          zpos = COM[2];          
          MPI_Send(&zpos, 1, MPI_DOUBLE_PRECISION, 0, tag, MPI_COMM_WORLD);       
          
          break;
             
        case EndOfRequest :
         
         done = true;
         break;
      }
      
    }
          
  }

  //Brocast the allRefZ to slave nodes;
  double* allRefZBuf;
  int nZConsMols;
  nZConsMols = indexOfAllZConsMols.size();
  
  allRefZBuf = new double[nZConsMols];
  
  if(worldRank == 0){

    for(int i = 0; i < nZConsMols; i++)
      allRefZBuf[i] = allRefZ[i];
  }    
   
    MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
  
  if(worldRank != 0){
    
    for(int i = 0; i < nZConsMols; i++)
      allRefZ.push_back(allRefZBuf[i]);  
  }
  
  delete[] allRefZBuf;
#endif

  
#ifdef IS_MPI
  update(); 
#else  
  int searchResult;
  
  refZ = allRefZ;

  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 
  
  fzOut = new ZConsWriter(zconsOutput.c_str());   
  
  if(!fzOut){
    sprintf( painCave.errMsg,
             "Memory allocation failure in class Zconstraint\n");
    painCave.isFatal = 1;
    simError();
  }
  
  fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
}

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();
  refZ.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);
      refZ.push_back(allRefZ[index]);      
    }
    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;
}

/** Function Name: integrateStep
 ** Parameter:
 **   int calcPot;
 **   int calcStress;
 ** Description:
 **  Advance One Step.
 ** Memo:
 **   The best way to implement z-constraint is to override integrateStep
 **   Overriding constrainB is not a good choice, since in integrateStep,
 **   constrainB is invoked by below line,
 **                  if(nConstrained) constrainB();
 **   For instance, we would like to apply z-constraint without bond contrain,
 **   In that case, if we override constrainB, Z-constrain method will never be executed;
 */
template<typename T> void ZConstraint<T>::integrateStep( int calcPot, int calcStress )
{
  T::integrateStep( calcPot, calcStress );
  resetZ();
  
  double currZConsTime = 0;
   
  //write out forces of z constraint
  if( info->getTime() >= currZConsTime){   
      fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); 
  }     
}

/** Function Name: resetZ
 ** Description:
 **  Reset the z coordinates
 */

template<typename T> void ZConstraint<T>::resetZ()
{
  double deltaZ;
  double mzOfZCons;   //total sum of m*z of z-constrain molecules
  double mzOfUncons; //total sum of m*z of unconstrain molecuels;
  double totalMZOfZCons;
  double totalMZOfUncons;
  double COM[3];
  double zsys;
  Atom** zconsAtoms;

  mzOfZCons = 0;
  mzOfUncons  = 0;
  
  for(int i = 0; i < zconsMols.size(); i++){
    mzOfZCons += massOfZConsMols[i] * refZ[i];    
  }

#ifdef IS_MPI
  MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#else
  totalMZOfZCons = mzOfZCons;
#endif

  for(int i = 0; i < unconsMols.size(); i++){
    unconsMols[i]->getCOM(COM);
    mzOfUncons += massOfUnconsMols[i] * COM[2];
  }
  
#ifdef IS_MPI
  MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#else
  totalMZOfUncons = mzOfUncons;
#endif  
  
  zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons; 

  cout << "current time: " << info->getTime() <<endl;  
  for(int i = 0; i < zconsMols.size(); i++){   
   
    zconsMols[i]->getCOM(COM);
    
    cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t";
    deltaZ = zsys + refZ[i] - COM[2];
    cout << "\tdistance: " << COM[2] +deltaZ - zsys;     
    //update z coordinate    
    zconsAtoms = zconsMols[i]->getMyAtoms();    
    for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
      zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ);  
    }    
    
    //calculate z constrain force
    fz[i] = massOfZConsMols[i]* deltaZ / dt2;
    
    cout << "\tforce: " << fz[i] << endl;
  }

      
}
Revision:	658
Committed:	Thu Jul 31 15:35:07 2003 UTC (20 years, 11 months ago) by tim
File size:	12159 byte(s)
Log Message:	Added Z constraint.
#	User	Rev	Content
1	tim	658	#include "Integrator.hpp"
2			#include "simError.h"
3
4			template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff)
5			: T(theInfo, the_ff), fz(NULL), indexOfZConsMols(NULL)
6			{
7
8			//get properties from SimInfo
9			GenericData* data;
10			IndexData* index;
11			DoubleData* sampleTime;
12			StringData* filename;
13
14
15			data = info->getProperty("zconsindex");
16			if(!data) {
17
18			sprintf( painCave.errMsg,
19			"ZConstraint error: If you use an ZConstraint\n"
20			" , you must set index of z-constraint molecules.\n");
21			painCave.isFatal = 1;
22			simError();
23			}
24			else{
25			index = dynamic_cast<IndexData*>(data);
26
27			if(!index){
28
29			sprintf( painCave.errMsg,
30			"ZConstraint error: Can not get property from SimInfo\n");
31			painCave.isFatal = 1;
32			simError();
33
34			}
35			else{
36			indexOfAllZConsMols = index->getIndexData();
37			}
38
39			}
40
41			//retrive sample time of z-contraint
42			data = info->getProperty("zconstime");
43
44			if(!data) {
45
46			sprintf( painCave.errMsg,
47			"ZConstraint error: If you use an ZConstraint\n"
48			" , you must set sample time.\n");
49			painCave.isFatal = 1;
50			simError();
51			}
52			else{
53
54			sampleTime = dynamic_cast<DoubleData*>(data);
55
56			if(!sampleTime){
57
58			sprintf( painCave.errMsg,
59			"ZConstraint error: Can not get property from SimInfo\n");
60			painCave.isFatal = 1;
61			simError();
62
63			}
64			else{
65			this->zconsTime = sampleTime->getData();
66			}
67
68			}
69
70
71			//retrive output filename of z force
72			data = info->getProperty("zconsfilename");
73			if(!data) {
74
75
76			sprintf( painCave.errMsg,
77			"ZConstraint error: If you use an ZConstraint\n"
78			" , you must set output filename of z-force.\n");
79			painCave.isFatal = 1;
80			simError();
81
82			}
83			else{
84
85			filename = dynamic_cast<StringData*>(data);
86
87			if(!filename){
88
89			sprintf( painCave.errMsg,
90			"ZConstraint error: Can not get property from SimInfo\n");
91			painCave.isFatal = 1;
92			simError();
93
94			}
95			else{
96			this->zconsOutput = filename->getData();
97			}
98
99
100			}
101
102
103			//calculate reference z coordinate for z-constraint molecules
104			double totalMass_local;
105			double totalMass;
106			double totalMZ_local;
107			double totalMZ;
108			double massOfUncons_local;
109			double massOfCurMol;
110			double COM[3];
111
112			totalMass_local = 0;
113			totalMass = 0;
114			totalMZ_local = 0;
115			totalMZ = 0;
116			massOfUncons_local = 0;
117
118
119			for(int i = 0; i < nMols; i++){
120			massOfCurMol = molecules[i].getTotalMass();
121			molecules[i].getCOM(COM);
122
123			totalMass_local += massOfCurMol;
124			totalMZ_local += massOfCurMol * COM[2];
125
126			if(isZConstraintMol(&molecules[i]) == -1){
127
128			massOfUncons_local += massOfCurMol;
129			}
130
131			}
132
133
134			#ifdef IS_MPI
135			MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
136			MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
137			MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
138			#else
139			totalMass = totalMass_local;
140			totalMZ = totalMZ_local;
141			totalMassOfUncons = massOfUncons_local;
142			#endif
143
144			double zsys;
145			zsys = totalMZ / totalMass;
146
147			#ifndef IS_MPI
148			for(int i = 0; i < nMols; i++){
149
150			if(isZConstraintMol(&molecules[i]) > -1 ){
151			molecules[i].getCOM(COM);
152			allRefZ.push_back(COM[2] - zsys);
153			}
154
155			}
156			#else
157
158			int whichNode;
159			enum CommType { RequestMolZPos, EndOfRequest} status;
160			//int status;
161			double zpos;
162			int localIndex;
163			MPI_Status ierr;
164			int tag = 0;
165
166			if(worldRank == 0){
167
168			int globalIndexOfCurMol;
169			int *MolToProcMap;
170			MolToProcMap = mpiSim->getMolToProcMap();
171
172			for(int i = 0; i < indexOfAllZConsMols.size(); i++){
173
174			whichNode = MolToProcMap[indexOfAllZConsMols[i]];
175			globalIndexOfCurMol = indexOfAllZConsMols[i];
176
177			if(whichNode == 0){
178
179			for(int j = 0; j < nMols; j++)
180			if(molecules[j].getGlobalIndex() == globalIndexOfCurMol){
181			localIndex = j;
182			break;
183			}
184
185			molecules[localIndex].getCOM(COM);
186			allRefZ.push_back(COM[2] - zsys);
187
188			}
189			else{
190			status = RequestMolZPos;
191			MPI_Send(&status, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
192			MPI_Send(&globalIndexOfCurMol, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
193			MPI_Recv(&zpos, 1, MPI_DOUBLE_PRECISION, whichNode, tag, MPI_COMM_WORLD, &ierr);
194
195			allRefZ.push_back(zpos - zsys);
196
197			}
198
199			} //End of Request Loop
200
201			//Send ending request message to slave nodes
202			status = EndOfRequest;
203			for(int i =1; i < mpiSim->getNumberProcessors(); i++)
204			MPI_Send(&status, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
205
206			}
207			else{
208
209			int whichMol;
210			bool done = false;
211
212			while (!done){
213
214			MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
215
216			switch (status){
217
218			case RequestMolZPos :
219
220			MPI_Recv(&whichMol, 1, MPI_INT, 0, tag, MPI_COMM_WORLD,&ierr);
221
222			for(int i = 0; i < nMols; i++)
223			if(molecules[i].getGlobalIndex() == whichMol){
224			localIndex = i;
225			break;
226			}
227
228			molecules[localIndex].getCOM(COM);
229			zpos = COM[2];
230			MPI_Send(&zpos, 1, MPI_DOUBLE_PRECISION, 0, tag, MPI_COMM_WORLD);
231
232			break;
233
234			case EndOfRequest :
235
236			done = true;
237			break;
238			}
239
240			}
241
242			}
243
244			//Brocast the allRefZ to slave nodes;
245			double* allRefZBuf;
246			int nZConsMols;
247			nZConsMols = indexOfAllZConsMols.size();
248
249			allRefZBuf = new double[nZConsMols];
250
251			if(worldRank == 0){
252
253			for(int i = 0; i < nZConsMols; i++)
254			allRefZBuf[i] = allRefZ[i];
255			}
256
257			MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
258
259			if(worldRank != 0){
260
261			for(int i = 0; i < nZConsMols; i++)
262			allRefZ.push_back(allRefZBuf[i]);
263			}
264
265			delete[] allRefZBuf;
266			#endif
267
268
269			#ifdef IS_MPI
270			update();
271			#else
272			int searchResult;
273
274			refZ = allRefZ;
275
276			for(int i = 0; i < nMols; i++){
277
278			searchResult = isZConstraintMol(&molecules[i]);
279
280			if(searchResult > -1){
281
282			zconsMols.push_back(&molecules[i]);
283			massOfZConsMols.push_back(molecules[i].getTotalMass());
284
285			molecules[i].getCOM(COM);
286			}
287			else
288			{
289
290			unconsMols.push_back(&molecules[i]);
291			massOfUnconsMols.push_back(molecules[i].getTotalMass());
292
293			}
294			}
295
296			fz = new double[zconsMols.size()];
297			indexOfZConsMols = new int [zconsMols.size()];
298
299			if(!fz \|\| !indexOfZConsMols){
300			sprintf( painCave.errMsg,
301			"Memory allocation failure in class Zconstraint\n");
302			painCave.isFatal = 1;
303			simError();
304			}
305
306			for(int i = 0; i < zconsMols.size(); i++)
307			indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
308
309			#endif
310
311			fzOut = new ZConsWriter(zconsOutput.c_str());
312
313			if(!fzOut){
314			sprintf( painCave.errMsg,
315			"Memory allocation failure in class Zconstraint\n");
316			painCave.isFatal = 1;
317			simError();
318			}
319
320			fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
321			}
322
323			template<typename T> ZConstraint<T>::~ZConstraint()
324			{
325			if(fz)
326			delete[] fz;
327
328			if(indexOfZConsMols)
329			delete[] indexOfZConsMols;
330
331			if(fzOut)
332			delete fzOut;
333			}
334
335			#ifdef IS_MPI
336			template<typename T> void ZConstraint<T>::update()
337			{
338			double COM[3];
339			int index;
340
341			zconsMols.clear();
342			massOfZConsMols.clear();
343			refZ.clear();
344
345			unconsMols.clear();
346			massOfUnconsMols.clear();
347
348
349			//creat zconsMol and unconsMol lists
350			for(int i = 0; i < nMols; i++){
351
352			index = isZConstraintMol(&molecules[i]);
353
354			if(index > -1){
355
356			zconsMols.push_back(&molecules[i]);
357			massOfZConsMols.push_back(molecules[i].getTotalMass());
358
359			molecules[i].getCOM(COM);
360			refZ.push_back(allRefZ[index]);
361			}
362			else
363			{
364
365			unconsMols.push_back(&molecules[i]);
366			massOfUnconsMols.push_back(molecules[i].getTotalMass());
367
368			}
369			}
370
371			//The reason to declare fz and indexOfZconsMols as pointer to array is
372			// that we want to make the MPI communication simple
373			if(fz)
374			delete[] fz;
375
376			if(indexOfZConsMols)
377			delete[] indexOfZConsMols;
378
379			if (zconsMols.size() > 0){
380			fz = new double[zconsMols.size()];
381			indexOfZConsMols = new int[zconsMols.size()];
382
383			if(!fz \|\| !indexOfZConsMols){
384			sprintf( painCave.errMsg,
385			"Memory allocation failure in class Zconstraint\n");
386			painCave.isFatal = 1;
387			simError();
388			}
389
390			for(int i = 0; i < zconsMols.size(); i++){
391			indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
392			}
393
394			}
395			else{
396			fz = NULL;
397			indexOfZConsMols = NULL;
398			}
399
400			}
401
402			#endif
403
404			/** Function Name: isZConstraintMol
405			** Parameter
406			** Molecule* mol
407			** Return value:
408			** -1, if the molecule is not z-constraint molecule,
409			** other non-negative values, its index in indexOfAllZConsMols vector
410			*/
411
412			template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
413			{
414			int index;
415			int low;
416			int high;
417			int mid;
418
419			index = mol->getGlobalIndex();
420
421			low = 0;
422			high = indexOfAllZConsMols.size() - 1;
423
424			//Binary Search (we have sorted the array)
425			while(low <= high){
426			mid = (low + high) /2;
427			if (indexOfAllZConsMols[mid] == index)
428			return mid;
429			else if (indexOfAllZConsMols[mid] > index )
430			high = mid -1;
431			else
432			low = mid + 1;
433			}
434
435			return -1;
436			}
437
438			/** Function Name: integrateStep
439			** Parameter:
440			** int calcPot;
441			** int calcStress;
442			** Description:
443			** Advance One Step.
444			** Memo:
445			** The best way to implement z-constraint is to override integrateStep
446			** Overriding constrainB is not a good choice, since in integrateStep,
447			** constrainB is invoked by below line,
448			** if(nConstrained) constrainB();
449			** For instance, we would like to apply z-constraint without bond contrain,
450			** In that case, if we override constrainB, Z-constrain method will never be executed;
451			*/
452			template<typename T> void ZConstraint<T>::integrateStep( int calcPot, int calcStress )
453			{
454			T::integrateStep( calcPot, calcStress );
455			resetZ();
456
457			double currZConsTime = 0;
458
459			//write out forces of z constraint
460			if( info->getTime() >= currZConsTime){
461			fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
462			}
463			}
464
465			/** Function Name: resetZ
466			** Description:
467			** Reset the z coordinates
468			*/
469
470			template<typename T> void ZConstraint<T>::resetZ()
471			{
472			double deltaZ;
473			double mzOfZCons; //total sum of m*z of z-constrain molecules
474			double mzOfUncons; //total sum of m*z of unconstrain molecuels;
475			double totalMZOfZCons;
476			double totalMZOfUncons;
477			double COM[3];
478			double zsys;
479			Atom** zconsAtoms;
480
481			mzOfZCons = 0;
482			mzOfUncons = 0;
483
484			for(int i = 0; i < zconsMols.size(); i++){
485			mzOfZCons += massOfZConsMols[i] * refZ[i];
486			}
487
488			#ifdef IS_MPI
489			MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
490			#else
491			totalMZOfZCons = mzOfZCons;
492			#endif
493
494			for(int i = 0; i < unconsMols.size(); i++){
495			unconsMols[i]->getCOM(COM);
496			mzOfUncons += massOfUnconsMols[i] * COM[2];
497			}
498
499			#ifdef IS_MPI
500			MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
501			#else
502			totalMZOfUncons = mzOfUncons;
503			#endif
504
505			zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons;
506
507			cout << "current time: " << info->getTime() <<endl;
508			for(int i = 0; i < zconsMols.size(); i++){
509
510			zconsMols[i]->getCOM(COM);
511
512			cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t";
513			deltaZ = zsys + refZ[i] - COM[2];
514			cout << "\tdistance: " << COM[2] +deltaZ - zsys;
515			//update z coordinate
516			zconsAtoms = zconsMols[i]->getMyAtoms();
517			for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
518			zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ);
519			}
520
521			//calculate z constrain force
522			fz[i] = massOfZConsMols[i]* deltaZ / dt2;
523
524			cout << "\tforce: " << fz[i] << endl;
525			}
526
527
528			}