ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE/libmdtools/ZConstraint.cpp

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 658 by tim, Thu Jul 31 15:35:07 2003 UTC vs.
Revision 726 by tim, Tue Aug 26 20:37:30 2003 UTC

#	Line 1 | Line 1
1		#include "Integrator.hpp"
2		#include "simError.h"
3	<
3	>	#include <cmath>
4		template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff)
5	<	: T(theInfo, the_ff), fz(NULL), indexOfZConsMols(NULL)
5	>	: T(theInfo, the_ff), fz(NULL), curZPos(NULL),
6	>	indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0)
7		{
8
9		//get properties from SimInfo
10		GenericData* data;
11	<	IndexData* index;
11	>	ZConsParaData* zConsParaData;
12		DoubleData* sampleTime;
13	+	DoubleData* tolerance;
14	+	StringData* policy;
15		StringData* filename;
16	+	double COM[3];
17	+
18	+	//by default, the direction of constraint is z
19	+	// 0 --> x
20	+	// 1 --> y
21	+	// 2 --> z
22	+	whichDirection = 2;
23	+
24	+	//estimate the force constant of harmonical potential
25	+	double Kb = 1.986E-3 ; //in kcal/K
26
27	<
28	<	data = info->getProperty("zconsindex");
16	<	if(!data) {
27	>	double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2;
28	>	zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox);
29
30	+	//creat force substraction policy
31	+	data = info->getProperty(ZCONSFORCEPOLICY_ID);
32	+	if(!data){
33		sprintf( painCave.errMsg,
34	<	"ZConstraint error: If you use an ZConstraint\n"
35	<	" , you must set index of z-constraint molecules.\n");
36	<	painCave.isFatal = 1;
37	<	simError();
34	>	"ZConstraint Warning: User does not set force substraction policy, "
35	>	"PolicyByMass is used\n");
36	>	painCave.isFatal = 0;
37	>	simError();
38	>
39	>	forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
40		}
41		else{
42	<	index = dynamic_cast<IndexData*>(data);
42	>	policy = dynamic_cast<StringData*>(data);
43
44	<	if(!index){
28	<
44	>	if(!policy){
45		sprintf( painCave.errMsg,
46	<	"ZConstraint error: Can not get property from SimInfo\n");
47	<	painCave.isFatal = 1;
48	<	simError();
49	<
46	>	"ZConstraint Error: Convertion from GenericData to StringData failure, "
47	>	"PolicyByMass is used\n");
48	>	painCave.isFatal = 0;
49	>	simError();
50	>
51	>	forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
52		}
53		else{
54	<	indexOfAllZConsMols = index->getIndexData();
54	>	if(policy->getData() == "BYNUMBER")
55	>	forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
56	>	else if(policy->getData() == "BYMASS")
57	>	forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
58	>	else{
59	>	sprintf( painCave.errMsg,
60	>	"ZConstraint Warning: unknown force substraction policy, "
61	>	"average force substraction policy is used\n");
62	>	painCave.isFatal = 0;
63	>	simError();
64	>	}
65		}
38	–
66		}
67
41	–	//retrive sample time of z-contraint
42	–	data = info->getProperty("zconstime");
68
69	+	//retrieve sample time of z-contraint
70	+	data = info->getProperty(ZCONSTIME_ID);
71	+
72		if(!data) {
73
74		sprintf( painCave.errMsg,
#	Line 67 | Line 95 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
95
96		}
97
98	<
99	<	//retrive output filename of z force
72	<	data = info->getProperty("zconsfilename");
98	>	//retrieve output filename of z force
99	>	data = info->getProperty(ZCONSFILENAME_ID);
100		if(!data) {
101
102
#	Line 82 | Line 109 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
109		}
110		else{
111
112	<	filename = dynamic_cast<StringData*>(data);
112	>	filename = dynamic_cast<StringData*>(data);
113
114		if(!filename){
115
#	Line 96 | Line 123 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
123		this->zconsOutput = filename->getData();
124		}
125
99	–
126		}
127
128	<
129	<	//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];
128	>	//retrieve tolerance for z-constraint molecuels
129	>	data = info->getProperty(ZCONSTOL_ID);
130
131	<	totalMass_local = 0;
132	<	totalMass = 0;
133	<	totalMZ_local = 0;
134	<	totalMZ = 0;
135	<	massOfUncons_local = 0;
136	<
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	>	if(!data) {
132	>
133	>	sprintf( painCave.errMsg,
134	>	"ZConstraint error: can not get tolerance \n");
135	>	painCave.isFatal = 1;
136	>	simError();
137		}
138	+	else{
139
140	<
141	<	#ifdef IS_MPI
142	<	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
140	>	tolerance = dynamic_cast<DoubleData*>(data);
141	>
142	>	if(!tolerance){
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);
144	>	sprintf( painCave.errMsg,
145	>	"ZConstraint error: Can not get property from SimInfo\n");
146	>	painCave.isFatal = 1;
147	>	simError();
148	>
149		}
150	<
150	>	else{
151	>	this->zconsTol = tolerance->getData();
152	>	}
153	>
154		}
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;
155
156	<	if(worldRank == 0){
157	<
158	<	int globalIndexOfCurMol;
159	<	int *MolToProcMap;
160	<	MolToProcMap = mpiSim->getMolToProcMap();
161	<
162	<	for(int i = 0; i < indexOfAllZConsMols.size(); i++){
163	<
164	<	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	<
156	>	//retrieve index of z-constraint molecules
157	>	data = info->getProperty(ZCONSPARADATA_ID);
158	>	if(!data) {
159	>
160	>	sprintf( painCave.errMsg,
161	>	"ZConstraint error: If you use an ZConstraint\n"
162	>	" , you must set index of z-constraint molecules.\n");
163	>	painCave.isFatal = 1;
164	>	simError();
165		}
166		else{
167
168	<	int whichMol;
169	<	bool done = false;
168	>	zConsParaData = dynamic_cast<ZConsParaData*>(data);
169	>
170	>	if(!zConsParaData){
171
172	<	while (!done){
173	<
174	<	MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
172	>	sprintf( painCave.errMsg,
173	>	"ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
174	>	painCave.isFatal = 1;
175	>	simError();
176
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	–
177		}
178	<
179	<	}
178	>	else{
179	>
180	>	parameters = zConsParaData->getData();
181
182	<	//Brocast the allRefZ to slave nodes;
183	<	double* allRefZBuf;
184	<	int nZConsMols;
247	<	nZConsMols = indexOfAllZConsMols.size();
248	<
249	<	allRefZBuf = new double[nZConsMols];
250	<
251	<	if(worldRank == 0){
182	>	//check the range of zconsIndex
183	>	//and the minimum value of index is the first one (we already sorted the data)
184	>	//the maximum value of index is the last one
185
186	<	for(int i = 0; i < nZConsMols; i++)
187	<	allRefZBuf[i] = allRefZ[i];
188	<	}
189	<
190	<	MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
191	<
192	<	if(worldRank != 0){
193	<
194	<	for(int i = 0; i < nZConsMols; i++)
195	<	allRefZ.push_back(allRefZBuf[i]);
196	<	}
197	<
198	<	delete[] allRefZBuf;
186	>	int maxIndex;
187	>	int minIndex;
188	>	int totalNumMol;
189	>
190	>	minIndex = (*parameters)[0].zconsIndex;
191	>	if(minIndex < 0){
192	>	sprintf( painCave.errMsg,
193	>	"ZConstraint error: index is out of range\n");
194	>	painCave.isFatal = 1;
195	>	simError();
196	>	}
197	>
198	>	maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
199	>
200	>	#ifndef IS_MPI
201	>	totalNumMol = nMols;
202	>	#else
203	>	totalNumMol = mpiSim->getTotNmol();
204	>	#endif
205	>
206	>	if(maxIndex > totalNumMol - 1){
207	>	sprintf( painCave.errMsg,
208	>	"ZConstraint error: index is out of range\n");
209	>	painCave.isFatal = 1;
210	>	simError();
211	>	}
212	>
213	>	//if user does not specify the zpos for the zconstraint molecule
214	>	//its initial z coordinate  will be used as default
215	>	for(int i = 0; i < parameters->size(); i++){
216	>
217	>	if(!(*parameters)[i].havingZPos){
218	>	#ifndef IS_MPI
219	>	for(int j = 0; j < nMols; j++){
220	>	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
221	>	molecules[j].getCOM(COM);
222	>	break;
223	>	}
224	>	}
225	>	#else
226	>	//query which processor current zconstraint molecule belongs to
227	>	int *MolToProcMap;
228	>	int whichNode;
229	>	double initZPos;
230	>	MolToProcMap = mpiSim->getMolToProcMap();
231	>	whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
232	>
233	>	//broadcast the zpos of current z-contraint molecule
234	>	//the node which contain this
235	>
236	>	if (worldRank == whichNode ){
237	>
238	>	for(int j = 0; j < nMols; j++)
239	>	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
240	>	molecules[j].getCOM(COM);
241	>	break;
242	>	}
243	>
244	>	}
245	>
246	>	MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);
247		#endif
248	+
249	+	(*parameters)[i].zPos = COM[whichDirection];
250
251	<
251	>	sprintf( painCave.errMsg,
252	>	"ZConstraint warningr: Does not specify zpos for z-constraint molecule "
253	>	"initial z coornidate will be used \n");
254	>	painCave.isFatal = 0;
255	>	simError();
256	>
257	>	}
258	>	}
259	>
260	>	}//end if (!zConsParaData)
261	>	}//end  if (!data)
262	>
263	>	//
264		#ifdef IS_MPI
265		update();
266		#else
267		int searchResult;
268	<
274	<	refZ = allRefZ;
275	<
268	>
269		for(int i = 0; i < nMols; i++){
270
271		searchResult = isZConstraintMol(&molecules[i]);
#	Line 281 | Line 274 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
274
275		zconsMols.push_back(&molecules[i]);
276		massOfZConsMols.push_back(molecules[i].getTotalMass());
277	<
277	>
278	>	zPos.push_back((*parameters)[searchResult].zPos);
279	>	//       cout << "index: "<< (*parameters)[searchResult].zconsIndex
280	>	//              <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
281	>	kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
282	>
283		molecules[i].getCOM(COM);
284		}
285		else
#	Line 294 | Line 292 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
292		}
293
294		fz = new double[zconsMols.size()];
295	+	curZPos = new double[zconsMols.size()];
296		indexOfZConsMols = new int [zconsMols.size()];
297
298	<	if(!fz || !indexOfZConsMols){
298	>	if(!fz || !curZPos || !indexOfZConsMols){
299		sprintf( painCave.errMsg,
300		"Memory allocation failure in class Zconstraint\n");
301		painCave.isFatal = 1;
302		simError();
303		}
304
305	<	for(int i = 0; i < zconsMols.size(); i++)
305	>	//determine the states of z-constraint molecules
306	>	for(int i = 0; i < zconsMols.size(); i++){
307		indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
308	+
309	+	zconsMols[i]->getCOM(COM);
310	+	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
311	+	states.push_back(zcsFixed);
312	+	else
313	+	states.push_back(zcsMoving);
314	+	}
315
316		#endif
317	+
318	+	//get total masss of unconstraint molecules
319	+	double totalMassOfUncons_local;
320	+	totalMassOfUncons_local = 0;
321
322	<	fzOut = new ZConsWriter(zconsOutput.c_str());
322	>	for(int i = 0; i < unconsMols.size(); i++)
323	>	totalMassOfUncons_local += unconsMols[i]->getTotalMass();
324	>
325	>	#ifndef IS_MPI
326	>	totalMassOfUncons = totalMassOfUncons_local;
327	>	#else
328	>	MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1,
329	>	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
330	>	#endif
331	>
332	>
333	>	//get total number of unconstrained atoms
334	>	int nUnconsAtoms_local;
335	>	nUnconsAtoms_local = 0;
336	>	for(int i = 0; i < unconsMols.size(); i++)
337	>	nUnconsAtoms_local += unconsMols[i]->getNAtoms();
338	>
339	>	#ifndef IS_MPI
340	>	totNumOfUnconsAtoms = nUnconsAtoms_local;
341	>	#else
342	>	MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1,
343	>	MPI_INT,MPI_SUM, MPI_COMM_WORLD);
344	>	#endif
345	>
346	>	// creat zconsWriter
347	>	fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);
348
349		if(!fzOut){
350		sprintf( painCave.errMsg,
#	Line 316 | Line 352 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
352		painCave.isFatal = 1;
353		simError();
354		}
355	<
356	<	fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
355	>
356	>	forcePolicy->update();
357		}
358
359		template<typename T> ZConstraint<T>::~ZConstraint()
360		{
361		if(fz)
362		delete[] fz;
363	+
364	+	if(curZPos)
365	+	delete[] curZPos;
366
367		if(indexOfZConsMols)
368		delete[] indexOfZConsMols;
369
370		if(fzOut)
371		delete fzOut;
372	+
373	+	if(forcePolicy)
374	+	delete forcePolicy;
375		}
376
377	+
378	+	/**
379	+	*
380	+	*/
381	+
382		#ifdef IS_MPI
383		template<typename T> void ZConstraint<T>::update()
384		{
#	Line 340 | Line 387 | template<typename T> void ZConstraint<T>::update()
387
388		zconsMols.clear();
389		massOfZConsMols.clear();
390	<	refZ.clear();
390	>	zPos.clear();
391	>	kz.clear();
392
393		unconsMols.clear();
394		massOfUnconsMols.clear();
#	Line 354 | Line 402 | template<typename T> void ZConstraint<T>::update()
402		if(index > -1){
403
404		zconsMols.push_back(&molecules[i]);
405	+	zPos.push_back((*parameters)[index].zPos);
406	+	kz.push_back((*parameters)[index].kRatio * zForceConst);
407	+
408		massOfZConsMols.push_back(molecules[i].getTotalMass());
409
410		molecules[i].getCOM(COM);
360	–	refZ.push_back(allRefZ[index]);
411		}
412		else
413		{
#	Line 367 | Line 417 | template<typename T> void ZConstraint<T>::update()
417
418		}
419		}
420	+
421	+	//determine the states of z-constraint molecules
422	+	for(int i = 0; i < zconsMols.size(); i++){
423	+	zconsMols[i]->getCOM(COM);
424	+	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
425	+	states.push_back(zcsFixed);
426	+	else
427	+	states.push_back(zcsMoving);
428	+	}
429	+
430
431		//The reason to declare fz and indexOfZconsMols as pointer to array is
432		// that we want to make the MPI communication simple
433		if(fz)
434		delete[] fz;
435	+
436	+	if(curZPos)
437	+	delete[] curZPos;
438
439		if(indexOfZConsMols)
440		delete[] indexOfZConsMols;
441
442		if (zconsMols.size() > 0){
443		fz = new double[zconsMols.size()];
444	+	curZPos = new double[zconsMols.size()];
445		indexOfZConsMols =  new int[zconsMols.size()];
446
447	<	if(!fz || !indexOfZConsMols){
447	>	if(!fz || !curZPos || !indexOfZConsMols){
448		sprintf( painCave.errMsg,
449		"Memory allocation failure in class Zconstraint\n");
450		painCave.isFatal = 1;
#	Line 394 | Line 458 | template<typename T> void ZConstraint<T>::update()
458		}
459		else{
460		fz = NULL;
461	+	curZPos = NULL;
462		indexOfZConsMols = NULL;
463		}
464
465	+	//
466	+	forcePolicy->update();
467	+
468		}
469
470		#endif
471
472	<	/**  Function Name: isZConstraintMol
473	<	**  Parameter
474	<	**    Molecule* mol
475	<	**  Return value:
476	<	**    -1, if the molecule is not z-constraint molecule,
477	<	**    other non-negative values, its index in indexOfAllZConsMols vector
472	>	/**
473	>	*  Function Name: isZConstraintMol
474	>	*  Parameter
475	>	*    Molecule* mol
476	>	*  Return value:
477	>	*    -1, if the molecule is not z-constraint molecule,
478	>	*    other non-negative values, its index in indexOfAllZConsMols vector
479		*/
480
481		template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
#	Line 419 | Line 488 | template<typename T> int ZConstraint<T>::isZConstraint
488		index = mol->getGlobalIndex();
489
490		low = 0;
491	<	high = indexOfAllZConsMols.size() - 1;
491	>	high = parameters->size() - 1;
492
493		//Binary Search (we have sorted the array)
494		while(low <= high){
495		mid = (low + high) /2;
496	<	if (indexOfAllZConsMols[mid] == index)
496	>	if ((*parameters)[mid].zconsIndex == index)
497		return mid;
498	<	else if (indexOfAllZConsMols[mid] > index )
498	>	else if ((*parameters)[mid].zconsIndex > index )
499		high = mid -1;
500		else
501		low = mid + 1;
#	Line 435 | Line 504 | template<typename T> int ZConstraint<T>::isZConstraint
504		return -1;
505		}
506
507	<	/** 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();
507	>	template<typename T> void ZConstraint<T>::integrate(){
508
509	<	double currZConsTime = 0;
510	<
511	<	//write out forces of z constraint
512	<	if( info->getTime() >= currZConsTime){
513	<	fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
514	<	}
509	>	//zero out the velocities of center of mass of unconstrained molecules
510	>	//and the velocities of center of mass of every single z-constrained molecueles
511	>	zeroOutVel();
512	>
513	>	curZconsTime = zconsTime + info->getTime();
514	>
515	>	T::integrate();
516	>
517		}
518	+
519
520	<	/** Function Name: resetZ
521	<	** Description:
522	<	**  Reset the z coordinates
523	<	*/
520	>	/**
521	>	*
522	>	*
523	>	*
524	>	*
525	>	*/
526	>	template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
527	>	double zsys;
528	>	double COM[3];
529	>	double force[3];
530	>	double zSysCOMVel;
531
532	<	template<typename T> void ZConstraint<T>::resetZ()
532	>	T::calcForce(calcPot, calcStress);
533	>
534	>	if (checkZConsState()){
535	>
536	>	#ifdef IS_MPI
537	>	if(worldRank == 0){
538	>	#endif
539	>	//       std::cerr << "\n"
540	>	//              << "*******************************************\n"
541	>	//              << " about to call zeroOutVel()\n"
542	>	//              << "*******************************************\n"
543	>	//              << "\n";
544	>	#ifdef IS_MPI
545	>	}
546	>	#endif
547	>	zeroOutVel();
548	>
549	>	#ifdef IS_MPI
550	>	if(worldRank == 0){
551	>	#endif
552	>	//       std::cerr << "\n"
553	>	//              << "*******************************************\n"
554	>	//              << " finished zeroOutVel()\n"
555	>	//              << "*******************************************\n"
556	>	//              << "\n";
557	>	#ifdef IS_MPI
558	>	}
559	>	#endif
560	>
561	>	forcePolicy->update();
562	>	}
563	>
564	>	zsys = calcZSys();
565	>	zSysCOMVel = calcSysCOMVel();
566	>	#ifdef IS_MPI
567	>	if(worldRank == 0){
568	>	#endif
569	>	//     cout << "---------------------------------------------------------------------" <<endl;
570	>	//     cout << "current time: " << info->getTime() << endl;
571	>	//     cout << "center of mass at z: " << zsys << endl;
572	>	//     cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
573	>
574	>	#ifdef IS_MPI
575	>	}
576	>	#endif
577	>
578	>	//do zconstraint force;
579	>	if (haveFixedZMols())
580	>	this->doZconstraintForce();
581	>
582	>	//use harmonical poteintial to move the molecules to the specified positions
583	>	if (haveMovingZMols())
584	>	this->doHarmonic();
585	>
586	>	//write out forces and current positions of z-constraint molecules
587	>	if(info->getTime() >= curZconsTime){
588	>	for(int i = 0; i < zconsMols.size(); i++){
589	>	zconsMols[i]->getCOM(COM);
590	>	curZPos[i] = COM[whichDirection];
591	>
592	>	//if the z-constraint molecule is still moving, just record its force
593	>	if(states[i] == zcsMoving){
594	>	fz[i] = 0;
595	>	Atom** movingZAtoms;
596	>	movingZAtoms = zconsMols[i]->getMyAtoms();
597	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
598	>	movingZAtoms[j]->getFrc(force);
599	>	fz[i] += force[whichDirection];
600	>	}
601	>	}
602	>	}
603	>	fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos);
604	>	curZconsTime += zconsTime;
605	>	}
606	>
607	>	zSysCOMVel = calcSysCOMVel();
608	>	#ifdef IS_MPI
609	>	if(worldRank == 0){
610	>	#endif
611	>	//    cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
612	>	#ifdef IS_MPI
613	>	}
614	>	#endif
615	>
616	>	}
617	>
618	>
619	>	/**
620	>	*
621	>	*/
622	>
623	>	template<typename T> double ZConstraint<T>::calcZSys()
624		{
625	<	double deltaZ;
626	<	double mzOfZCons;   //total sum of m*z of z-constrain molecules
627	<	double mzOfUncons; //total sum of m*z of unconstrain molecuels;
628	<	double totalMZOfZCons;
629	<	double totalMZOfUncons;
625	>	//calculate reference z coordinate for z-constraint molecules
626	>	double totalMass_local;
627	>	double totalMass;
628	>	double totalMZ_local;
629	>	double totalMZ;
630	>	double massOfCurMol;
631		double COM[3];
632	<	double zsys;
633	<	Atom** zconsAtoms;
632	>
633	>	totalMass_local = 0;
634	>	totalMZ_local = 0;
635	>
636	>	for(int i = 0; i < nMols; i++){
637	>	massOfCurMol = molecules[i].getTotalMass();
638	>	molecules[i].getCOM(COM);
639	>
640	>	totalMass_local += massOfCurMol;
641	>	totalMZ_local += massOfCurMol * COM[whichDirection];
642
643	<	mzOfZCons = 0;
644	<	mzOfUncons  = 0;
643	>	}
644	>
645
646	<	for(int i = 0; i < zconsMols.size(); i++){
647	<	mzOfZCons += massOfZConsMols[i] * refZ[i];
646	>	#ifdef IS_MPI
647	>	MPI_Allreduce(&totalMass_local, &totalMass, 1,
648	>	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
649	>	MPI_Allreduce(&totalMZ_local, &totalMZ, 1,
650	>	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
651	>	#else
652	>	totalMass = totalMass_local;
653	>	totalMZ = totalMZ_local;
654	>	#endif
655	>
656	>	double zsys;
657	>	zsys = totalMZ / totalMass;
658	>
659	>	return zsys;
660	>	}
661	>
662	>	/**
663	>	*
664	>	*/
665	>	template<typename T> void ZConstraint<T>::thermalize( void ){
666	>
667	>	T::thermalize();
668	>	zeroOutVel();
669	>	}
670	>
671	>	/**
672	>	*
673	>	*/
674	>
675	>	template<typename T> void ZConstraint<T>::zeroOutVel(){
676	>
677	>	Atom** fixedZAtoms;
678	>	double COMvel[3];
679	>	double vel[3];
680	>	double zSysCOMVel;
681	>
682	>	//zero out the velocities of center of mass of fixed z-constrained molecules
683	>
684	>	for(int i = 0; i < zconsMols.size(); i++){
685	>
686	>	if (states[i] == zcsFixed){
687	>
688	>	zconsMols[i]->getCOMvel(COMvel);
689	>	//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
690	>
691	>	fixedZAtoms = zconsMols[i]->getMyAtoms();
692	>
693	>	for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
694	>	fixedZAtoms[j]->getVel(vel);
695	>	vel[whichDirection] -= COMvel[whichDirection];
696	>	fixedZAtoms[j]->setVel(vel);
697	>	}
698	>
699	>	zconsMols[i]->getCOMvel(COMvel);
700	>	//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
701	>	}
702	>
703		}
704
705	+	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
706	+
707	+	zSysCOMVel = calcSysCOMVel();
708		#ifdef IS_MPI
709	<	MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
709	>	if(worldRank == 0){
710	>	#endif
711	>	//     cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;
712	>	#ifdef IS_MPI
713	>	}
714	>	#endif
715	>
716	>	// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
717	>	double MVzOfMovingMols_local;
718	>	double MVzOfMovingMols;
719	>	double totalMassOfMovingZMols_local;
720	>	double totalMassOfMovingZMols;
721	>
722	>	MVzOfMovingMols_local = 0;
723	>	totalMassOfMovingZMols_local = 0;
724	>
725	>	for(int i =0; i < unconsMols.size(); i++){
726	>	unconsMols[i]->getCOMvel(COMvel);
727	>	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
728	>	}
729	>
730	>	for(int i = 0; i < zconsMols.size(); i++){
731	>	if (states[i] == zcsMoving){
732	>	zconsMols[i]->getCOMvel(COMvel);
733	>	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
734	>	totalMassOfMovingZMols_local += massOfZConsMols[i];
735	>	}
736	>
737	>	}
738	>
739	>	#ifndef IS_MPI
740	>	MVzOfMovingMols = MVzOfMovingMols_local;
741	>	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
742		#else
743	<	totalMZOfZCons = mzOfZCons;
743	>	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
744	>	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
745		#endif
746
747	+	double vzOfMovingMols;
748	+	vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
749	+
750	+	//modify the velocites of unconstrained molecules
751	+	Atom** unconsAtoms;
752		for(int i = 0; i < unconsMols.size(); i++){
495	–	unconsMols[i]->getCOM(COM);
496	–	mzOfUncons += massOfUnconsMols[i] * COM[2];
497	–	}
753
754	+	unconsAtoms = unconsMols[i]->getMyAtoms();
755	+	for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
756	+	unconsAtoms[j]->getVel(vel);
757	+	vel[whichDirection] -= vzOfMovingMols;
758	+	unconsAtoms[j]->setVel(vel);
759	+	}
760	+
761	+	}
762	+
763	+	//modify the velocities of moving z-constrained molecuels
764	+	Atom** movingZAtoms;
765	+	for(int i = 0; i < zconsMols.size(); i++){
766	+
767	+	if (states[i] ==zcsMoving){
768	+
769	+	movingZAtoms = zconsMols[i]->getMyAtoms();
770	+	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
771	+	movingZAtoms[j]->getVel(vel);
772	+	vel[whichDirection] -= vzOfMovingMols;
773	+	movingZAtoms[j]->setVel(vel);
774	+	}
775	+
776	+	}
777	+
778	+	}
779	+
780	+
781	+	zSysCOMVel = calcSysCOMVel();
782		#ifdef IS_MPI
783	<	MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
783	>	if(worldRank == 0){
784	>	#endif
785	>	//     cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;
786	>	#ifdef IS_MPI
787	>	}
788	>	#endif
789	>
790	>	}
791	>
792	>	/**
793	>	*
794	>	*/
795	>
796	>	template<typename T> void ZConstraint<T>::doZconstraintForce(){
797	>
798	>	Atom** zconsAtoms;
799	>	double totalFZ;
800	>	double totalFZ_local;
801	>	double COMvel[3];
802	>	double COM[3];
803	>	double force[3];
804	>
805	>
806	>
807	>	//constrain the molecules which do not reach the specified positions
808	>
809	>	//Zero Out the force of z-contrained molecules
810	>	totalFZ_local = 0;
811	>
812	>	//calculate the total z-contrained force of fixed z-contrained molecules
813	>
814	>	for(int i = 0; i < zconsMols.size(); i++){
815	>
816	>	if (states[i] == zcsFixed){
817	>
818	>	zconsMols[i]->getCOM(COM);
819	>	zconsAtoms = zconsMols[i]->getMyAtoms();
820	>
821	>	fz[i] = 0;
822	>	for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
823	>	zconsAtoms[j]->getFrc(force);
824	>	fz[i] += force[whichDirection];
825	>	}
826	>	totalFZ_local += fz[i];
827	>
828	>	//cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
829	>	//      <<"\tcurrent zpos: " << COM[whichDirection]
830	>	//       << "\tcurrent fz: " <<fz[i] << endl;
831	>
832	>
833	>	}
834	>
835	>	}
836	>
837	>	//calculate total z-constraint force
838	>	#ifdef IS_MPI
839	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
840		#else
841	<	totalMZOfUncons = mzOfUncons;
842	<	#endif
841	>	totalFZ = totalFZ_local;
842	>	#endif
843	>
844
845	<	zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons;
845	>	// apply negative to fixed z-constrained molecues;
846	>	force[0]= 0;
847	>	force[1]= 0;
848	>	force[2]= 0;
849
850	<	cout << "current time: " << info->getTime() <<endl;
851	<	for(int i = 0; i < zconsMols.size(); i++){
852	<
853	<	zconsMols[i]->getCOM(COM);
850	>	for(int i = 0; i < zconsMols.size(); i++){
851	>
852	>	if (states[i] == zcsFixed){
853	>
854	>	int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
855	>	zconsAtoms = zconsMols[i]->getMyAtoms();
856
857	<	cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t";
858	<	deltaZ = zsys + refZ[i] - COM[2];
859	<	cout << "\tdistance: " << COM[2] +deltaZ - zsys;
860	<	//update z coordinate
861	<	zconsAtoms = zconsMols[i]->getMyAtoms();
517	<	for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
518	<	zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ);
519	<	}
857	>	for(int j =0; j < nAtomOfCurZConsMol; j++) {
858	>	//force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
859	>	force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
860	>	zconsAtoms[j]->addFrc(force);
861	>	}
862
863	<	//calculate z constrain force
864	<	fz[i] = massOfZConsMols[i]* deltaZ / dt2;
863	>	}
864	>
865	>	}
866	>
867	>	//   cout << "after zero out z-constraint force on fixed z-constraint molecuels "
868	>	//        << "total force is " << calcTotalForce() << endl;
869	>
870	>	//calculate the number of atoms of moving z-constrained molecules
871	>	int nMovingZAtoms_local;
872	>	int nMovingZAtoms;
873	>
874	>	nMovingZAtoms_local = 0;
875	>	for(int i = 0; i < zconsMols.size(); i++)
876	>	if(states[i] == zcsMoving)
877	>	nMovingZAtoms_local += zconsMols[i]->getNAtoms();
878	>
879	>	#ifdef IS_MPI
880	>	MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1,
881	>	MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
882	>	#else
883	>	nMovingZAtoms = nMovingZAtoms_local;
884	>	#endif
885	>
886	>	force[0]= 0;
887	>	force[1]= 0;
888	>	force[2]= 0;
889	>
890	>	//modify the forces of unconstrained molecules
891	>	for(int i = 0; i < unconsMols.size(); i++){
892	>
893	>	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
894	>
895	>	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
896	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
897	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
898	>	unconsAtoms[j]->addFrc(force);
899	>	}
900	>
901	>	}
902	>
903	>	//modify the forces of moving z-constrained molecules
904	>	for(int i = 0; i < zconsMols.size(); i++) {
905	>	if (states[i] == zcsMoving){
906
907	<	cout << "\tforce: " << fz[i] << endl;
907	>	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
908	>
909	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
910	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
911	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
912	>	movingZAtoms[j]->addFrc(force);
913	>	}
914	>	}
915		}
916
917	+	//cout << "after substracting z-constraint force from moving molecuels "
918	+	//      << "total force is " << calcTotalForce()  << endl;
919	+
920	+	}
921	+
922	+	/**
923	+	*
924	+	*
925	+	*/
926	+
927	+	template<typename T> void ZConstraint<T>::doHarmonic(){
928	+	double force[3];
929	+	double harmonicU;
930	+	double harmonicF;
931	+	double COM[3];
932	+	double diff;
933	+	double totalFZ_local;
934	+	double totalFZ;
935	+
936	+	force[0] = 0;
937	+	force[1] = 0;
938	+	force[2] = 0;
939	+
940	+	totalFZ_local = 0;
941	+
942	+	for(int i = 0; i < zconsMols.size(); i++) {
943	+
944	+	if (states[i] == zcsMoving){
945	+	zconsMols[i]->getCOM(COM);
946	+	//       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
947	+	//         << "\tcurrent zpos: " << COM[whichDirection] << endl;
948	+
949	+	diff = COM[whichDirection] -zPos[i];
950	+
951	+	harmonicU = 0.5 * kz[i] * diff * diff;
952	+	info->lrPot += harmonicU;
953	+
954	+	harmonicF =  - kz[i] * diff;
955	+	totalFZ_local += harmonicF;
956	+
957	+	//adjust force
958	+
959	+	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
960	+
961	+	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
962	+	//force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
963	+	force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
964	+	movingZAtoms[j]->addFrc(force);
965	+	}
966	+	}
967	+
968	+	}
969	+
970	+	#ifndef IS_MPI
971	+	totalFZ = totalFZ_local;
972	+	#else
973	+	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
974	+	#endif
975	+
976	+	force[0]= 0;
977	+	force[1]= 0;
978	+	force[2]= 0;
979	+
980	+	//modify the forces of unconstrained molecules
981	+	for(int i = 0; i < unconsMols.size(); i++){
982	+
983	+	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
984	+
985	+	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
986	+	//force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
987	+	force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
988	+	unconsAtoms[j]->addFrc(force);
989	+	}
990	+	}
991	+
992	+	}
993	+
994	+	/**
995	+	*
996	+	*/
997	+
998	+	template<typename T> bool ZConstraint<T>::checkZConsState(){
999	+	double COM[3];
1000	+	double diff;
1001	+
1002	+	int changed_local;
1003	+	int changed;
1004	+
1005	+	changed_local = 0;
1006	+
1007	+	for(int i =0; i < zconsMols.size(); i++){
1008	+
1009	+	zconsMols[i]->getCOM(COM);
1010	+	diff = fabs(COM[whichDirection] - zPos[i]);
1011	+	if (  diff <= zconsTol && states[i] == zcsMoving){
1012	+	states[i] = zcsFixed;
1013	+	changed_local = 1;
1014	+	}
1015	+	else if ( diff > zconsTol && states[i] == zcsFixed){
1016	+	states[i] = zcsMoving;
1017	+	changed_local = 1;
1018	+	}
1019	+
1020	+	}
1021	+
1022	+	#ifndef IS_MPI
1023	+	changed =changed_local;
1024	+	#else
1025	+	MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1026	+	#endif
1027	+
1028	+	return (changed > 0);
1029	+
1030	+	}
1031	+
1032	+	template<typename T> bool ZConstraint<T>::haveFixedZMols(){
1033	+
1034	+	int havingFixed_local;
1035	+	int havingFixed;
1036	+
1037	+	havingFixed_local = 0;
1038	+
1039	+	for(int i = 0; i < zconsMols.size(); i++)
1040	+	if (states[i] == zcsFixed){
1041	+	havingFixed_local = 1;
1042	+	break;
1043	+	}
1044	+
1045	+	#ifndef IS_MPI
1046	+	havingFixed = havingFixed_local;
1047	+	#else
1048	+	MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1049	+	#endif
1050	+
1051	+	return (havingFixed > 0);
1052	+	}
1053	+
1054	+
1055	+	/**
1056	+	*
1057	+	*/
1058	+	template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1059	+
1060	+	int havingMoving_local;
1061	+	int havingMoving;
1062	+
1063	+	havingMoving_local = 0;
1064	+
1065	+	for(int i = 0; i < zconsMols.size(); i++)
1066	+	if (states[i] == zcsMoving){
1067	+	havingMoving_local = 1;
1068	+	break;
1069	+	}
1070	+
1071	+	#ifndef IS_MPI
1072	+	havingMoving = havingMoving_local;
1073	+	#else
1074	+	MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1075	+	#endif
1076	+
1077	+	return (havingMoving > 0);
1078	+
1079	+	}
1080	+
1081	+	/**
1082	+	*
1083	+	*/
1084	+
1085	+	template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
1086	+	{
1087	+	double MVzOfMovingMols_local;
1088	+	double MVzOfMovingMols;
1089	+	double totalMassOfMovingZMols_local;
1090	+	double totalMassOfMovingZMols;
1091	+	double COMvel[3];
1092
1093	+	MVzOfMovingMols_local = 0;
1094	+	totalMassOfMovingZMols_local = 0;
1095	+
1096	+	for(int i =0; i < unconsMols.size(); i++){
1097	+	unconsMols[i]->getCOMvel(COMvel);
1098	+	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1099	+	}
1100	+
1101	+	for(int i = 0; i < zconsMols.size(); i++){
1102	+
1103	+	if (states[i] == zcsMoving){
1104	+	zconsMols[i]->getCOMvel(COMvel);
1105	+	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
1106	+	totalMassOfMovingZMols_local += massOfZConsMols[i];
1107	+	}
1108	+
1109	+	}
1110	+
1111	+	#ifndef IS_MPI
1112	+	MVzOfMovingMols = MVzOfMovingMols_local;
1113	+	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1114	+	#else
1115	+	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1116	+	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1117	+	#endif
1118	+
1119	+	double vzOfMovingMols;
1120	+	vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
1121	+
1122	+	return vzOfMovingMols;
1123		}
1124	+
1125	+	/**
1126	+	*
1127	+	*/
1128	+
1129	+	template<typename T> double ZConstraint<T>::calcSysCOMVel()
1130	+	{
1131	+	double COMvel[3];
1132	+	double tempMVz_local;
1133	+	double tempMVz;
1134	+	double massOfZCons_local;
1135	+	double massOfZCons;
1136	+
1137	+
1138	+	tempMVz_local = 0;
1139	+
1140	+	for(int i =0 ; i < nMols; i++){
1141	+	molecules[i].getCOMvel(COMvel);
1142	+	tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
1143	+	}
1144	+
1145	+	massOfZCons_local = 0;
1146	+
1147	+	for(int i = 0; i < massOfZConsMols.size(); i++){
1148	+	massOfZCons_local += massOfZConsMols[i];
1149	+	}
1150	+	#ifndef IS_MPI
1151	+	massOfZCons = massOfZCons_local;
1152	+	tempMVz = tempMVz_local;
1153	+	#else
1154	+	MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1155	+	MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1156	+	#endif
1157	+
1158	+	return tempMVz /(totalMassOfUncons + massOfZCons);
1159	+	}
1160	+
1161	+	/**
1162	+	*
1163	+	*/
1164	+
1165	+	template<typename T> double ZConstraint<T>::calcTotalForce(){
1166	+
1167	+	double force[3];
1168	+	double totalForce_local;
1169	+	double totalForce;
1170	+
1171	+	totalForce_local = 0;
1172	+
1173	+	for(int i = 0; i < nAtoms; i++){
1174	+	atoms[i]->getFrc(force);
1175	+	totalForce_local += force[whichDirection];
1176	+	}
1177	+
1178	+	#ifndef IS_MPI
1179	+	totalForce = totalForce_local;
1180	+	#else
1181	+	MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1182	+	#endif
1183	+
1184	+	return totalForce;
1185	+
1186	+	}
1187	+
1188	+	/**
1189	+	*
1190	+	*/
1191	+
1192	+	template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1193	+	//calculate the number of atoms of moving z-constrained molecules
1194	+	int nMovingZAtoms_local;
1195	+	int nMovingZAtoms;
1196	+
1197	+	nMovingZAtoms_local = 0;
1198	+	for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1199	+	if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1200	+	nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1201	+
1202	+	#ifdef IS_MPI
1203	+	MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1204	+	#else
1205	+	nMovingZAtoms = nMovingZAtoms_local;
1206	+	#endif
1207	+	totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1208	+
1209	+	#ifdef IS_MPI
1210	+	if(worldRank == 0){
1211	+	#endif
1212	+	//    std::cerr << "\n"
1213	+	//            << "*******************************************\n"
1214	+	//            << " fiished Policy by numbr()\n"
1215	+	//            << "*******************************************\n"
1216	+	//            << "\n";
1217	+	#ifdef IS_MPI
1218	+	}
1219	+	#endif
1220	+	}
1221	+
1222	+	template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1223	+	return totalForce / mol->getNAtoms();
1224	+	}
1225	+
1226	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, double totalForce){
1227	+	return totalForce / totNumOfMovingAtoms;
1228	+	}
1229	+
1230	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1231	+	return totalForce / mol->getNAtoms();
1232	+	}
1233	+
1234	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, double totalForce){
1235	+	return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1236	+	}
1237	+
1238	+	/**
1239	+	*
1240	+	*/
1241	+
1242	+	template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1243	+	//calculate the number of atoms of moving z-constrained molecules
1244	+	double massOfMovingZAtoms_local;
1245	+	double massOfMovingZAtoms;
1246	+
1247	+	massOfMovingZAtoms_local = 0;
1248	+	for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1249	+	if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1250	+	massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1251	+
1252	+	#ifdef IS_MPI
1253	+	MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1254	+	#else
1255	+	massOfMovingZAtoms = massOfMovingZAtoms_local;
1256	+	#endif
1257	+	totMassOfMovingAtoms = massOfMovingZAtoms_local + zconsIntegrator->totalMassOfUncons;
1258	+	}
1259	+
1260	+	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1261	+	return totalForce * atom->getMass() / mol->getTotalMass();
1262	+	}
1263	+
1264	+	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols( Atom* atom, double totalForce){
1265	+	return totalForce * atom->getMass() / totMassOfMovingAtoms;
1266	+	}
1267	+
1268	+	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1269	+	return totalForce * atom->getMass() / mol->getTotalMass();
1270	+	}
1271	+
1272	+	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, double totalForce){
1273	+	return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1274	+	}
1275	+

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines