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